Indirect methods for fuel-minimal rendezvous with a large population of temporarily captured orbiters

A main objective of this work is to assess the feasibility of space missions to a new population of near Earth asteroids which temporarily orbit Earth, called temporarily captured orbiter. We design rendezvous missions to a large random sample from a database of over 16,000 simulated temporarily captured orbiters using an indirect method based on the maximum principle. The main contribution of this paper is the development of techniques to overcome the difficulty in initializing the algorithm with the construction of the so-called cloud of extremals.


1.
Introduction. In the continued pursuit of knowledge regarding the origins of our solar system, scientists have been looking toward the minor bodies -asteroids and comets. These small yet abundant objects are distributed throughout the solar system and are pristine representatives of the history of our planet and the universe. Moreover, it has been recently suggested that missions to near Earth asteroids could be the key stepping stones for interplanetary human flight. The opportunity for research is huge, and potentially affordable, when considering those near Earth asteroids that come close to home.
Based on the asteroids humans have already visited, the target asteroids have been characteristically large, far away from Earth, and on elliptic, heliocentric orbits [13]. The so-called minimoons are different: they are temporarily captured near Earth asteroids and present intuitively many advantages for a space mission. There is, however, a trade-off since their capture time is on average only nine months and their orbits display a challenging complexity for an orbital transfer. The primary objective of this paper is to answer whether or not temporarily captured orbiters are suitable targets for spacecraft missions, and in so doing provide a strong push to develop detection programs for this new population of asteroids. More specifically, we assess the rendezvous suitability of a large sample of simulated temporarily captured orbiters in order to determine and characterize those that are better suited from those that are worse. We use indirect methods of optimal control theory and develop a foliation technique in order to successfully initialize the algorithms. We consider the Circular Restricted Four-Body Problem (CR4BP) with the Sun as a perturbation of the Circular Restricted Three-Body Problem (CR3BP). The main difficulty in utilizing indirect methods to numerically solve optimal control problems is providing the algorithm with a close enough initial guess so that the algorithm converges. To overcome this difficulty, a technique is developed to create a database of known extremals, which are then used with a continuation method to solve the rendezvous problem for large population of asteroids. We provide analysis of the results in terms of the application and discussion of the effectiveness of our techniques.
Our methods are applied to a large database of simulated asteroids and the results show that fuel efficient missions are indeed theoretically possible, while simultaneously taking into account the mission time constraints. More precisely, we demonstrate on a random sample of 3000 temporarily captured orbiters from the database produced in [11] that we can design a rendezvous transfer for 100% of them with a delta-v under 680 m/s for half of those and a duration under three months for most of the transfers. Note that this is a preliminary result and fuel consumption can be improved on a single-case basis once a target orbiter is selected for further study as it is done on a couple of temporarily captured orbiters at the end of this paper. In [7], the authors focus on asteroid 2006RH 120 and provide more detailed results for this specific target with realistic rendezvous and return transfers. In a broader context, the results suggest that the algorithm initialization problem may feasibly be overcome for a general optimal control problem following the techniques developed within this paper.
2. Temporarily captured orbiter. By definition a natural Earth satellite (NES) is a celestial body that orbits the Earth. A NES is defined as a temporarily captured orbiter (TCO) of the Earth if it simultaneously satisfies: 1. the geocentric Keplerian energy E geo < 0, 2. the geocentric distance is less than three Hill radii R H,⊕ , where 3R H,⊕ ≈ 0.03 AU for the Earth, 3. the object makes at least one full revolution around the Earth in the corotating frame while satisfying the first two criteria.
The number of revolutions is measured by recording the longitudinal angle traversed during capture in the co-rotating ecliptic coordinate system where the line connecting the Earth and the Sun is fixed, and forms a line of reference. This number is negative for retrograde orbits.
A few different names have been used to describe these objects, though there is yet to be a single accepted choice by the scientific community: natural Earth satellites (NES), minimoons, and finally, temporarily captured orbiter. As a convention for this paper, we choose to use temporarily captured orbiter (TCO), but the three are interchangeable.
Despite a large body of work on satellite capture by Saturn and Jupiter, there has been much less published about the natural satellites of the Earth other than the Moon. In [11] by Granvik et al., a study and characterization of the population statistics of the Earth's natural satellites is calculated for the first time. By generating, pruning, and integrating an enormous random sample of "test-particles" from the near Earth object (NEO) population, the authors of [11] calculate the size-frequency distribution (among other statistics) of the TCO population. To do so, an original random sample of 9,346,396,100 test-particles was generated over the 19-year Metonic Cycle (the time it takes for the Earth-Moon-Sun to approximately return to the same solar system configuration), so that after pruning there were 10,000,000 test-particles remaining for integration. The test particles were originally generated based on existing population statistics for near Earth objects, and the pruning was done according to the fundamental assumption that test particles which can get captured are those with heliocentric orbital elements similar to Earth's orbit. Using a high-precision n-body model, each test-particle was integrated over time, and classified as a TCO if for some portion of the integration the definition given above was satisfied. Of the 10 million integrated test-particles, the authors of [11] found that over 16,000 became temporarily captured orbiter. It is from this database of integrated temporarily captured orbiter that we select our targets for rendezvous missions.
As a convention for the rest of the paper we will refer to the 16,923 simulated TCOs as the TCO database or simply the database, and when we say TCO we are referring to an element in that database. To each TCO in the database, we assign an index number 1-16,923 so that we can refer to a specific simulated TCO by its index number. Each TCO in the database consists of a discretized time-history of the position and velocity of the asteroid during its capture. The Sun, Earth, and Moon ephemerides are also known at each time step, for every TCO. Figure 1 shows six TCO orbits from the database: TCO #6890, #17, #16, #110, #23, and #39. The distance units are normalized so that 1 Lunar Distance (LD) = 384,400 kilometers, an approximate average distance between the Earth and Moon. The diversity of the TCO's trajectories is striking. Indeed, these objects are very different than the types of asteroids for which prior missions have been developed. Instead of large asteroids, on far away, elliptic, heliocentric orbits, TCO are tiny asteroids and (during their capture) are on nearby, rapidly evolving, geocentric orbits. The distribution of capture times less than three years for the database of 16,923 temporarily captured orbiter is shown in Figure 2. A main objective is to characterize which TCOs are the best suited, since indeed it is expected that depending on their orbit characteristics some TCO are good targets while others would require too much fuel to be able to achieve rendezvous before the TCO escapes capture. Additionally, given a specific TCO we would like to determine the best feasible rendezvous location for that precise orbiter. The TCOs capture times vary from a few weeks to several years.
Statistically, the results from [11] show that at any given time there is at least one 1-meter-diameter TCO orbiting the Earth, although to date asteroid RH 120 has been the only documented TCO. The lack of more numerous discovery is likely due to the small size of the objects, which may go unnoticed or be dismissed as manmade debris, and due to the fact that the objects are only captured temporarily around Earth. In [3] theoretic work has been done to investigate the discoverability of TCOs by present or near-term ground-based and space-based facilities.
In addition to their statistical abundance and having been previously unexplored, TCOs have several characteristics which make them intuitively appealing targets for space missions. First, their closeness to Earth and the large amount of documentation already available regarding space travel in the Earth-Moon vicinity. Second, an orbiting object allows for more time for detection, planning, and execution of a space mission than an object which just flies by (i.e., does not complete a revolution of the Earth). Note, the mean capture time for the simulated TCOs is 286 days, with the mean number of revolutions 2.88 [11]. In particular, it is reasonable to think that TCOs are excellent asteroids to be redirected on to a stable orbit in the Earth-Moon system. This has been demonstrated by the work in [18] on RH 120 . A redirected TCO would then be an ideal target for a manned mission. Third, the small size of the TCOs allows us to envision returning, not only a sample, but the entire asteroid to Earth which would provide a breakthrough in the geophysical study of asteroid, see [7] for examples of return mission on RH 120 . Finally, the chaotic trajectories (non-planar, non-elliptic) is actually appealing in the sense that spacecraft maneuvering capabilities and mission planning strategies can truly be put to the test. Figure 3 shows the distribution of periapsis distances with respect to the Earth, Moon, L 1 , and L 2 . Recall that the definition of temporarily captured requires that the object be within 0.03 AU = 11.67 LD of the Earth, and in Figure 3 we see that all the temporarily captured orbiter come well-within this boundary. We mostly present this to observe that many TCOs come close to the EM Lagrangian points L 1 and L 2 .
3. Mission scenario. In order to assess the delta-v requirements necessary for a TCO rendezvous mission, fuel-minimal rendezvous transfers are computed for a chemical propulsion spacecraft with maximum thrust T max = 22 Newtons, specific impulse I sp = 230 seconds, and initial mass m 0 = 350 kilograms. Our choice of propulsion is motivated by the fact that electric-propulsion time-minimal transfers have very low transfer times, but their corresponding delta-v requirements were too high (over 2,000 m/s) [8,9,13], and the chosen specifics mimic the ones from GRAIL mission's spacecrafts. The spacecraft starts on a Halo orbit, with period 14.8 days and z-excursion of 5,000 kilometers, around the Earth-Moon L 2 pointwhich we denote by O H . Figure 3 shows that a majority of the TCOs come within 1 LD of L 2 , and moreover, the Artemis mission [15,19] has demonstrated that such an orbit can be maintained over long periods for very low delta-v. The influence of the Sun is taken into account to perturb the CR3BP, which gives the circular restricted four-body problem (CR4BP). The mass variation of the spacecraft is also included since it is a fuel minimal problem. Although we do not seek to minimize time, we are still constrained by the capture time of the TCOs. Therefore, we verify that the computed transfers take place entirely while their associated TCOs are in capture, meaning our optimization problem has a dual cost: the minimization of fuel while simultaneously keeping the transfer time low. 4. Circular restricted four-body problem. The circular restricted four-body problem is a perturbation of the CR3BP, see [12,16,17,20] for complete details. The idea is to begin with the same assumptions as for the CR3BP and then assume that a fourth body -the Sun -with mass µ s (in normalized CR3BP units) influences the motion of the spacecraft but not the Earth and Moon. The Sun is assumed to move uniformly on a circular orbit around the origin with angular momentum ω s , so that at time t its normalized position is (r s cos(θ(t)), r s sin(θ(t))), where r s is a constant representing the radial distance from the origin and θ(t) is the true anomaly of the Sun at time t with respect to the x-axis of the CR3BP. The specific parameter values for the Earth, Moon, Sun system are given in Table 1  Under these assumptions the uncontrolled equations of motion from the CR3BP are augmented to give the equations of motion for the CR4BP: where the potential energy V 4 = V + V s is the CR3BP potential energy function V : with the distances from the spacecraft to the Earth and Moon, respectively, plus the additional terms V s from the Sun's influence: and is the distance from the spacecraft to the Sun. Recall that the Earth and Moon are fixed points at (−µ, 0, 0) and (1 − µ, 0, 0) respectively; however, the system is nonautonomous since the Sun's position changes with time according to where θ 0 represents the initial Sun's true anomaly. Note that the system can be rewritten autonomously by considering θ(t) as a state variable which evolves according toθ = ω s . This system has no equilibrium points and no first integrals of motion. We assume the spacecraft starts on a Halo orbit around L 2 despite the fact that it is a periodic solution of the circular restricted three body problem, not four. Indeed Halo orbits have been utilized in actual missions including the Artemis mission [15,19], and since our computations have the spacecraft immediately departing from the Halo orbit, this is a reasonable approximation.
The controlled CR4BP is formulated by adding control terms u = (u 1 , u 2 , u 3 ) T , u ≤ 1, to the x, y, and z accelerations, giving where T max is the constant representing the maximum thrust of the spacecraft. The mass depletion is included in the model, so that m(t) decreases as the thrusters are fired according to the last equation of (5). The constant β = (g 0 I sp ) −1 is related to the Earth's gravity at sea level g 0 ≈ 9.81 m/s 2 and the spacecraft's specific impulse I sp which measures how quickly a spacecraft expends fuel.
The mathematical formulation of our problem is to compute solutions of the system (5) which minimize the expended fuel, which is equivalent to minimizing the delta-v. Delta-v, denoted ∆v, is a standard measure of the acceleration generated by the spacecraft in order to perform maneuvers, and is defined as The ∆v budgets for actual space missions can vary greatly depending on the specific mission objectives. For our purposes, we consider ∆v values less than 1,000 m/s as good, and below 500 m/s as great. From Equation (6) it's obvious that minimizing ∆v is also equivalent to maximizing the final mass m(t f ). As a convention, we will say "fuel-minimization" and "∆v-minimization" interchangeably. The relationship between ∆v and m(t) allows our cost functional to be stated as: where the control domain U = B R 3 (0, 1) is the unit ball in R 3 and the set of The initial time is assumed to be t 0 = 0, and the transfer time t f is free for the general problem, though is later considered fixed.

5.
Problem statements and the maximum principle. We present the general fuel-minimization problem and the corresponding application of the maximum principle. In order to overcome the difficulties inherent in initializing the indirect algorithms, we also present a restricted version of the problem which is computationally less sensitive to the initialization, and is a vital step toward computing solutions to the general problem.

General problem.
Let the control variables be given by u = (u 1 , u 2 , u 3 ) T and the state variables by X = (x, y, z,ẋ,ẏ,ż, m) T . It will sometimes be convenient to refer to the position s = (x, y, z) T and velocity v = (ẋ,ẏ,ż) T separately, or concatenated as q = (x, y, z,ẋ,ẏ,ż) T . With this notation the equations of motion can be written as the first order control system: where Let q rdvz be the position and velocity of a point along a chosen TCO's orbit, and let θ rdvz be the corresponding true anomaly of the Sun when the TCO is at q rdvz . Computing ∆v-minimal transfer from O H to q rdvz amounts to solving the optimal control problem Note the difference in expressing the initial and final position and velocity constraints: the TCO data is discretized and cannot be expressed as a manifold, hence a specific rendezvous location is chosen, namely q(t f ) = q rdvz . On the other hand, the Halo orbit O H is a continuous curve and so the departure point can be any point along that curve, i.e. q(0) ∈ O H .
Applying the Pontryagin maximum principle [2,14] to (9), we have the existence of a non-positive constant p 0 and an absolutely continuous function P : [t 0 , t f ] → R 7 , with (p 0 , P) = (0, 0), which is a solution of th pseudo-Hamiltonian equations (10). We denote the so-called costate variables by P = (p x , p y , p z , pẋ, pẏ, pż, p m ) T , the components of P which correspond to the state variables of X. We assume we are in the so-called normal case with p 0 = 0 so we can normalize p 0 = −1. We know that every solution X(t) to (9) is necessarily a projection of an extremal curve (X(t), P(t)) solution of the systeṁ where the pseudo-Hamiltonian function H is defined as with H i (P, X) = P, F i (X) , i = 0, ..., 4. Note in particular that H 4 (P, X) = p m , and also that The term u, p v is the only part of H which depends on the direction of u, and so by the maximization condition and from the Cauchy-Schwarz inequality we can determine that where α now denotes the yet to be determined magnitude of the control. Substituting and rewriting, we have the so called real Hamiltonian H r where Letting Φ = −1 + T max pv m − βp m denote the so-called switching function, the maximization condition further implies that The case when Φ(t) = 0 corresponds to a so-called bang arc where the control norm is either maximized or minimized. Isolated times where Φ(t) = 0 are the so-called switching times where the control switches between bang arcs. When Φ(t) ≡ 0 on a non-trivial time interval t ∈ (t 1 , t 2 ) we have a so-called singular arc, where more work must be done to determine the norm of the control. It is not in the scope of this work to analyze the singular arcs since we are focused primarily on the application. Moreover, in the next section restrictions will be made on the control structure which eliminate the possibility of singular arcs.
The transversality conditions of the maximum principle also give constraints on some components of the costate vector at the initial time and final times. Since m(t f ) is free, the maximum principle implies where p q are the components of P corresponding to the q components of X. Furthermore, as in [14], if the final time is free, the maximum principle implies that the real Hamiltonian is identically zero on [0, t f ], so we also have the constraint Denoting the pair Z = (X, P) T , the computation of solutions to (9) can be carried out by solving the shooting equation associated with the function Although m(0) = m 0 and θ(t f ) = θ rdvz are also boundary conditions, they can always be satisfied by construction and therefore are omitted from the shooting function. Note, equation (4) is equivalently rewritten θ(t) = θ rdvz + ω s (t − t f ).
A Newton-like algorithm, provided by the software Hampath [4,5], attempts to compute a zero of the shooting function. If successful, the corresponding triple (t f , X(t), P(t)) is an extremal and a candidate local minima. The difficulty of course is providing the solver with a good enough initial guess so that the algorithm converges.

5.2.
Restricted problem. We want to solve optimal control problems of the form (9) for a large variety of (q rdvz , θ rdvz ) from multiple points along over 16,000 TCOs, so there is a major need for computational efficiency. As a first approach we instead compute solutions to a restricted version of the problem. The main restriction we make is to reduce the set of admissible controls U by fixing the control structure. Additionally, to reduce the number of constraints in the shooting function we do not verify the transversality conditions associated with the transfer time (H r (P(t f ), X(t f )) = 0) or with the departure point on from the Halo orbit (p q ⊥ T q(0) O H ). Some explanation and justification for each choice is given below, followed by the statement of the restricted problem. We fix the control structure by assuming that the spacecraft performs three or less maximum-thrust boosts, so that the restricted control normᾱ(t) is given bȳ The times t i , i = 1, . . . , 4 are now called the switching times, and we will omit writing "i = 1, . . . , 4" whenever it is understood. The possible equalities allow for only one or two boosts instead of three. This choice is motivated by the physical application since utilizing less maneuvers is less risky in terms of the success of the mission. There is also some intuitive theoretical motivation: the first boost allows us to depart the initial orbit onto the orbit's unstable invariant manifold, the second boost puts us on course for rendezvous, and the third boost is used at arrival to actually achieve rendezvous. Moreover, tests showed that adding more boosts increases the computation time significantly without lowering the resulting delta-v values very much. Figure 4 illustrates an example of the control norm for the fixed control structure. In [6], the authors overcome the general fuel-minimal initialization for a specific mission scenario using a homotopic method which related the problem of minimizing the L 1 -norm (min t f 0 u dt) to the problem of minimizing the L 2 -norm (min t f 0 u 2 dt), rather than fixing the structure of the control. We explored this approach applied to the TCO rendezvous scenario, however determined it was too computationally demanding and therefore ill-suited for the large scale analysis we are aiming to conduct. It's worth mentioning however that the few solutions we did compute using this method consisted of either two or three boosts, which further supports the choice of restriction we make here.
Additionally, we now consider the transfer time t f to be fixed and therefore do not verify the transversality condition H r (X(t f ), P(t f )). As a consequence, the transfer times for computed transfers are not necessarily optimal. Similarly, the departure point from O H is now considered fixed, and denoted q dprt . We do not verify the transversality condition p q (0) ⊥ T q(0) O H , so the O H departure point q dprt for computed transfers is not necessarily optimal. We remark that with the fixed control structure, ∆v can be calculated from the switching times t i and final time t f . Recall ∆v = dt. Since u is now assumed to always be either 0 or 1, let τ be the total thrusting time (i.e. τ = t 1 +(t 3 − t 2 ) + (t f − t 4 )), then using equation (6) we obtain ∆v = − 1 β (ln(m(τ )) − ln(m 0 )) . Solving for τ , we get τ = m0 βTmax (1 − exp(−β∆v)). Therefore, for a fixed ∆v value, we know exactly the total thrusting time τ . The relationship between ∆v and τ is graphed in Figure 5 for the specific parameter values we use.
Note that these three restrictions only serve to make the restricted problem a special case of the general problem, and therefore solutions to the general problem will perform at least as well as solutions to the restricted problem. That is, the delta-v values we compute for the restricted problem are upperbounds for those of the general problem.
Let U = {u : R → B R 3 (0, 1) : u measurable, u(t) =ᾱ(t) for a sequence of times 0 ≤ t 1 ≤ t 2 ≤ t 3 ≤ t 4 ≤ t f }. Then the restricted optimal control problem is written We will attempt to solve ( ) for a variety of transfer times, departure points, and rendezvous points. As such, the parameters that define a specific instance of the problem ( ) are t f , q dprt , q rdvz and θ rdvz (m 0 is always assumed to be 350). For convenience we will refer to a tuple D = (t f , q dprt , q rdvz , θ rdvz ) as the problem data and the associated instance of the optimal control problem as ( ) D .
The derivation and expression of the real Hamiltonian from the maximum principle is the same as in the general case, except that the control norm is now fixed, yielding: The mass variation is then given bẏ That is, the mass is only depleting when the engines are firing (ᾱ = 1). Substituting (16) into H yields the expression of the real Hamiltonian H r : and we have the flow of the real Hamiltonian: By continuity of the Hamiltonian, we must have that the so-called switching conditions Φ(t i ) = 0 must hold at each switching time, as is shown in [10]. Again denoting the pair Z = (X, P) T , the computation of solutions to ( ) can be carried out by solving the shooting equation associated with the function S : R 11 → R 11 given by: As before, the software Hampath [4] attempts to compute a zero of the shooting function. If successful, the corresponding tuple (t f , t i , X(t), P(t)) is an extremal and a candidate local minima for ( ). In the next section, we address the major difficulty of algorithm initialization.
6. Cloud construction. The problem ( ) D , for D = (t f , q dprt , q rdvz , θ rdvz ), asks: what is the fuel-optimal way to get from q dprt to q rdvz in exactly t f units of time, and so that the Sun's true anomaly at arrival is θ rdvz ? We seek to answer this question for thousands of choices of (q rdvz , θ rdvz ) on TCO orbits, and a variety of choices for q dprt and t f . Figures 6 and 7 give depictions of (q rdvz , θ rdvz ). Solving ( ) D for any choice of problem data D amounts to properly initializing the (very sensitive) shooting algorithm. The so-called cloud we construct here provides an intuitive and effective way to initialize the shooting algorithm so that it converges to a solution for a large variety of problem data D.
Let t f be a fixed transfer time and let q dprt be a departure position and velocity on the halo orbit. If q rdvz represents a point along a TCO's orbit with associated Sun true anomaly θ rdvz , then we will call q rdvz a rendezvous point, and call an extremal for problem ( ) D , with D = (t f , q dprt , q rdvz , θ rdvz ), a rendezvous extremal.
If (t f , t i , X(0), P(0)) is an extremal for a problem ( ) D with D = (t f , q dprt , q 1 , θ 1 ), but q 1 is not on a TCO orbit, we will call q 1 a cloud point and (t f , t i , X(0), P(0)) a cloud extremal. These are locally optimal trajectories that go somewhere in space, not necessarily to a TCO. The idea is to compute rendezvous extremals by performing a continuation method from a nearby cloud extremal. We describe next how we can construct cloud extremals. We will refer to the collection of constructed cloud extremals as, simply, the cloud. With enough cloud extremals, we expect (and will show that indeed) rendezvous points will be very near cloud points, ensuring that the continuation method will have a high success rate. Thus, the first task is to construct a large pool of cloud extremals.
The construction of a cloud extremal relies on the fact that a random choice of (t i , P(0)) and of (t f , q(0), θ(0)) completely determines a trajectory (X(t), P(t)), thanks to (18). In particular, q(t f ) and θ(t f ) are determined. Note that (X(t), P(t)) is a trajectory from q(0) to q(t f ), but it is not necessarily an extremal of problem ( ) D , D = (t f , q(0), q(t f ), θ(t f )) since it might fail to satisfy the transversality  condition p m (t f ) = 0 and switching conditions φ(t i ) = 0. By attempting to solve the shooting function S D (·) = 0 for D = (t f , q dprt , q(t f ), θ(t f )) with initial guess (t i , P(0)), the first six components of S D (·) are by construction zero, allowing the algorithm to converge more easily. When we have convergence pf the shooting function, we obtain the so-called cloud extremals. The cloud points corresponding to those cloud extremals are compared to the rendezvous points as described below, and pairs which are close to each other are chosen. For each chosen pair, a continuation method is attempted from the cloud point to the rendezvous point. If successful, we've found a rendezvous extremal. We can that way construct a large database of cloud extremals. The process is summarized: 1. Choose values for (t f , t i , q dprt , θ 0 , P(0)). 2. Integrate according to (18) to determine q(t f ) and θ(t f ). 3. Define the problem ( ) D with probldata D = (t f , q dprt , q(t f ), θ(t f )). 4. Attempt to solve the shooting function (19) associated with ( ) D using the chosen values of (t i , P(0)) as the initial guess. 5. If successful, we have computed a cloud extremal, which is added as an element to the cloud. By repeating this process for varying choices of (t f , t i , q dprt , θ 0 , P(0)), we populate the desired large database of cloud extremals. Although it is unlikely that any of these extremals achieve rendezvous with any of the TCOs, by developing a dense enough cloud we expect that at least some of the extremals end near some of the TCO orbits, and that from there continuation techniques will prove successful. A two-dimensional projection of the complete trajectory is shown as the solid path in Figure 8, with q dprt marked as a triangle and q(t f ) marked as a square. The second boost after 10 days at q(t 2 ) is marked as a circle. By construction, all three boosts are one hour, so the total thrusting time τ = 3 hours which corresponds to ∆v = 807.6 m/s. 3. Use the choices for t f and q dprt with the integrated values for q(t f ) and θ(t f ) to define the problem data D = (t f , q dprt , q(t f ), θ(t f )) and corresponding optimal control problem ( ) D . 4. Attempt to solve ( ) D using (t i , P(0)) as the initial guess for the shooting method. Again, note that by construction now the first six component of S(t i , P(0)) are zero, and so our initial guess is in some sense close to a solution. The algorithm iteratively adjusts the values of (t i , P(0)) to satisfy the transversality and switching conditions. For this example, the algorithm converges to a solution, though convergence is in general not guaranteed. Our initial choices for the switching times t i and costate variables P(0) are compared to the the solution values (t * i , P(0) * ) in Tables 2 and 3 respectively. 5. The computed solution (t * f , t * i , X * (t), P * (t)) is an extremal to ( ) D . A twodimensional projection of the cloud extremal trajectory is shown as the dashed path in Figure 8, along with the initial guess trajectory in solid line. Note that both transfers take the same amount of time: t f = 20 days. Also notice     Table 3. Comparison of chosen values P(0) versus solution values P(0) * . that the two trajectories depart from and arrive at the same locations, but otherwise are different. The cloud extremal has a total thrusting time of τ * = 2.58 hours and ∆v = 674.8 m/s. 6.2. Populating the cloud. We repeat the above process 1,000,000 times, generating random initial values according to the following: • 0 ≤ t f ≤ 182.5 days. Recall that our transfer times are constrained by the capture time of the target minimoon. The average minimoon capture time is nine months, so as a first approach we'd like to have extremals taking between 0-6 months. • 0 ≤ τ ≤ 1.98 hours, which implies 0 ≤ ∆v ≤ 500 m/s. This should give us extremals with low delta-v values. As in the example, we always assume  three equal boosts of length τ /3 for the initial guess, with one boost at the start (t 1 = τ /3) and one boost at the end (t 4 = t f − τ /3). The middle boost from t 2 to t 2 + τ /3 is randomaly placed between boost one and boost three. Figure 9 illustrates the constructed initial guess control norm. Note, as in the example, the solver adjusts the values of t i so the resulting cloud extremals will not necessarily have three equal length boosts. • q dprt ∈ O H . We discretize the Halo orbit O H into 100 equally spaced (in terms of time, 3.55 hour steps) points, and randomally select one ( Figure 10). • 0 ≤ θ 0 ≤ 2π. The initial Sun true anomaly is randomly selected.
• P q (0) ∈ [−0.01, 0.01] 6 , −0.001 ≤ p m (0) < 0. These choices are made based on exploration which seem to give a higher rate of convergence to cloud extremals. We do not claim that these are the best choices, and admit that some extremals may be missed by not considering initial guesses outside of this range. To that remark, we reiterate that this is a first approach to generate a large pool of extremals and that it is not our intent to analyze the costate space. Figure 11. Delta-v distribution for the 12,558 cloud elements On a standard personal laptop, the 1,000,000 cloud construction iterations took approximately two weeks to run. Of the 1,000,000 attempts, the algorithm converged 12,558, yielding 12,558 cloud extremals. The initialization parameters were taken from uniform random distributions, but we may expect that the algorithm converges more easily for certain parameter values. The distribution of delta-v, transfer times t f , and initial Sun anomalies θ 0 for the 12,558 cloud elements are shown in Figures 11, 12, and 13, respectively. We observe a slight bell shaped distribution, with small tails for delta-v less than 100 m/s and also between 400 and 500 m/s. For the initial Sun true anomalies we see a fairly uniform distribution. The distribution is not uniform for transfer times, suggesting that the algorithm was more likely to converge for transfer times between approximately 10 and 30 days.
Twelve-thousand out of one million is a small percentage, which attests to the sensitivity of the shooting method. Albeit small, we are satisfied using the 12,558 cloud extremals as our reference database to initialize continuation methods to compute rendezvous transfers. As justification, we provide the following statistics to quantify the density of the cloud and its relationship to the TCO orbits.
For each cloud extremal (t * f , t * i , X * (t), P * (t)) in the database, let q cloud := q(t f ) be the associated cloud point and let θ cloud := θ(t f ) be the associated final Sun true anomaly. Of the 12,558 cloud extremals, 11,128 of the cloud points position s cloud are within 5 LD of the origin. For each of these cloud points, we compute its position distance to the closest neighboring cloud point. The maximum of this value over all 11,128 cloud points is 0.64 LD, and the mean and median distances are 0.05 LD and 0.03 LD respectively. Focusing closer to the origin, 6,771 cloud points are within 2 LD of the origin. The same maximum-minimum distance for these points is 0.37 LD, and the mean and median distances are 0.03 LD and 0.02 LD respectively. The small mean and median values indicate that our cloud highly populate the space within 5 LD of the origin. In Figure 14 the set of all positions of q cloud points are plotted, and the density closer to the origin is observable, and further away points begin to spread out.  To some degree the set of cloud points represents a portion of the accessibility set from the halo orbit (shown in Figure 14) for less than 500 m/s delta-v. The Earth is in the center and the Moon on the right where the halo orbit is. Note for different zoom levels, the halo orbit, Earth and Moon are not always visible. From the closer of the zoomed images, we see the "ghosts" of the energy level sets from the three-body problem, even though we are in the four-body problem.
We now briefly compare the cloud points to the TCO orbits. For every TCO in the database, we compare every point s rdvz along its orbit to every cloud point position s cloud , and compute the minimum Euclidean distance between any pair. In this sense, the furthest TCO from the cloud is TCO 13313 at 3.056 LD away, and  the closest is TCO 12752 at 0.00088 LD away. The mean distance over all TCOs is 0.055 LD and the median distance is 0.023 LD. The histogram in Figure 15 shows the distribution of Euclidean distances from every TCO to the cloud. It is important to consider the velocity components as well when discussing rendezvous. So, for a selected rendezvous point q rdvz on a TCO's orbit and a selected cloud point q cloud , we now define the distance between them to be the non-dimensional Euclidean norm of the difference of the position-velocity vectors: q cloud − q rdvz . For each TCO, we compute the minimum distance between any q rdvz on the TCO orbit and any cloud point q cloud . In this sense, the furthest TCO from the cloud is TCO # 3187, which is 3.196 units at its closest. The closest TCO to the cloud is TCO # 10425 which is 0.0279 units away. The mean value over all TCOs is 0.613 units, and the median is 0.546 units. The histogram in Figure 16 shows the distribution of non-dimensional distances from every TCO to the cloud. 7. Continuation methodology and results. With a database of cloud extremals established, we now attempt continuations to finally compute rendezvous transfers to TCOs. We present three different rounds of calculations. For the first initial results, we attempt 50 continuations to each of 100 TCOs in order to get a preliminary assessment of our methodology and to guide the forthcoming large scale analysis.
7.1. Initial results. For each TCO, we select the 50 pairs (q rdvz , q cloud ) with the smallest distance q cloud − q rdvz . Note in particular that the same q cloud may be paired with more than one q rdvz in the top 50, and similarly, the same q rdvz may be paired with multiple q cloud points. For each continuation, the cloud extremal associated with q cloud is used to initialize the continuation method. The continuation method, as described in [1,13], iteratively solves problems ( ) D λ for λ = 0, ..., 1, where the problem data is given by When λ = 0 we have the problem ( ) D0 for which the cloud extremal associated with q cloud is a solution. If the algorithm converges for λ = 1, we have successfully computed a rendezvous transfer to the TCO, and an extremal to problem ( ) D1 . For the 5,000 total attempted continuations, the algorithm succeeded for 1,338. At least one transfer was found for 93 of the 100 TCOs. Looking at the best computed transfer for each of the 93 TCOs, we get the mean and median ∆v values  Figure 18. The highest computed best ∆v for a TCO was 3,086.2 m/s for TCO #14 with a transfer time of 59.2 days. The transfer is shown in Figure 19, and notice the complexity of the spacecraft's trajectory compare to the transfer for TCO # 9 and high z-variation of the TCO. Figure 20 gives the ∆v distribution for best transfers for the 93 TCO where transfers were found.
We analyze characteristics of the cloud and rendezvous points which guides our forthcoming calculations. Figure 21 compares the position components of the q rdvz points versus ∆v. From the z rdvz plot, it is clear that q rdvz with low absolute z-coordinates gave lower ∆v values. Similarly, Figure 22 compares the velocity components of the q rdvz points versus ∆v which shows that q rdvz with low absolutė z-coordinates gave lower ∆v values.
Computing the CR3BP energy value E at each q rdvz and comparing it to the ∆v values, we get Figure 23 and see a positive trend as the rendezvous energy moves away from the departure energy (dashed line). Similarly, we can compare the CR3BP energy values of the q cloud points used for our calculations to their corresponding rendezvous extremal ∆v values, as shown in Figure 24.
We suggest one way to quantify the interaction of the spacecraft with the Earth and Moon, by computing the winding numbers of the spacecraft's trajectory around the Earth and Moon. More specifically, we look at the planar projection of q(t) and compute the winding numbers (w e , w m ) around the Earth at (−µ, 0) and Moon at (1 − µ, 0) respectively. Note that in order to have closed curves we connect q(0) to q(t f ) with a straight line. Table 4 lists all the unique pairs of winding numbers computed for the 1,338 transfers, along with some ∆v statistics. There is a lot of information to parse, but a few things are striking. First, that over 1,000 of the 1,388 have w e = w m . The ∆v ranges seem scattered, but compare those pairs (w e , w m ) where |w e − w m | ≤ 1 to those where |w e − w m | > 1. We notice the mean ∆v's for the first group are lower. In fact, the mean ∆v for all transfers where |w e − w m | ≤ 1 is 989.3 m/s, whereas for |w e − w m | > 1 the mean value is 1911.5 m/s -nearly 1,000 m/s higher! Figure 18. Initial results: Rendezvous with TCO #9. The Moon's orbit is shown as the ellipse around the Earth (gray dot). The thin grey path is the orbit of the TCO, starting from its capture point (triangle) to its escape point (square). The circle on the TCO orbit marks where the TCO was when the craft departed, and the star is the point of rendezvous. The dark gray path is the spacecraft's trajectory. It's three boosts are marked as white dots (including the final rendezvous boost) Figure 19. Initial results: Rendezvous with minimoon #14. The legend is the same as for Figure 18.
Moreover, by computing the winding numbers for the cloud extremals and comparing them to their paired rendezvous extremals, we see that for 1,275 of the 1,338 continuations the cloud extremal and rendezvous extremal have the same winding  It is not surprisingly actually that the continuation method seems to preserve the structure of the original trajectory. This, combined with the remark above, suggests we may improve our search for low ∆v transfers by only using cloud extremals where |w e − w m | ≤ 1. Figure 25 plots the transfer times versus ∆v, and we see that a majority of the transfers took between 20 and 60 days, with a small patch of very low ∆v values around 50 days. Note that even for transfer times of 150 days transfers were found with ∆v < 400 m/s. Figure 26 shows that the norm q cloud − q rdvz was correlated to the ∆v of the resulting rendezvous extremal. Figure 27 shows also a slight preference for continuation pairs where |θ cloud − θ rdvz | is small -the lower clusters near 0 and ±π can actually be identified together since the CR4BP uncontrolled dynamics are π-periodic.     From the above analysis, the following characteristics of q rdvz , q cloud and of the pair seem to be the best indicators of a successful continuation to a low ∆v transfer: • |z rdvz | < 0.5 and |ż rdvz | < 0.5 • Cloud extremal winding numbers (w e , w m ) such that |w e − w m | ≤ 1 Figure 28. Large scale results: Distribution of ∆v for the 3,000 TCOs.
• q cloud − q rdvz < 0.7 and |θ cloud − θ rdvz | < 0.05 We now run continuations on a larger sample of TCOs. Based on the criteria above, for each TCO we only consider rendezvous points q rdvz where |z rdvz | < 0.5 and |ż rdvz | < 0.5. Note that 2,238 of the TCOs have no points that satisfy this criteria. We also restrict our cloud to only use those extremals where |w e −w m | ≤ 1, which leaves 10,687 of the original 12,558 cloud extremals. Finally, for each TCO we select the 10 pairs (q rdvz , q cloud ) with smallest distance q cloud − q rdvz while simultaneously requiring |θ cloud − θ rdvz | < 0.05. The Sun restriction is chosen because the continuation runs significantly faster for very close θ cloud and θ rdvz . Note, that there are enough cloud extremals so that the Sun restriction does not eliminate any more of the TCOs from consideration. We randomly select 3,000 TCOs for our continuations, based solely on limited computational time. Future work will include analysis of the entire database.
Continuations are attempted starting with the closest pair. If a solution is found for any of the pairs before the tenth, the algorithm moves on to the next TCO. This allows the computer to get through more TCOs in less time, while still providing first results for each TCO. Note if none of the top 10 continuations succeed the algorithm also moves on to the next TCO. 7.3. Main Results. For each of the 3,000 attempted TCOs a transfer was found. The distribution of delta-v values is given in Figure 28. The mean and median ∆v values were 724.5 and 680 m/s, respectively. Our computations demonstrate that rendezvous transfers can be obtained for 100% of the considered TCOs, and we expect no difficulty in expanding to the entire database of TCOs aside this task being computationally very demanding since the 3000 TCOs have bee chosen randomly. This suggests that efforts in detection of this targets for future asteroid missions should be initiated. The delta-v is within a reasonable range, especially based on the fact that this is an upper bound and that refinement can be conducted on case by case as done in the next section. Notice, that the lowest ∆v was 88.4 m/s for TCO #9,046, with a transfer time of 41.4 days. The transfer is shown in Figure 29. The greatest ∆v was 4,825.8 m/s for TCO #16,844, with a transfer time of 99.6 days. The transfer is shown in Figure 30. Recall that the results discussed above are solutions to the restricted problem where the control structure is fixed, and the transfer time and Halo departure are not optimized. To approximate how much better solutions to the unrestricted problem could do, in terms of ∆v, we refine the initial results for the original 93 TCOs we computed transfers to. For each, we take the extremal with the lowest ∆v and initialize continuations methods on t f , q dprt and q rdvz . The refinement algorithm is designed as follows: 1. Start with a solution to ( ) D where D = (t f , q i dprt , q j rdvz , θ j rdvz ) 2. Use a continuation to attempt to solve ( ) Dn for n = 1, ..., 6, where the D n are defined as follows: D 1 = ( t f − 3 hours, q i dprt , q j rdvz , θ j rdvz ) D 2 = ( t f + 3 hours, q i dprt , q j rdvz , θ j where q i±1 dprt represents the next/previous discretized halo departure point, and (q j±1 rdvz , θ j±1 rdvz ) is the next/previous minimoon rendezvous point (with exceptional cases when we at the first or last rendezvous point for the minimoon). 3. If any of the new delta-v values ∆v i are better than the original, choose the problem ( ) Di with lowest ∆v i , and repeat. If all ∆v i are greater than the original ∆v, the algorithm terminates. Although this algorithm does not verify any transversality conditions, it provides a crude estimations of locally optimum t f , q dprt and q rdvz .
The mean and median improvements for the 93 TCOs was 104.25 and 70.769 m/s ∆v, respectively. The smallest improvement was 2.22 m/s ∆v for minimoon #20, which originally had a ∆v of 459.21 m/s. The greatest improvement was 831.35 m/s ∆v for TCO #5, which originally had a ∆v of 2,480.2 m/s. Figure 31 shows the original and refined transfers to TCO #5. Figure 32 shows the distribution of ∆v improvements for the 93 TCOs. An example that showed good improvement without having very large initial ∆v was TCO #40 -which improved by 473.42 m/s from 1055.5 m/s to 582.04 m/s. The original and refined transfers are shown in Figure 33.
Finally, we take the best transfer found in the large scale results, for TCO #9,046 shown in Figure 29 and apply the refinement algorithm. The resulting transfer has a ∆v of 65.3 m/s and a transfer time of 41.4 days. Both the original and refined transfers are shown in Figure 34 for comparison.