S Brelsford

M Chyba
email: chyba@hawaii.edu

T Haberkorn

G Patterson

Rendezvous Missions to Temporarily-Captured Near Earth Asteroids

Keywords: Temporarily Captured Objects, Three-body Problem, Optimal control, Indirect Numerical Methods 2010 MSC: 49M05, 70F07, 93C10

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In this paper we compute delta-v minimal spacecraft transfers from an Earth-Moon (EM) L 2 Halo orbit to rendezvous with temporarily captured Earth orbiters (TCO). All computed transfers are for a 350kg spacecraft with 22N maximum thrust and 230s specific impulse, and we impose that the trajectory utilize three or less max thrust boosts.

The only documented Earth TCO, known as 2006 RH 120 (or from now on simply RH 120 for brevity) serves as an important test target for our calculations. In addition, rendezvous missions are computed to several simulated TCOs from [START_REF] Granvik | The population of natural Earth satellites[END_REF].

All transfers are designed using one of two different gravitational models: 1. First, transfers are computed using the EM circular restricted threebody problem (CR3BP) for the gravitational dynamics, which is justified since Earth TCOs are naturally evolving near the Earth and Moon. For these calculations, the transfer time and precise Halo departure point are treated as free variables to be optimized. 2. Second, the influence of the Sun is included in the dynamics using the Sun-perturbed CR3BP, sometimes known as the Earth-Moon-Sun circular restricted four-body problem (CR4BP). In this case an exact location of the spacecraft is assumed at the moment of asteroid capture, so that a choice of rendezvous location and transfer time determines exactly the departure point from the Halo orbit (dubbed the synchronization problem and discussed in detail within).

Many methodologies have been developed over the past decades to design optimal transfers in various scenarios. Due to the complexity of the TCO orbits and the nature of the mission, techniques based on analytical solutions such as in [START_REF] Kluever | Using Edelbaum's Method to Compute Low-Thrust Transfers with Earth-Shadow Eclipses[END_REF] for circular Earth orbits are not suitable and we use a numerical approach. A survey on numerical methods can be found in [START_REF] Conway | A Survey of Methods Available for the Numerical Optimization of Continuous Dynamic Systems[END_REF], and for reasons related to the specifics of our problem we choose to use a deterministic approach based on tools from geometric optimal control versus an heuristic method such as in [START_REF] Besette | Optimal Space Trajectory Design: A Heuristic-Based Approach[END_REF][START_REF] Pontani | Optimal Finite-Thrust Rendezvous Trajectories Found via Particle Swarm Algorithm[END_REF] and [START_REF] Vaquero | Leveraging Resonant-Orbit Manifolds to Design Transfers Between Libration-Point Orbits[END_REF][START_REF] Zhu | Trajectory Optimization of Multi-Asteroids Exploration with Low Thrust[END_REF]. All computations are carried out using classical indirect methods based on the Pontryagin Maximum Principle, combined with sophisticated numeric methods and software. The well-known sensitivity to initialization for this type of approach is addressed via a combination of direct methods and continuation techniques. Validation of our approach can be seen by comparing our work to [START_REF] Dunham | Using lunar swingbys and libration-point orbits to extend human exploration to interplanetary destinations[END_REF], in which the authors develop a low delta-v asteroid rendezvous mission that makes use of a Halo orbit around Earth-Moon L 2 . Their situation is different from ours in that they have carefully chosen a idealized asteroid for rendezvous. With a one-year transfer time, the delta-v value they realize is 432 m/s, which is comparable to the delta-v values presented here in a less-ideal scenario.

Analysis of all computed transfers strongly suggests that the CR3BP energy may play a role in predicting suitable departure and rendezvous points for low delta-v TCO rendezvous missions. Moreover, we present evidence that TCOs with more planar and more circular orbits tend to yield lower delta-v transfers.

Temporarily Captured Orbiters (TCO), RH 120

The motivation for our work is to study asteroid capture missions for a specific population of near Earth objects. The targets, Temporarily Captured Orbiters (TCO), are small asteroids that become temporarily-captured on geocentric orbits in the Earth-Moon system. They are characterized as satisfying the following constraints:

• the geocentric Keplerian energy E Earth < 0;

• the geocentric distance is less than three Earth's Hill radii (e.g., 3R H,⊕ ∼ 0.03 AU);

• it makes at least one full revolution around the Earth in the Earth-Sun co-rotating frame, while satisfying the first two constraints.

In regard to the design of a round trip mission, the main advantage of the TCOs lies in the fact that those objects have been naturally redirected to orbit the Earth. This contrasts with recently proposed scenarios to design, for instance, a robotic capture mission for a small near-Earth asteroid and redirect it to a stable orbit in the EM-system, to allow for astronaut visits and exploration (e.g. the Asteroid Redirect Mission (ARM)).

RH 120 is a few meter diameter near Earth asteroid, officially classified as a TCO. RH 120 was discovered by the Catalina Sky Survey on September 2006. Its orbit from June 1 st 2006 to July 31st 2007 can be seen on Figures 1 and2, generated using the Jet Propulsion Laboratory's HORIZONS database which gives ephemerides for solar-system bodies. The period June 2006 to July 2007 was chosen to include the portion of the orbit within the Earth's Hill sphere with a margin of about one month. We can also observe that RH 120 comes as close as 0.72 lunar distance (LD) from Earth-Moon barycenter. In [START_REF] Granvik | The population of natural Earth satellites[END_REF], the authors investigate a population statistic for TCOs. Their work is centered on the integration of the trajectories for 10 million testparticles in space, in order to classify which of those become temporarilycaptured by the Earth's gravitational field -over eighteen-thousand of which do so. This suggests that RH 120 is not the only TCO and that it is relevant to compute a rendezvous mission to this specific asteroid to get insight whether TCOs can be regarded as possible targets for transfers with small fuel consumption, and thus cost.

The choice of targets for our rendezvous mission sets us apart from the existing literature where transfers are typically designed between elliptic orbits in the Earth-Moon or other systems [START_REF] Caillau | Minimum fuel control of the planar circular restricted three-body problem[END_REF][START_REF] Mingotti | Optimal Low-Thrust Invariant Manifold Trajectories via Attainable Sets[END_REF] and [START_REF] Mingotti | Mixed Low-Thrust Invariant-Manifold Transfers to Distant Prograde Orbits Around Mars[END_REF], or to a Libration point [START_REF] Folta | Earth-Moon libration point orbit station keeping: Theory, modeling, and operations[END_REF][START_REF] Picot | Shooting and numerical continuation method for computing time-minimal and energy-minimal trajectories in the Earth-Moon system using low propulsion[END_REF] and [START_REF] Ozimek | Low-Thrust Transfers in the Earth-Moon System Including Applications to Libration Point Orbits[END_REF][START_REF] Vaquero | Leveraging Resonant-Orbit Manifolds to Design Transfers Between Libration-Point Orbits[END_REF]. Rendezvous missions to asteroids in the inner solar system can be found in [START_REF] Dunham | Using lunar swingbys and libration-point orbits to extend human exploration to interplanetary destinations[END_REF][START_REF] Kuninaka | Asteroid Rendezvous of HAYABUSA Explorer Using Microwave Discharge Ion Engines IEPC-2005-10[END_REF] but they concern asteroids on elliptic orbits which is not the case for us since TCOs are present complex orbits and therefore require a different methodology. Our assumption on the hibernating location for the spacecraft, a Halo orbit around the Earth-Moon unstable Libration points L 2 , is motivated in part from the successful Artemis mission [START_REF] Russell | The ARTEMIS Mission (Google eBook)[END_REF][START_REF] Sweetser | ARTEMIS Mission Design[END_REF] and in part from the constraint on the duration of the mission, mostly impacted by the time of detection of the asteroid. Indeed, the Artemis mission demonstrated low delta-v station keeping on Halo orbits around L 1 and L 2 .

Optimal Control Problem and Numerical Algorithm

Equations of motion

In this paper, we introduce two models, the circular restricted three-body problem (CR3BP) [START_REF] Koon | Dynamical Systems, the Three-Body Problem and Space Mission Design[END_REF] is first used to approximate the spacecraft dynamics and then we refine our calculations with a Sun-perturbed model (CR4BP). The first approximation is justified by the fact that a TCO can be assumed of negligible mass, and that the spacecraft evolving in the TCO's temporary capture space is therefore attracted mainly by two primary bodies, the Earth and the Moon.

The CR3BP model is well known, and we briefly recall some basic properties and notation which is useful for the remainder of the paper. We denote by (x(t), y(t), z(t)) the spatial position of the spacecraft at time t. In the rotating coordinates system, and under proper normalization (see Table 1), the primary planet identified here to the Earth, has mass m 1 = 1 -µ and is located at the point (-µ, 0, 0); while the second primary, identified to the Moon, has a mass of m 2 = µ and is located at (1 -µ, 0, 0). The distances of the spacecraft with respect to the two primaries are given by 1 = (x + µ) 2 + y 2 + z 2 and 2 = (x -1 + µ) 2 + y 2 + z 2 respectively. The potential and kinetic energies, respectively V and K, of the system are given by

V = x 2 + y 2 2 + 1 -µ 1 + µ 2 + µ(1 -µ) 2 , K = 1 2 ( ẋ2 + ẏ2 + ż2 ). (1)
We assume a propulsion system for the spacecraft is modeled by adding terms to the equations of motion depending on the thrust magnitude and some parameters related to the spacecraft design. The mass of the spacecraft is denoted by m and the craft's maximum thrust by T max . Under those assumptions, we have the following equations of motion:

ẍ -2 ẏ = ∂V ∂x + T max m u 1 , ÿ + 2 ẋ = ∂V ∂y + T max m u 2 , z = ∂V ∂z + T max m u 3 (2)
where u(•) = (u 1 (•), u 2 (•), u 3 (•)) is the control, and satisfies the constraint u = u 2 1 + u 2 2 + u 2 3 ≤ 1, and with T max /m normalized. A first integral of the free motion is given by the energy of the system E = K -V . We will later use this energy value to analyze the choice of the rendezvous point and the parking orbit for the spacecraft. It is well known that the uncontrolled motion of the dynamical system has five equilibrium points defined as the critical points of the potential V . Three of them L 1 , L 2 and L 3 are aligned with the Earth-Moon axis and have been shown to be unstable, while the two others are stable and are positioned to form equilateral triangles in the plane of orbit with the two primaries. Since our goal is to maximize the final mass we must include the differential equation governing the variation of the mass along the transfer: ṁ = -βT max u

where the parameter β, the thruster characteristic of our spacecraft, is given by β = 1 Ispg 0 (it is the inverse of the ejection velocity v e ), with I sp the specific impulse of the thruster and g 0 the acceleration of gravity at Earth sea level.

As the TCO's distance from the Earth-Moon L 2 equilibrium point can go as far a 12 lunar distance in our transfer computations, it becomes necessary to take into account the potential force of the Sun. To this end, we will also consider CR4BP, a Sun-perturbed Earth-Moon CR3BP, as in [START_REF] Mingotti | A Method to Design Sun-Perturbed Earth-to-Moon Low-Thrust Transfers with Ballistic Capture[END_REF]. In this model, we assume that the Sun follows a circular orbit around the Earth-Moon barycenter and shares the same orbital plane as the Earth and the Moon. Denoting by µ S the normalized mass of the Sun and r S the normalized constant distance from the Sun to the Earth-Moon CR3BP's origin, the potential energy of the spacecraft becomes V + V S , where V S is given by

V S (x, y, z, θ) = µ S r S - µ S 2 S (x cos θ + y sin θ). (4)
In (4), θ(t) = θ 0 + tω S is the angular position of the Sun at a given time t, with ω S the angular velocity of the Sun and θ 0 the angular position of the Sun at initial time 0. The distance from the spacecraft to the Sun is

S (t) = (x -r S cos θ(t)) 2 + (y -r S sin θ(t)) 2 + z 2 .
The numerical values are given in Table 1. The equations of motion of this Sun-perturbed model take the same form as (2) but with V + V S replacing V .

Numeric Methods

We first introduce the necessary conditions for optimality in the Earth-Moon CR3BP model and expand to the CR4BP.

Let q = (q s , q v ) T where q s = (x, y, z) T represents the position variables and q v = ( ẋ, ẏ, ż) T the velocity ones. From section 3.1, our dynamical system without the Sun perturbation is an affine control system of the form:

q = F 0 (q) + T max m 3 i=1 F i (q)u i ( 5
)
where the drift, in R 6 , is given by

F 0 (q) =       ẋ, ẏ, ż 2 ẏ + x -(1-µ)(x+µ) 3 1 -µ(x-1+µ) 3 2 -2 ẋ + y -(1-µ)y 3 1 -µy 3 2 -(1-µ)z 3 1 -µz 3 2       , (6)
and the control vector fields are F i (q) = e 3+i with e i forming the orthonormal basis of IR 6 . When considering the CR4BP, the main difference is that the Sun's position depends explicitly on the time due to the angular position of the Sun and thus the dynamics is no longer autonomous, i.e. the fourth, fifth and sixth components of the drift vector field F 0 depend explicitly on t:

F S 0 (q, t) =       ẋ, ẏ, ż 2 ẏ + x -(1-µ)(x+µ) 3 1 -µ(x-1+µ) 3 2 -(x-r S cos θ)µ S 3 S (t) -µ S cos θ r 2 S -2 ẋ + y -(1-µ)y 3 1 -µy 3 2 -(y-r S sin θ)µ S 3 S (t) -µ S sin θ r 2 S -(1-µ)z 3 1 -µz 3 2 -zµ S 3 S (t)       (7)
This implies that we need to take into account the angular position of the Sun at the initial time, knowing that its angular position on June 1 st 2006 is θ 06-01-2014 = 1.10439 rad. Then, once we select a rendezvous point on RH 120 it imposes the time of rendezvous, that we denote by t rdv and based on the transfer duration t f we can compute the Sun's angular position at the initial time using θ(t rdv -t f ) = θ(t rdv ) -ω S t f . This allows us to overcome the complexity of free final time when the Sun is included in the model.

Let us denote by q rdv (•) the rendezvous transfer trajectory, we consider the rendezvous transfer from an initial point q rdv (t 0 ) on a parking orbit O 0 ∈ IR 6 to a final position and velocity q rdv (t f ) on the RH 120 orbit. Note that the initial and final positions and velocities are variables of the global optimization problem. The criterion to maximize is the final mass which is equivalent to minimizing the fuel consumption or the delta-v=

t f t 0 Tmax u(t) m(t)
dt. Since the mass evolves proportionally to the norm of the thrust, our criterion is equivalent to the minimization of the L 1 -norm of the control:

min u∈U t f t 0 u(t) dt, (8)
where U = {u(•); measurable bounded and u(t) ≤ 1 for almost all t}, t 0 and t f are respectively the initial and final time.

Remark that since we choose O 0 to be a Halo orbit around a libration point, it is uniquely determined by a single point of the orbit using the uncontrolled CR3BP dynamics. This fact will play an important role for one of the necessary optimality conditions below. Notice, that even though the Halo orbit is a periodic orbit in the CR3BP only and not in the CR4BP we chose this orbit to be our hibernating location for the spacecraft in both situations which is clearly justified by the results of the Artemis mission.

The large number of variables in our problem adds complexity to the search for a solution. In particular, in the case of free final time we expect an infinite time horizon with a control structure that mimics impulsive maneuvers. To simplify our optimal problem we have two options, either we fix the transfer time or we fix the structure of the control. If we fix the final time, the sensitivity of the shooting method can be addressed by using the solution of a smoother criterion than the L 1 -norm, for instance the L 2 -norm, and linking it to the target criterion by a continuation procedure, see [START_REF] Gergaud | Low thrust minimum-fuel orbital transfer: a homotopic appoach[END_REF]. In this paper, however, we take a different approach that can be applied to both our models and decide to fix the structure of the control. In the sequel we focus on designing transfers associated to controls with a piecewise constant norm with value in {0, 1} and four switchings. This choice is motivated by two reasons. First, as mentioned above from a purely computational point of view, leaving the control structure free adds up complexity to the optimization problem that is difficult to address. Second, our desire is to mimic impulse strategy justifies a piecewise constant control either at its maximum or zero and preliminary calculations with a free number of switchings with several rendezvous points demonstrated that statistically the optimal control naturally converges to one with at most three boosts. One to leave the departing orbit, one to redirect the spacecraft in direction of the rendezvous point and a final one to match the position and location of the asteroid at the encounter. More precisely, we consider control functions that are piecewise continuous such that there exists times

t 0 ≤ t 1 ≤ t 2 ≤ t 3 ≤ t 4 ≤ t f with u(t) = 1 if t ∈ (t 0 ; t 1 ) ∪ (t 2 ; t 3 ) ∪ (t 4 ; t f ) 0 if ∈ (t 1 ; t 2 ) ∪ (t 3 ; t 4 ). (9)
Here the final time t f is free in the CR3BP and fixed in the CR4BP, as explained in the next paragraph. Note that our numerical method will be able to select the best control strategy even if it has less than three boosts. When dealing with the CR4BP, we decide to add some practical constraints to the rendezvous problem. First, we add a synchronization constraint stating that we know the position of the spacecraft on O 0 at the capture time of the TCO, that is June 1 st 2006 for RH 120 . This initial position is denoted q HaloL2 , and is given in the numerical results section. Moreover, we impose that the departure date of the rendezvous transfer has to occur once the trajectory of the TCO is known, that is after the date it is detected to which we add 30 days to take into account the necessary computations to predict accurately it's trajectory. We call detection date this detection plus trajectory prediction date. This means that it is impossible for the transfer to aim at a target point that the TCO visited before this detection date or even a few days after this date. These two constraints are aimed at depicting more realistic transfers. The considered departure dates will be fixed with a discretization of 15 days. The rendezvous dates will also be fixed and use a one day discretization from RH 120 's detection to RH 120 's escape date.

Necessary Conditions for Optimality

The maximum principle provides first order necessary conditions for a trajectory to be optimal [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]. Details regarding the application of the maximum principle to orbital transfers can be found in many references, including [START_REF] Caillau | Minimum fuel control of the planar circular restricted three-body problem[END_REF][START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]. We denote by X(t) = (q(t), m(t)) ∈ IR 6+1 the state, where q = (x, y, z, ẋ, ẏ, ż) is the position and velocity of the spacecraft and m its mass. The conditions are mostly the same for both model, so we first give the conditions for the Earth-Moon CR3BP and then give the modifications to apply when dealing with the problem in the Sun-perturbed model. The maximum principle applied to our optimal control problem, in the Earth-Moon CR3BP, states that if (q(•), m(•), u(•)) = (X(•), u(•)) is an optimal solution defined on [t 0 , t f ], then there exists an absolutely continuous adjoint state (p 0 , p X (•)) = (p 0 , p q (•), p m (•)), defined on [t 0 , t f ] such that:

• (p 0 , p X (•)) = 0, ∀t ∈ [t 0 , t f ], and p 0 ≤ 0 is a constant. • Let H, the Hamiltonian, be H(t, X(t), p 0 , p X (t), u(t)) = p 0 u(t) + p X (t), Ẋ(t) , then Ẋ(t) = ∂H ∂p X (t, X(t), p 0 , p X (t), u(t)), for a.e. t ∈ [t 0 , t f ], (10)
ṗX (t) = - ∂H ∂X (t, X(t), p 0 , p X (t), u(t)), for a.e. t ∈ [t 0 , t f ], (11)
where , denotes the inner product.

• H(t, X(t), p 0 , p X (t), u(t)) = max ν ≤1 H(t, X(t), p 0 , p X (t), ν), ∀t s.t. u(t) = 1 (maximization condition).

•

Ψ(t i ) = 0 for i = 1, • • • , 4. • H(t f , X(t f ), p 0 , p X (t f ), u(t f )) = 0, if t f is free.
• p q (t 0 ), F 0 (q(t 0 )) = 0 (initial transversality condition).

• p m (t f ) = 0.

The function Ψ(•) is the so-called switching function corresponding to the problem with an unrestricted control strategy and we have Ψ(t

) = p 0 + T max pv(t) m(t) -p m (t)β . The maximization condition of the Hamiltonian H is used to compute the control on [t 0 ; t 1 ] ∪ [t 2 ; t 3 ] ∪ [t 4 ; t f ] and we have u(t) = pv(t) pv(t) for all t ∈ [t 0 ; t 1 ] ∪ [t 2 ; t 3 ] ∪ [t 4 ; t f ].
The initial transversality condition reflects the fact that the initial departing point is free on the Halo orbit O 0 . Remark that since the data for the TCO's trajectory are given as ephemerides, there are no dynamics equations to describe those orbits in the CR3BP or the Sunperturbed model. Thus, we cannot compute the tangent space to a TCO point and we cannot extract a transversality condition for p q (t f ) at the rendezvous point. Since we expect numerous local extrema for this optimal control problem, it is however preferable to solve the problem for fixed rendezvous points on a discretization of the TCO orbit.

When dealing with the problem in the CR4BP, the transfer duration is determined by the rendezvous point on the TCO orbit and the initial time of the transfer. Furthermore, as the position of the spacecraft at RH 120 's capture, on June 1 st 2006, is fixed to q HaloL2 , knowing the initial time gives the initial position of the spacecraft on the Earth-Moon Halo orbit. Thus, since the detection constraint leads us to consider a finite number of departure dates, each problem in the CR4BP has a fixed initial and final time. Thus, the necessary conditions of optimality of this problem do not include the initial transversality condition and the final Hamiltonian cancellation. In addition to these two changes, the Sun-perturbed model has a different expression of Ẋ.

Shooting Method

Our numerical method is based on the necessary conditions given in section 3.2.1. Let Z(t) = (X(t), p X (t)), t ∈ [t 0 ; t f ], and u(q, p) the feedback control expressed using the maximization condition. Then, we have Ż(t) = Φ(Z(t)) where Φ comes from equations [START_REF] Dunham | Using lunar swingbys and libration-point orbits to extend human exploration to interplanetary destinations[END_REF] and [START_REF] Folta | Earth-Moon libration point orbit station keeping: Theory, modeling, and operations[END_REF]. The goal is to find Z(t 0 ), t 1 , t 2 , t 3 , t 4 and t f (when we consider the free final time for the CR3BP) such that the following conditions are fulfilled:

1. Ψ(t i ) = 0 for i = 1, • • • 4; 2. X(t f
) is the prescribed rendezvous point; 3. X(t 0 ) ∈ O 0 , and the initial transversality condition is verified (only for the CR3BP); 4. H(t f ) = 0 in the CR3BP with free final time.

The problem has been transformed into solving a multiple points boundary value problem. More specifically, we must find the solution of a nonlinear equation S(Z(t 0 ), t 1 , t 2 , t 3 , t 4 , t f ) = 0, respectively S(p X (0), t 1 , t 2 , t 3 , t 4 ) = 0 in the CR4BP, where S is usually called the shooting function. When looking for Z(t 0 ), we are actually only looking for (t 0 , p X (0)) since X(t 0 ) is completely defined by one parameter.

The evaluation of the shooting function is performed using the high order numerical integrator DOP853, see [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF]. The search for a zero of the shooting function is done with the quasi-Newton solver HYBRD of the Fortran minpack package. Since S(•) is nonlinear, the Newton-like method is very sensitive to the initial guess. This leads us to consider heuristic initialization procedures. We combine two types of techniques, a direct method and a continuation method. The discretization of the TCO's orbit requires the study of thousands of transfers, we use a direct method for a dozen rendezvous and expand to other points on the orbit using a continuation scheme. The motivation is that direct methods are very robust but time consuming while continuation methods succeed for our problem in most cases very rapidly. The direct method uses the modeling language Ampl, see [START_REF] Fourer | AMPL: A Modeling Language for Mathematical Programming[END_REF], and the optimization solver IpOpt, see [START_REF] Waechter | On the Implementation of an Interior-Point Filter-Line Search Algorithm for Large-Scale Nonlinear Programming[END_REF], with a second order explicit Runge-Kutta scheme. Details on advanced continuation methods can be found in [START_REF] Gergaud | Low thrust minimum-fuel orbital transfer: a homotopic appoach[END_REF] and [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF], but we use a discrete continuation which is enough for our needs.

Moreover, note that in order to fulfill the initial transversality condition, we first prescribe X(t 0 ) on the parking orbit and find a zero of the shooting function satisfying all the other conditions. Afterward, we do a continuation on X(t 0 ) along the departing parking orbit in the direction that increases the final mass. Once a framing of the best X(t 0 ) with respect to the final mass is found, we perform a final single shooting to satisfy the initial transversality condition along with the other conditions. This decoupling in the search of an extremal is motivated by the fact that we could very likely find a local maximum on X(t 0 ) rather than a local minimum because of the periodicity of the initial parking orbit. We avoid this fact by first manually ensuring that the X(t 0 ) we find will be the one for the best final mass and not the worst. However, our continuation procedure on X(t 0 ) does not always succeed, mainly because of the high nonlinearity of the shooting function, as the trajectories we find can be very long. Even if some of the extremals we find do not satisfy the initial transversality condition with the aimed accuracy (typically a zero of the shooting function is deemed acceptable if S(Z(t 0 ), t 1,2,3,4,f ) ≤ 10 -8 ), they are still rather close to satisfy it (of the order of 10 -4 ). For theCR4BP, we took a different approach to limit the number of calculations and continuations. First as we mentioned previously, we still assume that the spacecraft is parked on an Halo orbit around the L 2 equilibrium point of the CR3BP but we now impose its location at the time of RH 120 capture. Its departure position from its hibernating orbit is then determined by the prescribed rendezvous point and the fixed duration of the transfer in that case. This permits to ignore the initial transversality condition and simplifies the algorithm.

Numerical Simulations and Results

This section is divided into subsections as follow. In section 4.1 we compute the best rendezvous transfer from the hibernating Halo orbit to RH 120 assuming a free final time. Since we assume a free final time the CR3BP model is used for our calculations. The main goal is to determine a lowest bound estimate for the delta-v with respect to the position of the rendezvous point on RH 120 orbit and the duration of the transfer, as well as to gain insight on specific characteristics of this rendezvous point. Our calculations are also expanded to four other TCOs to demonstrate the generality of our algorithm. In section 4.2, by fixing the final time we produce a plausible scenario with respect to detection time and transfer duration using the insights found in 4.1 and the CR4BP model. Section 4.3 expands on the analysis of the characteristics of the rendezvous points on RH 120 that produce the most efficient transfers, our observations are supplemented by a statistical study on a larger pool of TCOs. For all our calculations the spacecraft is assumed to be a monopropellant engine, and its characteristics are assumed to be an initial mass of 350 kg, a specific impulse I sp of 230 s. and a maximum thrust T max of 22 N . This choice is motivated by assuming a spacecraft with similar features than the Gravity Recovery and Interior Laboratory (GRAIL) spacecrafts. The Halo orbit from which the spacecraft is departing is chosen to have a z-excursion of 5,000 km around the EM libration point L 2 , see Figure 4. The point corresponding to the positive z-excursion is q HaloL 2 = (1.119, 0, 0.013, 0., 0.180, 0), and the period of this particular Halo orbit is t HaloL 2 = 3.413 in normalized time units or approximatively 14.84 days. During the period represented in Figure 1, asteroid RH 120 does 17 clockwise revolutions around the origin of the CR3BP frame, and 3.6 revolutions in Earth inertial reference frame. The evolution of the energy of RH 120 and its distance to the Earth-Moon libration point L 2 as the asteroid evolves on its orbit are given in Figure 5.

Rendezvous to RH 120 using CR3BP, free final time

The objective of this section is to determine the best rendezvous transfer to RH 120 assuming a free final time, and to present an analysis of the evolution of the fuel consumption with respect to the location of the rendezvous point on the RH 120 orbit. As mentioned in section 3.2, we restrict the study to transfers with at most three boosts.

To analyze the variations of the fuel consumption with respect to the rendezvous point on the orbit, we discretize uniformly the orbit of RH 120 using 6 hours steps. For each rendezvous point of this discretization, we compute an extremal transfer (i.e. a solution of the maximum principle) with free final time using the techniques explained in section 3.

Best Rendezvous Transfer

The best transfer is represented on Figure 6, and some relevant data is presented in Table 2. This transfer has a delta-v of 203.6 m/s. The rendezvous takes place on June 26 th 2006 and lasts 415.5 days which would require detection and launch about 14 months before June 1 st 2006. On Figure 7, we display the orbit of RH 120 in both the rotating and inertial frame with the rendezvous point for the best transfer. The best transfer exhibits 14 revolutions around the origin (in the rotating frame) and has a significant variation in the z-coordinate with respect to the EM plane. In particular, the z-coordinate of the rendezvous point is -1.04 normalized units, that is about 400,000 km, and the maximum z-coordinate along the trajectory is 5.34 normalized units, that is about 2 million km. The departure point on the parking orbit occurs 4.5 days after q HaloL 2 . The control strategy consists of a first boost of 19.7 s, a second boost starting after 154.7 days and lasting 51.1 min , finally, the last boost starting 261.1 days after the second boost and lasts 13.8 s, see Figure 8. The short initial boost could suggest an initial jump on an unstable invariant manifold but the boost direction does not match the eigenvector of the monodromy matrix associated to the unstable invariant manifold, see [START_REF] Koon | Dynamical Systems, the Three-Body Problem and Space Mission Design[END_REF]. However, this trajectory exploits the fact that a small initial boost leads to a far location where the gravity field of the two primaries is small and where the second boost can efficiently aim at the rendezvous point. In particular, we expect the existence of other local minima with a larger final time, going further away from the initial and final positions and providing an even better final mass. Also note that this kind of strategy could not have been obtained if we had restricted the control structure to have two boosts rather than three.

Fuel consumption with respect to rendezvous point

Figure 9 shows the evolution of the final mass, the delta-v and of the duration of the transfer with respect to the rendezvous point on the RH 120 orbit for a spacecraft departing from the Halo L 2 orbit and corresponding to the three boost control strategy. As explained in section 3.2.2, some of the departure points on the Halo orbit are not fully optimized -this is the case for about two thirds of them. Also note that a departure point different from q HaloL 2 implies a drift phase whose duration is not included in the transfer duration. It can be observed that the final mass has many local extrema and that the variation of the duration of the mission is not continuous (contrary to what we would expect). It is most likely due to local minima or to the fact that the value of the criterion (see Equation ( 8)) is discontinuous with respect to the rendezvous point. Unfortunately it is not possible to show a minimum is local without exhibiting a better minimum. As for the regularity of the criterion, its study would lead to the analysis of an ad hoc Hamilton Jacobi Bellmann equation that is out of scope of this paper as a numerical study is not possible due to the dimension of the problem.

Figure 10 shows the evolution of the departure point on the hibernating orbit for the spacecraft with respect to the rendezvous point on RH 120 (top) as well as the three most frequent departure points on the initial Halo orbit (bottom). More precisely, we represents the optimized argument of (y, z)(0) (up to the quadrant: arctan(z(0)/y(0))). Note that the initial position on the Halo orbit directly depends on the optimized initial drift time. Since this initial drift phase has not always been successfully optimized, this figure has to be interpreted with caution. However, we can see that the departure position on the initial periodic orbit seems to always be close to a multiple of π.

Comparing the evolution of the final mass from Figure 9 and the evolution of RH 120 energy from Figure 5, we can see that the best final masses are obtained on the first half of RH 120 trajectory, that is for rendezvous points with energies closer to the departing energy. For the best transfer, the energy difference between the rendezvous point on RH 120 and the departing point on the Halo orbit is about 0.046. The rendezvous point with the closest energy to the initial orbit occurs only slightly before the optimal rendezvous point. Its final mass is 319.3 kg which is only 0.5 kg worst than the best final mass. This remark suggests strongly that a small difference in energy between the rendezvous point on RH 120 orbit and the departing point for the spacecraft on the Halo orbit is advantageous.

Expansion to other TCOs

Further calculations on four synthetic TCOs obtained from the database produced in [START_REF] Granvik | The population of natural Earth satellites[END_REF] suggest that similar results can be expected on a large sample of TCOs. Indeed, in Table 3 we display data regarding the best three boost transfer for four other TCOs. These transfers have a duration of about one year each and produce delta-v values between 223.9 m/s and 344.2 m/s. The best transfer is for TCO 1 and occurs for an energy difference between terminal configurations of about 0.13, for TCO 16 the difference is about 0.6, 1.5 for TCO 19 and 2.3 for TCO 11 . This reinforces the relationship between the final mass of the transfer and the energy difference of the rendezvous point with respect to the one from the Halo orbit, see section 4.3 for more details. It can also be noted that the maximum distance of the spacecraft from L 2 during the transfer is similar for all four transfers which indicates that the long drift is used to pull away from the two primaries attraction fields to make the second boost more efficient.

Parameter

Symbol 4.2. Rendezvous to RH 120 using CR4BP, fixed final time In this section, we use the insights from the results obtained in section 4.1 to develop a more realistic scenario taking into account the Sun perturbation, the detection time of the asteroid, the time required to compute an accurate orbit and a transfer duration that will allow to reach RH 120 within its capture time. Asteroid RH 120 was actually detected 105 days after its capture by Earth gravity, on September 14th, 2006. Moreover if we add a window of 30 days to allow completion of the observations and calculations required to predict RH 120 orbit with enough precision, we obtain that the departure time must take place on or after 135 days after its capture on June 1 st 2006. For these simulations, the spacecraft is assumed to be at q HaloL 2 on June 1 st 2006.

Based on the calculations from section 4.1, we analyze all possible transfers departing from the Halo orbit 135 days after RH 120 capture by discretizing the transfer duration over the interval [100,290] where 290 represents the constraint that the spacecraft must reach RH 120 on or before it escapes Earth's gravity on July 31th 2007, 425 days after initial capture. We chose not to consider transfer durations that last less than 30 days since section 4.1 demonstrated that a shorter time implies a higher delta-v. We use a discretization of 15 days of the departure dates over the predefined interval [100,290]. The best rendezvous transfer using the Sun-perturbed model departing the Halo orbit exactly 30 days after RH 120 was detected has a duration of 255 days and the final mass is m f = 245.707kg, or equivalently delta-v=797.991 m/s, see Figure 11. Additional simulations on different starting dates and durations show that the best departure time that takes place at least 30 days after detection and before it escapes capture is 180 days after RH 120 capture. This scenario provides 75 days between the detection time and the departure of the spacecraft for the rendezvous mission to determine RH 120 orbit and design the transfer. It produces a final mass of 267.037 kg, or equivalently delta-v=610.224 m/s, and a rendezvous date 312 days after capture, that is April 9th 2007. Figure 12 shows the corresponding trajectory in an inertial frame and in the CR3BP frame, while Table 4 summarizes the main features of this transfer.

The three boosts of this transfer last respectively 4.1 mn, 0.89 h and 84.2 mn, while the two ballistic arcs durations are 70.6 and 61.3 days. Clearly from our simulations the Sun has an impact on the rendezvous transfers, in particular when considering the Sun the spacecraft does not drift as far from the Earth-Moon system since the Sun's gravity is now unavoidable. Notice also that for both CR4BP transfers the spacecraft shows a passage at close proximity to the Earth to quickly modify its energy level. More precisely, in Figure 11 the smallest distance between the spacecraft and the Earth is 105478 km (0.274397 LD) and in Figure 12 it is 35677.2 km (0.0928128 LD). Both of these are acceptable in terms of practicality with the second one at a distance about equivalent to the geostationary orbit.

Regression analysis

The goal of this section is to determine factors that can predict which location on a given TCO orbit is best suited for a fuel efficient rendezvous transfer. In section 4.1.2, we observed that for RH 120 a small energy difference between the departing point for the spacecraft and the rendezvous point provides a good candidate for an efficient transfer. Using rendezvous transfers on a larger pool of TCOs and transfers calculated using the Sun-perturbed model, we expand on this observation.

More precisely, we have that in the circular restricted three-body problem the energy E 3 is a first integral of motion, so any change in energy must be generated from the spacecraft's propulsion. Therefore, transfers between orbits of significantly different energy necessarily require more delta-v which is what we observed in section 4.1.2. In the circular restricted four-body problem however, energy is no longer a first integral, and theoretically the spacecraft can utilize the influence of the Sun to navigate between orbits of different energy for cheaper delta-v. For 1,000 selected rendezvous points on simulated TCOs from [START_REF] Granvik | The population of natural Earth satellites[END_REF], we computed a rendezvous transfer departing from the Earth-Moon libration point L 2 in the circular restricted four-body model, using the same features for the spacecraft as for our calculations on RH 120 , i.e. a 22N maximum thrust and 230 second I sp and transfer duration between 10 and 180 days. We see in Figure 13 a characterization of those rendezvous points based on the energy difference from departure |E 3 (rdvz) -E 3 (L 2 )| as well as the z and ż coordinates, colored by delta-v. The figure shows clearly that those rendezvous points with low absolute z and ż coordinates yield the lowest delta-v values. Moreover, although there are some low delta-v rendezvous with large energy difference, most of the low delta-v transfers have low energy differences, indicating the energy is still an important factor even in the four body model.

Further exploration to determine other predictors has been conducted as follows. We started from a random pool of 91 TCOs and their corresponding best rendezvous transfers. We discarded some extreme cases: transfers with inordinate expected fuel cost (delta-v> 1000) and TCOs with average energy vastly different from that of the departure point (absolute value of the differences greater than 4). Running a linear multiple regression using regstats in Matlab on the remaining 79 trajectories, we found that three predictors were statistically significant with 95% confidence. Listed below are the results. These delta-v values ranged from 130 to 1600 m/s.

Statistic Value

Adjusted R 2 0.397 f -value 18.106 p-value 6.17 × 10 -9

Table 5: Selected statistics for the multiple linear regression on the three significant variables of interest. This indicates that the three selected predictors explain roughly 40% of the variance in the data. The f -statistic and associated p-value indicate that this cannot be explained by random chance.

The regression shows that these three predictors account for approximately 40% of the variance in the transfer cost (adj-R 2 = 0.397), which implies that they may be useful in predicting future transfer costs, but that there may exist some other significant factors as well. Based on this data, we hypothesize the following practical trends from these findings:

• TCOs with average energy similar to the spacecraft's departure energy tend to have lower delta-v transfers.

Predictor Description p-value R 2

Energy Difference

The difference in the energy of the TCO and the spacecraft's departure point (equal in this case to the EM L 2 energy) 0.0169 0.163

Lunar Planarity

The average distance of the TCO from the plane in which the Moon orbits the Earth. 4.35 × 10 -6 0.216

Barycenter Variance

The variance in the distance of the TCO from the Earth-Moon barycenter. This is a rough indicator of how well the TCO adheres to a circular path around the earthmoon system.

3.10 × 10 -3 0.164 Table 6: Significant predictors of delta-v, the associated p-value from the multiple linear regression t-tests, and the associated R 2 value from the single linear regression. Subject to the inherent difficulties of regression analysis, this implies that Lunar planarity has the most impact on fuel costs while the other two predictors have roughly the same impact.

• TCOs which travel more aligned with the lunar plane tend to have lower delta-v transfers.

• TCOs which have roughly circular paths, or have perhaps even large sections of roughly circular paths, tend to have lower delta-v transfers.

It should be noted that these hypotheses come with fair caveat; we did not vary the departure energy, the sample size of TCOs was fairly low, and the pruning method for excluding extreme cases relies on calculating the transfer cost initially. However, the TCOs excluded due to extreme transfer costs alone (high delta-v and not high geocentric energy difference) were rare (3 out of 91), so we can probably safely ignore these as outliers.

Acknowledgements

We would like to specially thank Robert Jedicke and Mikael Granvik for their generous access to the database of synthetic TCOs, and in general, their help and support for this research. Geoff Patterson and Monique Chyba were partially supported by the National Science Foundation (NSF) Division of Mathematical Sciences, award #1109937.

Figure 1 :

 1 Figure 1: Orbit of RH 120 in the Earth-Moon CR3BP rotating reference frame.

Figure 2 :

 2 Figure 2: Orbit of RH 120 in the Earth-Sun CR3BP rotating reference frame.

Figure 3 :

 3 Figure 3: Control function modeling thrust impulses over time.

Figure 4 :

 4 Figure 4: Halo orbit from which the spacecraft is departing, z-excursion of 5,000 km around the EM libration point L 2 .

Figure 5 :

 5 Figure 5: Evolution of the energy (left) and of the Euclidean distance to the L 2 Libration point (right) for RH 120 . For energy, the horizontal line represents the energy of the Halo periodic orbit around L 2 which is about -1.58.

Figure 6 :

 6 Figure 6: Best 3 boost rendezvous transfer to RH 120 from a Halo orbit around L 2 : Inertial Frame (left). Rotating Frame (Right). Bottom: zoom on the start of the departure from the Halo orbit (rotating frame).

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Best 3 boost rendezvous transfer to RH 120 from a Halo orbit around L 2 : Inertial Frame (left). Rotating Frame (Right)

Figure 10 :

 10 Figure 10: The top picture represents the evolution of (y, z)(0) argument with respect to time of rendezvous on RH 120 . The bottom picture represents the three most frequent departure points on the initial Halo orbit.

Figure 11 :

 11 Figure 11: 3 boost rendezvous transfer to RH 120 from a Halo orbit around L 2 in inertial frame, using the CR4BP. Departure 135 days after June 1 st 2006, arrival 255 days later.

Figure 12 :

 12 Figure 12: Best 3 boost Sun-perturbed rendezvous transfer to RH 120 from a Halo orbit around L 2 in inertial frame, with transfer departure after RH 120 detection.

Figure 13 :

 13 Figure 13: Characterization of 1,000 simulated rendezvous transfers from L 2 .

Table 1 :

 1 Numerical values for the CR3BP and the CR4BP.

	CR3BP parameters		CR4BP parameters
	µ	1.2153 • 10 -1 µ S	3.289 • 10 5
	1 norm. dist. (LD) 384400 km r S	3.892 • 10 2
	1 norm. time	104.379 h	ω S -0.925 rad/norm. time

Table 2 :

 2 Best Transfer Data from Halo-L 2

	Parameter	Symbol Value
	Transfer Duration	t f	415.8 days
	Final Mass	m f	319.8 kg
	Delta-v	∆v	203.6 m/s
	Final Position	q rdv p	(2.25, 3.21, -1.04)
	Final Velocities	q rdv v	(2.92, -2.02, 0.46)
	Max Distance from L 2 d max L 2	10.63 LD

Table 2 :

 2 Data for the best transfer from q HaloL2 to asteroid RH 120 .

  TCO 1 TCO 11 TCO 16 TCO 19

	Transfer Duration (days) t f	362.0	386.6	362.2	364.9
	Final Mass (kg)	m f	316.9	300.5	311.0	307.1
	Delta-v (m/s)	∆v	223.9	344.2	266.1	294.6
	Max Distance from L 2 Time to d max L 2 (days)	d max L 2 t(d max L 2 )	12.7 232.9	11.5 196.6	11.5 197.3	12.8 232.3

Table 3 :

 3 Best transfer data from q HaloL2 to selected TCOs.

Table :

 : Best Sun-perturbed Transfer Data from Halo-L 2

	Parameter	Symbol Value
	Transfer Duration	t f	132 days
	Final Mass	m f	267 kg
	Delta-v	∆v	610.2 m/s
	Final Position	q rdv p	(2.27, 0.46, -0.37)
	Final Velocities	q rdv v	(0.79, -2.22, -0.58)
	Max Distance from L 2 d max L 2	4.81 LD

Table 4 :

 4 Data for the best Sun-perturbed transfer from q HaloL2 to asteroid RH 120 , with transfer departure after RH 120 detection.