# Systematic Design of Optimal Low-Thrust Transfers for the Three-Body Problem

###### Abstract

A computational approach is developed for the design of continuous low thrust transfers in the planar circular restricted three-body problem. The transfer design method of invariant manifolds is extended with the addition of continuous low thrust propulsion. A reachable region is generated and it is used to determine transfer opportunities, analogous to the intersection of invariant manifolds. The reachable set is developed on a lower dimensional Poincaré section and used to design transfer trajectories. This is solved numerically as a discrete optimal control problem using a variational integrator. This provides for a geometrically exact and numerically efficient method for the motion in the three-body problem. A numerical simulation is provided developing a transfer from a periodic orbit in the Earth-Moon system to a target orbit about the Moon.

15-757

## 1 Introduction

Designing spacecraft trajectories is a classic and ongoing topic of research. With the current fiscal constraints there is an increased focus on technologies with a critical impact on mass. Optimal expenditure of onboard propellant is critical to allowing a mission to continue for a longer period of time or to enable the launch of a less massive spacecraft. Electric propulsion systems offer a much greater specific impulse than chemical systems and are able to operate for extended periods of time. However, these electric propulsion systems typically have much less thrust than their chemical counterparts and therefore must operate over longer durations in order to impart the desired momentum change. Recent developments in miniature electric propulsion offer the potential for new research opportunities for small spacecraft [(1)]. With reduced development intervals and decreased launch costs, small satellites have gained increased attention as a cost effective means of scientific and technologic development. The merger of small satellites with miniature electric propulsion enables inexpensive and responsive missions requiring large changes in orbital /media/arxiv_projects/487592/energy or extended mission lifespan. With the potential for more demanding missions, even greater importance is placed on the mission design to ensure that optimal trajectories satisfy mission requirements. In addition, non-Keplerian orbits and multi-body dynamics have been shown to allow for a much greater range of potential missions at a reduced /media/arxiv_projects/487592/energy cost [(2)]. Future space missions are increasing in complexity and will require new classes of orbits that are not possible via the traditional conic approach [(3), (4)]. Optimally combining the dynamical structure of the three-body problem with low-thrust propulsion systems is vital for future mission success.

There has been extensive research focused on optimal control for spacecraft orbital transfers in the three-body problem [(5), (6)]. Typically, the optimal control problem is solved via direct methods, which approximate the continuous problem as a parameter optimization problem. The state and/or control trajectories are parameterized and solved in the form of a nonlinear optimization problem. References 5 and 6 use this direct approach in designing low-thrust transfers in the three-body problem. Alternatively, indirect methods apply calculus of variations to derive the necessary conditions for optimality. This yields a lower dimensioned problem than the direct approach and algebraic conditions that, when satisfied, guarantee local optimality in contrast to direct methods which result in sub-optimal solutions.

The application of optimal control methods for orbital /media/arxiv_projects/487592/trajectory design is nontrivial. The three-body system dynamics are nonlinear and exhibit chaotic behaviors. Small changes in initial conditions result in large variations of the resulting system /media/arxiv_projects/487592/trajectory. Therefore, any optimization routine is highly sensitive to the initial guess. In addition, insight into the problem or intuition on the part of the designer is often required to determine initial conditions that will converge which achieve satisfactory results. Efficient numerical implementation is dependent on correct initial conditions as well as accurate numerical integration.

Additionally, References 5 and 6 implement the solutions using conventional Runge-Kutta integration techniques. These techniques suffer from numerical instability and /media/arxiv_projects/487592/energy drift behaviors which make them ill-suited for long-term propagation. These dissipative effects are even more detrimental with the addition of low-thrust propulsion to the dynamic equations of motion. Conventional integration techniques fail to capture the physical laws and geometric properties of the dynamic system. As a result, the long term effects of low-thrust on the spacecraft /media/arxiv_projects/487592/trajectory are not accurately captured.

References 2 and 3 have illustrated the rich structure that exists in the three-body problem. Within the three-body problem, a spacecraft’s feasible region of motion is constrained by its /media/arxiv_projects/487592/energy, or Jacobi integral. It has been shown that there exist multi-dimensional tubes, or invariant manifolds, of constant /media/arxiv_projects/487592/energy trajectories that span the state space. Associated with periodic solutions of the three-body dynamics, these invariant manifolds allow for the spacecraft to traverse vast expanses of the state space with zero /media/arxiv_projects/487592/energy change. However, the results presented are highly case specific and difficult to generalize to arbitrary transfers. Also, these results are based on control-free trajectories which rely on the underlying structure of the three-body system. In addition, transfer orbits along an invariant manifold require a longer time of flight which may be undesirable for time critical missions. The addition of low-thrust propulsion offers the potential of reduced transit times and the ability to depart from the free motion /media/arxiv_projects/487592/trajectory to allow for increased transfer opportunities.

In order to address these issues, the authors propose an accurate and numerically stable method for long-term optimal orbital transfers. This approach will avoid the instability and dissipative effects of conventional integration schemes. In addition, the effects of low-thrust propulsion will be directly included in the integration method to ensure that extended duration optimal trajectories are computed correctly. This improved method will enable the derivation of a systematic method of generating optimal transfer orbits between arbitrary states. Indirect optimal control, based on the calculus of variations, will be implemented in order to generate optimal orbital transfers. This will avoid the approximation issues inherent in the previous work, which utilized direct optimal control methods. With this proposed method, the previous research on control-free trajectories will be generalized with the addition of low-thrust propulsion systems.

To achieve these objectives, computational geometric optimal control techniques are applied. The dynamics of the three-body system are derived from the discrete Lagrangian, which approximates the integral of the continuous time Lagrangian over a fixed discrete step. Application of the discrete Euler-Lagrange equations, or the discrete Legendre transform, results in the discrete equations of motion. This discrete update map, or variational integrator, shares the same geometric properties of the continuous time system and exhibits much better /media/arxiv_projects/487592/energy behavior than the traditional integration methods, especially over long time periods. A discrete optimal control problem is formulated from the discrete equations of motion. This approach, where explicit discretization occurs prior to optimization, is in contrast to the typical method, where the equations of motion are implicitly discretized during the optimization procedure. Formulating the problem in this manner results in more stable and accurate optimal solutions. In indirect methods the optimal control problem is expressed as a two-point boundary value problem. Optimal solutions are generally sensitive to small variations in the initial multipliers. As a result, the numerical stability of sensitivity derivatives is critical to accuracy and computational performance. The use of geometric integrators, which do not suffer the numerical dissipation of conventional integration methods, results in a more robust and efficient solution.

A discrete optimal control problem is formulated to determine the reachability set on a Poincaré section. Given an initial condition and fixed time horizon, the reachable set is the set of states attainable, subject to the operational constraints of the spacecraft. Generation of this reachability set allows for the extension of the previous control-free missions in the three-body problem. In addition, the generation of the reachable set allows for a more systematic method of determining initial conditions and eases the burden on the designer. The addition of low-thrust propulsion affords the ability to enlarge the reachable set as compared to the control-free case. Maximization of the reachability set, on an appropriately chosen Poincaré section, allows for a greater space of potential transfer trajectories. The use of the Poincaré section allows for design on a lower dimensional space and simplifies the design process. Once an intersection is determined on the Poincaré section a transfer is calculated.

In short, the authors present a systematic method of generating optimal transfer orbits in the three-body problem. Previous results in the design of optimal transfers have relied on suboptimal direct optimization methods and conventional integration techniques. This paper provides a discrete optimal control formulation to generate the reachability set on a Poincaré section. The use of a geometric integrators ensures numerical stability for long-duration orbit transfers. A numerical example is presented demonstrating a transfer /media/arxiv_projects/487592/trajectory from to the vicinity of the Moon and compared to the control-free design methodology.

## 2 System Model

The system model is based on the planar circular restricted three body problem (PCRTBP). The Earth is assumed to be the more massive primary, , while the Moon is the second, smaller primary . The equations of motion are developed in a rotating reference frame which allows for much greater insight into the structure of the dynamics. The axis is directed along the vector from the Earth to the Moon. The axis lies in the orbital plane and is orthogonal to . The rotating reference frame is centered at the system barycenter. It is assumed that the rotates with a constant angular velocity equal to the mean motion of the Moon. Following convention, the system is also non-dimensionalized by the characteristic mass, length, and time koon2000 (). As a result, the system can be characterized by a single mass parameter ,

(1) |

The larger primary, , is located at and the smaller is located at . In the rotating reference frame the Lagrangian is given by

(2) |

where the distance and of the spacecraft to each primary is defined in rotating coordinates as

(3) | ||||

(4) |

Application of the Euler-Lagrange equations results in following equations of motion defined in the rotating reference frame

(5) |

where the effective potential is defined as

(6) |

and the control inputs are defined as and . The state is defined as with and representing the position and velocity with respect to the system barycenter, respectively. The equations of motion may be rewritten in state space form as

(7) |

where the matrix and psuedo gravitational potential gradient are

(8) |

(9) |

### 2.1 Jacobi Integral

There exists a single integral, or constant of motion for the three-body problem lanczos1970 (); szebehely1967 (). This /media/arxiv_projects/487592/energy constant is analogous to the total mechanical /media/arxiv_projects/487592/energy, however is a non-physical quantity arising from the problem formulation szebehely1967 (). Also known as the Jacobi constant, it is defined as a function of the position and velocity in the rotating frame and given by

(10) |

Equation 10 divides the phase space into distinct regions based on the /media/arxiv_projects/487592/energy level of the spacecraft. Fixing the Jacobi integral to a constant defines zero velocity curves, which are the locus of points where the kinetic /media/arxiv_projects/487592/energy, and hence velocity vanishes. As seen in Figure 1, the phase space is divided into distinct realms based on the /media/arxiv_projects/487592/energy level. In the vicinity of or there exists a potential well. As the /media/arxiv_projects/487592/energy level increases there are five critical points of the effective potential Equation 6 where the slope is zero. Three collinear saddle points on the axis and two equilateral points. These equilibrium, or Lagrange points, are labeled and are shown in Figure 1. The Jacobi integral is a valuable invariant property of the three-body system that allows for greater insight into the motion of the spacecraft.

## 3 Variational Integrator

Geometric numerical integration deals with numerical integration methods which preserve the geometric properties of the flow of a differential equation, such as invaraint properties and symplecticity. Variational integrators are constructed by discretizing Hamilton’s principle rather than the continuous Euler-Lagrange equations marsden2001 (). As a result, integrators developed in this manner have the desirable properties that they are symplectic and momentum preserving. In addition, they exhibit improved /media/arxiv_projects/487592/energy behavior over long integration periods. A short background on the variational principle for mechanical systems is presented. A discrete approximation of the action integral is presented and allows for construction of a variational integrator for the PCRTBP.

### 3.1 Variational Principle

Consider a continuous mechanical system described by the Lagrangian, , for the generalized position, , and velocity, . In the standard approach of variational mechanics the action integral is formed by integrating the continuous Lagrangian along a path that the system follows from time to greenwood1988 (). In the continuous time the action integral is defined as

(11) |

Hamilton’s principle states that the actual path followed by a holonomic system results in a stationary action integral with respect to path variations for fixed endpoints. Taking the variation of Equation 11 gives

(12a) | ||||

(12b) | ||||

(12c) |

where we have used integration by parts and the conditions and . For Hamilton’s principle to be valid for all admissible variations , the integrand of Equation 12 must be zero for all , giving the continuous Euler-Lagrange Equations lanczos1970 ().

(13) |

Hamilton’s equations are derivable through the use of the Legendre transformation which is a mapping where is the generalized momenta.

(14) |

In the continuous time case the Hamiltonian is defined as

(15) |

Applying Equation 14 and taking the variation of Equation 15 allows us to derive the equations of motion in Hamiltonian form

(16a) | ||||

(16b) | ||||

(16c) |

Both Equations 16 and 13 result in equations of motion for the mechanical system and are equivalent via the Legendre transform. Equation 13 results in second order differential equations while Equation 16 results in first order differential equations.

### 3.2 Discrete Variational Mechanics

A discrete analogue of Hamilton’s principle and the action integral is formed. Rather than taking a position, , and velocity, , consider two positions and and a fixed time step . The two positions are points on the curve such that and . A discrete time Lagrangian is formed which approximates the action integral between and as

(17) |

Since Section 3.3 is calculated as a numerical integral, an appropriate quadrature rule is required. There are multiple possible methods one can use to approximate the integral in Section 3.3. An appropriate approximation rule is determined based on the ease of implementation and accuracy desired.

Rectangle | |
---|---|

Midpoint | |

Trapezoidal |

Table 1 shows several possible approximation rules that are typically applied. The rectangle rule is a first order accurate method and offers a straightforward implementation. The midpoint and trapezoidal rules are both second order accurate methods. However, the midpoint rule results in an implicit form which adds further complexity to the equations of motion. In this work, the trapezoidal approximation is applied to the PCRTBP.

Once an appropriate discrete Lagrangian is formed a discrete action sum is formed as the discrete analogue of Equation 11

(18) |

Once again a discrete version of Hamilton’s principle is applied to Equation 18. Applying summation by parts yields

(19a) | ||||

(19b) |

For the discrete action sum to be stationary with respect to all admissible path variations, with fixed endpoints, the discrete Euler-Lagrange equations must be satisfied for resulting in

(20) |

A discrete version of the Legendre transformation, referred to as a discrete fiber derivative, results in the equivalent Hamiltonian form expression. The discrete fiber derivative is given as

(21a) | ||||

(21b) |

This yields a discrete Hamiltonian map . A more extensive development of variational integrators can be found in Reference 9.

### 3.3 Discrete Equations of Motion

The discrete equations of motion for the PCRTBP are derived by choosing an appropriate quadrature rule to discretize the Lagrangian in Equation 2. In this work, the trapezoidal approximation is applied. The trapezoid rule allows for an explicit second order accurate approximation. The discrete Lagrangian is given by

(22) |

Applying a discrete version of the Lagrange-d’Alembert principle allows for inclusion of an external control force on the system marsden2001 (). Using Sections 3.3 and 21 and some manipulation, the equations of motion are given by

(23a) | ||||

(23b) | ||||

(23c) | ||||

(23d) |

The discrete equations of motion are given in the Lagrangian form after applying the discrete fiber derivative from Equation 21 as and . The state is defined as and the control input is . The discrete potential gradients are given by

(24a) | ||||

(24b) | ||||

(24c) | ||||

(24d) |

The distances to each primary are defined as

(25a) | ||||

(25b) | ||||

(25c) | ||||

(25d) |

Care must be taken during the implementation of Equation 23. As Equations 25 and 24 are defined at both step and they must be evaluated at both time steps. Equation 23 is implemented by first defining an initial state and control . The distances and gravitational potential at step are evaluated from Equations 24b, 24a, 25b and 25a. The discrete update steps in Equations 23b and 23a are evaluated to generate and . Next, the distances and gravitational potential at step are evaluated from Equations 24d, 24c, 25d and 25c. Finally, the update steps in Equations 23d and 23c are evaluated. This results in the complete discrete update map given .

A simulation comparing the variational integrator to a conventional Runge-Kutta method is given in Figure 2. A particle is simulated from an initial condition of for years in the Earth-Moon system. The variational integrator uses a step size of while the Runge-Kutta method uses a variable step size implemented via ODE45 in Matlab. creftypecap 1(a) shows the /media/arxiv_projects/487592/trajectory of the spacecraft in the rotating reference frame for both integration schemes. Both integration schemes result in trajectories that are initially nearly identical. The discrete equations of motion are an accurate approximation for the continuous dynamics as they closely match the solution of ODE45 over the initial portion of the simulation. However, as time progresses the trajectories begin to diverge due to the differences in system /media/arxiv_projects/487592/energy. creftypecap 1(c) shows the evolution of the Jacobi integral. The variational integrator exhibits a bounded behavior about the initial /media/arxiv_projects/487592/energy with a mean variation of . However, the conventional Runge-Kutta method demonstrates a clear /media/arxiv_projects/487592/energy drift of . Over long simulation horizons or with the addition of small control inputs this poor /media/arxiv_projects/487592/energy behavior limits the applicability of conventional techniques.

## 4 Invariant Manifolds

There exists a rich structure of tubes, or invariant manifolds, in the three-body problem koon2000 (); conley1968 (). The manifold structure associated with periodic orbits about the and Lagrange points are critical to the understanding of the motion of spacecraft as well as comets/asteroids. In addition, the stable and unstable manifolds serve as the boundaries of the phase space region that allow for the transport between realms in a single three-body system or between multiple three-body systems. These invariant manifolds only exist as a result of the dynamic formulation of the three-body problem in a rotating reference frame. In this way, it is possible to design trajectories that intersect the invariant manifolds and connect widely separated regions of the state space. Figure 3 shows some of the tubes, projected from the phase space onto the configuration space, associated with and periodic orbits in the Earth-Moon system. This technique has been heavily investigated in Reference 2, applied to operational missions in Reference 12, and applied to potential multi-moon orbiter missions in Reference 13. In addition, much recent work has focused on potential return missions to the Moon zanzottera2012 (); campagnola2012 (); mingotti2011 (); ozimek2010a (); mingotti2009 ().

Another useful technique in the analysis of the free motion of dynamical systems is that of the Poincaré section. The Poincaré section is the intersection of a periodic orbit of a dynamical system with a lower dimensional sub-space transverse to the flow. This allows for greater insight into the stability and dynamics of periodic solutions of dynamic system. Points are drawn as the periodic solution intersects the Poincaré section. Great insight and structure into the dynamical system is typically available through use of Poincaré section. For example, Figure 4 shows the intersections of the invariant manifolds of Figure 3 with the Poincaré sections defined by the black line segments. The Poincaré section and the Jacobi integral reduce the state space from four to two dimensions for the planar three-body problem koon2001 (). This greatly eases the design process and allows for a simple planar visualization of the intersection of the invariant manifolds. Intersecting states are easy to determine and allow for transfers between invariant manifolds.

The intersection of the invariant manifolds on a Poincaré section is central to determining the desired control-free /media/arxiv_projects/487592/trajectory. However, the results previously developed are highly case specific and difficult to generalize to arbitrary transfers. Also, these results are based on control-free trajectories which rely on the underlying structure of the three-body system. In addition, transfer orbits along an invariant manifold require large time of flights which may be undesirable for time critical missions. The addition of low-thrust propulsion offers the potential of reduced transit times and the ability to depart from the free motion /media/arxiv_projects/487592/trajectory to allow for increased transfer opportunities. In this work, an optimal control problem is formulated to generate a reachable set of the spacecraft. The reachable set is computed on an appropriate Poincaré section and is used to design an appropriate transfer /media/arxiv_projects/487592/trajectory.

## 5 Optimal Control Formulation for Reachability Set

In this section, an optimal control formulation is presented to determine and design transfers within the three-body problem. The application of variational integrators to optimal control problems is referred to as computational geometric optimal control. This formulation is based on the concept of the reachability set on a Poincaré section. This method allows for one to easily determine potential transfer opportunities by finding set intersections on a lower dimensional space and greatly reduces the design process. The addition of continuous low thrust propulsion extends the control free design process presented in Reference 2 and allows for a greater range of potential transfers with a reduced time of flight.

Reachability theory provides a framework to evaluate control capability and safety. The reachable set contains all possible trajectories that are achievable over a fixed time horizon from a defined initial condition, subject to the operational constraints of the system. Reachability theory has been applied to aerospace systems such as collision avoidance, safety planning, and performance characterization. The theory formally supporting reachability has been extensively developed and is directly derivable from optimal control theory varaiya2000 (); lygeros2002 (); lygeros2004 (). Computation of the reachable set for a system involves solving the Hamilton-Jacobi partial differential equation or satisfying a dynamic programming principle. Analytical computation of reachable sets is an ongoing problem and is only possible for certain classes of systems. Typically, numerical methods are used to generate approximations of the reachability set, but are generally limited by the dimensionality of the problem.

Reachability theory has recently been applied to space systems holzinger2009 (); komendera2012a (). Computation of reachable sets is critical to space situational awareness, rendezvous and proximity operations, and orbit determination operations. Specifically maintaining accurate estimates of a spacecraft state over extended periods is not trivial. The challenge is increased for multiple spacecraft operating in close proximity or when there are long periods of time between measurements. Coupling the ability for continuous low-thrust propulsion between measurements increases the measurement association complexity. Computing the reachability set given estimated states and control authorities allows one to better correlate subsequent measurements or determine sensor pointing regions in the event of a lost spacecraft.

The cost function is defined as

(26) |

The term is the final state of a control-free /media/arxiv_projects/487592/trajectory while the term is the final state under the influence of the control input. In this fashion, the aim is to maximize the distance of the final state from that of the control-free /media/arxiv_projects/487592/trajectory. A chosen Poincaré section is defined through the use of appropriate terminal constraints given by

(27a) | ||||

(27b) | ||||

(27c) |

where the angles and define the Poincaré section and a specific direction upon the section, respectively. The goal is to determine the control input such that the cost function Equation 26 is minimized subject to the state equations of motion Equation 23 and constraints Equation 27. Application of the Euler-Lagrange equations allows us to derive the necessary conditions for optimality bryson1975 (). The discrete variational integrator in Equation 23 is used rather than the continuous time counterpart. This results in a discrete optimal control problem and the discrete necessary conditions are given as

(28a) | ||||

(28b) | ||||

(28c) |

where the Hamiltonian is defined as

(29) |

and is the costate and are the additional Lagrange multipliers associated with the terminal constraints in Equation 27. This indirect optimal control formulation leads to a two point boundary value problem with split boundary conditions. By sweeping the angle one can approximate the reachable set on the Poincaré section subject to the bounded control input.

The costate equation of motion requires the Jacobian of Equation 23 and is given by

(30) |

The derivation of Equation 30 is given in Appendix A. In addition, the computation of Equation 30 requires inversion of the Jacobian matrix. This is a computationally expensive operation that is prone to numerical error and instability. A method is presented in Appendix B to avoid this inversion and determine an explicit update map .

The optimal control formulation presented in this section results in a two point boundary value problem (TPBVP). There exist many methods to solve TPBVPs such as gradient, quasilinearization, and shooting methods bryson1975 (). In this work, a multiple shooting method is implemented. Shooting methods are common in astrodynamic /media/arxiv_projects/487592/trajectory design problems and relatively simple to implement. In the shooting method, initial conditions are varied such that a terminal constraint is satisfied, similar to the way an archer modifies the bow in order to more accurately strike a target. Rather than numerical integration over the entire time interval, multiple shooting segments the interval into several smaller sub-arcs. Additional interior constraints are used to ensure a continuous solution. This has the benefit of decreasing the numerical sensitivity of the final states to changes in the initial conditions along each sub-arc.

## 6 Numerical Example

A numerical simulation is presented to demonstrate the transfer procedure. The goal is to design a transfer /media/arxiv_projects/487592/trajectory from a planar periodic orbit about the Lagrange point to a bounded orbit in the vicinity of the Moon. The target region is created by choosing an initial condition of with . The initial state is propagated over a period of in non-dimensional units which corresponds to approximately years.

Figure 5 shows that the desired trajectories remain in the vicinity of the Moon in the rotating reference frame. This type of orbit would be useful for a variety of mission scenarios. For example, a series of communication satellites could be placed in this type of orbit. The bounded trajectories of the vehicles and constant line of sight to both the Moon and the Earth would allow for constant communication for future manned missions and potential habitats.

As a baseline, the method introduced in Reference 2 is implemented. The method is based on the invariant manifolds associated with the periodic orbits of the three body system. These invariant manifolds are a set of trajectories that either asymptotically arrive or depart the periodic orbit. Linking the invariant manifolds of several periodic orbits, or coupled three body systems, has been used to generate orbital transfers. creftypecap 5(b) shows the unstable invariant manifold associated with the periodic orbit. The invariant manifolds are globalized using the eigenvectors of the Monodromy matrix along the periodic orbit. The eigenvectors serve as a local approximation of the stable and unstable directions. Perturbations along these directions serve to approximate the invariant manifolds.

The unstable invariant manifold is numerically propagated to the same Poincaré section defined on the axis. The intersection of the invariant manifold and the Poincaré section is denoted by the green markers in creftype 5(c). Only a single branch of the invariant manifold intersects with the ascending region of the target orbit. There are no intersections of the invariant manifold with the descending region of the target orbit. The numerical values associated with the green points denote the required time of flight along the invariant manifold in non-dimensional units. A transfer along the invariant manifold requires on average as compared to for a transfer using low thrust propulsion and the reachable set. Finally, it should be noted that an additional maneuver would be required for a transfer using the invariant manifold. The intersection on the Poincaré section only shows that the /media/arxiv_projects/487592/components intersect. An additional instantaneous would be required to transfer from the invariant manifold to the target.

A periodic orbit is generated about the Lagrange point. A Poincaré section is chosen to allow for design on a lower dimensional subspace. The section is defined along the axis, or defined with and intersects both the initial and target orbits. From the periodic orbit, a series of optimal trajectories are generated to approximate the reachable set. The trajectories are generated from a fixed initial state of over a fixed time span of . By varying the angle in Equation 27, a different direction along the Poincaré section is maximized. The intersection of the optimal trajectories as well as those of the target Moon orbit with the Poincaré section are shown in Figure 6.

The optimal trajectories, under the influence of the control input , are shown in creftype 5(a). Initially, the spacecraft is assumed to lie on the periodic orbit. As a result, the intersection of this periodic orbit with the Poincaré section are two points corresponding to the two crossing of the orbit. The use of the continuous low thrust propulsion expands the reachable set to region bounded by the red markers in creftype 5(c). The reachable set is an ellipsoidal region with a major axis aligned along .

The blue points in creftype 5(c) are the intersections of the target Moon orbit and the Poincaré section. The two circular regions are the ascending (right) and descending (left) intersections of the target orbit and Poincaré section. creftypecap 5(c) shows that the reachable set and those of the descending target region intersect. As both regions are discretely approximated a linear interpolation is used to determine the exact intersection state on the Poincaré section. This intersection generates a partial target state of . Using Equation 10 and the intersection state the final component is calculated. This results in a complete target state which lies on the reachable set and on the target orbit. A final optimal /media/arxiv_projects/487592/trajectory is generated such that the . This transfer /media/arxiv_projects/487592/trajectory is denoted by the green path in creftype 5(a).

## 7 Conclusions

In this paper, an optimal transfer process which combines concepts of reachability and Poincaré section is used to generate transfer between planar periodic orbits in the three-body problem. The Poincaré section allows for /media/arxiv_projects/487592/trajectory design on a lower dimensional phase space and simplifies the process. The indirect optimal control formulation enables straightforward method of incorporating additional path and control constraints. However, the use of optimal control techniques leads to open loop trajectories that are not robust to model uncertainties or disturbances. Lyapunov control theory, which has previously been applied to the two-body problem, is being investigated in the hope of designing closed loop control schemes for this three-body scenario chang2002 (). This analysis has also assumed perfect attitude control and the ability to orient thrust in any direction. The addition of attitude dynamics and realistic pointing constraints would significantly improve the applicability.

## References

- (1) S. E. Haque, M. Keidar, and T. Lee, “Low-Thrust Orbital Maneuver Analysis for Cubesat Spacecraft with a Micro-Cathode Arc Thruster Subsystem,” Proceedings of the thirty-third international electric propulsion conference, Electric Rocket Propulsion Society, Washington DC, USA, 2013.
- (2) W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, Dynamical Systems, the Three-Body Problem and Space Mission Design. World Scientific, 2000.
- (3) S. D. Ross, “The interplanetary transport network,” American scientist, Vol. 94, No. 3, 2006, p. 230.
- (4) G. Gómez, W. Koon, M. W. Lo, J. Marsden, J. Masdemont, and S. Ross, “Invariant Manifolds, the Spatial Three-Body Problem and Space Mission Design,” AAS/AIAA Astrodynamics Specialist Conference, Quebec City, Canada, 2001, American Astronautical Society, August 2001.
- (5) G. Mingotti, F. Topputo, and F. Bernelli-Zazzera, “Earth–Mars transfers with ballistic escape and low-thrust capture,” Celestial Mechanics and Dynamical Astronomy, Vol. 110, No. 2, 2011, pp. 169–188, 10.1007/s10569-011-9343-5.
- (6) D. J. Grebow, M. T. Ozimek, and K. C. Howell, “Design of Optimal Low-Thrust Lunar Pole-Sitter Missions,” The Journal of the Astronautical Sciences, Vol. 58, No. 1, 2011, pp. 55–79, 10.1007/BF03321159.
- (7) C. Lanczos, The Variational Principles of Mechanics, Vol. 4. Courier Corporation, 1970.
- (8) V. Szebehely, “Theory of orbits. The Restricted Problem of Three Bodies,” New York: Academic Press, Vol. 1, 1967.
- (9) J. E. Marsden and M. West, “Discrete Mechanics and Variational Integrators,” Acta Numerica 2001, Vol. 10, 2001, pp. 357–514.
- (10) D. T. Greenwood, Principles of Dynamics. Prentice-Hall Upper Saddle River, NJ, 1988.
- (11) C. C. Conley, “Low Energy Transit Orbits in the Restricted Three-Body Problem,” SIAM Journal on Applied Mathematics, Vol. 16, No. 4, 1968, pp. pp. 732–746.
- (12) W. S. Koon, M. W. Lo, J. Marsden, and S. Ross, “The Genesis Trajectory and Heteroclinic Cycles,” Astrodynamics 1999, Vol. 103, No. Part III, 1999, pp. 2327–2343.
- (13) K. Tanaka and J. Kawaguchi, “Low-Thrust Transfer between Jovian Moons using Manifolds,” Spaceflight Mechanics, Vol. 140, No. Part 3, 2011, pp. 1935–1942.
- (14) A. Zanzottera, G. Mingotti, R. Castelli, and M. Dellnitz, “Intersecting invariant manifolds in spatial restricted three-body problems: Design and optimization of Earth-to-halo transfers in the Sun–Earth–Moon scenario,” Communications in Nonlinear Science and Numerical Simulation, Vol. 17, 2 2012, pp. 832–843, http://dx.doi.org/10.1016/j.cnsns.2011.06.032.
- (15) S. Campagnola, P. Skerritt, and R. Russell, “Flybys in the planar, circular, restricted, three-body problem,” Vol. 113, No. 3, 2012, pp. 343–368, 10.1007/s10569-012-9427-x.
- (16) M. T. Ozimek and K. C. Howell, “Low-Thrust Transfers in the Earth-Moon System, Including Applications to Libration Point Orbits,” Journal of Guidance, Control, and Dynamics, Vol. 33, 2014/08/16 2010, pp. 533–549, 10.2514/1.43179.
- (17) G. Mingotti, F. Topputo, and F. Bernelli-Zazzera, “Low-/media/arxiv_projects/487592/energy, low-thrust transfers to the Moon,” Celestial Mechanics and Dynamical Astronomy, Vol. 105, No. 1-3, 2009, pp. 61–74.
- (18) W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, “Low /media/arxiv_projects/487592/energy transfer to the Moon,” Dynamics of Natural and Artificial Celestial Bodies, pp. 63–73, Springer, 2001.
- (19) P. Varaiya, “Reach set computation using optimal control,” Verification of Digital and Hybrid Systems, pp. 323–331, Springer, 2000.
- (20) J. Lygeros, “On the Relation of Reachability to Minimum Cost Optimal Control,” Decision and Control, 2002, Proceedings of the 41st IEEE Conference on, Vol. 2, Dec 2002, pp. 1910–1915 vol.2, 10.1109/CDC.2002.1184805.
- (21) J. Lygeros, “On reachability and minimum cost optimal control,” Automatica, Vol. 40, 6 2004, pp. 917–927, http://dx.doi.org/10.1016/j.automatica.2004.01.012.
- (22) M. Holzinger and D. Scheeres, “Reachability Analysis Applied to Space Situational Awareness,” Advanced Maui Optical and Space Surveillance Technologies Conference, 2009.
- (23) E. E. Komendera, D. J. Scheeres, and E. Bradley, “Intelligent Computation of Reachability Sets for Space Missions.,” IAAI, 2012.
- (24) A. E. Bryson and Y. C. Ho, Applied Optimal Control: Optimization, Estimation and Control. CRC Press, 1975.
- (25) D. E. Chang, D. F. Chichka, and J. E. Marsden, “Lyapunov-based Transfer between Elliptic Keplerian Orbits,” Discrete and Continuous Dynamical Systems Series B, Vol. 2, No. 1, 2002, pp. 57–68.

## Appendix A: Costate Equations of Motion

The development of the costate equations of motions begins with determining the second order partial derivatives of the gravitational potential. Due to the symmetry of partial derivatives only three terms are required and are given by

(31) | ||||

(32) | ||||

(33) |

The gradient of Equation 23a is given as

(34a) | ||||

(34b) | ||||

(34c) | ||||

(34d) |

The gradient of Equation 23b is given as

(35a) | ||||

(35b) | ||||

(35c) | ||||

(35d) |

The gradients of Equations 25d and 25c are given as

(36a) | ||||

(36b) |

The second order partial derivatives of the gravitational potential at are given as

(37a) | ||||

(37b) |

The gradient of Equations 23d and 23c are given as

(38a) | ||||

(38b) | ||||

(38c) | ||||

(38d) |

(39a) | ||||

(39b) | ||||

(39c) | ||||

(39d) |

These gradient equations are in a cascade type structure. Equations 38 and 39 are functions of Equations 23d and 23c. As a result, the accuracy of the Jacobian will tend to decrease as the first order approximation errors accumulate.

## Appendix B: Gauss Jordan Elimination

The costate equations of motion are given by Equation 30 and repeated here as

(40) |

To determine the discrete update map the inverse of the Jacobian matrix is required. In order to avoid the need of an explicit inversion a Gauss Jordan method is implemented. To begin, several terms are defined which are required to carry out the row operations and are defined as

(41a) | ||||

(41b) | ||||

(41c) | ||||

(41d) | ||||

(41e) | ||||

(41f) |

Equation 40 is transformed to row echelon form using elementary row operations and is defined as

(42) |

where the terms and are defined as follows

(43a) | ||||

(43b) | ||||

(43c) | ||||

(43d) | ||||

(43e) | ||||

(43f) | ||||

(43g) | ||||