DISCRETE AND DIFFERENTIAL HOMOTOPY IN CIRCULAR RESTRICTED THREE-BODY CONTROL

Abstract. The planar circular restricted three-body problem is considered. The control enters linearly in the equation of motion to model the thrust of the third body. The minimum time optimal control problem has two scalar parameters: The ratio of the primaries masses which embeds the two-body problem into the three-body one, and the upper bound on the control norm. Regular extremals of the maximum principle are computed by shooting thanks to continuations with respect to both parameters. Discrete and di↵erential homotopy are compared in connection with second order su cient conditions in optimal control. Homotopy with respect to control bound gives evidence of various topological structures of extremals.

Introduction. An important simplification of the general three-body problem in celestial mechanics is the planar circular restricted problem: Two primaries (e.g., Earth and Moon) rotate around the center of mass, while a third body evolves in the plane of their Keplerian motion under their attraction, but not influencing them. In a rotating frame centered on the barycenter of the primaries, the equations of the controlled motion arë q + 2@V µ (q) + 2iq = "u, |u| = q u 2 1 + u 2 2  1, where q and u 2 C ' R 2 are respectively the position and the control vector. The potential is primary having no mass); (ii) µ = 1/2 for primaries with equal masses defining a circular choreography.
In cartesian coordinates (q,q), the system is defined on the tangent bundle T Q µ where Q µ is the complex plane with two punctures at the primaries, µ and 1 µ, (singularities of the potential). The set Q µ has the topology of the eight figure, with a non-trivial fundamental group (Z ⇤ Z). One can alternatively use Jacobi energy first integral, J µ (q,q) = |q| 2 2 + V µ (q), to write the uncontrolled system in Hamiltonian form. Setting p =q + iq, J µ (q, p) = |p| 2 2 + p^q 1 µ |q + µ| µ |q 1 + µ| and the controlled dynamics iṡ q = @J µ @p ,ṗ = @J µ @q + "u, |u|  1.
In control form,ẋ = F 0 (x) + "(u 1 F 1 (x) + u 2 F 2 (x)), |u|  1, with x 2 X µ , and X µ = T Q µ or X µ = T ⇤ Q µ (cotangent bundle) depending on the choice of coordinates. In (q,q) variables, We consider minimum time control of the system. (See [5] for the case of energy). According to Pontryagin maximum principle, optimal trajectories must be projections on X µ of extremals, that is of integral curves of the following Hamiltonian: The Hamiltonian and the lifts H i are defined on T ⇤ X µ , and there are two possibilities because of homogeneity in (p 0 , p) (p 0 is a nonpositive scalar): (i) p 0 6 = 0 (normal case), (ii) p 0 = 0 (abnormal case). Restricting to normal extremals, we set p 0 = 1. Moreover, the Hamiltonian has to be maximized almost everywhere with respect to the control, so u = (H 1 , H 2 )/|(H 1 , H 2 )| and outside the codimension two switching surface ⌃ = {H 1 = H 2 = 0}. (See [8] for a detailed analysis). The flow of the smooth Hamiltonian (1) is studied using homotopies on the two parameters of the problem. First, "à la Poincaré", with a continuation on parameter µ. This amounts to connecting (controlled) two-and three-body problems, the former being well understood [3]. Discrete and di↵erential homotopies on µ are compared in section 1. A neat framework for the latter is provided by the theory of second order conditions in optimal control. In section 2, homotopy on the control bound " is considered. Using both discrete and di↵erential homotopies, insight into the corresponding paths is obtained despite di culties due in particular to the topology of Q µ as discussed in the last section.
1. Homotopy on µ. Continuation on µ allows to embed the two-body controlled problem into the three-body one. The simplest approach is then to devise heuristically a sequence (µ k ) k from µ = 0 to the targeted value µ, for instance µ = 1/2. At step k, the solution is used to initialize the (k + 1)-th computation. Good initializations are indeed crucial to ensure convergence of the so-called shooting method based on Pontryagin maximum principle.
Since the final time to be minimized is unknown, extremals are to be sought on the zero level set of the Hamiltonian. One defines the exponential mapping is the projection on the x-space of the extremal departing from (x 0 , p 0 ). Expliciting the dependence with the homotopic parameter = µ (" is fixed), the shooting function is , where x f is the target (a choice of coordinates in which to express the terminal conditions is assumed). Both the initial and final points may vary with µ, here, hence the dependence of x 0 and x f in the parameter. 2 Numerical results of discrete homotopy to compute time minimal trajectories towards the L 2 Lagrange libration point (equilibrium point of the uncontrolled system whose position depend on µ, see [12]) are displayed in Fig. 1. A remarkable feature is that, starting from µ = 0, one discrete step su ces to reach µ ' 1.21e 2 corresponding to the Earth-Moon system. Local optimality of the computed extremals is checked using the theory of conjugate points that we now briefly sketch [1]. 1. The first phase around the first primary is essentially a two-body control phase. This is typically the case for the Earth-Moon system as µ ' 1.21e 2 is even smaller. Right, µ = 0.5. The perturbation due to the second primary now with equal mass is clearly observed as the eccentricity of intermediate orbits is increased.
An extremal z = (x, p) is regular if it exists ↵ > 0 such that the strong Legendre condition holds along z (below m is the dimension of the control), Here, as u belongs to S 1 , one can set u = (cos ✓, sin ✓) (✓ 2 R), and extremals are regular outside the switching surface as The variational equation A normal regular extremal without conjugate point is C 0 -locally optimal among trajectories with same endpoints.
By restricting if necessary the open subset ⌦ of definition of the shooting function, we may assume that S is regular. Then, 0 is a regular value and each connected component of S 1 ({0}) (there are finitely many if ⌦ is bounded) is a one-dimensional embedded submanifold di↵eomorphic either to S 1 or to R. As each c = (⇠, ) 2 ⌦ is a regular point of S, it is possible to define ! T (c) as the unit vector in the kernel of S 0 (c). Such a vector is unique provided a choice of orientation is made according to det( t S 0 (c), ! T (c)) > 0 (see [2]). Starting from c 0 2 S 1 ({0}), the corresponding connected component or branch is parameterized by the di↵erential equation with respect to curvilinear abscissa dc ds At a point c(s) = (t f (s), p 0 (s), (s)), second order su cient conditions are checked using a simple numerical rank evaluation [7]. Indeed, at time t 2 (0, t f (s)], Jacobi fields span the image of the di↵erential with respect to (t, p 0 ) of (t, p 0 , ) 7 ! x(t, x 0 , p 0 , ).
So verifying that this mapping has a full rank derivative at all (t, p 0 (s)) for t 2 (0, t f (s)] ensures that we have a C 0 -local minimizer at t f (s). This eventually implies that @S/@⇠(c(s)) has full rank, so the path can locally be parameterized not only by s but also by (no turning point). The path computed thanks to [7] which, in contrast to standard di↵erential homotopy codes, integrates without correction step equation (2), is portrayed in Fig. 2. The strength of the di↵erential approach, when applicable, is that steps are automatically computed and tuned.
2. Homotopy on ". Having obtained three-body minimum time trajectories using an homotopy on µ, our second goal is to iterate these computations for smaller values of the control bound, 3 ". We are faced with several issues. Paths may have turning points associated with rank drops in the partial derivative of the shooting function that prevent the verification of second order su cient conditions (see the end of the section for an adaptation of such conditions to the case of terminal submanifold-a circular orbit around the Moon is targeted, here). See for instance Fig. 3-4. Moreover, to ensure global optimality, for a given " one has to compare costs on the di↵erent connected components of the set of zeros, and it turns that one has to jump from one component to another. While a given branch is followed using di↵erential homotopy (integration of 2), switching branches is obtained by discrete homotopy. Combining both approaches allow to obtain extremals whose local optimality is checked by second order conditions and whose global optimality is tested by comparing various branches. On both subgraphs, " is expressed in Newtons. In the Earth-Moon system, a 1 Newton thrust on a 1500 Kilogram spacecraft (typical mass for medium or low-thrust transfers) is equivalent to " ' 2.440e 1.
The same reference to Newtons is made in subsequent figures. On the left, the abscissa is the curvilinear abscissa on the path and two turning points are observed (see also Fig. 4). The corresponding cusp points on " ! t f (") are observed on the right. Optimality is necessarily lost past the self-intersection point near " ' 9.675e 1 (see also Fig. 9).
Before proposing a preliminary discussion of the topology of solutions observed on di↵erent branches (see Fig. 5), we conclude the section by recalling that, in the  case of a terminal condition defining a submanifold (here a circular orbit instead of a single point such as a Lagrangian point as in section 2), the conjugate point notion can be adapted as follows. When the target is a proper submanifold X f , a time t c 2 (0, t f ] is said to be focal (see [11] for a discussion in the Riemannian case, and [4] in optimal control) if there is a non-trivial Jacobi field z along z which is vertical at t = 0 and such that We are thus again able to check local optimality of transfers from circular to circular orbits, such as the one in Fig. 6. 3. Homological classification. Homotopy on " in the previous section suggests that to di↵erent connected components of the path are associated solutions of different homological nature: Either winding around the second primary changes from positive to negative, or the sign remain constant but the winding number changes. A coarse classification can be made using homology of (closed) curves in Q µ . (The classification is coarse as it takes into account neither theq coordinate, nor the adjoint p). First restricting to extremals with boundary conditions on q 2 = 0 (the line defined by the two punctures in Q µ ), one can associate to any such extremal a closed curve by patching to it the curve symmetric with respect to q 2 = 0 (and orientation reversed), see Fig. 7. The homology of this curve turns to be su cient to classify the di↵erent types of solutions observed on disjoint components of the homotopy path. Note nevertheless that in Fig. 5, the first two curves (upper part) belong to the same branch, have the same winding number around the first singularity but not around the second one (though the winding sign is unchanged). For arbitrary boundary conditions, the same symmetric closure has to be performed on the truncation of the q-projection of the extremal from its first intersection with q 2 = 0 until the last one. As a case study, the path between " = 2.441e 1 and " = 2.196e 1 is discussed in Fig. 8. The optimality status of each homology class is analyzed in  around the second primary, fixed frame, " = 2.220e 1. Controlled trajectory in blue, free trajectory (null control after the target is reached) green. The capture into a quasi-periodic motion around the second primary is clearly observed. Fig. 9. It is a matter of future investigation to integrate homological constraints [10] to cut out sub-optimal components from the path.  respectively. The two extremals belong to the same branch, labeled I (see Fig. 3) and have the same homology (positively winding up around both primaries). Middle pictures are transfers for " = 2.440e 1 (left) and " = 2.197e 1 (right), respectively. Extremals belong to the same branch, labeled II, and have the same homology (now negatively winding up around the second primary). Bottom pictures are transfers for " = 2.440e 1 (left) and " = 2.196e 1 (right), respectively. Extremals belong to the same branch, labeled III, and have the same homology (negatively winding up around the second primary with one more revolution around the first one that those on branches I and II).  Figure 9. Homotopy on ", Earth-Moon case (µ ' 1.21e 2), minimum time transfer towards a circular orbit around the second primary. Each curve represents the function " 7 ! t f (") evaluated along the three disjoint branches I (blue), II (red) and III (black) of the homotopy (see Fig. 8). Slightly below " ' 2.367e 1 (that is below 0.97 Newtons on the graph), branch II provides a better criterion than I, but both have already become worse than branch III. The infimum of the three curves provides an upper bound for the value function " 7 ! t f (").