A LOW-THRUST TRANSFER STRATEGY TO EARTH-MOON COLLINEAR LIBRATION POINT ORBITS
A Thesis Submitted to the Faculty of Purdue University by Martin T. Ozimek
In Partial Fulfillment of the Requirements for the Degree of Master of Science
December 2006 Purdue University West Lafayette, IN
ii
For Mom, Dad, Sarah, and Bops
iii
ACKNOWLEDGMENTS
I would like to thank my parents for their seemingly never ending support and confidence that I would persevere in my personal quest for an advanced degree. The decision to commit to a higher degree has been an adventure that I hope to continue along, and I can’t begin to explain the importance of that priceless feeling of simply knowing that someone is there when needed. Professor Kathleen Howell, my advisor, is also owed a great deal of gratitude, not only for posing the fateful words “low-thrust” one day in her office during a discussion about NASA’s potential Jupiter Icy Moons missions, but also for her personal standard for excellence that she instills in each of her many successful graduate students. I’ve always felt that Purdue University reached out to me and offered me that “extra” indefinable something from the moment I began seriously considering a graduate institution. In no other person is this ambiguous something “extra” exemplified than in Professor Howell, who has sought to ensure that my research efforts are guided and ultimately shared with others in the best possible way. This notion has also been exemplified by Professor James Longuski, whose door has always been open to me from day one, and whom I must also credit in heavily influencing my decision to attend Purdue University. On more than one occasion, Professor William Crossley has also had his door open to help along my path of solving what turned out to be a difficult optimization problem. I also owe many thanks to Daniel Grebow. Dan has been a close friend, colleague, and even roommate throughout my stay at Purdue, and this work is the direct continuation of a mission analysis that we worked on together. Often, many of the new ideas I have for current and future research are a result of simple dialogs that we frequently engage in.
iv
The idea to study mission applications toward lunar south pole coverage would never have originated had I not fortuitously been privileged to work at the NASA Goddard Spaceflight Center during the summers of 2005 (as a member of the NASA Academy by way of the Indiana Space Grant Consortium) and 2006. There, I benefited from the knowledge of some of the greatest libration point mission experts in the world, and owe particular thanks to my mentor David Folta.
Support from NASA under contract
numbers NNG05GM76G and NNX06AC22G is greatly appreciated. Finally, I would like to thank Purdue University for financial support, including the Andrews Fellowship, for the entirety of my M.S. program.
v
TABLE OF CONTENTS
Page LIST OF TABLES............................................................................................................ vii LIST OF FIGURES ......................................................................................................... viii ABSTRACT........................................................................................................................ x 1 INTRODUCTION ............................................................................................................1 1.1 Historical Overview of the Three-Body Problem.......................................................3 1.2 Developments in Low-Thrust Transfer Trajectories ..................................................5 1.2.1 Optimal Control .................................................................................................. 5 1.2.2 Application to Orbit Problems ............................................................................ 6 1.3 Focus of this Work .....................................................................................................7 2 BACKGROUND ............................................................................................................10 2.1 The Circular Restricted Three-Body Problem..........................................................10 2.1.1 Assumptions...................................................................................................... 11 2.1.2 Geometry........................................................................................................... 11 2.1.3 Equations of Motion.......................................................................................... 12 2.1.4 Libration Points................................................................................................. 15 2.1.5 Formulation Relative to P2 in the CR3BP......................................................... 17 2.2 Natural Periodic Orbits in the CR3BP......................................................................18 2.2.1 First-Order Variational Equations Relative to the Collinear Points.................. 19 2.2.2 The State Transition Matrix .............................................................................. 23 2.2.3 The Fundamental Targeting Relationships ....................................................... 25 2.2.4 Periodic Orbits .................................................................................................. 28 2.3 Invariant Manifolds ..................................................................................................31 2.3.1 Stable and Unstable Manifolds Associated with the Collinear Points .............. 32
vi
Page 2.3.2 Invariant Manifolds Relative to a Fixed Point .................................................. 35 2.3.3 Computation of Manifolds Corresponding to Fixed Points Along an Orbit..... 37 2.4 Optimal Control Theory ...........................................................................................39 2.4.1 Summary of the First Necessary Conditions for Optimal Control.................... 40 2.4.2 Tests for a Local Minimum Value of the Performance Index .......................... 43 3 LOW-THRUST TRANSFER ALGORITHM ................................................................45 3.1 Engine Model ...........................................................................................................46 3.2 Control Law Derivation............................................................................................48 3.3 Adjoint Control Transformation...............................................................................55 3.4 Numerical Solution via Direct Shooting: A Local Approach...................................59 3.5 Shotgun Method for Initial Conditions: A Global Approach ...................................63 4 MISSION APPLICATIONS...........................................................................................66 4.1 Orbits for Line-of-Sight Lunar South Pole Coverage (CR3BP) ..............................67 4.1.1 Three-Dimensional Periodic Orbits in the CR3BP ........................................... 67 4.1.2 Families of Orbits for Lunar South Pole Coverage ........................................ 68 4.1.3 Mission Orbit Selection Criteria ....................................................................... 73 4.2 Optimal Transfers to the Earth-Moon Stable Manifold............................................75 4.2.1 Transfers to a 12-Day L1 Halo Orbit ................................................................. 78 4.2.2 Transfer to a 14-Day L1 Vertical Orbit ............................................................. 87 4.2.3 Transfer to a 14-Day L2 Butterfly Orbit ............................................................ 91 5 SUMMARY AND RECOMMENDATIONS.................................................................96 5.1 Summary...................................................................................................................96 5.2 Recommendations for Future Work .........................................................................98 LIST OF REFERENCES...................................................................................................99
vii
LIST OF TABLES
Table
Page
Table 4.1 Dynamical and Propulsion Constants ................................................................77 Table 4.2 12-Day L1 Halo Orbit Transfer Data Summary .................................................82 Table 4.3 12-Day L1 Halo Orbit Long Transfer Data Summary........................................86 Table 4.4 14-Day L1 Vertical Orbit Transfer Data Summary............................................91 Table 4.5 14-Day L2 Butterfly Orbit Transfer Data Summary ..........................................95
viii
LIST OF FIGURES
Figure
Page
Figure 2.1 Geometry in the Restricted Three-Body Problem ............................................12 Figure 2.2 Equilibrium Point Locations for the CR3BP....................................................16 Figure 2.3 Geometry of P2-Centered Rotating Frame .......................................................18 Figure 2.4 Linearized L1 Periodic Orbit.............................................................................23 Figure 2.5 Basic Diagram for a Free Final Time Targeting Scheme.................................27 Figure 2.6 Targeting a Perpendicular X-axis Crossing in the CR3BP ..............................30 Figure 2.7 Several L1 Lyapunov Orbits Obtained Via Continuation .................................31 Figure 2.8 Stable and Unstable Manifold at X eq ...............................................................34 Figure 2.9 Global Manifolds for an Earth-Moon L1 Lyapunov Orbit, Ay = 23,700 km.....38 Figure 3.1 CSI Engine Example – Smart-1 Ion Engine.....................................................46 Figure 3.2 VSI Engine Example - VaSIMR Rocket ..........................................................47 Figure 3.3 Behavior of θΜ and τM Along the Stable Manifold Tube.................................51 Figure 3.4 Velocity Reference Frame................................................................................58 Figure 3.5 Numerical Algorithm – Direct Shooting Method via SQP ..............................62 Figure 3.6 Example Population Parameters for a 12-Day Halo Orbit ...............................65 Figure 4.1 Southern Halo Orbit Families: Earth-Moon L1 (Orange) and L2 (Blue); Moon Centered, Rotating Reference Frame ................................................................70 Figure 4.2 Vertical Orbit Family of Interest: Earth-Moon L1 (Magenta) and L2 (Cyan); Moon Centered, Rotating Reference Frame......................................................71 Figure 4.3 Southern L2 Butterfly Orbit Family; Moon Centered, Rotating Reference Frame.................................................................................................................72 Figure 4.4 Period versus Maximum x-Distance from the Moon (Left); Definition of Maximum x-Distance (Right) ...........................................................................74
ix
Figure
Page
Figure 4.5 Stability Index versus Maximum x-Distance from the Moon ..........................75 Figure 4.6 Optimal Orbit Raising from LEO.....................................................................77 Figure 4.7 Stable Manifold Tube for 12-Day L1 Halo Orbit (Green) and Target Reference Trajectory Along the Manifold (Blue) ..............................................................79 Figure 4.8 Low-Thrust Short Transfer to a 12-Day L1 Halo Orbit ....................................80 Figure 4.9 Position and Velocity Costate Time Histories for the 12-Day L1 Halo Orbit Transfer .............................................................................................................81 Figure 4.10 Time History of Propulsion Related Parameters for the 12-Day L1 Halo Orbit Transfer .............................................................................................................82 Figure 4.11 Stable Manifold Tube for 12-Day L1 Halo (Green) and Initial Target Reference Trajectory Along the Manifold (Blue) for Long Transfer ...............83 Figure 4.12 Low-Thrust Long Transfer to a 12-Day L1 Halo Orbit ..................................84 Figure 4.13 Position and Velocity Costate Time Histories for the 12-Day L1 Halo Orbit Transfer .............................................................................................................85 Figure 4.14 Time History of Propulsion Related Parameters for the 12-Day L1 Halo Orbit Long Transfer....................................................................................................86 Figure 4.15 Stable Manifold Tube for 14-Day L1 Vertical Orbit (Green) and Initial Target Reference Trajectory Along the Manifold (Blue).............................................88 Figure 4.16 Low-Thrust Transfer to a 14-Day L1 Vertical Orbit.......................................89 Figure 4.17 Position and Velocity Costate Time Histories for the 14-Day L1 Vertical Orbit Transfer....................................................................................................90 Figure 4.18 Time History of Propulsion Related Parameters for the ................................90 Figure 4.19 Stable Manifold Tube for 14-Day L2 Butterfly Orbit (Green) and Initial Target Reference Trajectory Along the Manifold (Blue) .................................92 Figure 4.20 Low-Thrust Transfer to a 14-Day L2 Butterfly Orbit .....................................93 Figure 4.21 Position and Velocity Costate Time Histories for the 14-Day L2 Butterfly Orbit Transfer....................................................................................................94 Figure 4.22 Time History of Propulsion Related Parameters for the 14-Day L2 Butterfly Orbit Transfer....................................................................................................94
x
ABSTRACT
Ozimek, Martin T. M.S.A.A., Purdue University, December, 2006. A Low-Thrust Transfer Strategy to Earth-Moon Collinear Libration Point Orbits. Major Professor: Kathleen Howell.
A strategy to compute low-thrust transfer trajectories in the Earth-moon circular restricted three-body problem is developed.
The dynamical model is formulated
assuming variable specific impulse engines, an advanced finite-thrust propulsion model. Originating in an Earth parking orbit, the spacecraft is delivered to a location along the stable manifold; the engines power off and the spacecraft asymptotically approaches the periodic libration point orbit of interest. Elements of optimal control theory are used to derive a primer vector control law as well as a set of additional dependent variables that characterize the solution to the corresponding two-point boundary-value problem (TPBVP).
A hybrid direct/indirect method results in transfer trajectories that are
associated with locally minimal propellant consumption. The generation of useful initial conditions is aided by an adjoint control transformation and a global “shotgun” method. The solution strategy is demonstrated in a detailed development of transfers to a 12-day L1 halo orbit, a 14-day L1 vertical orbit, and a 14-day L2 “butterfly” orbit. These target orbits are all selected from the families that meet line-of-sight coverage requirements in support of lunar south pole mission architecture.
1
1. INTRODUCTION
As understanding of the solar system and the dynamical structure of the space environment increases, ever more complex questions continue to emerge. In the past 60 years, access to space has opened to both human-crewed and robotic spacecraft. Both scientific interest and engineering capability have invariably played a vital role in this expansion of knowledge. New scientific demands often spur engineering advancements to accomplish a set of mission objectives; breakthroughs in engineering capability enlighten the scientific community and serve as a catalyst for original mission concepts. Perhaps, as a consequence, it is not surprising that libration point orbits have relatively recently risen as venues for robotic spaceflight. Beginning with NASA’s 1978 solar wind measuring satellite, the International Sun-Earth Explorer (ISEE-3) [1], libration point orbits are now considered viable options to meet a range of scientific goals. Although originally proposed for manned Apollo missions [1], such orbits were not exploited prior to ISEE-3. However, with the ever increasing speed of computers, trajectory design within the context of the n-body problem is now feasible. Although the n-body problem is unsolvable in closed form, certain simplifying assumptions expose equilibrium solutions, i.e., the libration points. Periodic orbits in their vicinity can be computed numerically. Not coincidentally, the success of the early missions like ISEE-3 and the continuing increase in computational capabilities, have led to more contemporary libration point missions such as WIND [2], SOHO [3], ACE [4], MAP [5], and Genesis [6]. Clearly, such trajectory designs fill a particular niche in mission applications where long-duration scientific observation is required. More recently, a geometrical approach in studies of the multi-body problem has also led to alternative strategies for transfer trajectory design, as well as stationkeeping maneuvers. This approach is based on a complete analysis of the phase space in the
2
neighborhood of the periodic libration point orbits. In-depth analysis of the phase space, as suggested by Poincaré [7] in 1892, has evolved into dynamical systems theory (DST). An important astrodynamics application from DST is the exploitation of invariant manifolds to design trajectory arcs that asymptotically arrive at, and depart from, the vicinity of the periodic libration point orbits. Propagation of these manifolds often yields very efficient transfers. In some cases, these manifolds even pass within the vicinity of a planet [6]. Typically, such transfer arcs are more applicable to robotic spaceflight, since an additional time-penalty is often incurred. This DST approach was a key component in the design of the Genesis low-energy trajectory.
The Genesis trajectory design
incorporated heteroclinic and homoclinic arcs to deliver a spacecraft to a Sun-Earth L1 libration point orbit with a subsequent return to Earth [6]. Such manifold transfer trajectories are also useful for applications in the Earth-moon system. Recent studies have identified libration point orbits as a potential component in the development of a communications relay between a manned facility at the lunar south pole and ground stations on the Earth. In the Earth-moon problem, however, many of the viable libration point orbits possess manifolds that pass no closer than 50,000 km to the Earth. These manifolds may still serve as transfers, but additional fuel is necessary to incorporate a leg from an Earth parking orbit to the manifold. One type of intermediate transfer from the parking orbit to the manifold involves the use of low-thrust propulsion. Low-thrust propulsion introduces a time penalty, but can yield lower fuel expenditure due to higher specific impulse engines. Incorporating lowthrust also adds dynamical complexity because a “steering” law for the thrust direction and magnitude must be determined and successfully implemented. The objective of the current work is a method to design low-thrust transfers that deliver a vehicle onto a manifold trajectory. The specific application of interest supports the establishment of lunar south pole communications relay infrastructure.
3
1.1 Historical Overview of the Three-Body Problem
The general problem of three bodies was first investigated by Isaac Newton in his 1687 landmark work, the Principia [8]. His successor, Leonhard Euler, receives much of the credit, however, for formulation of the restricted problem of three bodies. In 1765, Euler identified the equilibrium solutions in the restricted problem, i.e., the collinear libration points L1, L2, and L3 [9]. Euler also introduced a synodic reference frame in connection with the motion of the moon in 1772 [9]. Later, in 1772, the same year that Euler formulated the restricted problem, Lagrange determined the locations of two additional equilibrium points, the triangular libration points, L4 and L5 [9]. Euler’s formulation allows only one integral of motion in the restricted problem as determined by Jacobi in 1836, by balancing energy and angular momentum [10]. In 1878, George William Hill published his Researches in Lunar Theory [11], effectively modeling the motion of the moon as a satellite, exposed to the gravitational field of the Earth and the perturbing force of the Sun. In 1899, Henri Poincaré published his three-volume work, Les Méthodes Nouvelles de la Mécanique Celeste [7], the result of his unparalleled response to a contest in 1887. Participants were challenged to produce a definitive solution to the n-body problem. Ironically, Poincaré eventually won the prize by proving that the n-body problem cannot be solved in closed form. His work is highly regarded for the qualitative emphasis on behaviors in the n-body problem. In particular, Poincaré focused on the behavior of trajectories as time goes to infinity. Of course, the only trajectory that can be defined at infinite time is a periodic orbit. A detailed analysis of the phase space of a non-integrable system led Poincaré to invent a technique known as the “surface of section”. This work is considered the foundation of dynamical systems theory. In addition, Poincaré also proved that Jacobi’s Constant is the only integral of the motion in the restricted problem. Poincaré’s work generated great interest in the following decades. Periodic orbits are a common topic noted in the early 20th century work of Darwin [12], Plummer [13], and Moulton [14]. In the absence of extensive computational capabilities, these researchers exploited expansion procedures to construct analytical approximations.
Darwin and
4
Plummer are noted for approximating planar periodic orbits in the restricted problem. Beyond a general focus on planar periodic orbits, Moulton also studied in-plane and outof-plane orbits in the vicinity of the collinear libration points. The renewed attention on the restricted problem in the latter half of the 20th century is attributed primarily to the emergence of high-speed computing and the beginning of the space age. Several developments are notable. Szebehely’s 1967 book The Theory of Orbits [10], a comprehensive text on the three-body problem, is still regarded as one of the most thorough sources of information on the problem of three bodies. By the 1960’s, in support of the Apollo missions to the moon, trajectories computed in the restricted three-body problem were in development. Farquhar coined the term “halo” orbits and developed analytical approximations of these three-dimensional periodic orbits in the vicinity of L2 in the Earth-moon system. In addition to proposing orbits for use in the Apollo missions and ultimately for the unmanned ISEE-3 [1] spacecraft in the Sun-Earth system, Farquhar and Kamel [15] also developed third-order approximations for quasiperiodic orbits. Richardson and Cary [16] expanded these approximations to fourth order. Breakwell and Brown [17] are noted for a numerical study that generates families of periodic halo orbits. Howell [18] extended the analysis to include all collinear points and the families over all three-body systems. Hénon [19] has also produced thorough analyses of periodic orbits, including vertical orbits.
Approximations for nearly
rectilinear halo orbits at the collinear points were produced by Howell and Breakwell [20].
Perhaps the most rigorous numerical generation of periodic orbits is the
investigation by Dichmann, Doedel, and Paffenroth [21], that uses the software package AUTO to detail periodic orbits as well as their interrelated orbits via bifurcation.
5
1.2 Developments in Low-Thrust Transfer Trajectories
1.2.1 Optimal Control
The development of finite-burn trajectories (in particular, when the thrust level is low) is an application of optimal control theory and the calculus of variations. Optimization of curves and points can be traced to the 1600’s, when the calculus of variations was first introduced as an analysis tool for minimizing functions of functions, or “functionals”. It received particular attention with Johann Bernoulli’s proposal of the brachistochrone problem to the scientific community 1696 [22].
Generally, the focus is the set of
conditions on the functions that drive the functional to a maximum or minimum. Such functions are termed “extremals”. In 1755, Joseph Lagrange wrote a letter to Leonhard Euler in connection with their mutual interest in an analytical solution to the tautochrone problem [23]. The analytical development resulted in the formulation of the EulerLagrange Equations (and the corresponding transversality condition). These equations have since served as the basis of a widely known technique for determining extremals, and a foundation of the calculus of variations, a term created by Euler in 1766. Classification of these extremals is accomplished via the Legendre-Clebsch necessary condition, the Weierstrass Condition, and in the most general of terms, Pontryagin’s Minimum Principle. Details of the generalized theory in support of the applications of this methodology appear in Bryson and Ho [24], Hull [25], and Kirk [26]. Low levels of thrust in the computation of finite-burn spacecraft trajectories is achieved by the seemingly parallel availability of high-speed computing and the continuing development of advanced propulsion systems. In his 1963 book, Optimal Spacecraft Trajectories [27], Lawden used primer vector theory to outline a general procedure for determining optimal low-thrust trajectories. Primer vector theory blends a control law and switching structure common in many “indirect” optimization methods. Indirect, low-thrust, trajectory optimization methods are typically characterized by the two-point boundary-value problem formulation from optimal control theory, with a continuous parameterization of the thrust direction via the tangent to the primer vector.
6
Such approaches are termed “indirect” because once the two-point boundary-value problem is established, no minimization is required on the cost functional directly (although it can be included [28-30]). Alternatively, the solution involves a root-solving process on the kinematical as well as the natural boundary conditions as a result of the Euler-Lagrange equations, transversality condition, and a corresponding secondary test for a minimum. Conversely, “direct” methods involve an attempt to minimize the cost functional itself, and often use many variables to parameterize the thrust magnitude and direction. A common parameterization of the thrust direction is the development of a spline function [29,31]. While indirect methods typically require more precise initial conditions due to numerical sensitivities, the lower dimension on the search vector implies fewer computations. Marec [32] offers further mathematical analysis of the primer vector; examines high- and low-thrust propulsion systems and details the Contenson-Pontryagin Maximum Principle. Early applications to impulsive rendezvous problems are available in Jezewski [33] and to low-thrust rendezvous in Melbourne and Sauer [34]. Further increases in computing speed have allowed more sophisticated methods in numerical solutions to optimization problems. All trajectory optimization problems are typically solved with the following numerical methods:
direct shooting, indirect
shooting, multiple shooting, direct transcription, indirect transcription, dynamic programming, or genetic algorithms. For a detailed survey of all of these different methods, see Betts [35]. Thus, optimal control theory serves to set up the conditions and constraints that must be met to determine the existence of an optimal control, and the numerical optimization methods serve as the means to actually compute solutions that meet these exact conditions.
1.2.2 Application to Orbit Problems
A number of investigations are notable in examining optimal low-thrust trajectories. Many applications are formulated in the two-body problem [36-40]. In the three-body problem, there has been less attention. In one example, Herman and Conway [41]
7
compute optimal, low-thrust, Earth-moon transfers using equinoctial elements, with a direct collocation solution method. Kluever [29-30] also develops Earth-moon transfers, but utilizes a hybrid direct/indirect method. Golan and Breakwell [42] use matching of two trajectory segments to transfer into lunar orbit; transfers to L4 and L5 are also presented. In a Hill formulation of the three-body problem, Sukhanov and Eismont [43] establish a primer vector control law to develop a three-dimensional transfer into a SunEarth L1 halo orbit that includes a constrained thrust direction. Although Seywald, Roithmeyer, and Troutman [44] do not employ the circular restricted three-body equations of motion, they are notable for attempting to provide approximate analytical solutions of circle-to-circle orbit transfers with a variable specific impulse engine, and, furthermore, compare the engine model with results using constant specific impulse engines. Primer vector theory is employed by Russell [45] to develop a global search and local optimization method and establish a Pareto front on his resulting fixed time solutions. Russell subsequently applies his method to produce transfers to planar Earthmoon and Sun-Earth distant retrograde orbits. Senent, Ocampo and Capella [28] also use primer vector theory for free final time transfers to Sun-Earth libration point orbits via the stable manifold.
1.3 Focus of this Work
Optimal, free final time low-thrust transfer trajectory profiles to several Earth-moon libration point orbits is the objective of this work. The target orbits are selected based on lunar south pole coverage applications, including halo orbits, vertical orbits, and “butterfly” orbits. These transfer trajectories represent an extension of the exhaustive coverage analysis by Grebow et al. [46].
8
The work is organized as follows. Chapter 2: All background material is presented in Chapter 2. The equations of motion that govern the nondimensional, barycentric, cartesian, circular restricted three-body problem are developed.
A moon-centered model is also introduced. The linear variational
equations that result from using the collinear libration points as a reference are introduced, and, thus, form the basis for generating initial conditions in the nonlinear problem. Then, an alternate set of linear variational equations is developed, where the time-varying orbital trajectory is exploited as a reference. The second set of variational equations is useful for iterative orbit targeting. Invariant manifold theory is summarized for both libration points and for fixed points along a periodic libration point orbit. Finally, the necessary conditions for establishing a stationary value of a generalized performance index are introduced. These include both the Euler-Lagrange equations and Pontryagin’s Minimum Principle. Such conditions ultimately yield the full two-point boundary-value problem that may be solved via nonlinear programming methods. Chapter 3: Optimal control theory is applied to develop the well-known primer vector control law parameterization. Variable specific impulse engines (e.g., those in development for the VaSIMR rocket project) are modeled and produce notable improvement in numerical convergence. The problem is formulated as a free final time transfer from an Earth parking orbit to the stable manifold associated with the libration point orbit of interest. Using the stable manifold allows a stopping condition, and a secondary coast phase; the spacecraft asymptotically converges to the desired orbit. Once the control law and the elements of the two-point boundary-value problem are established, a solution method is presented. This scheme incorporates a “global” search method, and a subsequent local sequential quadratic programming (SQP) algorithm. Rather than solve the complete twopoint boundary-value problem, the local approach is set up as a direct shooting method, and is termed a hybrid direct/indirect approach.
The adjoint control transformation
(ACT) is used to map the initial costates into more meaningful physical quantities that are useful in the global method and the first step of the local, direct shooting method.
9
Chapter 4: Once the methodology is established, the criteria for target orbit selection is summarized and the orbit families of interest are introduced. A halo orbit, vertical orbit, and a “butterfly” orbit are selected and the solution procedure is applied to generate lowthrust transfer trajectories. The resulting thrust profiles, trajectory plots, and overall performance parameters are then discussed. Chapter 5: Conclusions
concerning
the
solution
methodology
are
presented.
Then,
recommendations on future work, including different solution methods, equations of motion, and higher fidelity models are detailed. Finally, potential future applications as a result of this work are presented.
10
2.
BACKGROUND
Some fundamental mathematical tools and concepts are necessary for the development of the trajectories and transfers in the current application. The cartesian, nondimensional, barycentric Circular Restricted Three-Body Problem (CR3BP) is initially formulated, and the corresponding dynamical equations of motion are derived. The five equilibrium points are the basis for the computation of special periodic orbits. The state transition matrix is introduced as a tool for predicting linear, and approximating nonlinear, motion. An alternate formulation of the equations of motion, centered at the smaller primary, is also presented. The fundamental structure underlying the dynamics in the restricted problem is analyzed with invariant manifold theory, yielding natural pathways to and from a periodic orbit. Basic numerical methods for a targeting procedure are detailed, forming the groundwork for determining periodic orbits. Finally, the basic concepts of optimal control theory are introduced, including the necessary conditions along a trajectory for the existence of a locally optimal transfer trajectory.
2.1 The Circular Restricted Three-Body Problem
The rapid growth of high-speed computing in the last three decades has helped spark the discovery of new and exciting trajectories. Many of these trajectories exist within the context of the Circular Restricted Three-Body Problem (CR3BP). An exact analytical solution to the CR3BP does not exist; thus, any solutions beyond the equilibrium points require numerical integration.
Nevertheless, at the expense of additional numerical
exploration, propagation of trajectories in this model result in non-Keplerian orbital motion, such as “figure-eight” orbits, “halo” orbits, and an infinite variety of other periodic orbits; quasi-periodic trajectories have also been identified.
11
2.1.1 Assumptions
Given an arbitrary inertial reference point, dynamical analysis indicates that 18 firstorder differential equations of motion are required to mathematically model the system comprised of the three bodies. This number, however, is reduced by considering the relative motion. An infinitesimally small point mass, P3 (of mass m3), moving with respect to two point masses, or primaries, P1 (of mass m1) and P2 (of mass m2), appears in Figure 2.1. The masses are defined such that m1 > m2 ≫ m3 , restricting the problem in the sense that all gravitational influence exerted by m3 is neglected.
With this
assumption, the motion of P1 and P2 is entirely Keplerian, and reduced to the solution of the two-body problem. Additionally, this two-body motion is constrained by assuming that the primaries move in a circular orbit about their common center of mass, or barycenter, B. As a result, the problem only requires 6 first-order differential equations.
2.1.2 Geometry An inertial reference frame, I, described in terms of unit vectors Xˆ − Yˆ − Zˆ , is centered at B such that the Xˆ − Yˆ plane is defined to be coincident with the orbital plane of the primaries. Since the primary motion is Keplerian, it is constrained to the Xˆ − Yˆ plane, however, the third body can move in any of the three spatial dimensions. A rotating frame, S, with coordinate axes xˆ − yˆ − zˆ is initially aligned with I, then rotates through the angle θ , such that the xˆ -axis is always directed from P1 toward P2. Both the zˆ direction and Zˆ -direction are parallel to the orbital angular velocity vector of the primaries, and, thus, the yˆ and Yˆ axes complete the respective right-handed systems. Due to the circular primary motion, the angular rate, θɺ , is constant and equal to the mean motion, n. The position of each body, Pi, with respect to the barycenter is defined by Ri , and the relative position of P3 with respect to P1 and P2 is defined by R13 and R 23 , respectively. Note that the overbars ( ) indicate vectors.
12
Yˆ yˆ
P3 ( m3 ) R 23 R3 R13
P1 ( m1 )
xˆ P2 ( m2 )
R2
θ
B R1
Xˆ
Figure 2.1 Geometry in the Restricted Three-Body Problem
2.1.3 Equations of Motion
One goal of this analysis into the restricted problem of three bodies is a description of the motion of the infinitesimal mass, P3, subject to the gravitational influence of the primaries. From Newton’s Second Law, the vector differential equation for motion of P3 is written
I 2 d R Gm m Gm m ∑ F = m3 dt 2 = − R13 3 R13 − R23 3 R 23 , 13 23
(2.1)
where the superscript I represents differentiation in the inertial frame. A standard nondimensionalization used in the CR3BP is employed here. Since the mass of the third body is negligible, the characteristic mass, m*, is the sum of the two primary masses, i.e.,
m∗ = m1 + m2 .
(2.2)
The characteristic length, l*, is then the constant distance between the primaries, i.e., the scalar distance,
l ∗ = R1 + R2 .
(2.3)
13 Finally, the characteristic time, τ*, is defined such that the nondimensonal gravitational constant, G, is unity, i.e., G = 1 . This property is accomplished through the use of Kepler’s third law, i.e.,
τ∗ =
l *3 , ɶ * Gm
(2.4)
where Gɶ represents the dimensional value of the gravitational constant for clarity. These newly defined natural units lead to the following nondimensional quantities, Ri Rij τ m ri = * , rij = * , µ = 2* , t = * , l l m τ
(2.5)
where µ is denoted as the mass ratio. The motion of the third mass is now expressed in terms of these quantities by dividing equation (2.1) by m3 and the appropriate characteristic units in equations (2.2)-(2.4), I 2 d r3 1− µ µ = − 3 r13 − 3 r23 . 2 dt r13 r23
(2.6)
The kinematical expansion of the (inertial) first and second derivatives on the left side of the expression in equation (2.6) exploits the well-known operator relationship, I S I ɺr = dr3 = dr3 + I ω S × r , (2.7) 3 dt dt S I d 2 r3 S d 2 r3 dr3 I S I S I ɺɺ I S (2.8) r3 = = + 2 ω × + ω × ω × r3 , 2 2 dt dt dt where I ω S is the angular velocity of the rotating frame, S, with respect to the inertial S 2 d r3 frame. The second derivative, , in equation (2.8) can also be expanded dt 2 kinematically in terms of the nondimensional, cartesian rotating frame, S, r = xxˆ + yyˆ + zzˆ , S dr3 ɺ ˆ + yy ɺ ˆ + zz ɺˆ , = xx dt S d S dr 3 xxˆ + ɺɺ yyˆ + ɺɺ zzˆ , = ɺɺ dt dt
(2.9) (2.10)
(2.11)
14
where dots indicate derivatives with respect to nondimensional time. In this case, I S ω = nzˆ = zˆ , since the nondimensional mean motion is equal to one. The inertial acceleration in the rotating frame is expressed by substituting equation (2.11) into equation (2.8), resulting in the kinematical expansion, I ɺɺ r3 = ( ɺɺ x − 2 yɺ − x ) xˆ + ( ɺɺ y + 2 xɺ − y ) yˆ + ɺɺ zzˆ .
(2.12)
The radius vectors of relative position can also be expanded in terms of the rotating coordinate frame, i.e.,
r13 = ( x − µ ) xˆ + yyˆ + zzˆ ,
(2.13)
r23 = ( x − (1 − µ ) ) xˆ + yyˆ + zzˆ .
(2.14)
Finally, the equations of motion in the rotating frame are derived by combining the kinematics (equations (2.12)-(2.14)) and the kinetics (equation (2.6)) associated with m3 to yield the following scalar, second-order differential equations, ɺɺ x − 2 yɺ = x −
(1 − µ )( x − µ ) − µ ( x + 1 − µ ) ,
(2.15)
(1 − µ ) y − µ y ,
(2.16)
r133
ɺɺ y + 2 xɺ = y − ɺɺ z=−
r233
r133
r233
(1 − µ ) z − µ z . r133
r233
(2.17)
Equations (2.15)-(2.17) are also written more compactly by introducing the pseudopotential function, U, U=
(1 − µ ) + µ r13
r23
+
1 2 x + y2 ) , ( 2
(2.18)
reducing equations (2.15)-(2.17), i.e., ɺɺ x − 2 yɺ = U x ,
(2.19)
ɺɺ y + 2 xɺ = U y ,
(2.20)
ɺɺ z = Uz .
(2.21)
15
where U j =
∂U . Equations (2.15)-(2.17) are particularly useful in numerical methods ∂x j
due to the inherent nondimensional scaling.
2.1.4 Libration Points
Since the equations of motion in the restricted problem do not possess time explicitly due to the formulation in a rotating frame, the possibility exists for time invariant equilibrium locations.
Such solutions are characterized by stationary position and
velocity in the synodic frame corresponding to the nonlinear system of differential equations.
These particular solutions are determined by nulling the velocity and
acceleration terms in equations (2.15)-(2.17), resulting in the scalar equations,
− xeq = −
(1 − µ ) ( xeq − µ ) µ ( xeq + 1 − µ ) −
r133 eq
− yeq = − 0=−
r233 eq
(1 − µ ) yeq µ yeq −
r133 eq
r233 eq
(1 − µ ) zeq µ zeq r133
−
r233
,
.
Equation (2.24) is readily solvable, that is, zeq = 0 .
,
(2.22)
(2.23)
(2.24) Substitution of this result into
equations (2.22) and (2.23) produces a coupled system of two equations and two unknowns, xeq and yeq . As discovered by Lagrange [9], if r13 = r23 = 1 , then equations (2.22) and (2.23) reduce to identity, implying that two of the equilibrium points are located at vertices of two unique equilateral triangles. Thus, in cartesian coordinates, the primaries comprise two of the common vertices of both triangles, with the remaining vertex defined by xeq =
1 3 − µ and yeq = ± . 2 2
16
Three other equilibrium points also exist along the x-axis. Discovered first by Euler [5], they are denoted the collinear points and can be computed by forcing yeq = zeq = 0 . Substitution into equation (2.22) yields,
xeq −
(1 − µ ) ( xeq − µ ) µ ( xeq + 1 − µ ) xeq − µ
3
−
xeq + 1 − µ
3
= 0.
(2.25)
Equation (2.25) is a quintic equation in xeq. These solutions require numerical rootsolving methods that ultimately yield three real solutions, labeled L1, L2, and L3. The L1 and L2 points are defined such that L1 is between the primaries, L2 is on the far side of the smaller mass, and L3 is nearly a unit distance from the larger primary. All five libration points appear in Figure 2.2.
yˆ L4
30
L3
60 L1
B
L2
60
30
L5 Figure 2.2 Equilibrium Point Locations for the CR3BP
xˆ
17
2.1.5 Formulation Relative to P2 in the CR3BP Using the same rotating frame and characteristic quantities as in the original development of the equations of motion, an alternate formulation in the CR3BP is employed when low-thrust terms are included; this alternative formulation can reduce numerical sensitivity when an additional force of very low magnitude is added to the model. The origin is shifted from the barycenter to the smaller primary, P2 (as defined in Figure 2.3), and the equations of motion are rewritten as a function of position, r , and velocity, v , relative to the rotating frame, ɺɺ r = g (r ) + h (v ) , (2.26) where
( x + 1 − µ ) + κ ( x + 1) + µ ρ g (r ) = y (1 + κ ) , κ r3
(2.27)
2 yɺ h ( v ) = −2 xɺ . (2.28) 0 Note that the kinematical terms in h ( v ) have been shifted to the right side of the equation. The intermediate term, κ , simply allows equation (2.27) to be written in a more compact form; it is defined,
κ =−
(1 − µ ) − µ ρ3
r3
.
(2.29)
The unit vectors associated with the P2-centered frame are defined parallel to xˆ − yˆ − zˆ , and are defined as rˆ1 − rˆ2 − rˆ3 . The radius vectors of relative position are expanded as follows, r = r1rˆ1 + r2 rˆ2 + r3rˆ3 ,
(2.30)
ρ = ( r1 − 1) rˆ1 + r2 rˆ2 + r3rˆ3 .
(2.31)
This formulation is preferred when attempting to accurately propagate finite-thrust transfers to libration point orbits at L1 and L2.
18
rˆ2
P3 ( m3 ) r
ρ
l* P1 ( m1 )
rˆ1
P2 ( m2 )
Figure 2.3 Geometry of P2-Centered Rotating Frame
2.2 Natural Periodic Orbits in the CR3BP
In his 1892 study of the n-body problem, “Méthodes Nouvelles de la Méchanique Celeste” [7], Poincaré focused on the behavior of the nonintegrable solutions as t → ∞ . There is no practical method of numerical integration to evaluate absolute behavior in this problem in terms of these conditions. Thus, Poincaré’s investigation focused on the behavior of periodic orbits – the only viable subset of solutions in the problem of three bodies for which motion can be predicted as t → ∞ .
For nonintegrable dynamical
systems, complete information on an orbit requires either an asymptotic, periodic, or almost periodic structure. In the CR3BP, the Hamiltonian consistent with a formulation relative to the rotating frame is time invariant, and an infinite number of periodic solutions are available. Initially, insight concerning the different types of orbital motion in the restricted problem is gained by linearizing relative to the libration points. A firstorder linearization process yields approximate gradient information for use in evaluating the stability of the equilibrium points, as well as applications to some nonlinear targeting algorithms. Also, the linear solutions are eventually extended into the actual nonlinear problem by numerically solving a two-point boundary-value problem that exploits symmetry across the x-axis.
19
2.2.1 First-Order Variational Equations Relative to the Collinear Points
Analyzing the linear equations in the neighborhood of a libration point offers information concerning potential bounded behavior.
Let the variables, ξ, η, and ζ
indicate the relative position of the third body, or spacecraft, with respect to the collinear libration point, i.e.,
ξ = x − xeq , η = y − yeq , ζ = z − zeq .
(2.32)
Linearizing the equations of motion in equations (2.15)-(2.17) relative to the libration point and using a first-order Taylor series expansion, results in an expression of the following form,
ξɺɺ ξɺ ξ ηɺɺ = B η + C ηɺ , ζɺɺ ζɺ ζ
(2.33)
where the matrices B and C are defined as follows,
U xx U xy U xz B = U yx U yy U yz , U zx U zy U zz X = X eq
(2.34)
0 2 0 C = −2 0 0 . 0 0 0
(2.35)
The subscripts “ij”, on U ij denote evaluation of the second partial derivative of the
∂ 2U . Note that the B matrix is evaluated on the reference, in ∂j∂i this case, the libration point, X = X eq . Thus, the differential equation in equation (2.33)
pseudopotential function,
is linear with constant coefficients. The relationship is rewritten in state-space form, (2.36) γɺ = Aγ , where
γɺ = ξ η ζ
T
ξɺ ηɺ ζɺ ,
(2.37)
20
0 I3 A= . B C
(2.38)
The A matrix is a 6x6 matrix composed of four 3x3 submatrices, the matrix 0 is simply a matrix of zeros, and I3 is defined as the 3x3 identity matrix. Evaluating the C-matrix at the collinear libration points, where yeq = zeq = 0, yields
U xz
X = X eq
U xy
= U yz X = X eq
X = X eq
= U yx
= U zx X = X eq
X = X eq
= U zy
= 0 , U xx
X = X eq
X = X eq
= 0 , U zz
> 0 , U yy
X = X eq
X = X eq
<0,
> 0.
(2.39) (2.40)
Equations (2.38)-(2.40) simplify the linearized equations of motion to the form,
ξɺɺ U xx 0 0 ηɺɺ = 0 U yy 0 ζɺɺ 0 0 U zz X = X eq
ξɺ ξ η + C ηɺ . ζɺ ζ
(2.41)
From inspection of the third row of equation (2.41), it is clear that the out-of-plane component, ζ , is completely decoupled from the independent variables, ξ and η . The characteristic equation for this out-of-plane motion in equation (2.41) also possesses purely imaginary roots. Thus, the solution is oscillatory,
ρ = C1 cos ω t + C2 sin ω t , where the frequency is ω = U zz
X = X eq
(2.42)
. The first two rows in equation (2.41) are coupled
through the matrix C and represent the in-plane motion, with a general solution of the following form, 4
ξ = ∑ Ai eλ t , i
(2.43)
i =1
4
η = ∑ Bi eλ t , i
(2.44)
i =1
where Ai and Bi are dependent. To determine the eigenvalues, λi , Szebehely [10] uses a special form of the characteristic equation (a quadratic), i.e., Λ 2 + 2 β1Λ − β 22 = 0 , where
(2.45)
21
β1 = 2 −
U xx
X = X eq
+ U yy
X = X eq
,
(2.46)
>0,
(2.47)
2
β 22 = − U xx
X = X eq
U yy
X = X eq
λ=± Λ.
(2.48)
The quadratic roots are determined first, i.e.,
Λ1 = − β1 + β12 + β 22 > 0 ,
(2.49)
Λ 2 = − β1 − β12 + β 22 < 0 .
(2.50)
Substituting equations (2.49)-(2.50) into equation (2.48) reveals the four characteristic roots,
λ1,2 = ± Λ1 (real),
(2.51)
λ3,4 = ± Λ 2 (imaginary).
(2.52)
Of course, the positive real root in equation (2.49) results in a positive real eigenvalue in equation (2.51) and unbounded behavior in the general solution in equations (2.43)-(2.44) as t → ∞ . The dependency between Ai and Bi is resolved by substituting equations (2.51)-(2.52) into the first two expressions in equations (2.41) and simultaneously solving the equations, thus,
βi =
λi − U xx
X = X eq
2λi
Ai = α i Ai .
(2.53)
Equation (2.43)-(2.44) and their derivatives are evaluated at the initial time to determine the unknown initial conditions. Equation (2.53) is substituted into equation (2.44) so that the results are entirely in terms of the independent variables Ai, i.e., 4
ξ0 = ∑ Ai eλ t , i 0
(2.54)
i =1
4
ξɺ0 = ∑ λi Ai eλ t , i 0
(2.55)
i =1
4
η0 = ∑ αi Ai eλ t , i 0
i =1
(2.56)
22
4
λi t0 ηɺ0 = ∑αiλi Ae . i
(2.57)
i =1
Since the coefficients A1 and A2 are associated with the real eigenvalues in equations (2.51)-(2.52), fixing the values of A1 and A2 to zero ensures that any exponential increase and decay is suppressed. Thus, a bounded planar solution emerges,
ξ = ξ0 cos s ( t − t0 ) +
η0 sin s ( t − t0 ) , β3
(2.58)
η = η0 cos s ( t − t0 ) − β 3ξ0 sin s ( t − t0 ) ,
(2.59)
λ3 = is ,
(2.60)
where
s = β1 + ( β 22 + β 32 )
12
,
(2.61)
* s 2 + U xx β3 = , 2s
(2.62)
α3 = iβ3 .
(2.63)
Once the initial conditions ξ0 and η0 are selected, ξɺ0 and ηɺ0 are predetermined to enforce A1 = A2 = 0. The resulting orbit is an ellipse with a collinear libration point at the center. The semimajor axis is parallel to the unit vector yˆ and the semiminor axis is in the x-direction. An example of such an orbit about L1 in the Earth-moon system appears in Figure 2.4. The period is IP =
2π for the planar motion, and the root, s, may be s
selected to achieve a specific value of IP.
In this linear model, the out-of-plane
frequency is not commensurate with the in-plane frequency. Nevertheless, these planar solutions form the basis for determining initial conditions consistent with planar nonlinear orbits; such planar orbits that are solutions in the nonlinear problem are then used in computing fully three-dimensional nonlinear orbits. Although the first-order linear set of initial conditions do not result in a closed periodic orbit in the nonlinear problem, initial conditions “sufficiently close” to the libration point still yield an initial guess that may be used in targeting algorithms.
23
Earth-Moon System
γ 0 = [0.1 0 0 0
- 0.837227214342055 0]
T
µ = 0.0121505649407351 , 1 unit = 385692.5 km
Moon
To Earth
L1
Figure 2.4 Linearized L1 Periodic Orbit
2.2.2 The State Transition Matrix
Besides stationary equilibrium points, dynamical information is also sought relative to time-varying solutions in the nonlinear problem. Numerical computation of trajectory arcs, such as periodic orbits and orbital transfers to a specific target in the CR3BP, requires the use of differential corrections procedures.
Before introducing these
corrections procedures, an important result from linearization of the equations of motion is required. Define X * ( t ) as a time-varying reference solution such that X ( t ) = X * ( t ) + δ X ( t ) . Let the perturbed state relative to the reference trajectory be defined as,
24
δ X ( t ) = {δ x δ y δ z δ xɺ δ yɺ δ zɺ} . T
(2.64)
Recall that a less general expansion was introduced in equation (2.36) with a constant reference. Linearizing the nonlinear equations of motion in equations (2.15)-(2.17) with a first-order Taylor series expansion results in the familiar equation δ Xɺ = A ( t ) δ X ,
(2.65)
where
I3 0 A (t ) = . B t C ( )
(2.66)
Since all linearization now employs a time-variable trajectory as the reference, the Amatrix is now time-dependent. This time dependency occurs specifically in the partial derivatives of the B-matrix (see equation (2.34)) where the second partials are evaluated along the reference trajectory path. The solution of the linear system of equations in equation (2.65) is expressed in terms of the state transition matrix, Φ , such that, δ X ( t + t0 ) = Φ ( t + t0 , t0 ) δ X ( t0 ) ,
(2.67)
where the matrix Φ is obtained via the matrix differential equation,
ɺ (t + t , t ) = A (t ) Φ (t + t , t ) . Φ 0 0 0 0
(2.68)
Numerical integration of equation (2.68) requires equation (2.66) to be evaluated along the trajectory at all times, and is, thus, simultaneously integrated along with equations (2.15)-(2.17). Once obtained, the state transition matrix is essentially a sensitivity matrix that is an approximate mapping of trajectories in the neighborhood of the reference. Decomposed into components, the 6x6 matrix is represented as,
25
∂x ∂x 0 ∂y ∂x 0 ∂z ∂x 0 Φ ( t + t0 , t0 ) = ∂xɺ ∂x 0 ∂yɺ ∂x 0 ∂zɺ ∂x 0
∂x ∂y0
∂x ∂z0
∂x ∂xɺ0
∂x ∂yɺ 0
∂y ∂y0
∂y ∂z0
∂y ∂xɺ0
∂y ∂yɺ 0
∂z ∂y0 ∂xɺ ∂y0 ∂yɺ ∂y0 ∂zɺ ∂y0
∂z ∂z0 ∂xɺ ∂z0 ∂yɺ ∂z0 ∂zɺ ∂z0
∂z ∂xɺ0
∂z ∂yɺ 0
∂xɺ ∂xɺ0
∂xɺ ∂yɺ 0
∂yɺ ∂xɺ0
∂yɺ ∂yɺ 0
∂zɺ ∂xɺ0
∂zɺ ∂yɺ 0
∂x ∂zɺ0 ∂y ∂zɺ0 ∂z ∂zɺ0 . ∂xɺ ∂zɺ0 ∂yɺ ∂zɺ0 ∂zɺ ∂zɺ0
(2.69)
If an exact solution to the equations of motion is available, it is possible that equation (2.68) can be integrated analytically. Otherwise, it may be necessary to numerically integrate equation (2.68) to generate a time-varying history for Φ .
2.2.3 The Fundamental Targeting Relationships
The availability of the STM is a critical component in any targeting algorithm. Consider the initial state of a spacecraft on a reference path, X * ( t0 ) ; then, the reference state at some future time can be denoted as X * ( t0 + ∆t ) . Note that any superscript ‘*’ refers to a condition on the reference path. Such points are represented as points A and C, respectively, in Figure 2.5. Point C, downstream from point A along the reference path, can be modeled as a numerical mapping of the initial state and the time interval, ∆t , X * ( t0 + ∆t ) = f X * ( t0 ) , ∆t . (2.70) The initial state at time t0 on some neighboring trajectory is represented by X ( t0 ) . A contemporaneous variation at point A shifts the spacecraft onto point B along X at time t0 , such that a new perturbed state is defined X ( t0 ) = X * ( t0 ) + δ X . After a specified time interval, ∆t , A is mapped to C and B is mapped to D′ . Thus, δ X ( t0 + ∆t )
(
)
26 represents the six-dimensional state at D′ with respect to C and is the contemporaneous variation at time t * = t0 + ∆t . Let point D be some other point along the neighboring path. Point D is achieved via an arbitrary time interval relative to B, i.e., ( t0 + ∆t + δ t ). Thus, comparing the state at D to the point C along the reference path defines a noncontemporaneous variation, δ X ( t0 + ∆t + δ t ) . The noncontemporaneous variation in the state between points C and D is mapped as a function of the reference state and any additional state and time perturbation, X * ( t 0 + ∆t ) + δ X f = f X * + δ X 0 , ∆ t + δ t .
(
)
(2.71)
For a general targeting algorithm, consider the governing differential equations from equations (2.15)-(2.17) as written in first-order form, ɺ X = f X ,t .
(
)
(2.72)
Let X = X * + δ X , and t = t * + δ t . The differential equations on the neighboring path BD are written as, ɺ ɺ X * + δ X = f X + δ X , t* + δ t .
(
)
(2.73)
Given equation (2.73), a first-order Taylor series expansion about the reference path A-C, yields, ɺ * ɺ * * ∂f X + δ X = f X ,t + ∂X
(
)
∂f δX + X =X* ∂t t =t
X =X* t =t
*
δt .
(2.74)
*
ɺ Equation (2.74) is reduced by noting that equation (2.72) eliminates X * and f X * , t * .
(
∂f Then, ∂X
X =X*
∂f is simply A ( t ) , and ∂t
t =t*
X =X *
)
is an additional matrix of partial time
t =t*
derivatives, denoted as K ( t ) , i.e.,
δ Xɺ = A ( t ) δ X + K ( t ) δ t .
(2.75)
27
The solution to equation (2.75) is, δ X ( t0 + ∆t + δ t ) = Φ ( t0 + ∆t , t0 ) δ X 0 + Xɺ
X =X * t =t
δt .
(2.76)
*
where the following definitions, apparent in Figure 2.5, apply, such that, X ( t0 + ∆t + δ t ) = X * ( t0 + ∆t ) + δ X f , δ X 0 = δ X ( t0 ) , δ X f = δ X ( t0 + ∆t + δ t ) .
X
Neighboring Path X B
A
X*
(2.78) (2.79)
Target D′
D
δ X ( t0 + ∆t )
δ X ( t0 )
(2.77)
δX( t0 +∆t +δτ) =δXf
C
Reference Path Initial State
t0
∆t
δt
t
Figure 2.5 Basic Diagram for a Free Final Time Targeting Scheme
In a more compact form, equation (2.76) becomes,
ɺ * δ X 0 X ( t0 + ∆t ) , (2.80) δ t ɺ where it is emphasized that Φ ( t0 + ∆t , t0 ) and X * ( t0 + ∆t ) are evaluated on the reference δ X f = Φ ( t0 + ∆t , t0 )
path. Equation (2.80) is the fundamental basis of most numerical targeting schemes. For time-fixed problems, equation (2.80) simplifies to
28
δ X f = Φ ( t0 + ∆t, t0 ) δ X 0 .
(2.81)
Thus, equations (2.80)-(2.81) comprise a process for approximating state sensitivities, given initial variations in the state vector. Note that the columns of Φ ( t0 + ∆t , t0 ) are associated with the control parameters, δ X 0 and δ t , and the rows correspond to the end point constraint parameters, δ X f .
2.2.4 Periodic Orbits
Simple targeting procedures to converge upon a periodic orbit often exploit symmetric behavior; additional conditions at various plane crossings can also benefit the process. Consider the planar problem of computing periodic Lyapunov orbits in the vicinity of L1 in the Earth-moon system. Due to the x-axis symmetry, only the first half-period of the orbit must be determined. Thus, for some given initial state on the x-axis, represented by X 0 in Figure 2.6, an arbitrary initial velocity, yɺ 0 , will produce the dashed curve in Figure 2.6. This initial state is specified by nonzero state components in only the xdirection and the yɺ -direction in terms of the linearized system, T X 0 = {x0 0 0 yɺ 0 } .
(2.82)
Correction of the nonlinear propagation results in a perpendicular crossing (the solid curve in Figure 2.6), but requires the solution of a two-point boundary-value problem. Satisfying the constraint is accomplished by adjusting one of the two nonzero initial states, as well as the propagation time. For purposes of demonstration, variations on the initial y-velocity, yɺ 0 , are permitted, i.e., δ yɺ 0 and δ t are controls, while all other initial variations remain fixed, i.e.,
δ x0 = 0 , δ y0 = 0 , δ xɺ0 = 0 .
(2.83)
29
When the trajectory terminates after a half-period, the final state will be of the following form,
X f = {x f
yɺ f } . T
0 0
(2.84)
Note that determination of the perpendicular x-axis crossing requires two constraints: one to enforce termination of the trajectory on the x-axis; and, one to enforce elimination of the final x-velocity, i.e.,
δ yf = 0 ,
(2.85)
δ xɺ f = 0 .
(2.86)
The control on the time variation, δ t , is only implicit since the x-axis crossing is always forced, and thus, equation (2.81) is always satisfied. The constraint in equation (2.86) remains active until the actual final x-velocity, xɺ fa , is equal to the desired final xvelocity, xɺ fd , i.e., zero. Therefore,
δ xɺ f = xɺ f − xɺ f = − xɺ f . d
a
a
(2.87)
The variational relationship between the initial and final states for a perpendicular crossing is the result of inserting equations (2.84)-(2.85) and equation (2.87) into equation (2.80),
∂y 0 ∂yɺ 0 ɺ = − x f ∂xɺ ∂yɺ 0
T
yɺ δ yɺ 0. δt ɺɺ x
(2.88)
Equation (2.88) is solved for the initial variation in the y-velocity,
∂xɺ
ɺɺ x ∂y
δ yɺ 0 = − − xɺ f . ɺ ɺ ɺ ∂ y y ∂ y 0 0 a
(2.89)
Equation (2.89) produces the estimate for the control necessary to meet the constraints. Since this update is based upon a linear approximation, an iterative process is employed. Thus, the initial y-velocity is updated until equation (2.86) is satisfied numerically, that is, until xɺ fd < ε , where ε is a user-defined error tolerance.
30
The computation of an entire family of periodic orbits requires a method of continuation to produce the new initial conditions. Given one periodic orbit, a second orbit is sought by generating a guess for initial conditions that are suitable to yield a second periodic orbit.
To produce the new guess, a new state is generated by
incrementing the value of x0 and using the previously successful value of yɺ 0 as an initial guess to initiate the new corrections process. An example result of such a process appears in Figure 2.7.
yˆ
X0
L1
X fd
X fa
xˆ
Figure 2.6 Targeting a Perpendicular xˆ -axis Crossing in the CR3BP
31
Earth-Moon System
µ = 0.0121505649407351 , 1 unit = 385692.5 km
L1
To Earth Moon
Figure 2.7 Several L1 Lyapunov Orbits Obtained Via Continuation
2.3 Invariant Manifolds
Henri Poincaré [7] determined that any dynamical system can be analyzed from a geometrical perspective. Knowledge of the phase-space of a dynamical system allows the decomposition of the flow into subspaces, thus characterizing the behavior of a system.
In obtaining transfer paths in the CR3BP, it is very useful to exploit the
knowledge of the flow in the vicinity of a reference solution.
Not only does this
dynamical systems analysis allow insight into the motion in the vicinity of the libration points, it is particularly significant for an understanding of the behavior near the periodic orbits.
32
2.3.1 Stable and Unstable Manifolds Associated with the Collinear Points
Analysis of dynamical systems from the perspective of manifolds is available in various references [47-51], but a summary is useful. Consider the first-order form of the equations of motion (equation (2.72)) that may be employed to equivalently represent a general nonlinear vector field. For such a general field, let the state vector X ( t ) ∈ ℝ n , and f : Ψ → ℝ n be a smooth function over the subset Ψ → ℝ n generating a flow, ΓT ( Ψ ) . The flow depends upon the initial condition X 0 , so it is also equivalently written as ΓT ( t0 ) , and may be propagated over time t. In the dynamical system of interest, one important aspect of the behavior of the flow can be evaluated via a linear stability analysis. To introduce this topic, use a time-invariant system as an example. In the Earth-moon system, the equilibrium point at Li is a solution to the nonlinear equations (2.15)-(2.17); these equations, of course, can be written in the form in equation (2.72). Recall that a Taylor series expansion results in the vector linear variational equation, equation (2.36), with constant Jacobian matrix, A. As in section 2.2.1, the reference solution is the fixed equilibrium point at Li, that is, the constant X ref = X eq . The solution to equation (2.36) appears in the following form [49], γ ( t ) = e A( t −t0 )γ ( t0 ) .
(2.90)
The eigenvalues, λ j , of the A-matrix are distinct [49] and thus, linearly independent. Equation (2.84) is further simplified as follows γ ( t ) = Ne Λ( t −t0 ) N -1γ ( t0 ) , (2.91) where the matrix of eigenvectors, N = N1 N 2 ... N n , corresponds to the diagonal matrix of eigenvalues, Λ , containing entries λ1 , λ2 ,..., λn . Thus, every N j corresponds to
λ j in equation (2.91).
The matrix e Λ (
Λ t − t0 )
is diagonal or block diagonal, and the
eigenvalues, λ j , are denoted as the characteristic exponents associated with the local flow. The linear system is decomposed into a sum that is a function of each mode,
33
n
γ (t ) = ∑ c je
λ j ( t − t0 )
j =1
Nj ,
(2.92)
where the coefficients cj are determined from γ ( t0 ) . The matrix A is composed of a total of n eigenvalues. Of this total set, nu eigenvalues are positive real or possess positive real components, ns are defined with negative real components, and nc eigenvalues are purely imaginary. Thus, n = nu + ns + nc. Due to the linear independence of the system, the solution space is comprised of the corresponding invariant subspaces Eu , E s , and Ec . Flow that is contained within one subspace at the initial time will remain in that subspace throughout the dynamical evolution of the system. Trajectories within the Eu subspace approach γ = 0 as t → −∞ , those within E s subspace approach γ = 0 as t → +∞ , and those within Ec neither grow or decay over time. Note that when all of the eigenvalues of A possess non-zero real parts, X eq is defined as a hyperbolic equilibrium point. These observations allow for the introduction of the Stable Manifold Theorem for Flows [49].
Theorem 2.1: (Stable Manifold Theorem for Flows). Suppose that Xɺ = f X
( )
has a
hyperbolic equilibrium point X eq . Then there exist local stable and unstable manifolds Wsloc X eq , Wuloc X eq , of the same dimensions ns, nu as those of the eigenspaces E s and
( )
( )
Eu of the linearized system, and tangent to E s and Eu at X eq .
Wsloc X eq
( )
and
Wuloc X eq are as smooth as the function f .
( )
A center subspace, Ec , associated with the linear system and corresponding to nc also exists, however, Wcloc X eq is not necessarily tangent to Ec [49]. As an example of
( )
Theorem 2.1, consider a constant equilibrium point. A conceptual representation of the Stable Manifold Theorem in a two-dimensional phase space appears in
34
Figure 2.8. Note the flow toward and away from the equilibrium point. The eigenvectors Vs and Vu correspond to the stable and unstable eigenvalues, respectively, and span the subspaces E s and Eu to form a vector basis in ℝ 2 . For the nonlinear system, Vs and Vu can be used to numerically approximate the manifolds in the nonlinear system ( Wsloc , Wuloc ) near X eq . The local stable manifold, Wsloc X eq , includes every initial
( )
condition that produces trajectories that asymptotically approach the equilibrium point, converging along the one-dimensional eigenvector tangent to E s , in both the positive and negative directions. Conversely, the local unstable manifold, Wuloc X eq , contains every
( )
initial condition that produces a trajectory that asymptotically departs the equilibrium point. Near X eq , it is tangent to Eu in the positive and negative directions along the corresponding one-dimensional eigenvector. The local manifolds, Wsloc and Wuloc , also possess global analogs Ws and Wu .
Actual computation of these global manifolds
requires numerical simulation forward and backward in time using initial conditions from
E s and Eu . Thus, the global flow is extended as a function of the localized behavior, Ws X ref = ∪ ΓT Wsloc X ref , (2.93) t ≤0
(
)
(
)
Wu X ref = ∪ ΓT Wuloc X ref . t ≤0
(
)
(
)
(2.94)
Note that these definitions can be extended to a generalized reference, X ref . Es
Vu
−Vs
Eu Wu +
Ws −
loc
loc
X eq Wu−
loc
−Vu
Vs
Figure 2.8 Stable and Unstable Manifold at X eq
Ws +
loc
35
2.3.2 Invariant Manifolds Relative to a Fixed Point
In addition to individual equilibrium points that serve as a reference solution, dynamical systems theory also supports periodic orbits as time-varying reference solutions. Consider the monodromy matrix, i.e., the state transition matrix at the end of one period, Φ ( IP,0 ) . The state transition matrix is a linear map, thus, the monodromy matrix defines a linear, stroboscopic map sampled over one period. Essentially, by creating a map, the continuous time system is transformed to a discrete time system. Any point along a periodic orbit can be used to create a map; this point is termed a “fixed point”. For planar periodic orbits in the vicinity of the collinear libration points, one of the two points on the x-axis is initially convenient as a fixed point, but not necessary. After one revolution at this discrete reference point, the eigenvalues and eigenvectors of the monodromy matrix characterize the local phase space similar to the subspaces in the analysis of an equilibrium point. However, in the case of a periodic orbit, the fixed point is redefined and the monodromy matrix must be recomputed at different points of interest along the orbit, as the behavior associated with the local phase space varies. This varying behavior along the orbit, characterized by different discrete points, results in unique eigenvectors. The corresponding eigenvalues, however, do not change since they are a property of the orbit [49]. For periodic orbits, recall that the linearized equations of motion (equation (2.65)) result in a time varying matrix A(t) to be evaluated at all points along the orbit. As a consequence of Floquet Theory [49], the STM can be rewritten in the general form,
Φ ( t ,0 ) = F ( t ) e Jt F-1 ( 0 ) ,
(2.95)
where F ( t ) is a periodic matrix, and J is a normal, block diagonal matrix whose elements are termed the Poincaré exponents.
Since F ( t ) is a periodic matrix,
F ( IP ) = F ( 0 ) . Thus, the monodromy matrix can be expressed,
Φ ( IP,0 ) = F ( 0 ) e JIP F-1 ( 0 ) . Rearranging equation (2.96),
(2.96)
36
e JIP = F-1 ( 0 ) Φ ( IP,0 ) F ( 0 ) .
(2.97)
Equation (2.97) is an explicit statement that F ( 0 ) is the eigenvector matrix associated with the monodromy matrix. The corresponding eigenvalues are contained in the matrix
e JIP . Defining the eigenvalues of Φ ( IP,0 ) as λi ; they are labeled the characteristic multipliers. Then, ρ i , the diagonal entries of J, are defined such that,
λi = e ρ IP , i
(2.98)
or equivalently,
ρi =
1 e ln λi . IP
(2.99)
The Poincaré exponents, ρ i , in equation (2.99) occur in positive/negative pairs according to Lyapunov’s theorem [51]. The stability associated with the fixed point of interest along the periodic orbit can now be summarized as follows. Assuming a hyperbolic system,
if λi < 1 for all λi , X = 0 as t → ∞ asymptotically stable, if λi > 1 for all λi , X = 0 as t → −∞ unstable. The stability of any fixed point, i.e., the behavior of λi , represents the stability of the periodic orbit. Moreover, due to the nature of the Poincaré exponents, along with an inspection of equation (2.98), it is apparent that the eigenvalues associated with a fixed point along the periodic orbit occur in complex conjugate pairs. Two of the eigenvalues will equal one, indicating a periodic orbit and, therefore, are subsequently used to reduce the dimension on the system and create the map. In Lyapunov orbits (in addition to all three-dimensional orbits considered later), the remaining real pair is used to specify the stable/unstable subspaces.
37
2.3.3 Computation of Manifolds Corresponding to Fixed Points Along an Orbit
Globalizing the stable and unstable manifold trajectories that correspond to fixed
E s and Eu .
points along a periodic orbit requires initial conditions in the subspaces
These conditions are estimated given the state and phase space information corresponding to any fixed point, X fp ( ti ) along the periodic orbit. Of course, any state along the orbit, X fp ( ti ) , will remain on the periodic orbit when propagated. To globalize the manifold trajectory and compute the flow toward and away from the periodic orbit, conditions on the manifold must be approximated near the fixed point. Thus, given X fp ( ti ) , a perturbation is added to shift the state into the desired subspace. Computation of this new state is accomplished by defining a perturbation in the direction of the stable or unstable eigenvector by some small distance d. If the eigenvectors associated with the stable and unstable mode are defined as Yˆ Ws ( ti ) = [ xs
Yˆ Wu ( ti ) = [ xu
yu
zu
xɺu
yɺ u
zɺu ]
T
ys
zs
xɺ s
yɺ s
zɺs ]
T
and
respectively, then a normalization in position
results in the definitions,
V Ws ( ti ) = V Wu ( ti ) =
Yˆ Ws ( ti ) xs2 + y s2 + z s2 Yˆ Wu ( ti ) xu2 + yu2 + zu2
,
(2.100)
.
(2.101)
The initial state vector to shift into E s or Eu is then represented by the full expressions X s ( ti ) = X fp ( ti ) ± d ⋅V Ws ( ti ) , (2.102) X u ( ti ) = X fp ( ti ) ± d ⋅V Wu ( ti ) . (2.103) The alternating signs on the displacement from X fp ( ti ) in equations (2.102)-(2.103) represent the fact that the trajectory may be perturbed in either direction in the stable or unstable subspace, as depicted in Figure 2.8. Propagating the states in equation (2.102) in positive time at fixed states along the entire orbit results in globalization of the stable manifold. Repeating the process on the states in equation (2.103) in negative time results
38
in globalization of the unstable manifold. Typically, in the Earth-moon system, a value of d = 50 km is sufficient to justify the linear approximation, yet still yield adequate integration time. An example Lyapunov orbit near L1 in the Earth-moon system with yamplitude Ay = 23,700 km appears in Figure 2.9. At various fixed points around the orbit, the eigenvectors of the monodromy matrix are computed and the initial states are shifted according to equations (2.102)-(2.103).
The initial conditions are then
numerically integrated in both directions. For the example in Figure 2.9, the number of fixed points examined is specified as nfp; for this case nfp = 50. The initial states are specified such that X fp ( ti ) = X fp ( iIP n fp ) , where IP is the orbital period and i = 1,2,3,…,nfp. Earth-Moon System
µ = 0.0121505649407351
Unstable Stable
Earth Unstable Stable Moon L1
Figure 2.9 Global Manifolds for an Earth-Moon L1 Lyapunov Orbit, Ay = 23,700 km
39
1,2,3,…,nfp. For each fixed point, the four sets of initial conditions on X s ( ti ) , X u ( ti ) are then numerically integrated in both directions until the Earth-relative flight path angle is zero for Earth-bound trajectories, or until the moon-relative flight path angle is zero for moon-bound trajectories. Thus, the red families of trajectories in Figure 2.9 represent the trajectories that asymptotically diverge from the orbit. Such trajectories are useful for transfer trajectories that begin at the orbit, and depart to another region of interest (for example, the surface of the moon). The green families of trajectories in Figure 2.9 represent the trajectories that asymptotically approach the orbit, and are exploited for purposes of transfers from a region of interest (for example, an Earth parking orbit) for delivery of a spacecraft into a libration point orbit.
2.4 Optimal Control Theory
Computation of manifold trajectories that depart from or arrive at libration point orbits comprises only part of the necessary background to develop a finite-thrust transfer. In the Earth-moon problem, for example, the manifolds do not pass very close to the Earth (see Figure 2.10). One design option to gain access to a manifold from the Earth, is a low-thrust transfer arc. Defining such an arc requires the determination of a thrust magnitude and directional history. Optimal control theory can be employed to solve this “steering” problem. The term “optimal control” is used to reference a branch of mathematics formally known as the calculus of variations, although the former term is commonly employed in dynamical applications [25]. The solution to a typical optimal control problem requires the determination of a set of conditions that minimize or maximize a scalar performance index with respect to admissible comparison paths, subject to differential constraints and boundary conditions.
To determine the solutions, i.e., “extremals”, a Taylor series
expansion leads to necessary and sufficient conditions for a local minimum.
As
expected, the first necessary condition requires that the differential of the performance index vanishes, and the second necessary condition requires non-negativity on the second
40
differential of the performance index. A stronger, and more general second necessary condition for a minimum is Pontryagin’s Minimum Principle. Note that the second necessary condition can also be viewed as a test to determine if the extremal (determined from the first necessary conditions) is actually a minimum or a maximum. Ultimately, for a solution to the optimal control problem, the determination of these extremals requires the solution of a two-point boundary-value problem (TPBVP) that may or may not be solvable in closed form. Because optimal control theory relies on comparing different neighboring paths, it is useful to note that a ‘**’ superscript indicates a condition on the optimal path.
2.4.1 Summary of the First Necessary Conditions for Optimal Control
Consider the free final time problem for the control history, u ( t ) , that minimizes the generalized scalar performance index, J, or cost function at a final time, tf ɺ J = φ X f , t f + ∫ H X , u , χ , t − χ T X dt , t
(
)
(
)
(2.104)
0
subject to the n differential constraints, ɺ X n = fn X , u, t ,
)
(2.105)
the n+1 prescribed initial constraints, ϕ n +1 X 0 , t0 = 0 ,
(2.106)
and the p prescribed final constraints ψp Xf =0,
(2.107)
(
(
( )
)
where p is some integer value such that 0 ≤ p ≤ n . Equation (2.107) is denoted the set of kinematic boundary conditions, and is a function of the final state. In equation (2.104), φ X f , t f is an endpoint function, and the second term is a path-wise function. The
(
)
scalar H is the Hamiltonian defined as,
41 H X , u, χ , t = L X , u, t + χ Τ fn X , u, t ,
(
)
(
)
(
)
(2.108)
where L X , u , t is a path-wise component to be minimized, and the second term is a set
(
)
of Lagrange multipliers or “costates”, χ X , u , t , adjoined to the differential constraints.
(
)
A minimization problem that only contains the endpoint function, φ X f , t f , is
(
)
commonly known as a Mayer problem. Those that only include path-wise components within the integral are Lagrange Problems. Finally, Bolza problems are defined as those that include both components. For the general Bolza problem stated in equation (2.104), the following assumptions are also imposed: 1. Define Ω as an arbitrary set in m-dimensional Euclidean space; the “measurable” control u ( t ) , defined over the closed time interval t0 , t f , possesses values lying
within Ω , i.e., u ( t ) ∈ Ω .
For all applications considered here, u ( t ) is a
bounded, piecewise continuous function. ∂f n X , u , t ∂f n X , u , t ∂L X , u , t 2. The functions f n X , u , t , , L X , u, t , , , and ∂X ∂t ∂X ∂L X , u , t ∂f n X , u , t are all continuous functions of X , u , t . Since and ∂t ∂u ∂L X , u , t do not need to exist, f n or L may contain instantaneous switching ∂u conditions on the control, i.e., u -type terms. 3. The endpoint component of the Bolza function, φ X f , t f ∈ C 1 , and the final
(
(
)
(
)
(
)
)
(
)
(
(
)
)
(
(
)
(
)
)
boundary conditions at ψ X f , t f ∈ C 1 , are continuously differentiable in all
( )
∂ψ i arguments; the vectors for i = 1,..., p are linearly independent in the region ∂X f of the optimal candidates X 0 , t0 , X f , t f .
(
)
42 Additional sets of Lagrange multipliers, ϖ and υ , are adjoined to the initial and final endpoint constraints, ϕ X 0 , t0 , and ψ X f , to construct the augmented performance
(
)
( )
index, tf
(
)
ɺ J ′ = Θ + ∫ H X , u , χ , t − χ T X dt , t0
(
)
(2.109)
where Θ = φ X f , t f + ϖ T ϕ X 0 , t0 + υ Tψ X f .
(
)
(
)
( )
(2.110)
Minimizing the performance index (equation (2.104)) is equivalent to minimizing the augmented performance index (equation (2.109)) if all constraints are satisfied.
Of
course, the minimization process initially requires that the first differential of the performance index, dJ ′ , vanishes. As a result, several necessary conditions arise. First, the well-known Euler-Lagrange equations must be satisfied to obtain an optimal state, X ** ( t ) , appropriate costates, χ ( t ) and optimal control history, u ** ( t ) [25], χɺ = − H xT X , u , χ , t (n-equations), (2.111)
(
0 = H u X , u , χ , t
(
)
)
(m-equations).
(2.112)
Because of the dependence on the states, the costate equations of motion (equations (2.111)) are integrated simultaneously with the state equations (equations (2.105)). Note that equations (2.105), (2.111), and (2.112) comprise the (2n + m) equations that are necessary to determine the n values of the state, X ( t ) , n costates χ ( t ) , and m controls
u ( t ) . To specify a complete two-point boundary-value problem (TPBVP), a total of (2n + 2) boundary conditions are required corresponding to the 2n initial and final states, as well as the initial and final times. From equations (2.106)-(2.107), (p + n + 1) boundary conditions are already available, leaving (n + 1 – p) boundary conditions that remain to be specified.
These remaining boundary conditions are determined from the
transversality condition (the final condition needed for dJ = 0 ), or from the natural boundary conditions (the final conditions needed for dJ ′ = 0 ). This development relies on the augmented performance index, however, further detail on the transversality
43
condition is available in Bryson and Ho [24]. Given the expression for the differential of the augmented performance index, the additional (n + 1 – p) boundary condition equations, that is, the natural boundary conditions arise, χ 0 = −ΘTX 0 , χ f = ΘTX f , H f = −Θt f .
(2.113)
When equations (2.113) are satisfied, the boundary conditions correspond to an identical extremal value of the performance index computed via the transversality conditions. Equations (2.111)-(2.113) are derived by equating the first differential of the augmented performance index to zero, and are, thus, only necessary conditions for a minimum value of J. But, for unbounded controls, equation (2.112) typically yields the control equation. For bounded controls, other approaches are often employed.
Of course, when the
performance index is to be maximized, J can simply be replaced by –J, and the same minimization process proceeds.
2.4.2 Tests for a Local Minimum Value of the Performance Index
The Euler-Lagrange equations result from the determination of stationary conditions on the first differential of the augmented performance index. Thus, the Euler-Lagrange equations are valid for either a maximum or minimum value. To ensure that the control actually results in the desired minimum value of the augmented performance index, a second necessary condition is required. The strongest and most general mathematical statement that guarantees minimal control is Pontryagin’s Minimum Principle [52], or more succinctly, the Minimum Principle. In a summary of the Minimum Principle, McShane [53] observes that the Hamiltonian, H , must be minimized over the set of all possible u . A rigorous proof of the theorem is offered by Pontryagin [54]. As a consequence, consider the problem posed in the previous section with the same assumptions. Then, the control u ** ( t ) ∈ Ω that yields a minimum value of J requires that the Hamiltonian H X ** ( t ) , u ** ( t ) , χ ( t ) , t to be continuous on t0 , t f and, H X ** ( t ) , u ** ( t ) , χ ( t ) , t ≤ H X ** ( t ) , u ( t ) , χ ( t ) , t . (2.114)
44 Moreover, a second-order Taylor series expansion about u demonstrates that if H is exists), the optimal control at t is interior, differentiable in u to second order (i.e., H uu
and u ( t ) is a weak variation, then, ( t ) ≥ 0 , H u ( t ) = 0, H uu
(2.115)
must be positive semi-definite. This condition equivalently ensures the noni.e., H uu
negativity of the second differential on J . Equation (2.115) is known as the LegendreClebsch necessary condition. Note that maximizing the performance index switches the inequality sign in equations (2.114)-(2.115). Satisfying the Euler-Lagrange equations and either the Minimum Principle or the Legendre-Clebsch condition forms two necessary conditions, and, thus, a necessary and sufficient condition for a minimum value of the augmented performance index.
Given that the Euler-Lagrange equations are
satisfied, if only the Legendre-Clebsch necessary condition is satisfied, then only a weak extremal [24] exists; thus, an unsatisfied Minimum Principle is still possible. In other ′ ′ . words, it is possible J strong < J weak An alternative test for a minimum is the Weirstrass condition [25]. The Weirstrass condition represents a more restricted version of the Minimum principle and is not discussed here. Regardless of the second necessary condition, however, optimal control yields only locally optimal solutions. Although different variations of the control result in different classifications of the extremals, variations of the state must still remain infinitesimally small.
45
3. LOW-THRUST TRANSFER ALGORITHM
In comparison to impulsive transfers, low-thrust transfers are typically characterized by increased fuel economy at the expense of increased time of flight. However, an additional difficulty is the determination of the control history, since an initially unknown engine thrust direction and magnitude history is required continuously along the path. For this investigation, invariant manifolds are used as an engine stopping condition, thus incorporating a thrust-free arc along the transfer.
Variable specific impulse (VSI)
engines are selected due to the advantage of a variable thrust profile. More robust numerical convergence is also observed. To derive a steering law, optimal control theory is applied, resulting in a primer vector law as detailed by Lawden [27] and Marec [32]. To meet the objectives of this analysis, it is not necessary to establish the complete twopoint boundary-value. (Not all natural boundary conditions are evaluated.) In fact, a hybrid direct/indirect optimization scheme is employed.
Nevertheless, many of the
benefits of optimal control theory are still employed. A direct minimization of the performance index is attempted, while components of optimal control theory are utilized to parameterize a control law. The successful numerical computation of low-thrust transfers begins with a direct shooting routine, using sequential quadratic programming (SQP) to update a search vector. The sensitivities associated with this search vector and the nonlinear constraints are evaluated numerically. Because the costates are so numerically sensitive, a useful mapping process, the adjoint control transformation (ACT), is employed, transforming several initial guess parameters into more physically realizable quantities. Despite this mapping, the initial guess parameters are still extremely sensitive. The state and the costate equations of motion are nonlinear and the propagation times are long (usually 20
46
days or more). To combat cumbersome blind initial guess strategies, a global set of initial conditions over a bounded range is undertaken first.
This “shotgun” step is
completed once to obtain a useful set of initial conditions for a desirable transfer. Once accomplished, iteration with the local SQP routine proceeds to satisfy the kinematic boundary conditions and achieve a stationary value, i.e., a local maximum, of the performance index.
3.1 Engine Model
In developing low-thrust transfers, the most commonly applied engine model is the constant specific impulse (CSI) engine. Such engines operate at maximum and minimum (typically zero) thrust levels, at a fixed engine impulse, and yield the control law and switching function detailed by Lawden [27]. These engine models have been exploited in the investigations by Russell [44], Coverstone-Carroll and Williams [40], Sukhanov and Eismont [41], Hiday [54], and many others. CSI engines find practical application in many current low-thrust missions, for example ESA’s Smart-1 mission, that uses a solar powered ion engine as shown in Figure 3.1 [55].
The determination of transfer
trajectories
Figure 3.1 CSI Engine Example – Smart-1 Ion Engine
47
transfer trajectories relies on the use of a switching function to govern engine on-and-off times, frequently resulting in “bang-bang” control. In practice, the switching behavior is difficult to predict given a set of initial conditions, and thus, a pre-determined switching structure is often specified. Forcing a switching structure, such as a “thrust-coast-thrust” profile implies that the two switching times must be incorporated as additional design variables. Other investigators [45] seek to remove any forced switching behavior, and predict the shape of the switching structure with the adjoint control transformation, and an intuitive selection of the initial conditions. Such cases yield trajectories that include several switching times, but are more sensitive to numerical convergence. There are also other studies that reformulate the optimal control problem as a minimum time problem, and assume an “always on” profile for the CSI engine. Variable specific impulse (VSI) engines, as the name implies, operate under a modulating specific impulse assumption. An example of such an engine is the Variable Specific Impulse Magnetoplasma Rocket (VaSIMR), seen in Figure 3.2 [52]. While the en
Figure 3.2 VSI Engine Example - VaSIMR Rocket
engine operates between a maximum and minimum power level, the varying impulse, Isp, allows a thrust magnitude that potentially varies (≥ zero thrust) along the entire transfer path. The assumption that the specific impulse, Isp, is unbounded allows the switching function to completely vanish from the derivation of the control law.
Without a
switching function present, additional variables in the search vector are not necessary.
48
Several applications that compare the two engine models, such as Ranieri and Ocampo [36] and Sakai [57], note that VSI engines are capable of producing more efficient transfers than CSI engines. In practice, the variable thrust profile also results in more rapid numerical convergence than a fixed thrust magnitude. Variable specific impulse engines are used here in response to the numerical sensitivity issues in the highly nonlinear CR3BP; a potential improvement in fuel economy is also noted.
3.2 Control Law Derivation
A finite-thrust transfer, one with the greatest economy of fuel, and hence a minimum expenditure of mass, is sought. To achieve this goal, a performance index is specified in terms of a Mayer problem in the calculus of variations. The first necessary conditions are derived using the Euler-Lagrange equations, as well as some of the natural boundary conditions assuming a force model consistent with the CR3BP and variable specific impulse (VSI) engines. Pontryagin’s Minimum Principle results in a steering law that is Lawden’s primer vector, to be implemented as a steering parameterization at all times on the extremals. The Minimum Principle, coupled with the Euler-Lagrange equations, comprise a partial representation of the necessary and sufficient condition for a weak relative minimum transfer trajectory. Due to the presence of VSI engines and a range on specific impulse that is assumed to be unbounded, the switching function to govern the thrusting profile vanishes upon inspection of the necessary conditions. The spacecraft is assumed to originate on a circular parking orbit of fixed radius with respect to the Earth, with the initial angular position, θ 0 , utilized as a free variable. The spacecraft then thrusts continuously until an insertion onto the stable manifold tube occurs. At this point, the engines turn off and the spacecraft asymptotically converges into the desired trajectory by means of the stable manifold associated with the target periodic orbit. The transfer time and the insertion point along the manifold tube are also utilized as free variables. Note that the optimal control formulation does not include a derivation of the natural boundary conditions corresponding to the variables associated with the initial
49
departure angle on the parking orbit, or the variables that specify the state location on the stable manifold. The derivation of the control law begins by revisiting the general form of the performance index, equation (2.104). Since the goal in the optimal control problem is to maximize the final mass of the spacecraft arriving on the stable manifold, the function is only dependent on endpoint conditions, and the resulting Mayer problem is max J = k ( m f ) or min J = − k ( m f ) ,
(3.1)
where k is currently an undefined constant that rescales the problem. This constant will serve a useful purpose later and does not alter the maximum mass goal. Six controls, i.e., the scalar elements of the vector uc , are used in the problem to ensure a stationary value of the performance index: three thrust directions, uˆT , the thrust magnitude, T, the engine power, P, and a slack variable, σ , associated with maintaining an engine power value within the prescribed bounds, uˆT T uc = . P σ
(3.2)
The performance index is subject to dynamical, control, and endpoint constraints. The dynamical constraints are comprised of the equations of motion in cartesian coordinates and conservation of mass along the entire trajectory, defined as, v ɺ X = f n X , uc , t = g ( r ) + h ( v ) + ( T m ) uˆT , -T 2 2 P
(
)
where X is the nx1 state vector (n = 7), defined by r X = v . m
(3.3)
(3.4)
The time invariant force field is defined consistent with the moon-centered CR3BP and the associated dynamical equations of motion (equation (2.26)). The expressions are
50 decomposed into functions of position, g ( r ) , and velocity, h ( v ) as defined in equations (2.27)-(2.28) . The control constraints require the thrust direction, uˆT , to be fixed on the unit sphere, and engine thrust power to be bounded, i.e., uˆT T uˆT = 1 ,
(3.5)
P = Pmax sin 2 σ ,
(3.6)
where σ is the slack variable used to ensure to that 0 ≤ P ≤ Pmax . Note that this model excludes the use of thrust and Isp constraints. Since the state vector is comprised of 7 elements (n = 7), a total of 2n + 2 = 16 boundary conditions are necessary to formulate the complete two-point boundary-value problem. The endpoint constraints are specified at the boundaries of the trajectory, with the spacecraft originating at a planar circular parking orbit (recall r0 is fixed, but θ 0 is free), and initial constraints (equation (2.106)) evaluated as,
ϕ n +1 ( t0 , X 0 )
x ( t0 ) − r0 cos θ 0 y ( t0 ) − r0 sin θ 0 z ( t0 ) I S v x ( t0 ) − − µ⊕ r0 sin θ 0 + ω r0 sin θ 0 = v y ( t0 ) − µ⊕ r0 cos θ 0 − I ω S r0 cos θ 0 vz ( t0 ) m ( t0 ) − m0 t0
( (
) )
= 0,
(3.7)
Equation (3.7) adds (n + 1) boundary conditions, with 8 yet required. The spacecraft terminates at a p x 1 target vector, ψ p X f , where p = 6, on the stable manifold,
( ) r ( t ) − r (θ ( X ) = v ( t ) − v (θ
,τ M ) = 0, (3.8) f , τ ) f M M M where the position and velocity states, rM (θ M ,τ M ) , and vM (θ M ,τ M ) are states along the
ψp
f
M
M
manifold tube parameterized by the free angle-like variable θ M , and the free time-like
51 variable τ M , as shown in Figure 3.3.
The angle-like variable specifies the stable
manifold trajectory given a fixed-point (red) along the libration point orbit (pink), and the time-like variable specifies the state at a given time along a specified stable manifold trajectory (blue). Since ψ p X f provides p = 6 boundary conditions, two more
( )
boundary conditions (one associated with the final mass, and one associated with the final time) are required for the TPBVP.
At this point, the augmented cost function is introduced by adjoining the two Lagrange multipliers vectors, ϖ T , and υ Τ to the end point constraints, a Lagrange multiplier vector η to the control constraints, and adding the result to the original performance index, tf
)
(
ɺ J ′ = Θ + ∫ Hˆ − χ T X dt , t0
(3.9)
where, Θ = km f + ϖ T ϕ X ( t0 ) , t0 + υ Τψ X ( t f ) ,
(3.10)
Hˆ = H + η1 ( uˆT uˆ − 1) + η2 ( P − Pmax sin 2 σ ) .
(3.11)
(
(
)
)
Recall that if all the constraints are satisfied, minimizing the augmented performance index is identical to minimizing the actual performance index.
The extended
Hamiltonian, Hˆ , is then the path-wise analog of adjoining Lagrange multipliers to the endpoint function, Θ . (See Hull [25] for further details.)
(
XM θMi ,τMj+2
(
XM θMi+1,τMj+2
) )
(
XM θMi ,τMj+1
(
)
XM θMi+1,τMj+1
(
XM θMi ,τMj
) (
XM θMi+1,τMj
)
Figure 3.3 Behavior of θΜ and τM Along the Stable Manifold Tube
)
52
To satisfy the first necessary condition for a maximum value of the performance index, it is required that the differential of the augmented cost function, dJ ′ , vanish. The expression dJ ′ = 0 leads to the Euler-Lagrange equations (equations (2.111)-(2.113)), applied to the extended Hamiltonian, and the natural boundary conditions (equations (2.115),
χɺ = − Hˆ XT ( X , uc , χ , t ) ,
(3.12)
Hˆ uTc X , uc , χ , t = 0 ,
(3.13)
(
)
( ) = Θ ( X ( t ) , t , X ( t ) , t ,ϖ ,υ ) , = −Θ ( X ( t ) , t , X ( t ) , t ,ϖ ,υ ) .
χ 0 = −ΘTX X ( t0 ) , t0 , X ( t f ) , t f ,ϖ ,υ ,
0
χf Hf
T Xf
0
T tf
0
0
f
0
f
f
f
(3.14) (3.15) (3.16)
The natural boundary conditions provide the additional 2 boundary conditions, since n +1− p = 2.
Solving the Euler-Lagrange equations requires the formation of the
Hamiltonian, as expressed in equation (2.108). Since the performance index is a Mayer problem, the Hamiltonian becomes, ɺ H = χ T X = χ rT v + χ vT g ( r ) + h ( v ) + (T m ) uˆT − χ m (T 2 2 P ) .
(
)
(3.17)
Equation (3.11) is then used to evaluate equation (3.12), yielding the costate equations of motion,
χɺ r = − B ( t ) χ v ,
T
(3.18)
χɺ v = − χ r − CT χ v ,
(3.19)
χɺ m = (T m 2 ) χ rT uˆT ,
(3.20)
where B(t) and C are defined in equations (2.34)-(2.35). Note that these matrices are ∂g ( r ) ∂h ( v ) simply and , respectively. Expanding equation (3.13) yields ∂r ∂v T
∂Hˆ = χ v T m + 2η1 uˆT = 0 , ∂uT
(3.21)
53
T
∂Hˆ T = χ v uˆT m − χ mT P = 0 , ∂T
(3.22)
T
∂Hˆ 2 2 = χ mT P + η2 = 0 , ∂P
(3.23)
T
∂Hˆ (3.24) = −2η2 Pmax sin σ cos σ = 0 . ∂σ Equation (3.21) implies that either χ v is parallel to uˆT , T and η1 are both zero, or χ v and η1 are both zero. Since the latter two cases will rarely apply in any practical problem, χ v is assumed to be always parallel to uˆT . As a result, two possible solutions for uˆT emerge, uˆT = ± χ v χ v . (3.25) Equation (3.22) can also be solved directly to yield, χ v P T= . χ mm
(3.26)
From equations (3.23)-(3.24), either η2 , cos σ , or sin σ must be zero. Inspection of equation (3.6) suggests the possibilities: (i) if cos σ = 0, then P = Pmax; (ii) if sin σ = 0, then P = 0; and, (iii) if η2 = 0 , then 0 ≤ P ≤ Pmax . Next, the remaining two natural boundary conditions (n + 1 – p = 2) are determined via equations (2.113), −∂Θ ∂X 0 = χ 0 → χ rT0 ∂Θ ∂X f = χ f → χ rTf
T
χ vT
T
χ m = − ϖ rT ϖ vT ϖ m ,
0
0
χ vT
f
χ m = υ rT υ vT
Hf = 0.
T
f
T
k ,
(3.27) (3.28) (3.29)
The constant, k, in equation (3.28) serves as the only non-trivial value determined from equations (3.27)-(3.28), and the constraint on the final value of the Hamiltonian in equation (3.29) provides the last boundary condition. The unknown initial mass costate,
χ m , may now be completely removed from the problem. Observing that the mass 0
costate monotonically increases in equation (3.20), that is, any initial value, χ m0 = ϖ m ,
54 allows the final costate value χ m f to approach the positive constant value k. Thus, the initial mass costate may be arbitrarily fixed as unity,
χm = 1 . 0
(3.30)
Although equation (3.29) is the final condition required for a fully specified TPBVP, it will later be ignored when a hybrid direct/indirect solution method is established. Pontryagin’s Minimum Principle is applied as a second necessary condition, and thus, the sufficient condition for a local maximum. For this specific problem, equation (2.114) is reduced to the following: χ v T (T ** m ) uˆT** − χ m (T **2 2 P** ) ≥ χ vT ( (T m ) uˆT ) − χ m (T 2 2 P ) .
(
)
(3.31)
It is clear from the above expression that a positive sign on uˆT in equation (3.25) must occur for the inequality to achieve a maximum value, resulting in the definition of p , the primer vector. The associated primer vector control law, as observed by Lawden [27], is expressed,
uˆT = λv λv ≡ p p .
(3.32)
The structure of this control law implies that the thrust is always in a direction along the primer vector. A physical interpretation of this optimal control law is that the thrust acceleration is always directed toward a neighboring point (also in motion) subject to the same gravitational field and thrust acceleration as the spacecraft. The observations (i)-(iii) on the engine power, P, resulting from equations (3.23) and (3.24) are further reduced by substituting equation (3.26) and equation (3.32) into equation (3.31),
χ v2 P** χ v2 P ≥ , 2 χ m m 2 2 χ mm 2
(3.33)
that yields two possible values of P, P = Pmax, χ m ≥ 0 ,
(3.34)
P = 0, χ m < 0 ,
(3.35)
But, since the first necessary conditions identified a value χ m that begins at one and monotonically increases, equation (3.35) is never possible for this problem.
Thus,
55
equation (3.34) is exclusively employed, so the engine will always operate at maximum power, Pmax. The control law for the power, equation (3.34), automatically satisfies equations (3.23)-(3.24) when cos σ is always zero.
This requirement also reduces
equation (3.26) such that, T=
χ v Pmax . χ mm
(3.36)
Four results that are that necessary to ensure a local maximum value of the performance index emerge, and are summarized as follows:
1. The thrust direction is always tangential to the primer vector, i.e., uˆT = χ v χ v ≡ p p . 2. There is no requirement for a switching function on any of the controls, the engine always operates at maximum power, i.e., P = Pmax . 3. The initial mass multiplier always equals unity, χ m0 = 1, and monotonically increases to reach the arbitrary final mass multiplier χ m f = k. 4. The thrust magnitude is always defined as T=
χ v Pmax . χ mm
These four conditions ultimately comprise the “indirect” components in the formulation of a hybrid direct/indirect numerical solution.
3.3 Adjoint Control Transformation
Not surprisingly, one of the most difficult aspects in obtaining a solution to a trajectory optimization problem that uses elements of optimal control theory and the associated TPBVP, is generating an accurate initial guess for the costates. Low-thrust problems are also typically characterized by long propagation times, rendering the problem even more sensitive to initial conditions. Typical strategies to address the initial guess dilemma
56
often solve several smaller sub-problems [29-30] or use analytical results [44]. An alternate approach, first investigated by Dixon and Biggs [58], introduces physical control variables and their derivatives as an estimate of the initial costates. These control variables allow exploitation of physical intuition to produce a guess for the values of the initial costate variables. Such insight can reduce the problem sensitivity. Consider a reference frame centered at the spacecraft, defined by the unit vectors vˆ − wˆ − hˆ . The vˆ -axis of this frame is aligned with the relative velocity vector, v . The hˆ -axis is aligned with the instantaneous angular momentum vector, h . Finally, the wˆ axis is defined to complete a right-handed system. These unit vectors, and associated time derivatives that create this frame, are defined as, r ×v v ˆv = , hˆ = , wˆ = hˆ × vˆ , r ×v v
(3.37)
vɺˆ = vɺ v − vvɺ v 2 , ɺ ɺ hˆ = h h − hhɺ h 2 , ɺ wɺˆ = hˆ × vˆ + hˆ × vɺˆ.
(3.38)
Given a vector and its time derivative, the following relationships are used to fully determine equation (3.38),
vɺ = v ⋅ vɺ v , hɺ = h ⋅ hɺ h.
(3.39)
As is apparent in Figure 3.4, two spherical angles, α , β , and their time derivatives αɺ and βɺ specify the orientation of the thrust direction relative to this frame, uˆTvwh , and also the time derivative of the thrust direction, uɺˆTvwh , that is, uˆTvwh = [cos α cos β uˆɺTvwh
sin α cos β
sin β ] , T
−αɺ sin α cos β − βɺ cos α sin β = αɺ cos α cos β − βɺ sin α sin β . βɺ cos β
(3.40)
(3.41)
57
However, since the equations of motion are integrated in the cartesian, moon-centered rotating frame (with unit vectors iˆ − ˆj − kˆ ), a rotation matrix, R, is required to transform the thrust direction, uˆTvwh (and uɺˆTvwh ), iˆ ⋅ vˆ iˆ ⋅ wˆ iˆ ⋅ hˆ R = ˆj ⋅ vˆ ˆj ⋅ wˆ ˆj ⋅ hˆ , ˆ ˆ ˆ ˆ k ⋅ vˆ k ⋅ wˆ k ⋅ h
(3.42)
iˆ ⋅ vɺˆ iˆ ⋅ wɺˆ iˆ ⋅ hɺˆ ɺˆ ɺ ɺ ɺ ˆ ˆ ˆ ˆ ˆ R = j ⋅v j ⋅w j ⋅h , kˆ ⋅ vɺˆ kˆ ⋅ wɺˆ kˆ ⋅ hˆɺ
(3.43)
uˆTijk = RuˆTvwh ,
(3.44)
ɺ uˆ + Ruɺˆ . uɺˆTijk = R Tvwh Tvwh
(3.45)
The direction, uˆTijk , denotes the thrust direction as expressed in terms of unit vectors in the moon-centered, rotating frame. The definition of the primer vector, p , from equation (3.32), is employed to parameterize the velocity costate vector, p = χ v = χ v uTijk ,
(3.46)
where χ v is the magnitude of the velocity costate, χ v = χ v . The equation of motion for the velocity costate, equation (3.19), is directly involved in parameterizing the position costate vector,
χ r = − χɺ v − CT χ v .
(3.47)
The derivative of the velocity costate vector, χɺ v , is available by differentiating equation (3.46), and substituting the result into equation (3.47), to yield an expanded version of χ r ,
χ r = χɺ v uT + χ v uɺT − CT χ v . ijk
ijk
(3.48)
58 The magnitude of the velocity costate time derivative vector, χɺ v , is approximated by assuming an initial value of the Hamiltonian, H0 = 0, substituting equation (3.36) into equation (3.17), and rearranging,
(
))
1 χ v uɺTijk v + χ vT Cv − χ vT g ( r ) + h ( v ) . uTijk
χɺ v =
(
(3.49)
Additionally, equation (3.36) is used to parameterize χ v in terms of the thrust, T. Thus,
(
)
an important mapping sequence M α , αɺ , β , βɺ , T → ( χ r , χ v ) is now available, where it is noted that the dimension of the initial conditions in the mapping are reduced by one through the assumption on H0. (A similar implementation of this step is provided by Senent et al. [28].) Rather than an initial guess that requires explicit values of the position and velocity costates, the physically meaningful quantities α , αɺ , β , βɺ , and T , determined via equations (3.36)-(3.49), are a practical alternative.
hˆ
βɺ β
α
vˆ
uˆT
wˆ
αɺ
Figure 3.4 Velocity Reference Frame
59
3.4 Numerical Solution via Direct Shooting: A Local Approach
An algorithm is constructed to iteratively solve the optimization problem by direct shooting, while incorporating components of the indirect TPBVP. For the actual numerical process, a new state vector, Y ( t ) is defined with elements that correspond to the states and the costates. The governing differential equations then appear in the following form, v rɺ ɺ g ( r ) + h ( v ) + (T m ) uˆT v −T 2 2 P mɺ ɺ . Y (t ) = = f Y = (3.50) T ɺ − B t χ () v χr χɺ − χ r − CT χ v v 2 T χɺ m (T m ) χ r uˆT Note that the new vector Y , is comprised of 14 states. However, for the transfer problem, an initial numerical search vector, S ′ is composed of 9 design variables: 5
( )
variables for the adjoint control transformation; one variable for the initial parking orbit angular position, θ 0 ; one variable corresponding to the final time, tf ; the angle-like variable, θ M , along the manifold tube; and, the time-like variable, τ M , along the manifold tube. The sensitivities of θ 0 , θ M , and τ M are acquired completely numerically, and are not available in closed form. The inclusion of these three additional search variables in S ′ is not formulated in the model of the indirect TPBVP. If these additional variables ( θ 0 , θ M , τ M ) are to be incorporated in the formulation, additional criteria are required to produce additional natural boundary conditions. In practice, the absence of the other natural boundary conditions in the TPBVP (including equation (3.29)) is offset by direct iteration on the performance index. The initial search vector, S ′ is defined,
60
θ0 α αɺ β S ′ = βɺ , (9 x 1), T 0 tf θ M τ M
(3.51)
where the mapping, M (α , αɺ , β , βɺ , T0 ) → ( χr , χv ) is supplied only once to determine the initial value of the elements in the actual (10 x 1) search vector, S , 0
0
θ0 χ x0 χy 0 χ z0 χ v x0 S= , (10 x 1). χ v y0 χ vz0 t f θM τ M
(3.52)
Subsequently, an iteration process is incorporated to update the search vector such that the kinematic boundary conditions are all satisfied. In the numerical algorithm, the kinematic boundary conditions specify the entire constraint vector, c , i.e., r ( t f ) − rM (θ M ,τ M ) = 0 , (6 x 1). c =ψ X f ,t f = (3.53) v ( t f ) − vM (θ M ,τ M )
(
)
A nonlinear programming algorithm, based on medium-scale Sequential Quadratic Programming (SQP), using “fmincon” in MATLAB®, updates the values of the 10 design variables in S to resolve any potential constraint violations in c . In this particular scheme, the optimizer solves a quadratic programming (QP) sub-problem every iteration,
61
and computes a quasi-Newton approximation of the Hessian of the Lagrangian using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) formula [59]. A medium-scale algorithm is used due to the nonlinear constraints. numerically.
All gradient information is approximated
This problem might benefit from analytical gradients, although
recalculation is required any time the problem is slightly reformulated. This solution scheme is termed a hybrid direct/indirect method since the performance index is minimized directly, while elements of the indirect optimal control problem are still employed. Thus, formulation of all the natural boundary conditions is ignored, and a stationary value of the performance index is detected numerically, effectively replacing a condition that could only be satisfied through the solution of the complete TPBVP. An objective function, Q, is supplied to the numerical optimizer to be the equal to the performance index, – J, i.e., Q = − J = −m f .
(3.54)
The flow chart in Figure 3.5 represents the details of the numerical algorithm that is implemented to determine the locally optimal solutions.
62
{
Guess S ′ = θ 0
α 0 αɺ 0 β 0 β 0 T0 t f
(
) (
θM τ M }
T
Mapping M θ 0 , α 0 , αɺ 0 , β 0 , βɺ0 , T0 → χ r , χ v 0 0
)
using
equations (3.35)-(3.34)
SQP ROUTINE
SQP ROUTINE
{
S = θ0
χx
0
χy
0
χz
0
χv
x0
χv
y0
χv
z0
tf
θM τ M
}
T
Initialize Y0
Propagate equations
(3.50)
forward to time tf, objective function,
Using equation (3.8),
Q, (equation(3.54)),
compute
terminal
is computed,
constraint
values,
ψ
Q = − J = −m f
f
Optimizer performs Is Q locally minimal and ψ f ≤ tolerance ?
several small updates on S to estimate numerical gradients,
YES
NO
and then SQP
algorithm updates S
Optimal S is stored for plotting and future continuation processes
Figure 3.5 Numerical Algorithm – Direct Shooting Method via SQP
63
3.5 Shotgun Method for Initial Conditions: A Global Approach
Despite the added benefits of the adjoint control transformation, many initial guesses are typically required to determine a vector S within the convergence radius of the local SQP optimizer. To avoid a cumbersome manual search for appropriate initial conditions, a numerical, low-fidelity, automated process is developed to establish a “global” set of initial conditions. Such a step is intended to precede any SQP optimization attempt. Upper and lower bounds on each design variable in the initial search vector, S ′ , are selected, except the variables t f , θ M , and τ M . These three design variables are selected differently. For the first six elements of S ′ , the variables are each within the bounds corresponding to the individual global ranges, i.e., θ 0LB , θ 0UB α LB , αUB αɺ LB , αɺUB β LB , βUB βɺLB , βɺUB S g′ = T , T , 0LB 0UB t fh h =1,...,d θ M i i =1,...,a τ M j j =1,...,b
(3.55)
where the subscript “LB” denotes the lowest possible value on design variable and “UB” denotes the highest possible value. Thus, the first six design variables may fluctuate. The propagation time is always initially selected as a maximum value, t f fixed = td (typically 50 days), although samples over the entire range are subsequently tested. The full range of the angle-like parameters, θ M i =1,...,a (of length a), and time-like parameters,
′ (of length b) are always initialized. (See Figure 3.3.) The first 7 variables in S j =1,..., b allow Y to be propagated forward for the fixed length specified for the maximum time
τM
64
interval, t f fixed . Note that any discretized point along the stable manifold is a potential “match” point, since every combination of the angle-like and time-like variables θ M i =1,...,n and τ M j =1,...,m is utilized. From this pool of starting values, representative combinations are sampled. All elements in S g′ are first fixed; t fh is sampled through each time element, h = 0,..., d , over a propagation of Y . Along each time step, t fh , each possible combination of the current kinematic boundary condition errors along the stable manifold tube, ψ hij , are computed,
ψ hij
( ) ( )
( (
r t f − rM θ M ,τ M h i j = v t − v θ ,τ M Mi Mj fh
) . )
(3.56)
For each initial state, the terminal error vector along the trajectory, ψ hij , with the lowest
norm is then stored. Another candidate S ′ within the desired bounds of S g′ is then propagated, and the process is repeated. Once a pre-determined number of S ′ combinations is propagated, only the five search vectors with the lowest error norm, ψ hij , are selected for attempted convergence in the local SQP algorithm represented in Figure 3.5. Note that once the indices, h, i, j, that identify the lowest error vector, ψ hij ,
are isolated, a single S ′ vector is available for use in the local method, (i.e., tf = th ,
θ M = θ M , τ M = τ M ). This process is similar to the first level of a global parameter i
j
optimization scheme that establishes a “population”, such as a genetic algorithm. It is also analogous to propagating a representative set of “open loop” initial conditions since no iterative procedures are occurring, merely propagation. A random number seed is incorporated to propagate a specific population among the bounded variables that is repeatable. Thus, this step is loosely regarded as the first level of a “hybrid” parameter optimization sequence (not to be confused with the hybrid direct/indirect trajectory optimization method), where the best condition(s) are supplied to the “closed-loop” local SQP method. For an example involving a 12-day halo orbit (section 4.2.1), sample
65
populations for α , β , αɺ , and βɺ appear in Figure 3.6.
Note that a uniform random
distribution is employed.
Bounds: α 0 = [0, 0.2] , β 0 = [ −0.3, 0.3] , αɺ 0 = [ −0.02, 0.02] , βɺ0 = [ −0.3 0.3]
Population Size: 500
Figure 3.6 Example Population Parameters for a 12-Day Halo Orbit
66
4. MISSION APPLICATIONS
The free final time targeter in equation (2.80) is used to generate families of periodic orbits. These families include L1 and L2 halo orbits, L1 and L2 vertical orbits, and L2 “butterfly” orbits. All of the families of orbits are generated to meet a basic altitude constraint. The trajectory design process involves multiple spacecraft in combinations of potentially different orbits to meet the design objective, that is, nearly continuous line-ofsight coverage of the lunar south pole. Analysis of periodicity and orbital stability yields individual orbit selection criteria. Many of these orbits have already been incorporated into a more rigorous coverage and stationkeeping analysis in a full ephemeris model [46]. Once identified as a viable option, transfers to these orbits are required. Investigations of transfers from Earth into such orbits are limited. Low-thrust transfer trajectories to these orbits are the focus here. The algorithm for the design of low-thrust transfers is applied to a mission scenario that requires the delivery of a spacecraft into an orbit for potential lunar south pole coverage.
All transfers are propagated within the P2-centered Earth-moon CR3BP
(equations (3.3)) to reduce numerical sensitivities in the vicinity of the insertion point. Since the invariant manifolds, like the associated orbits, are generated in the barycentric frame associated with equations (2.15)-(2.17), a coordinate transformation aligns the axes to the P2-centered frame.
As assumed in Chapter 3, a spacecraft originates at an
unspecified departure angle in a circular parking orbit about the Earth. It then reaches an unspecified location along the stable manifold tube, and inserts into a manifold trajectory. The invariant manifold coast includes a significant length of time that the engine power remains off.
The local hybrid direct/indirect optimization process results in locally
optimal transfer trajectories. When necessary, the global “shotgun” method is used to produce initial conditions within the convergence radius of the direct shooting routine.
67
4.1 Orbits for Line-of-Sight Lunar South Pole Coverage (CR3BP)
Due to the potential existence of frozen volatiles [60,61], one current location of interest for future space exploration is the region near the lunar south pole. This goal has been identified in the President’s Vision for Space Exploration announcement in January 2004, as NASA indicates that water ice at the lunar poles may help facilitate exploration of the solar system [62]. NASA’s Exploration Communication and Navigation Systems (ECANS) Team, specifically the Lunar Communications and Navigation Systems (LCNS) group, is interested in spacecraft architectures for communications with ground stations on the lunar surface. Such a ground station on the moon would benefit from a system of satellites that are always within direct view of the Earth and that provide constant communications between the lunar surface and the Earth. Various CR3BP orbits are potentially applicable in mission design of lunar relay communication satellites for lunar coverage due to the fixed geometry in the rotating frame and line-of-sight capability. For example, L1 and L2 southern halo orbits possess a line-of-sight with the lunar south pole over the majority of the orbital period, and a line-of-sight with the Earth for the entire orbital period. Additional orbital information on the periods and stability indices aids in the selection of specific orbits.
4.1.1 Three-Dimensional Periodic Orbits in the CR3BP
Similar to the two-dimensional, Lyapunov example from Section 2.2.4, periodic orbits are determined with the free final time targeting algorithm (equation (2.80)). Depending on the specific type of orbit, different control parameters are employed, and different states serve as targets to successfully converge on an initial condition that yields a periodic orbit. (Recall that the term “control”, used in this context, simply refers to parameters that may vary in targeting periodic orbits.) All orbits generated share a symmetry across the x-z plane. Thus, an initial state vector is aligned on the x-y plane, with nonzero elements in x , z , and yɺ only, such that,
68
X 0 = {x0
0 z0
0
yɺ 0
0} . T
(4.1)
The variational targeting equation (equation (2.80)), is now reduced to,
δx 0 ɺ δ z0 ∂X f ∂X f ∂X f δXf = X f , δ yɺ 0 ∂x0 ∂z0 ∂yɺ 0 M δτ
(4.2)
where M is a 3 x 4 matrix. The three elements of the target vector, δ X f , are selected depending on the type of orbit. From observation of equation (4.2), there are four possible controls. In a process similar to that for Lyapunov orbits, one nonzero initial parameter is be fixed, while the other two (and time) are allowed to vary as control parameters. Once a periodic solution is available, a method of continuation is applied to create a new orbit. The targeting process is repeated to generate successive orbits in the family. A description of the full targeting scheme is presented in Grebow [63].
4.1.2 Families of Orbits for Lunar South Pole Coverage
Families of orbits for potential application in the problem of lunar south pole coverage are obtained with lunar altitudes between 50 km and 100,000 km. The maximum bound is assumed as a communications instrument constraint; the minimum bound is selected arbitrarily to avoid a subsurface arc. Orbits within the acceptable range from the L1 and L2 southern halo orbit families appear in the moon-centered frame in Figure 4.1 [46,63].
The halo orbits (a term first used by Farquhar [64]), bifurcate from both the L1 and L2 Lyapunov family of orbits, and resemble a halo-shape about the moon when viewed from the Earth in the rotating frame. As previously stated, the orbits are particularly effective in the lunar south pole coverage problem since the motion is almost always within lineof-of sight to the Earth.
The family is composed of halo orbits that resemble the
traditional “halo” shape in addition to highly “elliptic”, near-rectilinear orbits with passage very close to the moon’s surface. For almost the entire period of motion, a spacecraft in any near-rectilinear halo orbit possesses a line-of-sight to the lunar south
69
pole.
Most recently, the halo orbit families have been thoroughly investigated by
Farquhar [64], Breakwell and Brown [17], Howell [18], and Gómez et al. [65]. Members of the southern L1 and L2 vertical orbit family are depicted in Figure 4.2 [46,63]. The motion consists of a doubly symmetric, “figure-8” shaped pattern when viewed in the y-z plane. These orbits occur near the libration points. The existence of these orbits was predicted by Moulton in 1920 [14], and have also been studied recently by Dichman et al. [21]. Large amplitude L1 vertical orbits terminate when they become exactly vertical, while large amplitude L2 vertical orbits encompass both primaries (although these trajectories are not included due to the mission constraints). The orbits also possess the characteristic of bending toward both the north and south poles of the moon, a favorable trait for maintaining line-of-sight over a pole. An additional family also includes orbits that remain in view of the lunar south pole for significant intervals of time. Some of these orbits possess characteristics similar to the near-rectilinear halo orbits. The orbits bifurcate from a 6-day near-rectilinear L2 halo orbit and might be described as a “butterfly” shape. (See Figure 4.3 [46,53]). Comparable motions around the smaller primary have been documented by Robin and Markellos [66]. Similar to vertical orbits, the motion in a butterfly orbit resembles a “figure-8” shape, however, these orbits wrap around both the near and far side of the moon, such that a direct line-of-sight to the lunar south pole exists for nearly the entire orbital period.
70
Moon To Earth
To Earth
Moon
Figure 4.1 Southern Halo Orbit Families: Earth-Moon L1 (Orange) and L2 (Blue); Moon Centered, Rotating Reference Frame
71
To Earth
Moon
Moon
To Earth
Figure 4.2 Vertical Orbit Family of Interest: Earth-Moon L1 (Magenta) and L2 (Cyan); Moon Centered, Rotating Reference Frame
72
Moon
To Earth
Moon
To Earth
Figure 4.3 Southern L2 Butterfly Orbit Family; Moon Centered, Rotating Reference Frame
73
4.1.3 Mission Orbit Selection Criteria
The time to complete one full period is used as a design parameter for orbit selection to be applied in the coverage problem. Let the maximum excursion distance identify a particular orbit, as indicated using a halo orbit in Figure 4.4 [46,63] (right). Maximum excursion distance is defined as the maximum x-distance for each orbit in the moon centered, rotating frame. In Figure 4.4 [46,63] (left), orbital periods are plotted against maximum excursion distance during initial design selection. Commensurate orbits are sought to phase multiple spacecraft for complete line-of-sight coverage. One such region might consist of orbits in L1 and L2 halo families sharing periods between 7.9 and 12.2 days. An example that exhibits feasible south pole coverage is a 12-day L1 and 12-day L2 halo orbit combination, illustrated by the black dashed line in Figure 4.4. Another region with commensurate combinations consists of orbits with a ratio of periods equal to 2:1, that is, one period is exactly twice that of the other. Note that L2 halo orbits with periods between 6.0 and 7.2 days exhibit this behavior with the entire L2 butterfly orbit family. This is not actually surprising when the shapes of the orbits are viewed in Figure 4.1 and Figure 4.3. An example from this region consists of a 14-day L2 butterfly orbit and a 7day L2 halo orbit combination, as noted by the two red dashed lines in Figure 4.4. The information in Figure 4.4 serves as a basis for the determination of many other commensurate orbit combinations that lead to complete south pole coverage. Also useful for design purposes is the stability index, S . I .
The stability index,
corresponding to one orbit period, IP, is defined as,
S.I . =
1 1 λmax + , λmax 2
(4.3)
where λmax is the maximum eigenvalue from the monodromy matrix, Φ ( IP + t , t ) , computed at the end of one revolution. A stability index of one indicates a stable orbit, whereas stability indices greater than one reflect instability. Of course, a large stability index indicates a divergent mode that departs from the vicinity of the orbit very quickly. Generally, the stability index is directly correlated to the station-keeping costs and is inversely related to transfer costs. The stability indices for orbits from the various
74
families appear in Figure 4.5 as functions of the maximum excursion distance from the moon. In general, the stability index increases with maximum excursion distance from the moon. Once the orbits of interest have been selected to yield the appropriate coverage of the lunar south pole, transfers to deliver the spacecraft into such orbits must also be available. An analysis of the coverage schemes is discussed in references [42,59]. Within the context of the multi-body problem, the stability index must be sufficiently large to produce stable manifold trajectories that arrive at the orbits during numerical simulations. As stability indices approach a value of one (a stable orbit), the manifolds become more difficult to produce numerically due to the increasingly stable behavior within the phase space. In some situations, this complexity is offset by increasing the initial orbital displacement distance, d, described in Section 2.3.3. In general, such orbits are excluded from this investigation. As a result, because the 6-day and 7-day halo orbits possess stability indices very close to one, they are not considered for use in this transfer scheme.
14-day L2 butterfly and 7-day L2 halo orbit combination
L1 Halo
z
Orbit Example x Moon 12-day L1 halo and 12-day L2 halo orbit combination
Max |x|
Figure 4.4 Period versus Maximum x-Distance from the Moon (Left); Definition of Maximum x-Distance (Right)
75
Figure 4.5 Stability Index versus Maximum x-Distance from the Moon
4.2 Optimal Transfers to the Earth-Moon Stable Manifold
For the design of low-thrust transfers to orbits for lunar south pole coverage, a number of specific scenarios are identified.
Due to the numerical sensitivities in the
determination of three-dimensional transfers in the CR3BP model, the design process for the low-thrust transfers includes the global and local solution method detailed in Chapter 3. At the initiation of the local SQP routine, a hybrid direct/indirect method reduces the objective function and satisfies the kinematic boundary conditions. The fixed dynamical and propulsion constants for all scenarios are listed in Table 4.1. An initial spacecraft mass is assumed at 1,500 kg, with VSI engines capable of delivering a maximum engine power of 10 kw. For all simulations with low-thrust acceleration terms, the differential equations that model the system are the P2-centered, i.e., moon-centered equations of motion, as detailed in equations (2.26)-(2.31), and in equation (3.3) with the associated
76
thrust terms. Additionally, these nondimensional equations include a scale factor on the initial spacecraft mass such that it is equal to 1 in the numerical process, thus avoiding scaling issues that arise from the magnitude of the characteristic Earth-moon mass, m*. Any simulations using the barycenteric equations of motion, e.g., computation of the orbits and the stable manifolds, are shifted into moon-centered coordinates as well. A fixed, circular orbit parking radius of 20,000 km is established for all low-thrust transfer examples. Of course, a wide range of Earth orbits may serve as the departure orbit. For convenience, it is assumed that the departure orbit is circular. Parking orbits at radii lower than 20,000 km result in long integration times and numerical sensitivity issues. Of course, these numerical computations can be offset by establishing a parking orbit at LEO and utilizing higher initial thrust values. But the larger radius was selected to demonstrate the capability of engine parameters previously examined [28]. Nevertheless, using the same solution method, a fuel optimal, planar, circle-to-circle transfer is determined, as seen in Figure 4.6, to achieve optimal orbit raising from LEO (200 km altitude; dotted red line in Figure 4.6) to the nominal departure orbit at radius 20,000 km (dotted green line in Figure 4.6). Obviously, the ∆v required to reach the 20,000 km radius departure is dependent on the initial parking orbit supplied from a launch vehicle. The value of 3.21 km/s for ∆vVSI represents the cost in terms of an equivalent maneuver magnitude using VSI engines to reach the nominal departure orbit. For comparison, a transfer arc modeled as a two-burn Hohmann transfer requires a ∆vH value of 3.09 km/s.
77
Table 4.1 Dynamical and Propulsion Constants Parameter
Value
Units
Circular Parking Orbit Radius
20000
km
m0
1500
kg
Pmax
10
kw
m*
6.04680403834987 x 1015
kg
l*
385692.5
km
t*
377084.152667039
sec
Earth GM
398600.432896939
km3/s2
Moon GM
4902.80058214776
km3/s2
Earth Radius
6378.14
km
Moon Radius ∆vVSI From LEO to Parking Orbit
1737.4
km
3.206418
km/s
3.088713
km/s
∆vH From LEO to Parking Orbit
Figure 4.6 Optimal Orbit Raising from LEO
78
4.2.1 Transfers to a 12-Day L1 Halo Orbit
4.2.1.1 “Short” Transfer
For the first demonstration of a free final time, optimal, low-thrust transfer, a 12-day L1 halo orbit is selected. As noted in Figure 4.1, there are two different 12-day L1 halo orbits that may be used. The lower z-amplitude orbit (Az = 13,200 km) is selected for this example to reduce the sensitivity issues on the out-of-plane costates. The orbit, and the associated stable manifold tube, are depicted in Figure 4.7. The stable manifold (green) is propagated (backwards in time) for two successive Earth flybys; each flyby is recognized during the simulation by the point at which the flight-path-angle changes sign. The manifold tube is then parameterized in terms of 50 state vectors, X s ( ti ) . Each X s ( ti ) (nfp = 50), calculated from equation (2.102), is propagated backwards as detailed in Section 2.3.3. Recall that each fixed point M is associated with one trajectory along the manifold. So, the angle-like parameterization, θM, tags a particular fixed point and, thus, a specific manifold trajectory. The time-like parameter, τM, corresponds to a timeindex along the specified manifold trajectory. In this case, 0 ≤ θ M ≤ 50 (50 fixed-points equally spaced in time), and 0 ≤ τ M ≤ 2000 (for each fixed point, a trajectory is composed of 2000 time elements). Once the optimization routine is initialized, a two-dimensional cubic spline is employed to represent the corresponding state along the tube; interpolation for states at any point along the tube is then available. The global method, as detailed in Section 3.5, initially generates a population of 100 uniformly distributed random initial search vectors, S ′ . In this case, only the initial condition with the lowest kinematic boundary condition error norm, ψ ( t f ) is retained. The manifold trajectory that is targeted as a result of this process is highlighted in blue in Figure 4.7. The local, hybrid direct/indirect routine completed 90 iterations for convergence to an insertion point with tolerances of 1 x 10-7 on position and velocity.
79
Earth
Initial Guess for Insertion Point Moon
Figure 4.7 Stable Manifold Tube for 12-Day L1 Halo Orbit (Green) and Target Reference Trajectory Along the Manifold (Blue)
The final low-thrust solution appears in Figure 4.8. The solid blue trajectory indicates periods of engine on-time, and the dotted line represents the fact that the trajectory has arrived on the stable manifold. The engines are off during an asymptotic approach to the orbit. During the powered phase, the spacecraft engine is “on” for a total of 25.96 days, and the remaining translunar coast on the stable manifold is 10.85 days, for a total transfer time of 36.82 days.
80
Earth Moon
Figure 4.8 Low-Thrust Short Transfer to a 12-Day L1 Halo Orbit
The trajectory exhibits the “spiral” structure that is commonly observed in low-thrust applications, as the spacecraft gradually builds up momentum to insert into the target trajectory. The position and velocity costate histories appear in Figure 4.9. Clearly, an initial out-of-plane component on the thrust direction is required to enable the planar
81
orbit to gradually shift to the inclined insertion point. This out-of-plane component also oscillates however, and all of the velocity costates increase in magnitude during the approach. In comparison, the position costates diminish with increasing time. It is also apparent from Figure 4.9 that the position and velocity costates maintain a relatively well-behaved and organized structure during the initial spiraling of the powered phase. But, further along the path, a distinctly nonlinear behavior is clear in the vicinity of the insertion point. Finally, note that the magnitudes of the costates along the converged solution are not intuitive, demonstrating the usefulness of the adjoint control transformation (ACT).
Figure 4.9 Position and Velocity Costate Time Histories for the 12-Day L1 Halo Orbit Transfer The performance of the mass costate, χm , appears in Figure 4.10. The spacecraft mass and thrust histories are also plotted in the figure. The monotonic increase in χm, observed upon inspecting equation (3.19), is clear in the plot. The spacecraft mass ultimately reaches 1139.52 kg, for a total fuel mass consumption of 360.47 kg. The thrust profile maintains an oscillatory structure due to its dependence on the magnitude of the primer vector in equation (3.36). The thrust magnitude initially peaks at 2.32 N, but equation
82
(3.36) also indicates that, due to the dependence on mass, a decrease occurs as propellant is gradually expelled, reaching values as low as 0.94 N.
These observations are
summarized in Table 4.2.
Figure 4.10 Time History of Propulsion Related Parameters for the 12-Day L1 Halo Orbit Transfer
Table 4.2 12-Day L1 Halo Orbit Transfer Data Summary Parameter
Value
Units
Powered Time
25.965818
day
Coast Time
10.849792
day
Total Transfer Time
36.815610
day
Final Spacecraft Mass
1139.524870
kg
Propellant Mass Consumed
360.4751292
kg
Tmax
2.317819
N
Tmin ∆v
0.936037
N
3.020521
km/s
83
4.2.1.2 “Long” Transfer
Using the same 12-day L1 halo orbit, a different locally optimal trajectory is also generated. Incorporating an initial search vector that targets the vicinity of the second closest approach (in backwards time) results in a “long” transfer scenario, such that a larger percentage of the transfer time occurs when the vehicle moves on the stable manifold, e.g., with the engines off. The target point is illustrated in Figure 4.11 along the same stable manifold tube that previously appeared in Figure 4.7.
Initial Guess for Insertion Point
Moon
Earth
Figure 4.11 Stable Manifold Tube for 12-Day L1 Halo (Green) and Initial Target Reference Trajectory Along the Manifold (Blue) for Long Transfer
The final solution for the transfer as plotted in Figure 4.12 reflects the increased time on the stable manifold. For this “long” transfer, the powered phase lasts 21.53 days and the coast phase 21.98 days for a total time-of-flight, TOF = 43.5 days. This represents a decrease of 4.43 over the powered arc when compared to the “short” transfer. But, an additional 11.13 days is added to the coast time. The total increase in transfer time compared to the “short” transfer is 6.7 days. Note also that the low-thrust transfer arc clearly departs from the original spiral structure once the spacecraft reaches the stable
84
manifold. The shape of the manifold trajectory also differs significantly from that along the “short” transfer.
Note also that a different manifold trajectory is incorporated
compared to that of the “short” transfer.
Earth
Moon
Figure 4.12 Low-Thrust Long Transfer to a 12-Day L1 Halo Orbit
85
The corresponding time histories of the costate variables corresponding to the “long” transfer appear in Figure 4.13 and Figure 4.14.
The position and velocity costates
confirm the oscillatory, nonlinear behavior observed previously. The final spacecraft mass is 1083.87 kg, yielding a total fuel consumption of 416.13 kg. Thus, even though the engine is “on” for a shorter duration, the long transfer consumes 55.83 kg more fuel than the short transfer. Some information that contributes to an understanding of this mass penalty is apparent in Figure 4.14. As the transfer approaches the periodic orbit, the engine supplies a thrust level that oscillates with increasing amplitude. In this case, the thrust magnitude reaches a maximum value of 3.66 N and a minimum level equal to 0.72 N. These thrust magnitudes obviously correspond to lower specific impulse ranges, and less efficient thrusting. These results are summarized in Table 4.3.
Figure 4.13 Position and Velocity Costate Time Histories for the 12-Day L1 Halo Orbit Transfer
86
Figure 4.14 Time History of Propulsion Related Parameters for the 12-Day L1 Halo Orbit Long Transfer
Table 4.3 12-Day L1 Halo Orbit Long Transfer Data Summary Parameter
Value
Units
Powered Time
21.982297
day
Coast Time
21.529817
day
Total Transfer Time
43.512114
day
Final Spacecraft Mass
1083.871837
kg
Propellant Mass Consumed
416.128162
kg
Tmax
3.657608
N
Tmin ∆v
0.715391
N
2.913369
km/s
87
4.2.2 Transfer to a 14-Day L1 Vertical Orbit
A 14-day L1 vertical orbit is selected from Figure 4.2. The associated manifold and the target transfer trajectory appear in Figure 4.15. The initially targeted insertion point is over 150,000 km from Earth, that is, approximately 50,000 km further than the insertion points on the stable manifold corresponding to the 12-day halo orbit.
Due to the
maximum z-amplitude of the orbit, Az = 57,000 km, the manifold trajectories are significantly out-of-plane at several points along the path. A global population of 500 uniformly distributed initial search vectors, S ′ , are created to initially propagate a more “dense” set during the “shotgun” process. The local optimization scheme required 200 iterations to determine a solution. The larger number of iterations is not unexpected since the thrust profile is required to “lift” the initial state along the planar parking significantly out-of-plane. The final transfer arc from the Earth parking orbit to vertical orbit insertion appears in Figure 4.16. The spacecraft achieves insertion onto the stable manifold after 18.26 days, and then coasts for 27.81 days. Thus, the total transfer time is 46.06 days. Time histories of the costates and the propulsion parameters are plotted in Figure 4.17 and Figure 4.18. Note that thrust amplitudes are bounded between 0.75 N and 2.3 N. This thrust range ultimately drives the spacecraft mass to a final value of 1128.4 kg; the total mass of fuel that is consumed is 371.64 kg. Despite the out-of-plane requirement and an insertion point that is located at a distance farther from the Earth, this transfer still requires less fuel than the long transfer to the 12-day L1 halo orbit.
88
Earth
Initial Guess for Insertion Point
Moon
Figure 4.15 Stable Manifold Tube for 14-Day L1 Vertical Orbit (Green) and Initial Target Reference Trajectory Along the Manifold (Blue)
89
Moon
Earth
Figure 4.16 Low-Thrust Transfer to a 14-Day L1 Vertical Orbit
90
Figure 4.17 Position and Velocity Costate Time Histories for the 14-Day L1 Vertical Orbit Transfer
Figure 4.18 Time History of Propulsion Related Parameters for the 14-Day L1 Vertical Orbit Transfer
91
Table 4.4 14-Day L1 Vertical Orbit Transfer Data Summary Parameter
Value
Units
Powered Time
18.257074
day
Coast Time
27.807328
day
Total Transfer Time
46.064402
day
Final Spacecraft Mass
1128.364331
kg
Propellant Mass Consumed
371.635669
kg
Tmax
2.307188
N
Tmin ∆v
0.748454
N
3.855672
km/s
4.2.3 Transfer to a 14-Day L2 Butterfly Orbit
The final orbit selected from Figure 4.2 for demonstration of the design of low-thrust transfers in this problem is a 14-day L2 butterfly orbit. The associated manifold and first guess for a manifold insertion trajectory appears in Figure 4.19. The initially targeted insertion state along the stable manifold, like the vertical orbit, is over 150,000 km (actually 162,000 km) from Earth.
This orbit is also unique because the manifold
trajectories exhibit unreliable extreme sensitivities to the angle-like parameter, θM. This behavior is observed in several L2 orbits as well as orbits with very low stability indices. When the stable manifolds associated with these orbits pass within the vicinity of the Earth, the erratic behavior is apparent in Figure 4.19. As a result, the variable θM is fixed (θM = 8.0 always), that is, the initial condition associated with the target stable manifold trajectory, X s ( ti ) , does not change. (A lower cost might be achieved by selecting new values θM and solving the problem again.) A global population of 500 uniformly distributed initial search vectors, S ′ , is established in the “shotgun” phase. From the global method, a transfer trajectory resembling the “long” transfer seen in Figure 4.12 is generated.
92
The final converged trajectory appears in Figure 4.20.
The spacecraft achieves
insertion into the stable manifold after 32.95 days, and coasts for 49.18 days. Clearly, this transfer is longer than all other examples, with a total transfer time of 82.13 days. Figure 4.21 and Figure 4.22 include the time history of the costates and propulsion parameters. The thrust amplitudes are bounded between the ranges of 0.59N and 1.93 N, the lowest values in the four examples here. These thrust-ranges ultimately drive the spacecraft to a final mass of 1166.84 kg, thus yielding a propellant cost of 333.16 kg.
Initial Guess for Insertion Point Earth Moon
Figure 4.19 Stable Manifold Tube for 14-Day L2 Butterfly Orbit (Green) and Initial Target Reference Trajectory Along the Manifold (Blue)
93
Earth Moon
Figure 4.20 Low-Thrust Transfer to a 14-Day L2 Butterfly Orbit
94
Figure 4.21 Position and Velocity Costate Time Histories for the 14-Day L2 Butterfly Orbit Transfer
Figure 4.22 Time History of Propulsion Related Parameters for the 14-Day L2 Butterfly Orbit Transfer
95
Table 4.5 14-Day L2 Butterfly Orbit Transfer Data Summary Parameter
Value
Units
Powered Time
32.950829
day
Coast Time
49.175507
day
Total Transfer Time
82.126336
day
Final Spacecraft Mass
1166.835826
kg
Propellant Mass Consumed
333.164174
kg
Tmax
1.928668
N
Tmin ∆v
0.596077
N
3.226374
km/s
96
5. SUMMARY AND RECOMMENDATIONS
5.1 Summary
A strategy for the determination of finite-thrust transfer trajectories that deliver a spacecraft into collinear libration point orbits is developed. The computational scheme involves a systematic sequence of steps to generate these transfers. The resulting control law requires that the engine always operates at maximum power. These transfers can be decomposed into two phases – a fully powered phase followed by the arrival at a trajectory along the stable manifold tube, i.e., an un-powered phase yielding an asymptotic arrival to the vicinity of the target orbit.
Clearly, transfers to collinear
libration point orbits are unique due to the availability of the stable manifold that allows for natural transfer arcs. The manifolds do not always pass conveniently close to the Earth, however. Low-thrust options serve as a bridge between an Earth parking orbit and the stable manifold. The initial powered phase, however, requires the determination of an unknown thrust magnitude and direction history, i.e., the primary controls in the “steering” problem. Variable specific impulse engines, an advanced prototype low-thrust engine, allow the thrust profile to vary via throttling of the specific impulse. To resolve these unknown controls, optimal control theory is used to derive the two-point boundaryvalue problem with free final time that corresponds to a necessary and sufficient condition for a maximal value of final spacecraft mass. Thus, optimal control is useful not only for determining a control law, but also for obtaining locally minimal transfer costs. Solving the hybrid direct/indirect local method via direct shooting requires the determination of extremely sensitive initial costates, with magnitudes that supply little or no physical meaning. The adjoint control transformation yields a mapping process to
97
resolve these costates into meaningful quantities such as initial angles, angular rates, and thrust magnitude. Even though this mapping does not completely remove the numerical sensitivities, it is used to determine bounds on the initial conditions.
Randomized
collections of initial conditions within these bounded regions are then propagated in a global “shotgun” approach. The terminal errors in the elements of the state vector for each trajectory, at each instant in time, are then computed for any potential insertion point along the stable manifold tube, resulting in a “best” initial guess, i.e., the initial search vector that provides information such as the “starter” state vector. The initial search vector is then used to initiate a local hybrid direct/indirect optimization problem.
Rather than solve the complete optimal control, two-point
boundary-value problem (an indirect problem), direct shooting is utilized to directly reduce the performance index. Additional design variables are incorporated into the search vector, such as the initial departure angle along the parking orbit, and the anglelike and time-like variables that parameterize the current state along the manifold tube. Components of the indirect problem, such as parameterizations for the controls, characterize the nature of the hybrid direct/indirect shooting algorithm. Direct shooting is accomplished via a medium-scale SQP algorithm, where the only active constraints are the kinematic boundary conditions. To demonstrate this methodology for the determination of low-thrust transfers, the problem of lunar south pole coverage in the Earth-moon problem is introduced. Families of halo, vertical, and “butterfly” orbits are generated using simple targeting procedures. The orbits generated in these families meet the specified criteria and comprise potential orbits in a relay satellite configuration between a facility at the lunar south pole and Earth-based ground stations. Ultimately, information such as period (used in satellite phasing), and the stability index (used in estimating potential stationkeeping costs) determine the orbits that are selected for further analysis. A sample orbit is selected from each of the different families to validate the low-thrust transfer concept.
98
5.2 Recommendations for Future Work
These transfer trajectories are preliminary in nature. Design within the context of the circular restricted three-body problem is only an initial step in the design of actual trajectory arcs in the full ephemeris model. Further development in this problem includes an extension of the control law and solutions using a higher fidelity gravity model. Numerical determination of a single transfer trajectory is nontrivial, thus, an exhaustive trade study is not currently practical. However, this preliminary analysis is a crucial step towards balancing the lowest possible transfer cost with the lowest possible orbital maintenance costs for long duration coverage missions.
Other, more sophisticated
transfer schemes incorporating multiple manifolds, and potential lunar flybys are not included at this time. Examination of a wider range of orbits is also likely to yield alternatives not yet apparent. Future investigations of transfers to these orbits, including highly elliptical orbits, would also benefit from the mission analysis results of lunar south pole coverage. As mentioned previously, computational issues significantly impact the capability to produce results. Integration time is significant, and hundreds of iterations are necessary, with typical run-times of 2 or more hours. Developing an algorithm in a programming language optimized for speed will allow more efficient numerical processes and expand analysis options. Finally, a full range of local methods, or more highly sophisticated global methods, to identify potential solutions were not attempted. Other approaches, mentioned in Chapter 1, that use control parameterizations for multiple shooting, collocation, and dynamic programming may be useful.
While some of these methods may require the
discretization of the continuous trajectory into many impulsive burns, it may still result in more rapid convergence to a solution. Clearly, the low-thrust problem opens many avenues for future research efforts in numerous problems, including the lunar south pole coverage problem.
99
LIST OF REFERENCES
[1] R. Farquhar, “The Flight of ISEE-3/ICE: Origins, Mission History, and a Legacy.” Journal of the Astronautical Sciences, Vol. 49, No. 1, January-March 2001, pp. 23-73. [2] H. Franz et al., "WIND Nominal Mission Performance and Extended Mission Design." Journal of the Astronautical Sciences, Vol. 49, No. 1, pp. 145-167. [3] D. Dunham et al., “Transfer Trajectory Design for the SOHO Libration-Point Mission.” IAF-92-0066, 43rd Congress of the International Astronautical Federation, Washington, D.C., USA, August/September, 1992. [4] P. Sharer and T. Harrington, “Trajectory Optimization for the ACE Halo Orbit Mission, Paper AAS 96-3601.” AAS/AIAA Astrodynamics Specialist Conference, San Diego, California, July, 1996. [5] M. Mesarch and S. Andrews, “The Maneuver Planning Process for the Microwave Anisotropy Probe (MAP) Mission.” AIAA/AAS Astrodynamics Specialists Conference and Exhibit, Monterey, California, August 5-8, 2002. [6] M. Lo, B. Williams, W. Bollman, D. Han, Y. Hahn, J. Bell, E. Hirst, R. Corwin, P. Hong, K. Howell, B. Barden, and R. Wilson, “Genesis Mission Design.” Paper No. AIAA 1998-4468, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Boston, Massachusetts, August 10-12, 1998. [7] H. Poincaré, “Les Méthodes Nouvelles de la Mécanique Celeste.” Vol II, Paris, Gauthier-Villars, 1893. [8] A. Whitehead and B. Russell. Principia Mathematica, 3 Vols., University Press, Cambridge, 1910, 1912, and 1913. [9] W. Wiesel. Spaceflight Dynamics. Mcgraw-Hill Book Company, New York, 1962. [10] V. Szebehely, Theory of Orbits. Academic Press, New York, 1967.
100
[11] G. Hill, “Mr. Hill’s Researches in the Lunar Theory.” Monthly Notices of the Royal Astronautical Society, Vol. 39, November 1879, pp. 258-261. [12] G. Darwin, “Periodic Orbits.” Acta Mathematica, Vol. 21, 1897, pp. 99-242. [13] H. Plummer, “On Oscillating Satellites.” Monthly Notices of the Royal Astronautical Society, Vol. 63, June 1903, pp. 436-443. [14] F. Moulton (in collaboration with D. Buchanan, T. Buck, F. Griffin, W. Longley, and W. MacMillan), “Periodic Orbits.” Carnegie Institution of Washington, Washington, 1920. [15] R. Farquhar and A. Kamel, “Quasi-Periodic Orbits about the Translunar Libration Point.” Celestial Mechanics, Vol. 7, 1973, pp. 458-473. [16] D. Richardson and N. Cary, “A Uniformly Valid Soluton for Motion about the Interior Libration Point of the Perturbed Elliptic-Restricted Problem.” Paper No. AAS 75-021, AAIA/AAS Astrodynamics Conference, Nassau, Bahamas, July 28-30, 1975. [17] J. Breakwell and J. Brown, “The ‘Halo’ Family of 3-Dimensional Periodic Orbits in the Earth-Moon Restricted 3-Body Problem.” Celestial Mechanics, Vol. 20, 1979, pp. 389-404. [18] K. Howell, “Three-Dimensional, Periodic, ‘Halo’ Orbits.” Celestial Mechanics, Vol. 32, No. 1, Janurary 1984, pp. 29-52. [19] M. Hénon, “Vertical Stability of Periodic Orbits in the Restricted Problem.” Celestial Mechanics, Vol. 8, 1973, pp. 269-272. [20] K. Howell and J. Breakwell, “Almost Rectilinear Halo Orbits.” Celestial Mechanics, Vol. 32, No. 1, January 1984, p. 29-52. [21] D. Dichmann, E. Doedel, and R. Paffenroth, “The Computation of Periodic Solutions of the 3-Body Problem Using the Numerical Continuation Software AUTO.” Libration Point Orbits and Applications, G. Gómez, M. Lo, J. Masdemont (Editors), Hong Kong, World Scientific, 2003. [22] H. Goldstine, A History of the Calculus of Variations from the 17th Through the 19th Century, Springer Publishing, Berlin, 1980.
101
[23] C. Padra, “The Beginnings of Variational Calculus, and its Early Relation with Numerical Methods.” Variational Formulations in Mechanics: Theory and Applications, E. Taroco, E. de Souxa Neto, and A. Novotny (Editors), The International Center for Numerical Methods in Engineering (CIMNE), Barcelona, Spain, 2006. [24] A. Bryson and Y. Ho, Applied Optimal Control. Blaisdell Publishing Company, Waltham, Massachusetts, 1969. [25] D. Hull, Optimal Control Theory for Applications. Springer Publishing, New York, 2003. [26] D. Kirk, Optimal Control Theory: An Introduction. Prentice-Hall Publishing Company, Englewood Cliffs, New Jersey, 1970. [27] D. Lawden, Optimal Trajectories For Space Navigation. Butterworths Publishers, London, 1963. [28] J. Senent, C. Ocampo, and A. Capella, “Transfers to and Station-keeping for Unstable Periodic Orbits with Low Thrust Propulsion.” Paper No. AAS 03-633, AAS/AIAA Astrodynamics Specialists Conference, Big Sky, Montana, August 3-7, 2003. [29] C. Kluever and B. Pierson, “Three-State Approach to Optimal Low-Thrust Earth-Moon Trajectories,” Journal of Guidance, Control, and Dynamics, Vol.17, No. 6, 1994, pp. 1275-1282. [30] C. Kluever and B. Pierson, “Optimal Low-Thrust Three-Dimensional EarthMoon Trajectories,” Journal of Guidance, Control, and Dynamics, Vol.18, No. 4, 1995, pp. 830-837. [31] J. Désidéri, S. Peigin, and S. Timchenko, Application of Genetic Algorithms to the Solution of the Space Vehicle Reentry Trajectory Optimization Problem. Institut National de Recherche en Informatique et en Automatique (INRIA) Sophia Antipolis, Research Report No. 3843, 1999. [32] J. Marec, Optimal Space Trajectories. Elsevier Scientific Publishing Company, Amsterdam, 1979. [33] W. Melbourne and C. Sauer, “Optimum interplanetary rendezvous trajectories with power limited vehicles.” AIAA Journal, Vol. 1, January 1963, pp. 54-60. [34] D. Jezewski, “Primer Vector Theory and Applications.” Technical Report NASA TR R-454, Lyndon B. Johnson Space Center, Houston, Texas, November 1975.
102
[35] J. Betts, “Survey of Numerical Methods for Trajectory Optimization.” Journal of Guidance, Control, and Dynamics, Vol. 21, No. 2, March-April 1998, pp. 193-207. [36] C. Ranieri and C. Ocampo, “Optimization of Round-Trip, Time-Constrained, Finite-Burn Trajectories via an Indirect Method.” Journal of Guidance, Control, and Dynamics, Vol. 28, No. 2, 2005, pp. 306-314. [37] S. Zimmer and C. Ocampo, “Analytical Gradients for Gravity Assist Trajectories Using Constant Specific Impulse Engines,” Journal of Guidance, Control, and Dynamics, Vol. 28, No. 4, 2005, pp. 753-760. [38] K. Chen, T. McConaghy, M. Okutsu, and J. Longuski, “A Low-Thrust Version of the Aldrin Cycler.” Paper No. AIAA 2002-4421, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Monterey, California, August 2-8, 2002. [39] T. Debban, T. McConaghy, and J. Longuski, “Design and Optimization of LowThrust Gravity-Assist Trajectories to Selected Planets.” Paper No. AIAA 20024729, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Monterey, California, August 2-8, 2002. [40] V. Coverstone-Caroll and S. Williams, “Optimal Low Thrust Trajectories Using Differential Inclusion Concepts.” Journal of the Astronautical Sciences, Vol. 42, No. 4, 1994, pp. 379–393. [41] A. Herman and B. Conway, “Optimal, Low-Thrust, Earth-Moon, Orbit Transfers.” Journal of the Guidance, Control, and Dynamics, Vol. 21, No. 1, Janurary-February, 1998, pp. 141-147. [42] O. Golan and J. Breakwell, “Minimum Fuel Lunar Trajectories for a Low-Thrust Power-Limited Spacecraft,” Paper No. AIAA 90-2975, AIAA/AAS Astrodynamics Conference, Portland, Oregon, August 20-22, 1990. [43] A. Sukhanov and N. Eismont, “Low Thrust Transfer To Sun-Earth L1 and L2 Points With A Constraint On The Thrust Direction.” Libration Point Orbits and Applications Conference, Parador d'Aiguablava, Girona, Spain. June 10-14, 2002. [44] H. Seywald, C. Roithmayr, and P. Troutman,, “Fuel-Optimal Orbital Transfers For Variable Specific Impulse Power Spacecraft.” Paper No. AAS 03-122, AIAA/AAS Space Flight Mechanics Meeting, Puerto Rico, February 9-13, 2003.
103
[45] R. Russell, “Primer Vector Theory Applied to Global Low-Thrust Trade Studies.” Paper No. AAS 06-156, AIAA/AAS Space Flight Mechanics Meeting, Tampa, Florida, January 22-26, 2006. [46] D. Grebow, M. Ozimek, K. Howell, and D. Folta, “Multi-Body Orbit Architectures for Lunar South Pole Coverage.” Paper No. AAS 06-179, AIAA/AAS Space Flight Mechanics Meeting, Tampa, Florida, January 22-26, 2006. [47] T. Parker and L. Chua, Practical Numerical Algorithms for Chaotic Systems. Springer-Verlag Publishing, New York, 1989. [48] J. Guckenheimer and P. Holmes, Nonlinear Oscillations. John Wiley and Sons, New York, 1995. [49] L. Perko, Differential Equations and Dynamics Systems. Springer-Verlag Publishing, New York, 1991. [50] B. Marchand, “Temporary Satellite Capture of Short-Period Jupiter Family Comets From the Perspective of Dynamical Systems.” M.S.A.A. Thesis, School of Aeronautics and Astronautics, Purdue University, 2000. [51] V. Yakubovich and V. Starzhinskii, Linear Differential Equations with Periodic Coefficients. John Wiley and Sons, New York, 1975. [52] L. Pontryagin, “The mathematical theory of optimal processes.” Interscience Publishers, New York, 1962. [53] E. McShane, “On Multipliers for Lagrange Problems.” American Journal of Mathematics, Vol. 61, New York, 1939, pp. 809-819. [54] L. Hiday, “Optimal Transfers Between Libration-Point Orbits in the Elliptic Restricted Three-Body Problem.” Ph.D. Thesis, School of Aeronautics and Astronautics, Purdue University, 1992. [55] Smart-1 Engine Image, http://esamultimedia.esa.int/images/smart_1/. European Space Agency website, Copyright 2000-2006. [56] VaSIMR Rocket Image, http://www.nasaexplores.com/. NASA Explores website, 2006. [57] T. Sakai. “A Study of Variable Thrust, Variable Isp Trajectories for Solar Exploration.” PhD Thesis, School of Aerospace Engineering, Georgia Institute of Technology, 2004.
104
[58] L. Dixon and M. Biggs, “The Advantages of Adjoint-Control Transformations when Determining Optimal Trajectories by Pontryagin’s Maximum Principle,” Aeronautical Journal, Vol. 76, No. 735, 1972, pp. 169-174. [59] G. Vanderplaats, Numerical Optimization Techniques for Engineering Design. Vanderplaats Research & Development, Inc., Colorado Springs, 2001, pp. 250256. [60] S. Nozette, C. Lichtenberg, P. Spudis, R. Bonner, W. Ort, E. Maleret, M. Robinson, and E. Shoemaker, “The Clementine Bistatic Radar Experiment.” Science, Vol. 274, Issue 5292, 1996, pp. 1495-1498. [61] E. Kozlova, “The Presence of Volatiles in the Polar Regions of the Moon.” 35th Committee on Space Research (COSPAR) Scientific Assembly, Paris, France, July 18-25, 2004. [62] “The Vision for Space Exploration,” National Aeronautics and Space Administration Publication, NP-2004-01-334-HQ, February 2004. [63] D. Grebow, “Generating Periodic Orbits in the Circular Restricted Three-Body Problem with Applications to Lunar South Pole Coverage.” MSAA Thesis, School of Aeronautics and Astronautics, Purdue University, 2006. [64] R. Farquhar, “The Utilization of Halo Orbits in Advanced Lunar Operations.” NASA TND-365, GSFC, Greenbelt, Maryland, 1971. [65] G. Gómez, J. Libre, R. Martínez, and C. Simó, Dynamics and Mission Design Near Libration Points; Vol. 1, Fundamentals, The Case of Collinear Libration Points, World Scientific Publishing Company, Singapore, 2001. [66] I. Robin and V. Markellos, “Numerical Determinations of Three-Dimensional Orbits Generated From Vertical Self-Resonant Satellite Orbits.” Celestial Mechanics, Vol. 21, 1980, pp. 395-434.
LIST OF REFERENCES
