Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Autonomous interplanetary constellation design
(USC Thesis Other)
Autonomous interplanetary constellation design
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
AUTONOMOUS INTERPLANETARY CONSTELLATION DESIGN by Cornelius Channing Chow II A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ASTRONAUTICAL ENGINEERING) May 2012 Copyright 2012 Cornelius Channing Chow II Dedication To my wonderful wife Irina ii Acknowledgements First and foremost, I would like to express my sincerest gratitude to Dr. Benjamin F. Villac for giving me the rare opportunity to join his research team and participate in the most valuable learning experience of my life. His guidance, enthusiasm, and immense knowledge have nourished me throughout my research career. I could not have imagined having a better advisor and mentor for my Ph.D. studies. Special thanks are due to Dr. Martin W. Lo for his unyielding support and dedication, without which, my involvement with this work would not have been possible. It is truly an honor to work with both of these esteemed scholars. I am also grateful to Dr. Gerald R. Hintz and James K. Miller for sharing their wisdom and oering their advice on writing technical papers and giving presentations. And, of course, I am thankful to the faculty and sta of USC's Department of Astronautical Engineering who have supported me all of these years. Finally, I wish to thank my friends for their patience and humor throughout my time in graduate school. And most endearing of all, I thank my family for their unconditional love, encouragement, and belief in me since day one. iii Table of Contents Dedication ii Acknowledgements iii List Of Tables vii List Of Figures viii Abstract xi Chapter 1 : Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Previous Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Problem Denition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Dissertation Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Chapter 2 : Background 9 2.1 Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 Circular Restricted 3-Body Problem . . . . . . . . . . . . . . . . . 10 2.1.2 Ephemeris Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Dynamical Systems Concepts . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.1 Invariant Manifold Structures . . . . . . . . . . . . . . . . . . . . . 16 2.2.2 Periodic Orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.3 Bifurcation Diagrams . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3 Continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3.1 Numerical Continuation . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3.2 The Implicit Function Theorem . . . . . . . . . . . . . . . . . . . . 23 2.3.3 Keller's Pseudo-Arclength Technique . . . . . . . . . . . . . . . . . 24 2.4 Spacecraft Navigation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.4.1 Statistical Orbit Determination . . . . . . . . . . . . . . . . . . . . 27 2.4.2 Observability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.4.3 The LiAISON Strategy . . . . . . . . . . . . . . . . . . . . . . . . 31 iv Chapter 3 : Constellation Design Methodology 34 3.1 General Optimization Problem . . . . . . . . . . . . . . . . . . . . . . . . 35 3.1.1 Parametrization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.1.2 Formulating the Minimization Problem . . . . . . . . . . . . . . . 39 3.1.3 Existence of Local Extrema . . . . . . . . . . . . . . . . . . . . . . 40 3.1.4 Scope of the Solution Space . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Solution Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.2.1 Heuristic Approach and Local Optimizers . . . . . . . . . . . . . . 42 3.2.2 Mapping the Design Space using Continuation . . . . . . . . . . . 44 Chapter 4 : Tailoring the Constellation Design Framework 49 4.1 Modeling the Dynamical System . . . . . . . . . . . . . . . . . . . . . . . 49 4.1.1 Restricted Earth-Moon System . . . . . . . . . . . . . . . . . . . . 50 4.1.2 Subsets of Periodic Orbit Families . . . . . . . . . . . . . . . . . . 52 4.2 Choosing the Cost Function . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2.1 Geometric Dilution of Precision . . . . . . . . . . . . . . . . . . . . 57 4.2.2 Estimation Error Covariance . . . . . . . . . . . . . . . . . . . . . 61 Chapter 5 : Results for the Earth-Moon System 65 5.1 Multi-Spacecraft Arrangements . . . . . . . . . . . . . . . . . . . . . . . . 65 5.1.1 L1 Halo Case for 2,3, and 4 Satellites . . . . . . . . . . . . . . . . 67 5.2 Design Space Maps and Optimal Congurations . . . . . . . . . . . . . . 70 5.2.1 L1 Halo Family . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2.2 L1 Axial Family . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.2.3 L2 Vertical Family . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.2.4 L5 Vertical Family . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.2.5 Ephemeris Simulations . . . . . . . . . . . . . . . . . . . . . . . . . 78 Chapter 6 : Performance Characterization 81 6.1 Computational Eciency . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.2 Continuation of Stochastic Systems . . . . . . . . . . . . . . . . . . . . . . 83 6.3 Isolated Islands of Local Extrema . . . . . . . . . . . . . . . . . . . . . . . 84 Chapter 7 : Conclusion 85 7.1 Summary of Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.2 Insights Gained . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 7.3 Recommendations for Future Study . . . . . . . . . . . . . . . . . . . . . 89 References 90 Appendix A Coordinate Transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Appendix B Performance Index Partials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 v Appendix C Initial Conditions for Periodic Orbits . . . . . . . . . . . . . . . . . . . . . . . . 106 Appendix D GMAT Conguration Settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 vi List Of Tables 3.1 Number of Computations . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1 Sun-Earth-Moon Neighborhood Physical Constants . . . . . . . . . . . . . 51 4.2 Sun-Earth-Moon Neighborhood CR3BP Constants . . . . . . . . . . . . . 51 4.3 Distances of Collinear Libration Points from Nearest Primary . . . . . . . 51 5.1 Initial Conditions for the Representative L1 Halo Orbit . . . . . . . . . . 67 5.2 Initial Conditions for Optimal Constellations (Synodic Frame) . . . . . . 78 5.3 Initial Conditions for Optimal Orbits (Sidereal Frame) . . . . . . . . . . . 79 C.1 L1 Halo Family . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 C.2 L1 Axial Family . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 C.3 L2 Vertical Family . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 C.4 L5 Vertical Family . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 vii List Of Figures 2.1 Schematic representation of the Earth-Moon CR3BP geometry. Five La- grange points are shown, along with the reference unit circle and equilateral triangles. The third-body restricted mass (i.e. spacecraft) is shown in the upper right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Collinear libration points for = 0:2. Distances from each libration point to the nearest primary ( i ) are colored light green. Note that the distance between the primaries is normalized to 1 dimensionless length unit. . . . . 13 2.3 The Earth-Moon L2 vertical Lyapunov family projected in the synodic position space as (a) a set of discrete orbits and (b) a smooth continuous family; both colored by the continuum index. . . . . . . . . . . . . . . . . 18 2.4 Bifurcation diagram of several sets of L1 and L2 periodic orbit families associated with the Earth-Moon CR3BP. The color dimension shows the orbital period per branch, with red/purple indicating high/low values, re- spectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5 Example of the implicit function theorem nding a local solution manifold M M in the neighborhood of (a 0 ;b 0 ). The unit circle,M, and subsets,A andB, containinga 0 andb 0 respectively, are also shown. . . . 24 2.6 Graphical representation of pseudo-arclength continuation. Depicted are the initial solution, z 0 , the associated path tangent, _ z 0 , and several pre- dictions during Newton's iteration, z () 1 , as bordered white dots. The converged solution,z 1 , is shown as a bordered white square. . . . . . . . . 26 3.1 Phase angle parametrization, i , along an L1 halo orbit projection. . . . . 38 4.1 A selected subset of the L1 axial family, plotted in the synodic frame. . . 53 4.2 A selected subset of the L1 halo family, plotted in the synodic frame. . . . 53 4.3 A selected subset of the L2 vertical family, plotted in the synodic frame. . 54 viii 4.4 A selected subset of the L5 vertical family, plotted in the synodic frame. . 54 4.5 The admissible Earth-Moon periodic orbit subspace, plotted in the synodic frame. The Earth and Moon are shown as the large and small black spheres (exaggerated in scale for visualization), respectively. . . . . . . . . . . . . 55 4.6 Example congurations for (a) \good" and (b) \bad" 4-spacecraft geome- tries (red dots), with respect to an end-user (central black dot), as seen mapped onto the unit sphere. . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.7 Volumetric plots of GDOP values in a region about a sample Earth-Moon L1 halo orbit, where the white dots show the locations of eight satellites at a given instance in time. Shown in the color dimension are the accuracies seen by an end-user using a constellation that has spacecraft (a) evenly spaced in phase angle, and (b) unevenly clustered towards one side. . . . . 60 5.1 Cost function values, , shown on a log scale to accentuate details, for a sample halo constellation of (a) two, (b) three, and (c-f) four spacecraft. . 69 5.2 (a) Continuation-generated bifurcation diagram of the design space for the L1 halo family, overlaid on top of a background surface map compiled from an exhaustive grid-search. The \optimal" orbit is highlighted by the star marker, whose projection in three-dimensions is plotted separately (b) in context along the manifold and (c) with orthogonal projections and spacecraft at one instance in time. . . . . . . . . . . . . . . . . . . . . . . 74 5.3 (a) Continuation-generated bifurcation diagram of the design space for the L1 axial family, overlaid on top of a background surface map compiled from an exhaustive grid-search. The \optimal" orbit is highlighted by the star marker, whose projection in three-dimensions is plotted separately (b) in context along the manifold and (c) with orthogonal projections and spacecraft at one instance in time. . . . . . . . . . . . . . . . . . . . . . . 75 5.4 (a) Continuation-generated bifurcation diagram of the design space for the L2 vertical family, overlaid on top of a background surface map compiled from an exhaustive grid-search. The \optimal" orbit is highlighted by the star marker, whose projection in three-dimensions is plotted separately (b) in context along the manifold and (c) with orthogonal projections and spacecraft at one instance in time. . . . . . . . . . . . . . . . . . . . . . . 76 ix 5.5 (a) Continuation-generated bifurcation diagram of the design space for the L5 vertical family, overlaid on top of a background surface map compiled from an exhaustive grid-search. The \optimal" orbit is highlighted by the star marker, whose projection in three-dimensions is plotted separately (b) in context along the manifold and (c) with orthogonal projections and spacecraft at one instance in time. . . . . . . . . . . . . . . . . . . . . . . 77 5.6 Optimal orbits corresponding to (a) L1 halo, (b) L1 axial, (c) L2 vertical, and (d) L5 vertical cases, with red/black lines indicating CR3BP/ephemeris solutions, respectively, shown in the non-dimensional synodic frame. . . . 80 A.1 A sample Earth-Moon L2 vertical periodic orbit plotted in black for both the (a) synodic and (b) sidereal frames. Orthogonal projections are shown in light gray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 x Abstract According to NASA's integrated space technology roadmaps, space-based infrastruc- tures are envisioned as necessary ingredients to a sustained eort in continuing space exploration. Whether it be for extra-terrestrial habitats, roving/cargo vehicles, or space tourism, autonomous space networks will provide a vital communications lifeline for both future robotic and human missions alike. Projecting that the Moon will be a bustling hub of activity within a few decades, a near-term opportunity for in-situ infrastructure development is within reach. This dissertation addresses the anticipated need for in-space infrastructure by inves- tigating a general design methodology for autonomous interplanetary constellations; to illustrate the theory, this manuscript presents results from an application to the Earth- Moon neighborhood. The constellation design methodology is formulated as an optimization problem, in- volving a trajectory design step followed by a spacecraft placement sequence. Modeling the dynamics as a restricted 3-body problem, the investigated design space consists of fam- ilies of periodic orbits which play host to the constellations, punctuated by arrangements of spacecraft autonomously guided by a navigation strategy called LiAISON (Linked xi Autonomous Interplanetary Satellite Orbit Navigation). Instead of more traditional ex- haustive search methods, a numerical continuation approach is implemented to map the admissible conguration space. In particular, Keller's pseudo-arclength technique is used to follow folding/bifurcating solution manifolds, which are otherwise inaccessible with other parameter continuation schemes. A succinct characterization of the underlying structure of the local, as well as global, extrema is thus achievable with little a priori intuition of the solution space. Furthermore, the proposed design methodology oers benets in computation speed plus the ability to handle mildly stochastic systems. An application of the constellation design methodology to the restricted Earth-Moon system, reveals optimal pairwise congurations for various L1, L2, and L5 (halo, ax- ial, and vertical) periodic orbit families. Navigation accuracies, ranging fromO(10 1 ) meters in position space, are obtained for the optimal Earth-Moon constellations, given measurement noise on the order of 1 meter. xii Chapter 1 Introduction In early October 1957, the former Soviet Union grabbed the world's attention with the successful launch of Sputnik 1: the rst articial satellite put into orbit about the Earth. Since then, the space age has ushered in thousands of spacecraft into Earth's lo- cal neighborhood; counting communications, navigation, weather, and scientic research satellites among the majority. As an undeniably important part of our daily lives, satel- lite constellations continue to expand and advance, oering greater service coverage and enhanced capabilities. This chapter heralds the motivation for this research with a spe- cic outgrowth in constellation design, geared towards the goal of lunar colonization. Since interplanetary constellations do not yet exist, some past studies on the pre-cursory eorts of satellite constellation designs in the 2-body approximation are highlighted. A generalized problem statement is then posed, and an overview of the remainder of the dissertation is outlined. 1 1.1 Motivation According to NASA's integrated space technology roadmaps, 30 space-based infras- tructures are envisioned as necessary ingredients to a sustained eort in continuing space exploration. Whether it be for extra-terrestrial habitats, roving/cargo vehicles, or space tourism, autonomous space networks will provide a vital communications lifeline for both future robotic and human missions alike. The Moon is a critical destination for the development of space. While it certainly is a stepping stone on the way to the rest of the solar system, it is itself a world of wonder and resources as yet almost entirely unknown and unexplored. One can envision that, within a few decades, it will be a busy place of human activities, involving both private and public ventures. As such, infrastructures are needed to aid in the transport of both humans and cargo between the Earth and the Moon. While ground support can provide the necessary tracking/monitoring infrastructure for the rst return missions to the Moon, this resource may not be a viable option for a long-term enterprise that would involve repeated access to the lunar surface. 31 Experience with conducting operations in space (e.g. International Space Station ), indeed reveals an important reliance on having adequate navigation and communication capabilities. For near Earth missions, these capabilities are provided by a combination of ground-based and space-based infrastructures, as exemplied with the Global Positioning System 20 (GPS) and the Tracking and Data Relay Satellite System y (TDRSS). Together, these in-space assets provide a robust network capable of handling a wide array of missions http://nasa.gov/station y http://tdrs.gsfc.nasa.gov/ 2 in the immediate vicinity of Earth. As no such support structure exists for the larger Earth-Moon neighborhood, further innovation is required to extend this capability beyond low-Earth orbit. One idea that was proposed to address this deciency in cis-lunar coverage blended the notions of periodicity together with spacecraft autonomy in a concept called Inter- planetary Navigation and Communications Network 40 (INCnet), that featured a novel autonomous spacecraft navigation strategy known as LiAISON 17,18 . Functionally, this concept is equivalent to the GPS and TDRSS networks combined, but designed instead for use in the multi-body gravitational dynamics of the Earth-Moon neighborhood, at a reduced overall systems cost. This research originally aimed to contribute to the INCnet concept, by developing a design process to systematically evaluate possible multi-spacecraft congurations for the Earth-Moon system. However, in tackling this specic problem, a more general constellation design framework was born; both of which are captured in this dissertation. But before delving into the details of the proposed design methodology, a brief literature review is given to place this work into a more precise context. 1.2 Previous Research Formally studied by the likes of Kepler (c. 1609), Newton (c. 1687), Euler (c. 1767), Lagrange (c. 1772), and Poincar e (c. 1887), the motion of celestial bodies was, and still is, a scholarly pursuit that has garnered widespread attention, especially over the past Linked Autonomous Interplanetary Satellite Orbit Navigation. 3 several centuries. Originally built on the astronomical theories of Kepler in his study of planetary motion, classical mechanics was later revised by Newton to incorporate his universal law of gravitation and the three fundamental laws of generalized motion. While the problem of two bodies was completely solved by Bernoulli, Newton was among the rst to consider the motion of three bodies under the in uence of their mutual gravi- tational attraction. Euler made important strides in nding particular solutions to the restricted problem and discovered the existence of three collinear equilibrium points, all by re-formulating the problem in a rotating reference frame. Lagrange subsequently found two more equilateral triangular libration points (though all ve typically bear his name- sake), as well as simplied the description of the mechanics by using Hamilton's principle of stationary action. Poincar e then radically transformed conventional thinking by in- troducing the idea of observing the qualitative behavior of solutions rather than always seeking a quantitative analysis. These selected historical contributions are, however, only mere snippets among a vast profusion of discoveries/advancements that paved the way for modern day space ight. Today, solving for the motion of spacecraft orbiting a celestial body, classically in- volves approximating the problem as a restricted 2-body system; that is, with one central gravitational body and a negligible point mass within the primary's sphere of in uence. 28 Devoid of any other perturbations, the only admissible conic section that satises an \orbit " is an ellipse (of which the circle is a special case y ). This traditional modeling reduction works remarkably well in a close enough proximity to the central body and An orbit, from the study of dynamical systems, is the image of a function; i.e. a subset in phase space. y Making the distinction that an orbit is a closed path, eliminates counting parabolic/hyperbolic arcs. 4 on suciently small timescales. Even for discussing rudimentary interplanetary trajecto- ries, the patched-conics approach resorts to using decoupled 2-body systems. 4 However, as the need arises for higher delity models, tertiary perturbations (e.g. spherical har- monics, atmospheric drag, and solar radiation pressure) can be successively added to the equations of motion to bring the model closer to representing the actual dynamics. Often, a dierential corrections process is added to adjust a prototypical conic section to achieve continuity in a perturbed environment. 26 With regards to the design of an entire constellation or a cluster of formation- ying spacecraft , much of the available research literature focuses on Earth orbiting congu- rations, where indeed using the previous approach can safely ignore third-body eects. In some recent work, for example, Mortari et al. 29 promotes the concept of a \Flower Constellation Set" that is based on phasing a set of design parameters (of which a few are classical orbital elements) as a systematic way to orient various orbital planes and place spacecraft among them. On the other hand, Schaub et al. 35 developed a dynamics based J2-invariant model that exploits the Earth's oblateness eects on orbital drift rates in order to properly congure the chief and deputy spacecraft of a cluster formation. Before fully venturing into the 3-body problem, an interesting scheme suggested by Chow et al., 9 bridges the two models by launching a eet of microsatellites to the Earth-Moon L1 libration point and utilizing the instabilities in the dynamics to provide dierent dis- persion routes to target 2-body orbits about the Earth. Although representing vastly dierent strategies, ultimately, all of these techniques at their core still rely on a very Essentially a satellite constellation where the relative distances between the spacecraft are usually much smaller as compared to the orbital radii, and are kept xed over time. 25 5 specic dynamical system for the nal constellation design: the perturbed, restricted, 2-body system. Among the rst to notably mention interplanetary orbits for mission design (i.e. using 3-body dynamics), Farquhar advocates the utility of the Earth-Moon L2 halo orbits for telecommunications relays for an Apollo mission. 14 Howell 21 has since then sparked in- tense interest in this subject, beginning with her demonstration of nding precise periodic orbits using numerical routines. More recently, Doedel et al. 13 conducted an in-depth numerical continuation survey on several families of periodic orbits bifurcating from each equilibrium point in the Earth-Moon system. Casoliva et al. 5 categorized dierent classes of periodic families called resonant cyclers, that create repeated access gateways between the Earth and Moon using high-energy orbits and heteroclinic connections. As evidenced by this sample of recent advances, great lengths have been taken to characterize trajecto- ries and classify orbit spaces for use in Earth-Moon mission concepts. There is, however, little emphasis placed on the follow-on problem of arranging multi-spacecraft congu- rations along these interplanetary manifolds to achieve a particular constellation design objective. One such example, however, would be the aforementioned INCnet concept that ad- dresses the spacecraft placement problem in conjunction with interplanetary orbit design. In particular, leveraging the LiAISON technique (further detail given in Chapter 2) to guide spacecraft arrangements along non-Keplerian orbits in an asymmetrical dynamical system, represents a link between traditional 2-body constellation design techniques and trajectories in the 3-body problem, thus providing the impetus for this thesis work. 6 1.3 Problem Denition Constellation design indeed encompasses many challenges, ranging from conguring the orbital geometry to devising the actual functionality of the satellites as a collective. While some parameters are discrete (e.g. the number of spacecraft), others fall into a more complicated spectrum of values (e.g. surface area coverage or navigation performance indices). In any case, as the number of parameters taken into account increases, this problem represents a combinatorial growth in complexity. 24 At its most basic level, however, a satellite constellation can be composed from a host orbit(s) punctuated by two or more spacecraft. The problem of constellation design therefore involves at least the survey of a trajectory continuum (dictated by a chosen dy- namical system) over which every combination of spacecraft arrangements can be tested. Yielding potentially a very large design space, the task of this research is thus to nd a mechanism with which to eciently map this solution space and down-select the optimal congurations contained within. 1.4 Dissertation Outline Chapter 2 continues this dissertation by providing the mathematical setting and lan- guage used throughout this manuscript. Following a review of the relevant background topics, Chapter 3 builds the theoretical framework for a general constellation design methodology as a parametrized optimization problem, complete with the identication of a solution strategy. Chapter 4 then tailors the general design concept for a specic application to autonomous interplanetary constellations. Simulation results for optimal 7 congurations from four periodic orbit families in the restricted Earth-Moon system are collected and showcased in Chapter 5. Next, Chapter 6 discusses the performance of such simulations and their utility in space ight. Lastly, a summary of the results from this research, including conclusions drawn from numerical experimentation, and directions of possible future research are oered in Chapter 7. Several appendices are also provided at the end to elaborate on certain details covered in the text. 8 Chapter 2 Background The body of research contained herein this dissertation, largely draws from several elds of study, such as space ight mechanics and applied mathematics. It is therefore instrumental in understanding this work to review the underlying concepts and provide the appropriate context and sucient background. This chapter prepares the mathe- matical framework and tools encountered in the design of autonomous interplanetary constellations. Such material can be found in standard references as Arnold, 3 Danby, 11 H enon, 16 Hirsch et al., 19 Keller, 22 Koon et al., 23 Spivak, 36 Szebehely, 38 Tapley et al., 39 and Wiggins, 42 to name a few. 2.1 Models Mathematical models are used to describe a system with a particular set of governing laws and physical properties; in eect, simulating reality up to some desired delity. A model can be viewed as an abstract object whose structure is made precise by qualifying its behavior (e.g. rates of change, ow, etc.) and quantifying its characteristics (e.g. physical constants, forces, etc.). By comparing the dierences between mathematical 9 attributes, certain strengths and limitations of models can be rigorously identied. This distinction is particularly useful when dealing with models that describe the same system but in varying levels of detail. For example, this study is concerned with the problem of spacecraft motion over time, using: (1) an approximated dynamical system and (2) a fuller \unabridged" version based on planetary ephemerides. This section reviews both models, given in order of increasing complexity. 2.1.1 Circular Restricted 3-Body Problem The restricted 3-body problem considers the motion of a particle of negligible mass under the gravitational in uence of two bodies (i.e. primaries) with nite masses denoted m 1 and m 2 (here, m 1 m 2 and assumed to be point masses). In the context of space- ight, this model represents the rst-order approximation to the motion of a spacecraft in the vicinity of two celestial bodies. The motion of the primaries about their mutual barycenter denes the x-y plane. Thez-direction is along the angular momentum vector and completes the right-handed coordinate system with the origin at the barycenter. This approximate model represents a time-independent, conservative, dynamical system . The specic case, where the orbits of the primaries are conned to circular orbits, is considered here and is appropriately termed the circular restricted 3-body problem (CR3BP). 38 A synodic (i.e. rotating) coordinate system is used such that the x-axis is always aligned to and rotates with the syzygy axis (i.e. line extending through both primaries). From this reference frame, the primaries appear stationary, as seen in Fig. 2.1 for the Earth-Moon system. The force eld is derived from the negative Jacobian of a scalar potential eld. 10 Figure 2.1 Schematic representation of the Earth-Moon CR3BP geometry. Five Lagrange points are shown, along with the reference unit circle and equilateral tri- angles. The third-body restricted mass (i.e. spacecraft) is shown in the upper right. Since the general formulation of the CR3BP contains several variables (masses, po- sitions, and velocities) that are inter-related, it is possible to describe the system with fewer parameters than state elements. Conveniently, the system can be parametrized by a single scalar parameter, 2R. This dimensionless mass ratio is dened as, = m 2 (m 1 +m 2 ) 82 (0; 0:5] (2.1) Following this denition three specic quantities are normalized to unity: (1) the sum of the two masses, m 1 +m 2 = 1, (2) the distance between the primaries,kr 2 r 1 k = 1, and (3) the mean motion of the primaries about their barycenter, ! = 1. The latter can also be thought of as the angular velocity of the line of syzygy about the z-axis. The normalizations of (2) and (3) thus dene the dimensionless length and time units, LU 11 and TU, respectively. With this parametrization the gravitational constant, G, is unity by construction. 11 Using the convention of placing the smaller primary on the positive x-axis, the equa- tions of motion in the synodic frame are given by the time-invariant dynamics: x 2 _ y =x (1) (x +) d 3 1 (x 1 +) d 3 2 = U x y + 2 _ x =y (1)y d 3 1 y d 3 2 = U y z = (1)z d 3 1 z d 3 2 = U z (2.2) where U represents the eective potential, and the Euclidean distances (from the third- body to each of the primaries) dened as, d 1 = p (x +) 2 +y 2 +z 2 ; d 2 = p (x 1 +) 2 +y 2 +z 2 Equation (2.2) describes a system of second-order dierential constraints on a three- dimensional vector space, with a state given by X = [x; y; z] T . By re-writing these equations as a system of rst-order dierential constraints on a six-dimensional phase space, 23 with an associated state of X = [x; y; z; _ x; _ y; _ z] T , the vector eld can be written concisely as, _ X = F (X)2 R 6 (an elaboration on dynamical concepts is given later in this chapter). While this transformation is not unique, it is this formulation of the dynamics that is most useful when using numerical integrators to solve the system. For the equations of motion in the sidereal (i.e. inertial) coordinate system, refer to Appendix A for the appropriate transformations. 12 A particular feature of interest to this study is the presence of libration (i.e. stationary) points. Often referred to as Lagrange points, the locations of these equilibria are obtained by nding the critical points of the eective potential eld (r U = 0). The ve known solutions are classied into two categories: collinear and equilateral (see Fig. 2.1). Figure 2.2 Collinear libration points for = 0:2. Distances from each libration point to the nearest primary ( i ) are colored light green. Note that the distance between the primaries is normalized to 1 dimensionless length unit. Determining the locations of the collinear points (y = 0) reduces to nding the value of their abscissas along the x-axis. These solutions are the positive real roots of the quintic polynomials for each pointL 1 (between the primaries),L 2 (far-side of the smaller primary), and L 3 (far-side of the larger primary). For each collinear libration point, a quintic polynomial is obtained via the partials ( U x = 0) of the eective potential. 23 Though ultimately the search is for thex-coordinate of the collinear points, the quintic equations can be simplied when changing coordinates with respect to the nearest libra- tion point ( i ), whose relationship is shown in Fig. 2.2. After solving each corresponding quintic polynomial, the locations of the collinear points are given by the coordinates, The values and i for specic system are provided in Chapter 4. 13 L 1 = (1 1 ; 0; 0) ; L 2 = (1 + 2 ; 0; 0) ; L 3 = ( 3 ; 0; 0) (2.3) The location of the triangular libration pointsL 4 (leading Greek camp) andL 5 (trail- ing Trojan camp) fall out naturally from the dynamics. Appropriately named, these solutions exist for d 1 =d 2 = 1, which furnish the case where y6= 0. From geometry, the locations of the equilateral points are given by the coordinates, L 4 = 1 2 ; p 3 2 ; 0 ! ; L 5 = 1 2 ; p 3 2 ; 0 ! (2.4) With the CR3BP capturing the salient phenomena in a basic model, the next level of delity comes from analyzing tertiary perturbations (e.g. obtained from planetary ephemerides) that augment these dynamics in linear combinations. 2.1.2 Ephemeris Model Often during preliminary studies, ideal models (such as the CR3BP) are chosen as starting points to prove the existence of certain phenomena before treating the full prob- lem. In practice, these restricted models are usually insucient in representing the actual dynamics as they largely ignore tertiary perturbations which, albeit small at rst, can become appreciable when integrated over time. Once the phenomena are conrmed using low-delity models, further investigation is required to extend the work to apply in higher 14 delity models. For space ight, more realistic models of the solar system use the JPL DE series of planetary ephemerides. Considering the same local environment as in the CR3BP (e.g. Earth-Moon neigh- borhood), additional solar system perturbations are now introduced to more accurately re ect the physical system. These perturbations include forces such as third-body grav- itational eects, nutations, and librations. Note that in the presence of these additional perturbations, the approximated time-invariant dynamics of the CR3BP now become time-dependent; that is, the dynamics are relevant to a particular epoch. Refer to Ap- pendix D for conguration settings and an abridged lunar ephemeris gleaned from the JPL DE405 37 library. In order to use the ephemeris model as provided by JPL, the non-dimensional CR3BP coordinates of a given orbit must be mapped from the barycentered synodic frame to the geocentric sidereal frame, and then re-dimensionalized into physical units (see Appendix A for transformations). Once a trajectory is specied in the Earth's inertial coordinate system at a particular epoch, a numerical integration is performed using the transformed CR3BP solution as the initial guess for the perturbed n-body dierential equations of motion. The software package provided by NASA, known as the General Mission Analysis Tool y (GMAT), is used to integrate a desired trajectory using the JPL DE405 ephemeris (see Chapter 5 for simulation results). For the purposes of this study, no further attempt at dierential correction is made to adjust an orbit in the ephemeris model to match a similar orbit in the CR3BP setting. Jet Propulsion Laboratory Development Ephemerides. y http://gmat.gsfc.nasa.gov/ 15 2.2 Dynamical Systems Concepts With the mathematical models establishing the framework for the dynamics, the next task is to dene some useful concepts with which to describe additional features of the dynamics. In part, the theory of dynamical systems seeks to qualitatively describe the behavior of a process as it evolves over time. In other words, the task is to reconstruct the past and predict the future based on the present knowledge of the local laws of evolution of a system. 3 These laws dene the transformations of the state and can be discussed as either continuous (dierential equations or ows, _ X(t) = F (X(t);t)) or discrete (dierence equations or maps,X(t) ='(t; t 0 ;X(t 0 ))). In either case, the state of a process can be represented as a nite tuple of numbers, whose entire collection over time is viewed as the phase space of the system. This section highlights some special structure of these phase spaces and underlines their geometric interpretations as related to space ight. 2.2.1 Invariant Manifold Structures When discussing solutions as belonging to a continuous family of curves, rather than as isolated trajectories, it is convenient to describe the collection using the topological concept of a manifold. The structure of manifolds captures the behavior of parametrized sets of trajectories and species characteristics associated with the family (as opposed to a discrete individual). In this manner, the global properties of an entire solution set can be explored. 16 Along with the representation of solutions as points in phase space 36 (i.e. ordered tuples of numbers), the notion of dimension follows naturally. A manifold has an asso- ciated dimension that is dened by the dimension of a point contained within it. Using terms familiar to geometry, the manifold can therefore be thought of as ann-dimensional surface . Although the entirety of a manifold may be seen as a complicated object, a manifold locally has the well-behaved structure of Euclidean space. For example, a textbook illustration of a generalized n-dimensional manifold is the surface of the unit n-sphere of radius r, embedded in (n + 1)-dimensional Euclidean space, dened as the setS n =fr2R (n+1) :krk = 1g. A 2-sphere is thus the ordinary unit sphere embedded in three-dimensional space. Of particular importance to this study are manifolds that are invariant under the ow of the dynamics. 42 In other words, a point on an invariant manifold will remain on the manifold under the forward (or backward) propagation of time. A glimpse of the deterministic nature of the dynamics is thus imprinted onto the structure of such manifolds. In addition to the convenient visualization provided by the identication of a manifold, its local Euclidean structure allows familiar branches of mathematics to be imported for their study. Operations of calculus, such as dierentiation and integration, now have meaning when applied on a manifold. The use of manifolds in astrodynamics helps in understanding dierent phenomena ranging from foliations of constant energy surfaces that bound spacecraft motion, to collections of periodic orbits that form the admissible space for satellite constellations. Each lower dimension has corresponding projections of a manifold, which are themselves surfaces. 17 2.2.2 Periodic Orbits Among the many interesting features of dynamical systems, Poincar e found periodic solutions (i.e. orbits) to play a critical role in understanding the ow of the dynamics. 16 In general, periodic orbits reveal regimes of motion where initial solutions return to their original state after propagating for some nite amount of time; i.e. a closed orbit . More rigorously, these orbits satisfy the relation,X(t 0 +T ) =X(t 0 ), for some periodT , where X(t 0 +T i )6=X(t 0 ) forT i <T . Note that as the period approaches zero, the periodic orbit degenerates into a stationary point. Thus, all periodic orbits emanate from equilibrium points as the period is varied. Shown in Fig. 2.3 is such a class of orbits known as the vertical Lyapunov family emanating from the Earth-Moon L2 libration point. A good review of periodic orbits in the CR3BP can be found in H enon, 16 with computational techniques supplemented by Romanov. 33 a) Sample subset of twelve orbits b) Smooth continuum Figure 2.3 The Earth-Moon L2 vertical Lyapunov family projected in the synodic position space as (a) a set of discrete orbits and (b) a smooth continuous family; both colored by the continuum index. If the orbit does not strictly repeat exactly after one period, then it is said to be nearly-periodic. 18 A periodic orbit can be succinctly represented by a tuple consisting of an initial position-velocity 6-state and its scalar orbital period, such that (X(t 0 );T )2 R 7 . Therefore, an entire orbit is characterized by a zero-dimensional point inR 7 . A continuous collection or family of periodic orbits (i.e. points) can thus be represented as a one- dimensional curve (i.e. manifold) embedded in a seven-dimensional phase space. Furthermore, since a related set of periodic orbits is naturally formed by a one- parameter variation (e.g. T ), the one-dimensional manifold can be locally parametrized by various other scalar \in-family" parameters, 2 R, as well. Though there are no unique denitions for these in-family parameters, they are typically taken to be the orbital period, Jacobi constant, or continuum index (i.e. any number between 0 and 1). In this fashion, families of periodic orbits can be represented by parametric curves, S() =f ()2R 7 :X(t 0 +T ();) =X(t 0 ;)g. 2.2.3 Bifurcation Diagrams With the existence of several equilibria and associated periodic orbits of dynamical systems, in general, there is a desire for a mechanism with which to organize all of the solutions together. Closely related to the idea of phase portraits, this need is satised with bifurcation diagrams. This map, or graph, allows multiple solution manifolds to be linked as a system evolving from one branch to the next. A bifurcation occurs when small variations in a bifurcation parameter causes sudden changes in the behavior of solution manifolds. The variations cause a forking of the solution manifold, giving rise to new families (i.e. branches). Thus, distinct families that 19 Figure 2.4 Bifurcation diagram of several sets of L1 and L2 periodic orbit families associated with the Earth-Moon CR3BP. The color dimension shows the orbital period per branch, with red/purple indicating high/low values, respectively. intersect in phase space are said to bifurcate from each other. The intersection point is referred to as a branch point or bifurcating orbit. To brie y illustrate when a bifurcation occurs, the stability of solution points is in- vestigated. A branch point is detected along a solution manifold when its corresponding vector eld Jacobian yields an eigenvalue (i.e. Floquet multiplier) with zero real part. A Floquet multiplier equal to exactly zero is termed steady-state, whereas a purely imagi- nary eigenvalue gives rise to Hopf bifurcations. 12,42 Other than mentioning the existence of other avors of bifurcations, such as transcritical, pitchfork, and period-doubling, their details will not be elaborated here and the reader is referred to Wiggins 42 for deeper insight. Bifurcation diagrams oer the convenient ability to highlight particular details by choosing to display appropriate projections of the phase space. That is, this selective 20 characterization allows the summary of dierent information to be captured in a one- dimensional graph. For example, Fig. 2.4 shows several sets of periodic orbit families that are related through bifurcations (two-level deep ) associated with the L1 and L2 equilibrium points. 13 Each point along the curves of this bifurcation diagram represents a periodic orbit, colored according to its orbital period (where the full spectrum is spanned per branch), and projected onto the x-average and z-amplitude space. Equipped with the knowledge of how to detect bifurcations, the systematic procedure of following a branch, detecting bifurcations, switching branches, and continuing along the new branches, enables the discovery and construction of an intricate web of solution manifolds. This process is made possible by a technique known as continuation. 2.3 Continuation In an abstract sense, a continuation principle describes the notion of nding subse- quent solutions based on information gained from an initial solution. On a more concrete level, a practical implementation of interest is captured in a computational technique known as numerical continuation. This section outlines the framework for this scheme and its motivating theorem, as well as focuses on a specic method called pseudo-arclength continuation. In principle, the number of branch switches can continue indenitely, so the depth of bifurcation usually needs to be constrained in order to restrict the space to encompass a manageable set. 21 2.3.1 Numerical Continuation Numerical continuation is a method of approximating solution manifolds of parametrized nonlinear equations, given an initial solution. Following Doedel and Keller, 12,22 in gen- eral, a problem is suitable for one-parameter continuation if it can be written in the form: :R p+1 !R p (u;) = 0; where (;);u2R p and 2R This formulation considers a system of p equations and p unknowns, with the vector of unknowns denoted as u and a free parameter denoted as (i.e. the continuation parameter). Here, another parametrization is made by dening a single point in phase space to bez (u;)2R p+1 . The Jacobian of the mapping, D (z), may now be dened at the initial point as 0 z z (z 0 ). A solution, z 0 , is called regular if the corresponding Jacobian is full rank , Rank( 0 z ) =p, and singular if the Jacobian is rank decient. Furthermore, if the Jacobian, 0 z , is indeed of maximal rank, then by the implicit function theorem, there also exists a unique local one-dimensional continuum of solutions nearz 0 and passes through it; i.e. a solution branch. 12,42 This regularity (or more strongly the guarantee of a locally continuously dierentiable manifold) allows continuation methods to be applied. The task of numerical continuation is then to take an initial solution, z 0 , and nd the next solution,z 1 , along the branch, by iteratively solving (z 1 ) = 0, using informa- tion obtained from the initial Jacobian, 0 z . To better understand the feature enabling continuation, a digression is made to describe this fundamental theorem. The Jacobian has dimension p (p + 1), and thus has at most p linearly independent columns. 22 2.3.2 The Implicit Function Theorem In multivariable calculus, the implicit function theorem gives a sucient condition to determine the existence and uniqueness of local solution manifolds. 36 Consider a continuously dierentiable implicit function, , in an open set containing the point (a 0 ;b 0 ), :R p R q !R p ; with (a 0 ;b 0 )2R p R q If the function satises, (i) (a 0 ;b 0 ) = 0 (ii)Rank(D (a 0 ;b 0 )) =p (2.5) then there exists open sets AR p containinga 0 and BR q containingb 0 , such that, 9! :B!A ( (b);b) = 0; where (b)2A; 8 b2B The injective mapping, , is therefore an explicit function dened implicitly by the equation, (a;b) = 0. Since the manifold M =f( (b);b)2 AB : ( (b);b) = 0g locally overlaps the entire solution manifold dened byM =f(a;b)2R p R q : (a;b) = 0g, it is in principle possible to recover this part of the manifold in the vicinity of (a 0 ;b 0 ), using only the two conditions given in Eq. (2.5). Note thatM (and any MM) is a q-dimensional surface 42 embedded inR p+q . This notion of overlapping sets is shown with an example of the unit circle in Fig. 2.5. Of dierentiability class C 1 . 23 Figure 2.5 Example of the implicit function theorem nding a local solution mani- foldMM in the neighborhood of (a 0 ;b 0 ). The unit circle,M, and subsets,A and B, containing a 0 and b 0 respectively, are also shown. Essentially, the criteria given in Eq. (2.5) ensures the existence of a local tangent space, thereby consequently inferring the presence of an underlying continuum. A powerful conclusion of this theorem is therefore the guarantee that there exists a unique localized portion of the solution manifold, M, that is recoverable without having to know the topology of the rest of the manifold,M. 2.3.3 Keller's Pseudo-Arclength Technique Returning to the discourse on continuation, a particularly robust algorithm is ex- pounded. Although many continuation algorithms exist, the pseudo-arclength technique proposed by H. Keller 22 is advantageous because of its ability to follow folds along so- lution branches. In the following paragraphs, an abridged overview of the generalized concept of pseudo-arclength is provided. 12 24 Keller's pseudo-arclength technique involves parametrizing the solution branch as a function of the arclength,z(s). That is, applying the standard predictor-corrector notion, pseudo-arclength approximates a solution by predicting along the path tangent to the solution branch and then recursively correcting orthogonally to the tangent predictor until convergence. This iterative scheme is illustrated in Fig. 2.6. Since the arclength, s, is used as the continuation parameter, the notation dening a path tangent is thus _ z dz=ds. The initial path tangent can be computed from the null space of the Jacobian , _ z 0 =N ( 0 z ). With the addition of the pseudo-arclength constraint, Keller's method seeks to nd the solution to, (u 1 ; 1 ) = 0 (u 1 u 0 ) T _ u 0 + ( 1 0 ) _ 0 s = 0 (2.6) Using a modied form of Newton's method y as the corrector, pseudo-arclength con- tinuation amounts to solving the following augmented square system of equations: 2 6 6 4 1 u () 1 () _ u T 0 _ 0 3 7 7 5 2 6 6 4 u 1 () 1 () 3 7 7 5 = 2 6 6 4 (u 1 () ; 1 () ) (u 1 () u 0 ) T _ u 0 + ( 1 () 0 ) _ 0 s 3 7 7 5 (2.7) where () is the iteration counter. The initial solution point,z 0 = (u 0 ; 0 ), is given and the initial guess (predictor) for the Newton iteration can be taken asz 1 (0) =z 0 + s _ z 0 . The dierential correction terms represent the changes between iterations, z 1 () The null space is non-trivial because the rank-nullity theorem guarantees the kernel has, at least, dimension 1. y In general, solving (z) = 0 by dierential correction: z (+1) =z () D 1 (z () ) (z () ). 25 Figure 2.6 Graphical representation of pseudo-arclength continuation. Depicted are the initial solution, z 0 , the associated path tangent, _ z 0 , and several predictions during Newton's iteration,z () 1 , as bordered white dots. The converged solution, z 1 , is shown as a bordered white square. z 1 (+1) z 1 () . Lastly, once Newton's iteration has converged, the next (normalized) path tangent can be computed from, 2 6 6 4 1 u 1 _ u 0 _ 0 3 7 7 5 2 6 6 4 _ u 1 _ 1 3 7 7 5 = 2 6 6 4 0 1 3 7 7 5 (2.8) With the core concepts establishing the dynamics (and the tools to compute orbits therein), the next layer of sophistication comes from adding spacecraft into the mix; that is, addressing the nautical part of astronautics. 26 2.4 Spacecraft Navigation The eld of spacecraft navigation typically encompasses two branches of study: (1) orbit determination and (2) ight path control. In this dissertation, the meaning of spacecraft navigation will focus only on the former. This section thus encapsulates the aspects of space ight dealing specically with the reconstruction and prediction of a satellite's position and velocity in a given dynamical system. Furthermore, the act of determining the state of the spacecraft on-board, distinguishes this theory as autonomous navigation. One particular application of note, which enables constellation autonomy, is called LiAISON 17,18 (Linked Autonomous Interplanetary Satellite Orbit Navigation). 2.4.1 Statistical Orbit Determination The problem of estimating a spacecraft's unknown motion, from observations in u- enced by random and systematic errors, using inexact mathematical models, is termed statistical orbit determination. 39 This process can be modeled as a constrained optimiza- tion problem and is generally treated with tools from estimation theory, linear systems theory, and linear algebra. The following discussion formulates the problem of batch orbit determination to provide the context in which the navigation strategies in this dissertation are based. Consider a satellite that is subject to nonlinear dynamics and given a set of observa- tions along with an associated nonlinear measurement model. The optimal solution (i.e. best estimate of the state) of the system must satisfy: 27 _ X =F (X;t); _ X;X2R n and t2R (continuous) Y i =G(X i ;t i ) + i ; Y i ; i 2R d for i = 1;:::;` (discrete) (2.9) where X X(t) denotes the state of the spacecraft, Y i Y (t i ) is an observation set with i corresponding errors, and ` is the number of observation sets (the total number of observations are m =d`). In general, orbit determination problems are overdeter- mined, such that the number of observations far exceed the number of unknown state elements, nm. The function F :R n+1 !R n represents the continuous dynamics and G:R n+1 !R d represents the discrete measurement model. Note, for time-independent systems, the functions are no longer explicitly dependent on time and thus the dimension of the domains are correspondingly reduced by one. However, for the sake of generality, time is kept as an explicit parameter. Traditionally, the system of nonlinear equations given in Eq. (2.9) is solved by lin- earizing about a reference trajectory and then minimizing over the linearized deviations. By dening the state and observation residuals to be x = XX 0 and y i = Y i Y 0 i , respectively, the original nonlinear equations can be expanded in a Taylor series about the nominal states,X 0 andY 0 i respectively, and truncated after the linear term, leading to the linearized dynamics and linearized measurement model, _ x =Ax y i = ~ H i x i + i ; for i = 1;:::;` (2.10) where the corresponding Jacobians, evaluated on the reference trajectory, are dened as, 28 A = @F @X 0 ; ~ H i = @G @X 0 i Through the aid of the state transition matrix, (t i ;t k ), between the measurement time (t i ) and a xed epoch (t k ) at which the state is to be estimated, the general solution to the linear dynamical system can be found using thet-advance mapping,x i = (t i ;t k )x k . Conveniently, this dierence equation is used to express the linearized measurement model in terms of a single unknown state-deviation vector at a common epoch, thus reducing the number of unknowns from n` to n. Using the relation H i = ~ H i (t i ;t k ) and due to the independence of the observations, the measurement model given in Eq. (2.10) may be accumulated to produce, y =Hx k +; where x k 2R n and y;2R m (2.11) Using a technique introduced by Gauss in the late eighteenth century, the method of nonlinear least-squares formulates an optimization problem to approximate a solution to Eq. (2.11). The procedure consists of minimizing a cost function , J, dened as, J(x k ) = 1 2 (x k x k ) T W k (x k x k ) + 1 2 (yHx k ) T W k (yHx k ) where W k is the a priori information matrix of the a priori state deviations x k and W k is the weighting matrix for the observation errors. Further details regarding the origin of The cost function represents the square of the observation errors. 29 this particular cost function and the relation between the weighted least-squares process to that of the minimum variance estimate can be found in Tapley, et al. 39 Satisfying the rst-order necessary conditions to determine extrema (@J=@x k = 0), the so-called normal equations are formed after simplication, x k =Z (2.12) where = H T W k H + W k is the information matrix and Z = H T W k y + W k x k is the data vector. Thus, the optimal solution, that provides the best-estimate of the state deviations, is achieved using the pseudoinverse relation , ^ x k =arg min x k J(x k ) = (H T W k H + W k ) 1 (H T W k y + W k x k ) (2.13) 2.4.2 Observability The solvability of the orbit determination problem relates directly to the notion of observability. That is, a solution in the form of a state vector is said to be fully observable y if each element can be uniquely determined with a given observation set. 39 Conversely, if one or more elements are not observable, then the problem has no unique solution. Since the problem is nonlinear, the process is iterated until a pre-dened convergence criteria is met. y As opposed to the meaning of physically observed. 30 Using the linearized dynamics and measurement model of Eq. (2.10), the observability can be determined from the rank of the information matrix in Eq. (2.12) which can also be expressed as, = W k + ` X i=1 T (t i ;t k ) ~ H T i W k ~ H i (t i ;t k ) (2.14) That is, the linear system is observable if and only if the square matrix, in Eq. (2.14), is nonsingular (i.e. full rank); therefore, this matrix is equivalent to the observability grammian 27 used in linear control theory. Note, though this discrete time observability grammian is accumulated over the t span of interest dened by the measurement arc, the measurements themselves do not directly factor into the information matrix. Since the observability of the system depends on the reference trajectory and the measurement type (e.g. range), the observation set only indirectly aects observability from its contribution to the least-squares solution given in Eq. (2.13), after at least one iteration. In other words, the linear independence of the vectors comprising the observability grammian ensures that the sequence of observations is produced by a distinguishable set of initial conditions. This concept of observability plays an important role in the following navigation strategy, in determining the feasibility of using candidate orbits for autonomous navigation. 2.4.3 The LiAISON Strategy Proposed by Hill, et al., 17,18 LiAISON is an orbit determination strategy that is capable of recovering the absolute states of a pair (or more) of spacecraft using only satellite-to-satellite tracking. As the spacecraft move along their orbits, relative range 31 measurements are processed in a nonlinear lter designed to estimate the orbital states of all participating spacecraft simultaneously, without any supplemental input from the ground (i.e. autonomous navigation). Of course, this technique only works for the limiting cases where the underlying trajectories are suciently distinguishable so that the problem maintains full observability without a priori information; herein lies the novelty of this strategy. Since the LiAISON strategy uses only scalar intersatellite range measurements, it is fundamentally oriented towards pairs of spacecraft. The notion of using pairs also satises the simplest conguration of what can be construed as a constellation . Throughout this dissertation, the reduced set of a pair of spacecraft will be used to illustrate the theory, while tacitly acknowledging that an N-spacecraft conguration is possible. The state space for a pair of spacecraft is thus the span of their positions and velocities, X2R 12 . The associated vector eld can now be described by _ X = F (X). Noting the observation type to be a single scalar (i.e. range), the measurement model is formed as Y i = G(X i ) + i where Y i ; i 2 R, for i = 1;:::;m. Here, G is simply the Euclidean norm between the position vectors, G(X i ) = p (x 1 x 2 ) 2 + (y 1 y 2 ) 2 + (z 1 z 2 ) 2 , and is therefore not an explicit function of time. The LiAISON problem thus takes the form of a nonlinear statistical orbit determination problem. As discussed previously, traditional linearization techniques are employed and a nonlinear regression t is applied to estimate the unknown state deviations of both spacecraft. Although the general LiAISON framework does not restrict spacecraft to share the same orbit, the meaning of LiAISON in this dissertation will be used solely in the context of single-orbit constellations. 32 If the observation errors are modeled well, in a least-squares sense (i.e. the distribution is Gaussian), the convergence criteria becomes easier to compute when handling only one data type and is based on the weighted root-mean square of the estimated errors of Eq. (2.11), WRMS = " 1 m m X i=1 ^ T i W ^ i # 1=2 (2.15) Furthermore, in the case of one data-type, ^ reduces to an m-vector and W becomes a scalar, thus ultimately reducing Eq. (2.15) to WRMS =k^ k p W=m. This simplication is particularly useful in this study for processing orbits with many observations. The success of the LiAISON strategy ultimately depends on the spacecraft states being fully observable. Thus, care must be taken to choose host orbits and to arrange spacecraft along those orbits such that a distinguishable set of intersatellite range measurements can be produced. In this dissertation, this process is referred to as autonomous constellation design. 33 Chapter 3 Constellation Design Methodology A constellation consists of multi-spacecraft congurations evolving in time and per- forming a distributed task, as exemplied by the GPS space segment that provides nav- igation beacons to end-users. Designing such a constellation infrastructure, at its most basic level, involves selecting the host orbits and choosing the appropriate number and placement of satellites along these orbits. All other requirements, such as the purpose and functionality of the constellation, are thought of as constraints used to bound admis- sible congurations or as performance criteria to balance. The challenge is in deciding which conguration suits best; that is, exploring the dierent possible arrangements and isolating several preferred candidates, altogether in an ecient and organized manner. The contribution of this work is, indeed, in the systematic methodology of characteriz- ing the constellation design space; i.e. formulating the constellation design problem as a low-dimensional, nested, optimization problem whose solution space is mapped out via numerical continuation. A particular novelty of this work is that it provides a exible framework that can be applied in general, as opposed to traditional design methodologies that rely on 2-body dynamics. 34 3.1 General Optimization Problem In the general setting, formulating an optimization problem involves two fundamental operations: (1) modeling and parametrizing the system in terms of design variables, = ( 1 ; 2 ;:::; n ), and (2) down-selecting a particular set of design parameters, ^ , to achieve a particular design objective. Often, the second task is formulated as a constrained optimization problem whose solution provides the desired set of design parameters, ^ = arg min J(), with the cost function J representing the design objective. Due to the large number of parameters associated with a multi-body constellation design problem, the diculty in formulating an optimization problem in this fashion, is in nding the appropriate parametrization that yields a manageable set of variables, such that the necessary dependencies of the problem are adequately captured and the problem remains solvable in a reasonable amount of time. This section reviews these two general operations as they pertain to constellation design, highlighting the assumptions and simplications made. 3.1.1 Parametrization The rst step in formulating an optimization problem is determining the appropriate set of parametrization variables. What may seem like obvious choices, a priori, frequently require some modication to properly capture the salient features of the problem while maintaining a manageable sized set. In the following paragraphs, this evolving process is described for the selection of constellation design parameters. 35 Since the design of a constellation involves modeling the motion of the satellites over time, in eect, observing the behavior of solutions of the dynamical system, _ X(t) = F (X(t);t), a natural parametrization is to consider the states of each spacecraft as de- sign variables, X i (t) for i = 0;:::;N 1 where N denotes the number of spacecraft forming the constellation . Fortunately, due to the deterministic nature of the dynamics, an immediate reduction of this variable space can be made, such that the description of an entire trajectory can be reduced to a phase point at a given epoch, t 0 , asX i (t 0 )2R 6 . That is, this simplication conveniently allows the specication of a single spacecraft tra- jectory to be captured in one six-dimensional state vector. In this fashion, the collection of N satellites can be characterized by 6N +1 variables (accounting for the scalar epoch as a possible optimization variable). With an iteration of the possible set of design variables at hand, the next step is to ensure that the parameters satisfy the conditions of the problem. In this case, such a collection of states must be further constrained to satisfy the notion of periodicity; that is, a constellation is generally understood as a formation of satellites with a near repeating pattern, rather than an unstructured cluster. Looking at existing constella- tions as benchmarks, each spacecraft orbit is kept nearly periodic in order to minimize stationkeeping requirements and maintain a desired closeness to the nominal orbit. Thus, in keeping with the same paradigm, a quasi-periodicity constraint is imposed such that, kX i (t 0 +T i )X i (t 0 )k < , for some tolerance, , and period, T i . In turn, these con- straints add N more parameters (one for each spacecraft); i.e. the quasi-periods, T i for i = 0;:::;N1, resulting in 7N+1 design variables. While this parametrization represents The rst point is indexed with 0. 36 a signicant improvement over the starting design space, there remains the possibility of reducing this set even further by considering a dierent set of parameters (other than the state coordinates) and still capture the main dependencies of the problem without loss of generality. From the standpoint of preliminary design, an approximate model of the dynamics, that faithfully captures the essence of the dynamics, is preferred over a more complicated, higher delity model. As such, the proposed interplanetary constellation design strategy initially leverages the periodic orbit structure of the familiar CR3BP . Besides the avail- ability of a large class of periodic orbits, another advantage of using the CR3BP is the simplication aorded by the time-invariant dynamics. This property removes the need to optimize over the epoch, t 0 , which also shows that t 0 has, a priori, little in uence on the design problem in the actual dynamics. In fact, over the course of the constellation life-time, any relative orientation is expected to be encountered, thus further reinforcing the removal of the unnecessary optimization variable with respect to a specic epoch; leaving 7N design variables. The next reduction in the variable space is not strictly essential, but rather is seen as a necessary step before proceeding to a more complicated case. More specically, the optimization variables are further reduced by restricting allN spacecraft to lie on a single periodic orbit y ; i.e. only considering single-orbit constellations z . This constraint allows the location of each spacecraft to be represented by a simple phasing relative to the \lead" For constellations about one central body, the well-known Keplerian structure of conic sections suces. y For more general Lissajous or quasi-periodic orbits, two phase parameters are required instead of one. z Indeed, this assumption represents a limiting case, but is seen as a fundamental building block for other studies involving multi-orbit constellations. Additionally, using a single orbit reduces the specication of dierent periods, Ti, to just one period, T . 37 spacecraft, as illustrated in Fig. 3.1. The phase angle, 17 analogous to the Keplerian mean anomaly, is dened similarly as, i = 2 (t i t 0 ) T ; i 2 [0; 2] (3.1) Figure 3.1 Phase angle parametrization, i , along an L1 halo orbit projection. Although the designation of the lead spacecraft is somewhat arbitrary, it can be made precise by specifying its state to coincide with a subset of the same phase point that de- nes the host periodic orbit itself, namelyX 0 (t 0 ) (). The states of the remaining (N 1) spacecraft are thus given by thet-advance mappings,X i (t 0 ) ='( i (T=2); t 0 ;X 0 (t 0 )) fori = 1;:::;N1, using the time displacements implied in Eq. (3.1). That is, the state at any time t i on the orbit is found via a xed propagation of the state at time t 0 under the action of the ow. The periodic orbit, ()2R 7 , minus the specication of the period, T , leavesX(t0)2R 6 . 38 This \time-phasing" parametrization of spacecraft locations, along with the local parametrization of periodic orbit families, nally reduces the number of optimization variables to N parameters; i.e. changing the state coordinate parametrization to N1 phase angles ( 2 R N1 ) plus one scalar in-family parameter (2 R). In other words, given a periodic orbit family , this nal parametrization scheme allows to select a specic orbit from the continuum 13 and to specify the placements of the satellites along the orbit relative to a lead spacecraft 17 (whose location is simultaneously determined by ). The nal set of parameters can be succinctly described as the N-tuple, = (;). 3.1.2 Formulating the Minimization Problem With the establishment of the appropriate parametrization of the design space, the next task is to formulate the optimization problem. The idea of using the structure of an optimization problem in the rst place, comes from the notion of wanting to select the \best" conguration. Though the meaning of best has not been dened, the problem can still be formulated, in general, even in the absence of an actual concrete performance measure. In particular, the constellation design problem can be cast as a minimization problem of the form, (^ ( ^ ); ^ ) =arg min ; J N (;); J N :R N !R (3.2) where the solution lies in the domain of the generic cost function, J N (still yet to be determined but should at least be a function of the design variables). While, in principle, In principle, the subspace of admissible orbits can be computed in advance. 39 the positive integerN could be included as an optimization variables, it is instead consid- ered here to dene a discrete parametric class of optimization problems; that is, one for each N2Z + , with N = 2 being the minimum required to be considered a constellation. Thus, for a given number of constellation spacecraft,N, thisN-dimensional optimization problem seeks to nd the most eective periodic orbit (orbit down-selection) on which the best performing conguration can be arranged (spacecraft placement). The coupled optimization problem, given in Eq. (3.2), can be thought of as being composed from a sequence of nested optimization problems as follows: (^ ( ^ ); ^ ) =arg min farg min J N (;)g (3.3) This formulation makes it more apparent that the problem is to nd the optimal place- ments of N satellites for each periodic orbit, (^ ();), and then nding the optimal of this set, (^ ( ^ ); ^ ), by performing a one-dimensional line search over , such that J N (^ ( ^ ); ^ )J N (;) for all (;)2R N . 3.1.3 Existence of Local Extrema Before attempting to solve the problem, a check should be performed to determine whether or not one or more solutions are available. Using Eq. (3.1) that denes a phase angle as belonging to a closed circle, 2 [0; 2], the space of all phase angles can be generalized to form a compact set described by the topology of the (N1)-torus. Fur- thermore, by setting the periodic orbit in-family parameter to be the continuum index, this scalar variable can also be mapped onto the closed unit interval, namely 2 [0; 1]. 40 Therefore, by applying the extreme value theorem, the existence of at least one minimum (and a corresponding maximum) is guaranteed over the compact space of the \thick" (N1)-torus formed by the direct product of both sets . That is, there must exist at least one solution to the minimization problem posed by Eq. (3.2). Though the existence of a pair of extrema can be proven using the properties of compactness, there is no such guarantee that these are unique. In fact, multiple local extrema are expected to occur, and as such, the proposed design methodology should not preclude the exploration of such additional solutions. 3.1.4 Scope of the Solution Space Given the current set-up of the problem and the assurance of at least one solution, some attention is now necessary to scope the level of eort required to solve the problem. Without imposing any intuition regarding the structure of the solution space, this space will a priori span the subspace ofR N associated with the thick (N1)-torus that forms the domain of J N . By itself, this is a continuous space that, even when discretized for processing, can yield an extremely large y solution space to survey. The most straightforward way to solve this problem, is to evaluate each permutation of and with the cost function J N (;) to test for minimality; i.e. an exhaustive search z . To get a rough sense of the number of computations involved, suppose the phasing intervals are divided into one degree increments, and the number of periodic orbits to investigate is 100. Thus, for a pair of spacecraft, the number of permutations The union of the (N1)-torus with the unit interval, results in a collection of (N1)-tori, or for short, a \thick" (N1)-torusR N . y The exact size of the solution space is a function of the increment of discretization. z While the most cumbersome, a \brute-force" grid search gives an upper bound on the level of eort. 41 is 36; 000. For three spacecraft, the required computations reach 12:96 million. This number, of course, increases exponentially with the addition of more spacecraft and is directly proportional to the number of periodic orbits considered within a family. With this \brute force" search technique quickly becoming insucient to handle the exponentially growing problem, alternate solution methods are desired to circumvent these inevitable computational burdens. 3.2 Solution Methodology While anN-dimensional optimization problem admits a potentially huge search space, certain innovations are leveraged to avoid direct computation of the entire set. Among the many techniques available, the selected methodologies addressed in this dissertation are all numerically based. This section covers several techniques applied to solve the full optimization problem, starting with the attempted trials and ending with the nal implemented strategy. 3.2.1 Heuristic Approach and Local Optimizers Initially, the strategy applied to solve the optimization problem involved using heuris- tic arguments to abbreviate the problem, such that the full problem need not be solved explicitly. The heuristic approach decouples the problem into an orbit down-selection sub-problem followed by a single spacecraft placement step. This separation of the prob- lem is, in principle, possible if there exists a metric that can dependably approximate the results of the spacecraft placement problem. If such a metric is used, the full problem Subsequently discarded due to the discovery of a better alternative. 42 breaks into two sub-problems: a one-dimensional and an (N1)-dimensional optimization problem, ^ =arg min J N () followed by ^ ( ^ ) =arg min J N (; ^ ) whereJ N () represents an easily computed function that mimics optimal satellite con- gurations (i.e. the heuristic metric). Although this method has the potential to reduce the number of computations to the order of the number of periodic orbits, the particular diculty with this approach is the non-trivial task of actually nding an appropriate \heuristic indicator" y . At each stage of the decoupled problem, numerical optimizers (e.g. fmincon z ) were used in an attempt to solve for the respective global minimums. However, as is well known, numerical optimizers suer from a strict dependence on initial guesses. This reliance on the initial estimates restricts the utility of this technique to nding only local extrema in the neighborhood of convex basins of attraction. These local optimizers, therefore, make unsuitable tools for nding global extrema, unless the ability to guess suciently close to the global solution is available. The unsuccessful tactic of dodging the responsibility of explicitly solving the optimiza- tion problem, by way of an approximate heuristic metric supplemented with sensitive local optimizers, warranted a complete re-thinking of the solution strategy. Can be theoretically computed alongside the dynamics at minimal extra cost. y Though never found, a simple test metric was instead introduced to validate the methodology. z http://www.mathworks.com/help/toolbox/optim/ug/fmincon.html 43 3.2.2 Mapping the Design Space using Continuation With the previous methods proving ineective for solving the full N-dimensional problem, a more global approach was adopted to survey the structures of entire solution manifolds rather than just at critical points. This way, less focus is placed on the nal solution, and more attention is drawn towards the overall behavior of the solutions as they evolve over the phase space. Though producing ultimately the same result, the added benet of solving the problem in this manner, is in gaining a better understanding of the topology of the whole solution space along the way. To accomplish this global mapping approach, the solution strategy hinges on the techniques of numerical continuation; in particular, Keller's pseudo-arclength technique. An advantage of writing the full optimization problem as the nested equation in Eq. (3.3), is that the innermost optimization problem has the structure of the root-nding problem: @J N (;) @ = 0 2R N1 (3.4) In other words, the roots of this optimization sub-problem are the critical points of the underlying solution space for a particular choice of (i.e. the spacecraft placement problem for a particular periodic orbit). Over a given set of (i.e. for a periodic orbit family), the innermost optimization problem thus represents the space of interest, making it the primary focus of this method. Admitting N1 equations and N1 unknowns plus one known parameter , this system is well suited for a one-parameter continuation scheme. Note, though there is some freedom in choosing an appropriate cost function to model the desired objective, 44 care should be taken to select a suciently well-behaved function,J N (;)2R, such that @J N (;)=@ is smooth in (the dependence on is assumed to be known). In the spirit of generality, the preparation of the optimization problem, for use with pseudo-arclength continuation, leaves J N temporarily undened. Appending the pseudo-arclength constraint to the system of dierential equations, as outlined in Eq. (2.6), gives the augmented set of equations, @J N ( 1 ; 1 ) @ 1 = 0 ( 1 0 ) T _ 0 + ( 1 0 ) _ 0 s = 0 where the vector of unknowns is the (N 1)-tuple of phase angles, 1 , and the free parameter is the in-family parameter, 1 . The step-size along the pseudo-arclength (i.e. path tangent) is denoted as s. The framework of the modied Newton's method, given in Eq. (2.7), can now be populated as, 2 6 6 4 1 () 1 () _ T 0 _ 0 3 7 7 5 2 6 6 4 1 () 1 () 3 7 7 5 = 2 6 6 4 ( 1 () ; 1 () ) ( 1 () 0 ) T _ 0 + ( 1 () 0 ) _ 0 s 3 7 7 5 (3.5) where, ( 1 () ; 1 () ) = @J N ( 1 () ; () 1 ) @ 1 () (3.6) 1 () 1 () = " @ 2 J N ( 1 () ; () 1 ) @ 2 1 () @ 2 J N ( 1 () ; () 1 ) @ 1 () @ 1 () # (3.7) 45 Recall, since the pseudo-arclength technique requires an initial guess for the solution of the root-nding problem given in Eq. (3.4), a single exhaustive search is performed at 0 to nd the initial set of 0 corresponding to specic local extrema . Additionally, the initial path tangent to the solution manifold at ( 0 ; 0 ) is required and can be computed from the null space of the Jacobian y given in Eq. (3.7). Note that in this formulation, if the analytical derivatives are not available, then a discrete nite dierence approximation may be substituted for each unknown partial derivative. The continuation process proceeds by iteratively solving Eq. (3.5). The new solution is then set as the initial guess for the next run, and the associated new path tangent is computed via, 2 6 6 4 1 1 _ 0 _ 0 3 7 7 5 2 6 6 4 _ 1 _ 1 3 7 7 5 = 2 6 6 4 0 1 3 7 7 5 In this fashion, an entire solution manifold is recovered from a single seed solution and the knowledge of the local curvature at that point. This continuation process can, a priori, be repeated for as many local extrema as desired. In fact, separate continuation runs are necessary if solution manifolds encounter bifurcating solutions (folds in the branches are continued through unhindered). In the event of a branch termination or in ection z , another exhaustive search can be performed to retrieve the associated local extrema for the current orbit, and thus provide more potential avenues to follow. Altogether, the collection of solution branches obtained This continuation process is required to run separately for an initial guess at each local extremum. y Though second-order derivatives are present, Eq. (3.7) is considered the Jacobian of Eq. (3.6). z The transition of a local minimum to a local maximum, or vice versa. 46 by this method can be succinctly represented in a bifurcation diagram, with as the bifurcation parameter. Mapping the design space using continuation algorithms signicantly reduces the num- ber of computations involved in solving N-dimensional optimization problems. Due to the exponential growth in number of permutations asN increases, the savings in compu- tation using continuation is much more pronounced at larger N values, as exemplied in Table 3.1. For example, using the approximation that it takes roughly 10 iterations for Newton's method to converge, at one degree increments in , for 100 periodic orbits, the comparison gives: Exhaustive Search: 360 (N1) 100 Continuation: 360 (N1) + 10 100 Table 3.1: Number of Computations N Exhaustive Search Continuation 2 36:00 10 3 1:360 10 3 3 12:96 10 6 130:6 10 3 4 4:666 10 9 46:66 10 6 Moreover, the intricate branch structure of local extrema is revealed as a tangled web of one-dimensional manifolds embedded in N-dimensional space; that is, a branch can be represented by a QN matrix, where Q is the nite quantization of a manifold. Although the focus has been directed towards the structure of the solution space spanned by the innermost optimization problem, if the global optimums are still desired, they can No matter the size of N, the structure of the problem will always yield one-dimensional manifolds. 47 be easily found by a one-dimensional line search over each solution branch computed in the phase space; i.e. completing the outermost optimization step of Eq. (3.3). With the general framework established for solving the constellation design problem as an optimization problem, the only task left is to determine the specic cost function, J N (;), that captures the desired notion of what is considered \optimal". 48 Chapter 4 Tailoring the Constellation Design Framework Since the previous chapter formulated the constellation design problem as a general optimization problem, specic quantities related to a particular problem were left to be populated. The features remaining to be rigorously dened are the vector eld associated with a particular planetary system (i.e. the dynamics) and the cost function (i.e. perfor- mance index) to quantify the notion of optimality. This way the problem can be exibly solved for individually tailored cases. With the motivation of this work stemming from the goal of lunar colonization, the interplanetary constellations considered in this disser- tation fall mainly under the regime of the restricted Earth-Moon system. This chapter focuses on tailoring the optimization framework for this interplanetary system, using a cost function that measures navigational autonomy; hence, an example of autonomous interplanetary constellation design. 4.1 Modeling the Dynamical System In the context of this research, the preliminary design eort models the spacecraft dynamics as a CR3BP. This section gives specic details related to the Earth-Moon 49 system, including the identication of selected subsets of periodic orbits that form the admissible constellation design space; i.e. the invariant manifold structures reviewed in Chapter 2, here form a non-contiguous subspace of periodic orbit families. 4.1.1 Restricted Earth-Moon System While actually a multi-body dynamical system with a variety of perturbations, the Earth-Moon neighborhood is approximated here as an unperturbed restricted 3-body system for the sake of prototyping. In this fashion, the familiar CR3BP model can be imposed to simplify the design space, while still capturing the essential characteristics of the dynamics. Also, this way, a diverse selection of truly periodic orbits are available to study. In order to standardize the reproducibility of the results simulated in this research, Table 4.1 and Table 4.2 list the physical and CR3BP constants relevant to Sun-Earth- Moon neighborhood (with the Sun related parameters given only for comparison). Since the length scales dier so drastically between the Sun-Earth and Earth-Moon systems, a dimensionless distance error on the order ofO(10 6 ) yields a physical error of100 km and1 km, in each system respectively. Likewise, with such large time dierences, periods for periodic orbits can vary tenfold between the two systems (making observation t spans realistically dicult to maintain over larger intervals). So, although certain Sun-Earth periodic orbits may encompass a greater Earth-Moon neighborhood, the large length/time scales make them impractical for consideration as cis-lunar constellations. With the addition of any perturbations, these \true" periodic orbits are only nearly-periodic. 50 The masses and radii of the primaries are obtained from the online Astronomical Almanac 2011 . The gravitational parameters, , are computed from Eq. (2.1), and the orbital periods and distances between primaries are read directly from JPL's Solar System Dynamics website y . Using the dimensionless conformity of the CR3BP, the period/2 denes the time unit (TU) and the distance denes the length unit (LU). The universal gravitational constant is G = 6:67428 10 20 km 3 /kg/s 2 . Table 4.1: Sun-Earth-Moon Neighborhood Physical Constants Body Mass (kg) Equatorial Radius (km) Sun 1.9884 10 30 696000.00 Earth 5.9722 10 24 6378.14 Moon 7.3458 10 22 1737.40 Table 4.2: Sun-Earth-Moon Neighborhood CR3BP Constants System Period (days) Mean Distance (km) Sun-Earth 3.04045 10 6 365:25 1.49598 10 8 Earth-Moon 1.21505 10 2 27:32 3.84400 10 5 Every CR3BP system admits ve equilibrium points. The coordinates of each libration point is given by Eq. (2.4) and Eq. (2.3). Table 4.3 lists the distances of the collinear points from their nearest primary, in non-dimensional units. Table 4.3: Distances of Collinear Libration Points from Nearest Primary System 1 2 3 Sun-Earth 0:0100110 0:0100783 0:9999982 Earth-Moon 0:1509340 0:1678323 0:9929121 Stemming from these libration points, emanate families of periodic orbits. Of particular use for cis-lunar constellations are periodic orbits that are associated with the L1, L2, http://asa.usno.navy.mil/ y http://ssd.jpl.nasa.gov/ 51 and L4/L5 points (with L3 orbits omitted due to being too distant from the Earth-Moon corridor). 4.1.2 Subsets of Periodic Orbit Families Due to the nature of the proposed optimization structure, the construction of the orbit space (to host the constellation satellites) is conveniently decoupled from the satellite placement problem. Subsets of trajectory manifolds can thus be computed prior to any assignment of number and/or arrangements of spacecraft. For the purposes of this work, the space of interest is set to be a subspace of Earth-Moon periodic orbits. Each periodic orbit family is numerically continued from an associated libration point by variations in the periods. AUTO y (a continuation and bifurcation software) is used to compute the initial conditions for the periodic orbits (see Appendix C for a list of initial conditions), and then the post-processing is performed in Matlab z to complete the manifolds using the built-in ordinary dierential equations integrator ode113, with associated relative/absolute tolerances of 10 9 and 10 12 , respectively. Since AUTO uses a collocation method to dierentially correct periodic orbits rather than a multiple- shooting algorithm, slight deviations from nominal periodic orbits are observed in the initial conditions when re-integrated in Matlab. As such, residual imperfections from the AUTO-generated initial conditions are rst polished via a two-level, asymmetric, constrained dierential corrections process 26 in Matlab, before any further simulations. The L4 and L5 vertical families are actually the same family. 13 y http://indy.cs.concordia.ca/auto/ z http://www.mathworks.com/ 52 a) Sample subset of twelve orbits b) Smooth continuum Figure 4.1 A selected subset of the L1 axial family, plotted in the synodic frame. a) Sample subset of twelve orbits b) Smooth continuum Figure 4.2 A selected subset of the L1 halo family, plotted in the synodic frame. 53 a) Sample subset of twelve orbits b) Smooth continuum Figure 4.3 A selected subset of the L2 vertical family, plotted in the synodic frame. a) Sample subset of twelve orbits b) Smooth continuum Figure 4.4 A selected subset of the L5 vertical family, plotted in the synodic frame. 54 a) Sample subsets of the non-contiguous design space b) Smooth continuums of the non-contiguous design space Figure 4.5 The admissible Earth-Moon periodic orbit subspace, plotted in the syn- odic frame. The Earth and Moon are shown as the large and small black spheres (exaggerated in scale for visualization), respectively. 55 Four classes (i.e. families) of periodic orbits are chosen to represent the admissible subspace for cis-lunar constellations: (1) L1 axial { Fig. 4.1 and Table C.2, (2) L1 halo { Fig. 4.2 and Table C.1, (3) L2 vertical { Fig. 4.3 and Table C.3, and (4) L5 vertical { Fig. 4.4 and Table C.4. The subsets of periodic orbits are truncated in order to avoid diculties with the LiAISON strategy (to be discussed in the Chapter 5). The plots represent three-dimensional synodic projections of periodic orbit invariant manifolds. In each gure, a sample subset of twelve orbits is paired with its corresponding continuum to highlight a particular family's constituent elements. Each family is showcased separately, concluding with a composition containing the entire orbit subspace of all four classes together (see Fig. 4.5). The color gradient of each plot visually represents the unit interval. That is, each periodic orbit manifold is parametrized by the scalar in-family parameter 2 [0; 1], which in turn is mapped onto a full chromatic spectrum. The extremities correspond to the bounds of the unit interval; i.e. red indicates 0 and purple indicates 1, and with each number in between matching a unique color within the gradient. 4.2 Choosing the Cost Function Once the orbit subspace is constructed, the last remaining detail to be specied to fully tailor the general optimization problem, is to choose the performance index as a function of the parametrized variables, J N (;). This cost function will determine in what sense a particular conguration is optimal. Since the notion of optimality will vary Arbitrarily shortened at each end of the continuum; more noticeable in the larger period orbits. 56 per mission objective, the denition of this cost function is not unique. This section considers two possible cost functions which are relevant for constellations used to support lunar colonization eorts: (1) a factor that captures the end-user navigation performance, and (2) a metric that measures navigational autonomy between constellation spacecraft. 4.2.1 Geometric Dilution of Precision Assuming, for now, that the primary objective of the nal constellation is to provide navigational support for cis-lunar trac (i.e. an extended GPS system for the Earth- Moon neighborhood), a reasonable notion of an optimal constellation would be one whose conguration produces the most accurate point-positioning within its serviceable region. Borrowing from the current Earth-based GPS infrastructure, a common metric used to quantify the accuracy of position xes, obtained via guide-beacon navigation, is the scalar factor called DOP 20 (dilution of precision). The core foundation of this concept is based on calculating the unknown positions of receivers with respect to several disparately placed transmitters with (relatively well) known locations. In eect, this type of point-positioning is an abstraction of the famil- iar geometric process of triangulation, using pseudo-range measurements (together with receiver clock-error corrections ) instead of actual ranges. Though, in principle, the position of an observer can be determined using any subset of four spacecraft, the quality of the measurement is not the same for each combination. That is, the accuracy of the user location estimate depends on the geometric placement of the constellation spacecraft with respect to the observer. Furthermore, this estimate With the addition of clock-error, the minimum number of satellites required for a position x is four. 57 is also a function of the precision with which the locations of the beacons are known, as well as the precision of the range measuring devices themselves . In fact, beacon location uncertainties can be considered to be pseudo-range measurement errors, whose variances can be denoted generally as 2R + . In this context, the covariance matrix is given by, P = 1 , where is the informa- tion matrix formed by the outer product of the Jacobian of the measurement model with itself. All of the dilution of precision metrics are garnered from this covariance matrix. Part of the general DOP family y , GDOP z specically measures the overall eect of the geometric placement of satellites on the navigational precision attainable by the end-users, who are using these spacecraft as guide-beacons. This dimensionless factor is captured by the trace of the covariance matrix obtained from the linearized observation- state relationship, GDOP = p Trace ( 1 ) where = @ 0 @x T i @ 0 @x i where the state deviation associated with the end-user is dened as,x (X; Y; Z; i ), and the measurement model is given by the Euclidean distance, 0 i =kX i X 0 k, between the end-user's reference position,X 0 2R 3 , and each satellite,X i 2R 3 , for i = 1;:::;N. Therefore, in order for a solution,x2R 4 , to exist, the number of spacecraft in view must match at least the number of unknowns, (i.e. N 4). A priori this hardware-based error is largely neglected. y The DOP family consists of ve components: position (PDOP), time (TDOP), horizontal (HDOP), vertical (VDOP), and geometric (GDOP). z GDOP is the most versatile because it encompasses the full spectrum of the covariance matrix and is invariant under similarity transformations. 58 a) Good geometry b) Bad geometry Figure 4.6 Example congurations for (a) \good" and (b) \bad" 4-spacecraft ge- ometries (red dots), with respect to an end-user (central black dot), as seen mapped onto the unit sphere. A convenient way to visualize GDOP is to consider the geometry described by the Jacobian of the measurement model (@ 0 i =@x). The rst three columns of this Jacobian dene the unit vectors that point from the beacon satellites to the end-user's reference position. These unit vectors, in turn, describe a polyhedral surface (with the number of sides equal to the number of satellites, as shown in Fig. 4.6 for four satellites). Lever- aging the properties of Euclidean geometry and linear algebra, the volume enclosed by this surface is measured by one-half of the determinant of the matrix containing the unit vectors. The insight gained from this geometrical abstraction of the Jacobian, is that large serviceable volumes (corresponding to \good" spacecraft geometries), mirror \low" GDOP values. In other words, to achieve better measurement accuracy, the ideal place- ment of satellites, relative to the observer, should be widespread in the sky with large angular separations between them. 59 a) Good service volume b) Bad service volume Figure 4.7 Volumetric plots of GDOP values in a region about a sample Earth- Moon L1 halo orbit, where the white dots show the locations of eight satellites at a given instance in time. Shown in the color dimension are the accuracies seen by an end-user using a constellation that has spacecraft (a) evenly spaced in phase angle, and (b) unevenly clustered towards one side. While this metric captures the main factors in uencing the end-user positioning ac- curacies, it is not the best constellation design cost function as it strongly depends on end-user location and time ; i.e. it is arbitrary (see Fig. 4.7). That is, a priori the user location is unknown and the serviceable volume is undened in the design phase, so the cost function cannot strictly depend on these additional parameters. This observation leads to the conclusion that GDOP-based metrics should only be used as a posteriori characterizations of constellations that match a discrete set of potential solutions to a particular problem. Instead, attention should be redirected towards the uncertainties in Though a four-dimensional integral over time and volume averages this metric for a particular con- stellation, the limits of integration are, a priori, also design parameters. 60 the locations, , of the constellation spacecraft themselves, as this parameter is more intrinsic to the problem. 4.2.2 Estimation Error Covariance Digressing from the utility of the constellation to end-users, focus is now shifted to the navigational performance of the constellation. Addressing yet another, arguably more precursory, goal of lunar colonization, comes the notion of low-cost space mission designs. From the standpoint of infrastructure planning and development, the feasibility of the outcome greatly depends on the resources required to build and maintain the desired architectures. With respect to a constellation, a leading budgetary concern stems from the recurring costs associated with ground-support and in-space maintenance. 32 One strategy, proposed to mitigate these persistent expenditures, is to utilize some form of spacecraft autonomy. With subsets of periodic orbits already chosen as baseline trajectories for the con- stellation design space, the problem is well primed to accept the novel mechanism for achieving navigational autonomy called LiAISON. 18 Since this method, which is based on intersatellite range measurements, shows that it is possible to achieve a mutual under- standing of orbital states independently from third-party involvement, it should a priori palliate scal concerns regarding the expenses associated with monitoring/tracking. In other words, autonomous constellations would require less ground-based resources to op- erate, and given a highly accurate means of orbit determination, would also require a Assuming at least the structure of non-Keplerian orbits oered by the CR3BP. 61 smaller fuel budget for stationkeeping. As such, it is this navigation strategy that will provide the basis from which a cost function is derived. As noted in a previous chapter, the fundamental orientation of LiAISON towards pairs of spacecraft, allows the discussion of a constellation as being composed from multiple satellite pairs. Though several pairwise connectivities are possible, such as pairs formed by nearest neighbors (sequentially connected) or independent oset pairs (disjoint), the simplest case is indeed that of a single pair. The case exemplied here, for the purposes of verifying the methodology, is thus the base set of a 2-spacecraft constellation. Recall, that the state of a 2-spacecraft conguration is the union of the positions and velocities of both satellites,X2R 12 . Therefore, in the least-squares formulation of the orbit determination problem, the linearized system yields a square-symmetric information matrix of dimension, 2R 12 R 12 . And, if the system is fully observable, the associated estimation error covariance matrix is given by, P = 1 . Originally described by Hill et al., 17,18 one measure of navigation performance of a particular constellation is proportional to the sum of the eigenvalues of the batch weighted least-squares 39 position covariance sub-matrix associated with the lead spacecraft. In this vein, estimation errors are summarized by taking successive averages over the t span and then over combinations of congurations for a particular group of spacecraft, i = q Trace [P i ] 33 ; j = 1 m m X i i ; ave = 1 N N X j j Performance is quantied by the magnitude of the state estimation errors. 62 for m number of observations in the tting interval (with i indexing time) and N space- craft. Without loss of generality, a slight adjustment is made to these denitions; namely, the entire spectrum of the estimation error covariance is used instead of just a subset of eigenvalues and the square-root is omitted , leaving a more universal performance index, = 1 mN N X j " m X i Trace(P i ) # j (4.1) This scalar,2R, is indeed the same variable used for the spacecraft state uncertain- ties in the GDOP analysis. In a sense, this performance metric provides a conservative upperbound to both the position and velocity uncertainties. Equation (4.1) thus cre- ates a viable constellation design cost function by dening a dimensionless performance characteristic for an entire constellation, over one full period. Furthermore, because the LiAISON concept is a function of the constellation parameters, the phase space is also the space of the \thick" (N1)-torus spanned by the set = (;). The nested optimization framework of Eq. (3.3), can now be rigorously described by setting, J N (;) =(;) (4.2) With this choice, the best performing LiAISON-guided constellation satises the no- tion of optimality; that is, the \best" constellation is one that can autonomously deter- mine the absolute states of all participating spacecraft, with the least amount of uncer- tainty. Trading the description of standard deviations for variances. 63 Following the design approach of Chapter 3, mapping the constellation design space thus involves rst nding several \seed" solutions (corresponding to local optimums) and then \growing" entire solution manifolds from them using methods of numerical continuation. 64 Chapter 5 Results for the Earth-Moon System With the dynamical system chosen as the Earth-Moon CR3BP and the cost function measuring LiAISON navigation performance, four case studies are examined to illustrate the design of autonomous interplanetary constellations. Each simulation represents actual navigation simulations, complete with noise and initial state uncertainties. This chapter begins with an investigation of multi-spacecraft arrangements for a sample L1 halo orbit, in order to exemplify a standalone spacecraft placement problem before moving on to the continuation of optimal solutions over the prescribed orbit manifolds. The next section captures the results for a pair of satellites evolving separately over the L1 halo, L1 axial, L2 vertical, and L5 vertical periodic orbit families. Four corresponding optimal congurations are then subsequently identied via one-dimensional line searches over each branch contained in the bifurcation diagrams (i.e. design space maps). 5.1 Multi-Spacecraft Arrangements In solving the proposed optimization problem over and , recall that a decoupling is noted between the construction of the host orbits and the arrangement of satellites. As 65 mentioned in Chapter 3, this independence allows the general problem to be decoupled into a periodic orbit continuation step separated from a spacecraft placement sequence. In other words, given a periodic orbit, the spacecraft placement problem can be solved locally for an optimizedN-spacecraft conguration (^ on any). The spacecraft placement sub- problem is thus the innermost optimization problem of Eq. (3.3). By separately analyzing the spacecraft placement problem on a xed baseline orbit, a better understanding of the structure of local extrema for that orbit is gained. To illustrate this optimization sub-problem, pairs of spacecraft on a representative L1 halo orbit have been analyzed and the solution spaces are plotted on a log scale to accentuate the details (see Fig. 5.1). Note, for the case of three or more spacecraft, the constellations are a synthesis of coupled pairs of satellites on the same periodic orbit. Key results from this study reveal the presence of many local extrema as well as expose the computational burden of such simulations. Although the existence of a solution is guaranteed by the extreme value theorem, there is no such assurance that this solution is unique. In fact, the presence of several local extrema indicates that there exist congurations of spacecraft that alternate between \good" and \bad" placements; i.e. the topology of the solution space of Eq. (4.1) is not convex. Since one local extrema may oer dierent properties than another, it is of value to focus not only on the global optimum, but rather include all local extrema in the analysis. As such, an exhaustive search approach is conducted to fully survey the solution space of a particular constellation instead of using numerical optimizers to nd Where the joining spacecraft is both the last member of the previous pair and the rst member of the following pair. 66 just the global extremum. This method is, however, computationally intensive; requiring, for example, on the order ofO(10 2 ) navigation simulations for a pair of spacecraft and increasing to roughlyO(10 5 ) for three coupled spacecraft, at a 1 resolution grid spacing. 5.1.1 L1 Halo Case for 2,3, and 4 Satellites A representative periodic orbit is considered, and a baseline trajectory is established by integrating the initial conditions (given in Table 5.1) forward for one full period in the synodic system. To populate the constellation, the host orbit is punctuated with N satellites at varying degrees of separation. The rst spacecraft is placed on the x-axis (i.e. line of symmetry), thereby dening 0 in -space. Every subsequent phase shift is taken in increments of = 10 (implicitly constraining the minimum separation between neighboring spacecraft as well). Sequential loops are then constructed to perform an ex- haustive search over for each combination of spacecraft taken two at a time. Table 5.1: Initial Conditions for the Representative L1 Halo Orbit T = 11:3507 (days) r (km) v (km/s) x 325:2942 10 +03 _ x 37:2198 10 15 y 3:2401 10 09 _ y 270:6137 10 03 z 65:5020 10 +03 _ z 53:0241 10 15 The two spacecraft are modeled using an observation t span of one full orbital period. Perturbations with standard deviations of 0:1 meters and 0:1 millimeters per second (for position and velocity components, respectively) are applied to the initial conditions of the While this increment is arbitrary, a resolution of 10 in phase was found to be sucient to capture the salient features of the navigation performance; interpolation is used to achieve a ner resolution if desired. 67 periodic orbit to emulate an initial guess for the state of the lead spacecraft. To obtain the state for the second spacecraft, the lead spacecraft is integrated under the ow of the dynamics up to a time corresponding to some of interest, and then similarly perturbed with noise. Numerical integration is performed using Matlab's ordinary dierential equa- tions solver (e.g. ode113) with relative and absolute tolerances set to 10 9 and 10 12 , respectively. As the spacecraft move along the orbit, simulated range measurements are taken every 0:0005 time units (187:6 seconds) with a 1 meter random white noise added. For simplicity, occultations are not modeled. A standard weighted batch least-squares lter 39 is then applied to process the series of measurements gathered over the period interval. This statistical orbit determination procedure concludes with the estimate of the states of both spacecraft. Estimation errors obtained during this ltering process are captured in the performance metric given in Eq. (4.1). As part of the exhaustive search, each com- bination of spacecraft pairing represents a dierent dimension for the error metric. The set of spacecraft combinations is evaluated and Fig. 5.1 plots are generated to graphically represent the solution space, thus allowing the visual location of the global optimum (up to N = 4). Simulation results from two, three, and four spacecraft congurations all indicate that LiAISON navigation exhibits a preference towards spacecraft arrangements with large relative separations, and in fact, fails altogether when separations drop below about = 10 . Though globally optimal congurations are found in these cases to resemble a regular N-gon placement of satellites, this convenient structure is just a coincidence. 8 Locally optimal congurations more generally favored scenarios where optimally spaced 68 a) Two spacecraft b) Three spacecraft c) Four spacecraft { yz-slices d) Four spacecraft { zx-slices e) Four spacecraft { xy-slices f) Four spacecraft { tangent slice Figure 5.1 Cost function values,, shown on a log scale to accentuate details, for a sample halo constellation of (a) two, (b) three, and (c-f) four spacecraft. 69 pairs of satellites were used; reinforcing a reasonable notion that having good performers (in a navigation sense) as part of an ensemble will be the leading driver of the overall performance. For each additional spacecraft beyond a pair, the number of points required in this exhaustive search increases exponentially as the power of N1. Furthermore, such a simulation investigates only one orbit at a time (albeit for N spacecraft). Thus, this full search approach presents a great computational diculty when more than one orbit is of interest, and especially so, for more than a pair of spacecraft. This numerical limitation elicits the consideration of alternate means of surveying the solution space for a variety of orbits, with a potentially varying number of spacecraft. 5.2 Design Space Maps and Optimal Congurations The single-orbit spacecraft placement problem described above, shows how increasing the dimension of the problem (i.e. the search space) compounds the numerical dicul- ties encountered. Since a constellation design study will inevitably involve analyzing a multitude of orbits and the associated arrangements of satellites among them, numerical continuation techniques are introduced here to eciently map the design space. Pseudo-arclength continuation involves repeatedly solving Eq. (3.5) using a modied Newton iteration scheme. This, in turn, means that the implicit function given in Eq. (3.6) is also solved numerous times; even more so due to using nite dierence approximations . Though the analytical derivatives for the partials of the cost function given by Eq. (3.6) and Eq. (3.7), can be derived for this performance index, some deciencies were observed while using them in the simulations (see Appendix B for details), and so nite dierences are used instead. 70 As such, each attempt to solve Eq. (3.6) requires exactly ve executions of Eq. (4.2), each of which, is referred to as a LiAISON simulation. 6 At the beginning of each continuation simulation over a specic family, an initial spacecraft placement problem is solved via an exhaustive search (LiAISON simulations at every 1 increment in ) to compute the full navigation prole of the \starting" peri- odic orbit. This step is necessary to obtain the initial guesses for the pseudo-arclength process; i.e. the seed solutions are the local extrema of the rst orbit. Once the seeds are found, separate continuation runs follow each local extrema branch until either the branch stops or bifurcates. In the event of a branch termination, another full search is performed on the periodic orbit corresponding to the stopping point. A new set of local extrema is found and then the continuation process proceeds with evolving the new seeds. Naturally, a manual restart is required to follow bifurcating solutions and similarly with continuing in both directions when starting in the middle of a solution manifold. Thus, the bifurcation diagram is gradually built one branch at a time. Note that since the continuation process could theoretically trace back and forth over one solution branch indenitely, some intuition is required to know when to consider a simulation nished. A constraint on the periodic orbits themselves is imposed by noting the inconsisten- cies in estimation errors due to noise. Though a crucial part in modeling the LiAISON strategy, the measurement noise incorporated into these simulations adversely aects the stability of convergence in the estimator. That is, due to the inherent chaos introduced with statistical aberrations, LiAISON simulations are not exactly deterministic . Experi- ments show this uctuation in performance to be the most prominent at the periphery of In fact, the corresponding system of dierential equations models an intrinsically stochastic process. 71 a particular family. For this reason, the continuum of a periodic orbit family is truncated to investigate only those orbits that lie between 0:2 0:9, where the performance is more well-behaved. This phenomena is consistent with observations made by Hill and Born, whereby LiAISON performance deterioration is seen to fall proportionally to congurations that progressively become more planar. 17 The initial conditions for the periodic orbits themselves are extracted from solutions generated with AUTO and then subsequently dierentially corrected to obtain true pe- riodic orbits, as discussed in Chapter 4. Since AUTO uses a variable step size in its numerical continuation routine, the number of computed orbits that comprise a family is not necessarily the same. In order to consistently span the data, the boundaries of a particular family are mapped onto the closed unit interval [0; 1]. This mapping allows an interpolation within the interval such that each orbit can be uniquely identied by a scalar between 0 and 1. This scalar is, indeed, the local in-family parameter 2R. At the conclusion of a pseudo-arclength simulation, the results can be summarized and succinctly represented in bifurcation diagrams, shown in Figs. 5.2a, 5.3a, 5.4a, and 5.5a, for the case studies of the L1 halo, L1 axial, L2 vertical, and the L5 vertical periodic orbit families, respectively. These curves, in particular, are overlaid on top of surface maps corresponding to whole design spaces that are generated via exhaustive searches, with each point corresponding to a full navigation simulation (in practice, these maps are unnecessary but are used here for context). Branches of local extrema can be seen to turn and fold, develop in isolated regions, and roam unhindered across the whole solution space. Some manifolds are even seen to connect branches of local minima to branches of local maxima. 72 For validation purposes, the coarse grid is scanned for crude outlines of local extrema, with which the continuation-generated branches are compared against as rough bench- marks. These \skeleton paths" are shown with square markers in the design maps of Figs. 5.2a { 5.5a, and the continuation-generated branches are represented with circular markers. Although some idea of where local extrema exist can be guided from analyz- ing just these skeleton paths, understanding the ner contours and the interconnections between branches is missing without a more detailed search. That is, the accuracy of- fered by exhaustive searches is limited by the resolution of the grid. On the other hand, pseudo-arclength continuation is able to recover the optimal congurations alluded to by the brute force approach, but at a much ner quality than attainable from a crude grid search. In addition to being able to identify the global optimum from continuing local ex- trema, the insight gained from analyzing the continuation maps aords a more unied understanding of the overall structure of the solution space. 7 By surveying the bifurca- tion diagrams obtained via the continuation process, the global optimum of each family are identied using one-dimensional line searches over each branch (shown in the de- sign maps with a star marker). Three-dimensional projections of the optimal orbits are shown in parts (b) and (c) in Figs. 5.2 { 5.5. Also depicted in the gures, are the opti- mal placements of the pair of spacecraft, at one instant in time. To recapitulate, these globally optimal congurations represent the best performing LiAISON constellation ar- rangements from the standpoint of constellation-spacecraft navigation (as opposed to the end-user). The initial conditions corresponding to the down-selected, optimal periodic orbits, are provided in Table 5.2. 73 5.2.1 L1 Halo Family a) Bifurcation diagram b) Optimum along manifold c) Optimum with projections Figure 5.2 (a) Continuation-generated bifurcation diagram of the design space for the L1 halo family, overlaid on top of a background surface map compiled from an exhaustive grid-search. The \optimal" orbit is highlighted by the star marker, whose projection in three-dimensions is plotted separately (b) in context along the manifold and (c) with orthogonal projections and spacecraft at one instance in time. 74 5.2.2 L1 Axial Family a) Bifurcation diagram b) Optimum along manifold c) Optimum with projections Figure 5.3 (a) Continuation-generated bifurcation diagram of the design space for the L1 axial family, overlaid on top of a background surface map compiled from an exhaustive grid-search. The \optimal" orbit is highlighted by the star marker, whose projection in three-dimensions is plotted separately (b) in context along the manifold and (c) with orthogonal projections and spacecraft at one instance in time. 75 5.2.3 L2 Vertical Family a) Bifurcation diagram b) Optimum along manifold c) Optimum with projections Figure 5.4 (a) Continuation-generated bifurcation diagram of the design space for the L2 vertical family, overlaid on top of a background surface map compiled from an exhaustive grid-search. The \optimal" orbit is highlighted by the star marker, whose projection in three-dimensions is plotted separately (b) in context along the manifold and (c) with orthogonal projections and spacecraft at one instance in time. 76 5.2.4 L5 Vertical Family a) Bifurcation diagram b) Optimum along manifold c) Optimum with projections Figure 5.5 (a) Continuation-generated bifurcation diagram of the design space for the L5 vertical family, overlaid on top of a background surface map compiled from an exhaustive grid-search. The \optimal" orbit is highlighted by the star marker, whose projection in three-dimensions is plotted separately (b) in context along the manifold and (c) with orthogonal projections and spacecraft at one instance in time. 77 5.2.5 Ephemeris Simulations The previous results represent a prototypical study of autonomous constellations in the Earth-Moon CR3BP, using the general constellation design framework addressed in this dissertation. While interesting in their own right, these results using the CR3BP may not necessarily be the end goal, but rather serve as initial guesses for higher delity simulations. To demonstrate the next step, the four optimal orbits identied in this study are mapped from the CR3BP frame (see Table 5.2) to the Earth Mean Equator and Equinox of Reference Date coordinate system and reference frame (see Table 5.3), and then integrated using Earth-Moon-Sun ephemerides. For the purposes of investigating the Earth-Moon system, the Sun's presence is added as a third-body perturbation. The NASA GMAT software package is used to simulate the resulting optimal trajectories using the JPL DE405 ephemeris library (see Appendix D for conguration settings). Table 5.2: Initial Conditions for Optimal Constellations (Synodic Frame) Parameters L1 Halo L1 Axial L2 Vertical L5 Vertical ^ (deg) 198.74 52.644 145.26 217.45 ^ (none) 0.64762 0.37204 0.55096 0.73051 ln() (none) -17.4426 -18.9731 -19.4822 -15.7267 ave (m) 3.1716 0.7279 0.5551 25.1978 T (days) 10.5798 17.4809 22.5624 27.3949 x (km) 3.34643e5 3.61135e5 3.82028e5 2.67651e5 y (km) -2.83358e4 3.48972e4 3.26388e4 3.76588e4 z (km) 6.21431e4 -6.16999e4 -9.94797e4 -2.72554e5 _ x (km/s) -8.56135e-2 -2.64092e-2 2.05111e-1 7.43409e-1 _ y (km/s) 2.18020e-1 -1.60989e-1 2.94640e-2 -1.26624e0 _ z (km/s) 1.38248e-1 1.76318e-1 2.77457e-1 5.64709e-1 78 Table 5.3: Initial Conditions for Optimal Orbits (Sidereal Frame) Epoch: 22 Feb 2011 00:00:00.0000 (UTC) Parameters L1 Halo L1 Axial L2 Vertical L5 Vertical T (days) 10.5798 17.4809 22.5624 27.3949 x (km) -3.25755e5 -3.13315e5 -3.29904e5 -2.08494e5 y (km) -1.16139e5 -1.29928e5 -1.20135e5 -1.67694e4 z (km) -1.39141e4 -1.54244e5 -1.92945e5 -3.25748e5 _ x (km/s) 4.51531e-1 4.22765e-1 2.97343e-1 -8.53943e-1 _ y (km/s) -1.00806e0 -7.23048e-1 -1.05391e0 1.74235e-2 _ z (km/s) -2.42756e-1 -8.01115e-2 -1.25215e-1 5.43062e-1 Figure 5.6 shows the initial CR3BP trajectories (red) superimposed over the ephemeris simulations (black), for each of the four cases separately. As is immediately obvious, the addition of perturbations to the approximated dynamics of the CR3BP quickly distorts the previously well-behaved orbits into paths only vaguely resembling periodic orbits. Traditionally, in order to fully extend the capabilities of a trajectory design technique into the ephemeris model, a dierential correction routine must be employed to nd the corresponding periodic orbits in that dynamical system. 26 But since the proposed con- stellation design methodology has already been demonstrated using a textbook dynamical system, no further attempt is made in this work to tailor the framework to another dy- namical system. The ephemeris simulations shown here, are for context only and are used to reveal some limitations of the CR3BP. 79 a) Optimal L1 halo orbit b) Optimal L1 axial orbit c) Optimal L2 vertical orbit d) Optimal L5 vertical orbit Figure 5.6 Optimal orbits corresponding to (a) L1 halo, (b) L1 axial, (c) L2 ver- tical, and (d) L5 vertical cases, with red/black lines indicating CR3BP/ephemeris solutions, respectively, shown in the non-dimensional synodic frame. 80 Chapter 6 Performance Characterization Following the successful demonstration of pseudo-arclength in nding optimal constel- lations, certain important characteristics in its performance are noted. In the presence of techniques such as exhaustive searches or numerical optimizers, comparisons are made to draw out some advantages of numerical continuation. This chapter highlights both notable merits and demerits of using pseudo-arclength continuation to solve the proposed constellation design problem. 6.1 Computational Eciency Although it may seem attractive to simply generate a full surface map (via the lengthy \brute force" method) to survey the solution space, the limitation is that this surface is only visualizable forN parameters, ifN 4. Conversely, as a general technique numerical continuation is well suited to eciently extend solutions over an arbitrary N dimensional design space, using little a priori intuition about the overarching topology. That is, from the standpoint of parametrizing a solution point with the scalar pseudo-arclength, z(s), 81 each solution branch maintains the topology of a one-dimensional curve irrespective of the number of dimensions in the design space. So even when the design parameters exceed the comfortable limit of visualization, pseudo-arclength can still follow the solution branch. Numerical continuation thus provides a concise framework for following the behavior of solution points embedded in any N-dimensional space, without having to compute the entire space. In an optimization setting, the solution branches correspond to the null level set. Though this is the natural goal of an optimization problem, in principle, continuation processes can recover any level set; i.e. any contourC2R which satises the root-nding problem (u;)C = 0. One advantage of being able to select specic contours to evolve, is in the ability to investigate regions nearby optimal curves. This discrimination is particularly benecial from the perspectives of providing initial guesses for convex basins of attraction, as well as revealing how the local neighborhood around optimal solutions deforms. Each level set generation is comparable in speed to recovering the optimal branches. Comparing the run-times for the exhaustive search routines and the continuation schemes, shows a 2{4 times improvement in computational speed favoring continuation. This improvement, however, is observed here for a low-dimensional continuation problem (N = 2), as compared to a coarse resolution grid search y . This trend is expected to increase exponentially as the dimensionality of the problem increases and geometrically An N-tuple is viewed as a coordinate, or \point", in N-space, whose collection is thus a one- dimensional curve. 36 y The grid is divided by 10 in and 0:02 dimensionless units in , yielding 1260 points. 82 with the resolution of the grid. Thus, continuation methods become more ecient at solving \larger" problems. 6.2 Continuation of Stochastic Systems Another merit of continuation oddly comes from observing its diculty with recover- ing certain solution branches. Given a suciently small step size, the modied Newton method is guaranteed to converge, however, not always within the same number of itera- tions or at the same resolution 12 (in fact, some contours may be missed if the step size is too large). In the case of using pseudo-arclength continuation in a spacecraft navigation setting, a non-deterministic implicit function dened by Eq. (3.6), incorporates statisti- cal noise that perturbs the stability of the state estimation to varying degrees, causing Newton's method some level of diculty. Congurations of LiAISON constellations that are robust to noise tend to produce consistent estimation errors and hence are easily con- verged upon using pseudo-arclength. The opposite is true with congurations that are sensitive to noise, taking considerably longer to converge and often at tenfold smaller step sizes than their counterparts. This eect is demonstrated in Fig. 5.5a, where the den- sity of points near the fold in the branch is far greater than elsewhere along the branch; corresponding to a strenuous simulation likely caused by an unstable branch of local ex- trema. Therefore, in some sense, the number of iterations required plus the density of solution points can be viewed as indicators of the sensitivity of particular congurations towards perturbations. This insight is notably absent when performing the \brute force" exhaustive searches. 83 This surprising ability to handle stochastic systems of nonlinear dierential equations seems to stem from the numerics in the Newton iteration process. Since Newton's method is based on recursively rening an initial guess, a radius of convergence must be dened to alert the program to stop; this convergence criterion is numerically set as a tolerance in the conditional statement of the iteration loop. Adjusting this \dial" thus allows the routine to accommodate dierent level sets of noise . The trade-o comes, however, at a detriment to the accuracy of the solution obtained during the root-nding process. To achieve a generally well-behaved simulation of stochastic systems, as a rule the conver- gence tolerance for Newton's method should not be set lower than the standard deviation of noise. 6.3 Isolated Islands of Local Extrema Notwithstanding its eciency in extending solutions overN-space, and its robustness to noise, continuation methods in general have a minor demerit. In the case of studying extrema (always occurring in pairs) there exists a situation where a previously unknown min/max pair could develop between a set of known local extrema. Without sucient a priori information, numerical continuation schemes have no guarantee of recovering these isolated islands of local extrema. As an example, these islands appear towards the center portion of Fig. 5.5a undetected by the standard continuation routine. In these scenarios, extra logic must be applied to detect potential gaps in coverage; such as using dierent contour levels as mentioned above. Up to some limit, yet without claiming that this process can absorb highly chaotic systems. 84 Chapter 7 Conclusion In this dissertation, a generalized constellation design methodology is proposed. Addi- tionally, a particular implementation of this theory is shown with examples of designing autonomous interplanetary constellations in the restricted Earth-Moon system. This chapter concludes this dissertation by summarizing the results of this research and draw- ing insights related to its application in space ight. Lastly, some thoughts on possible directions for future study are oered. 7.1 Summary of Contributions As seen in previous chapters, the basic notion of a constellation consists of periodic host orbits punctuated by N spacecraft. The art of designing such a constellation thus seeks to nd appropriate combinations of such orbits and ensembles of satellites placed among them. But since the orbits are local structures of a particular set of dynamics, and what is considered an optimal placement of satellites varies per mission objective, a general design strategy must incorporate some exibility to accept these variations. The main contribution of this work is indeed in developing a tailorable constellation 85 design methodology that can, in principle, be adapted to handle arbitrary dynamical systems along with a variety of continuously dierentiable scalar cost functions that dene optimality. This framework diers from conventional constellation design techniques in that it is extensible beyond the regime of standard 2-body Keplerian mechanics (i.e. beyond conic sections). Formally, the constellation design process is approached as an optimization problem. The construction of a constellation is viewed as a sequence of decoupled sub-problems involving a trajectory design step followed by a spacecraft placement step; i.e. forming a nested optimization problem. This two-step process aims to nd the \best" (dened by the cost function) arrangement of spacecraft along the most eective periodic orbits. A noteworthy feature is the parametrization scheme that allows the solution of the nested optimization problem to be found via numerical continuation techniques (e.g. pseudo-arclength continuation). While the specic variables used for the parametrization are not unique, in general, it is shown that an N-dimensional problem can be reduced to a one-dimensional problem by ensuring the problem contains a transformation that maps R N toR N1 , and that the Jacobian of such a mapping is full rank. In this fashion, the inner optimization problem yields a one-dimensional solution manifold (by the implicit function theorem), over which the outer optimization problem can further rene via a simple line search. By using the pseudo-arclength continuation strategy to map the solution space of the optimization problem, a succinct characterization of the underlying structure is achievable with little a priori intuition of the map's overarching topology. Furthermore, continuation techniques reduce computational overhead as compared to more traditional sequential 86 exhaustive search methods. While an achievement in its own right, providing a dramatic savings in terms of processing speed also opens the door to analyzing a much larger class of problems that are otherwise intractable using other search techniques (e.g. global/local optimizers and evolutionary algorithms). 7.2 Insights Gained Using the optimization framework developed in this dissertation, a proof of concept study is conducted to design autonomous Earth-Moon constellations. Here, the dynamical system is chosen to be the CR3BP and the cost function that measures optimality is set as the navigation performance index obtained from the LiAISON technique. The results of the simulations revealed optimal interplanetary constellations in four dierent periodic orbit families. Moreover, besides demonstrating the proposed constella- tion design methodology, this work serves as verication of the eectiveness of the LiAI- SON technique as applied to invariant manifolds of periodic orbits other than the initially investigated halo families (e.g. extending to the vertical and axial families). Navigation performance errors, obtained from the position covariance sub-matrices, yield approxi- mately 0:6; 0:7; 3; and 25 meters, for the L2 vertical, L1 axial, L1 halo, and L5 vertical optimal congurations, respectively; that is, the batch lter converges to the order of 10 1 times the measurement noise (here, 1 meter standard deviation). The congurations that produce these estimation accuracies seem to indicate that the LiAISON strategy exhibits a preference for arrangements that have large angular separations between the spacecraft, and in fact fails altogether if the relative spacing falls below roughly 10 degrees in phase. 87 There also seems to be some correlation between LiAISON accuracy and the dynamic stability of the periodic orbits, however, this connection is not well understood. Although illustrated using only the base set of a pair of spacecraft on the same orbit, the design approach is well suited to accommodate any number of spacecraft per constellation. While this research is indeed focused on \constellation" design, there is actually no inherent restriction of the developed framework to work only with closed orbits. In fact, since preliminary studies did not incorporate a dierential correction routine, the man- ifolds originally studied were actually unconnected. This ability to accept non-compact manifolds insinuates a wider eld of application, reaching beyond conventional constella- tion design to more generalized trajectory continuums. Surprisingly, the problem of interest also need not be fully deterministic. The suc- cess of the continuation strategy as applied to a statistical orbit determination setting, uncovered the unexpected ability of the pseudo-arclength technique to handle stochastic systems of nonlinear dierential equations. Thus, the process of numerically continuing solution branches of an intrinsically non-deterministic problem could potentially act as an indicator of a system's sensitivity to noise. Despite its diculty in recovering remote islands of local extrema, numerical continuation oers a natural means to eciently ex- plore the interconnectedness of solutions over arbitrarily large design spaces, making it a robust option for trade/survey studies. 88 7.3 Recommendations for Future Study The interplanetary constellations studied throughout this research considered only basic periodic orbits. Clearly, a much richer orbital environment is available, and as such, a natural extension of this work would be to apply the constellation design methodology to other classes of host orbits: several interesting candidates include resonant periodic orbits (e.g. multi-revolution cyclers 34 orp:q resonances 2 ), KAM tori, 41 homoclinic/heteroclinic connections, 15 and nearly-symmetric Keplerian orbits. Presenting perhaps the most immediate in uence in modern constellation design would be the extension of LiAISON to 2-body orbits. This addendum to the exist- ing arsenal of autonomous navigation strategies used in nearly-symmetric gravity elds, could potentially lead to a novel cost-eective mechanism for achieving spacecraft au- tonomy, by employing only relative intersatellite range measurements bolstered by infre- quent updates. Termed \LiAISON bootstrapping", this technique could improve orbit determination functionality onboard spacecraft, as well as oer a means to alleviate the dependence on ground-based resources. Since the application of this work focuses on single-orbit constellations (and by ex- tension, over single orbit manifolds), it would be interesting to further investigate how to adequately parametrize the union of multiple orbit manifolds, so that a \fuller" con- stellation design methodology could be developed. One thought would be to interleave other levels of nested optimization problems into the current framework to overcome the a priori construction of admissible orbit spaces. Quasi-periodic structures attributed to small nonlinear perturbations of a dynamical system as pre- dicted by the Kolmogorov-Arnold-Moser (KAM) theory. 89 As noted earlier, limitations of this methodology in using closed orbits can be removed by considering intentionally unbounded trajectories. This modication leads to the ex- tension of this framework to the design of other trajectory manifolds (e.g. stable/unstable invariant manifolds associated with the arrival/departure to and from libration point or- bits 42 ). Although this path of research represents a deviation from constellation design, it is still a worthwhile avenue of exploration. Several enhancements that are desirable, of course, involve processing improvements through software, and by extension, to computer hardware. Most notably absent from the current algorithms is the usage of parallelism. Due to the independent grid-like nature of survey problems, large design studies would benet from porting the computations from multi-core CPU's over to the massively parallel architectures aorded by GPU's y . Enabling parallel processing further broadens the scope of surmountable problems in astrodynamics, as well as ameliorates the delity of current problems with unprecedented levels of detail. 10 Together with eorts to ensure automation (e.g. automatic bifurcation- detection and branch-switching) in the numerical continuation routines, these renements to the existing design framework could culminate in a software tool that would be useful to satellite network designers and mission planners, as an ecient means to prototype nominal designs for beginning-phase research. Central Processing Units. y Graphical Processing Units. 90 References [1] R. Anderson. Low Thrust Trajectory Design for Resonant Flybys and Captures using Invariant Manifolds. PhD thesis, Univesity of Colorado, Boulder, Boulder, CO, 2005. [2] R. Anderson and M. Lo. A Dynamical Systems Analysis of Resonant Flybys: Bal- listic Case. Journal of Astronautical Sciences, 58(2):167{194, 2011. [3] V.I. Arnold. Ordinary Dierential Equations. The MIT Press, 1973. [4] R. Bate, D. Mueller, and J. White. Fundamentals of Astrodynamics. Dover, 1971. [5] J. Casoliva, J.M. Mondelo, B.F. Villac, K.D. Mease, E. Barrabes, and M. Olle. Two Classes of Cycler Trajectories in the Earth-Moon System. Journal of Guidance, Control, and Dynamics, 33(5):1623{1640, 2010. [6] C. Chow and B. Villac. Numerical Continuation of Optimal Spacecraft Place- ments for LiAISON Constellations. AAS/AIAA Astrodynamics Specialist Confer- ence, Girdwood, AK, AAS:11{542, Jul 31 { Aug 04, 2011. [7] C. Chow and B. Villac. Mapping Autonomous Constellation Design Spaces Using Numerical Continuation. Journal of Guidance, Control, and Dynamics (accepted), 2012. [8] C. Chow, B. Villac, and M. Lo. Optimizing Spacecraft Placement for LiAISON Con- stellations. AAS/AIAA Space ight Mechanics Meeting, New Orleans, LA, AAS:11{ 229, Feb 13{17, 2011. [9] N. Chow, E. Gralla, J. Chase, and N.J. Kasdin. Low Earth Orbit Constellation Design using the Earth-Moon L1 Point. AAS/AIAA Space ight Mechanics Meeting, Maui, HI, AAS:04{248, Feb 8{12, 2004. [10] A.C. Colombi. Quick Evaluation of Small Body Gravitation. PhD thesis, University of Illinois at Urbana-Champaign, Urbana-Champaign, IL, 2008. [11] J. Danby. Fundamentals of Celestial Mechanics. Willmann{Bell, Inc., 2nd edition, 2003. [12] E. Doedel. Numerical Analysis and Bifurcation Problems. Lecture Notes, Montreal, Canada, May 2000. 91 [13] E. Doedel, V. Romanov, R. Paenroth, H. Keller, D. Dichmann, J. Galan-Vioque, and A. Vanderbauwhede. Elemental Periodic Orbits associated with the Libration Points in the Circular Restricted 3-Body Problem. International Journal of Bifur- cation and Chaos, 17(8):2625{2677, 2007. [14] R. Farquhar. The Control and Use of Libration Point Satellites. PhD thesis, Stan- ford University, Stanford, CA, 1968. [15] G. G omez, W.S. Koon, M.W. Lo, J.E. Marsden, J. Masdemont, and S.D. Ross. Connecting Orbits and Invariant Manifolds in the Spatial Restricted Three-Body Problem. Nonlinearity, 17(5):1571{1606, 2004. [16] M. H enon. Generating Families in the Restricted Three-Body Problem. Lecture Notes in Physics: Monographs, volume 1. Springer{Verlag, 1997. [17] K. Hill and G.H. Born. Autonomous Interplanetary Orbit Determination Us- ing Satellite-to-Satellite Tracking. Journal of Guidance, Control, and Dynamics, 30(3):679{686, 2007. [18] K. Hill, M.W. Lo, and G.H. Born. Linked, Autonomous, Interplanetary Satellite Orbit Navigation (LiAISON). Advances in the Astronautical Sciences, 123, Part III:2353{2367, 2006. [19] M. Hirsch, S. Smale, and R. Devaney. Dierential Equations, Dynamical Systems & An Introduction to Chaos. Elsevier, 2nd edition, 2004. [20] B. Hofmann-Wellenhof, H. Lichtenegger, and J. Collins. GPS Theory and Practice. Springer{Verlag, 2nd edition, 1992. [21] K.C. Howell. Three-Dimensional, Periodic, Halo Orbits. Celestial Mechanics, 32:53{ 71, 1984. [22] H. Keller. Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems. Applications of Bifurcation Theory, pages 359{384, 1977. [23] W. Koon, M. Lo, J. Marsden, and S. Ross. Dynamical Systems and the Three- Body Problem and Space Mission Design. Marsden Books, ISBN:978{0{615{24095{ 4, 2008. [24] M. Lo. Satellite-Constellation Design. Computing in Science and Engineering, pages 58{67, Jan{Feb 1999. [25] B. Marchand and K. Howell. Control Strategies for Formation Flight in the Vicinity of the Libration Points. Journal of Guidance, Control, and Dynamics, 28(6):1210{ 1219, 2005. [26] B. Marchand, K. Howell, and R. Wilson. An Improved Corrections Process for Constrained Trajectory Design in the n-Body Problem. Journal of Spacecraft and Rockets, 44(4):884{897, 2007. 92 [27] F.L. Markley and J.R. Carpenter. Generalized Linear Covariance Analysis. Journal of Astronautical Sciences, 57(1{2):233{260, 2009. [28] O. Montenbruck and E. Gill. Satellite Orbits: Models, Methods, Applications. Springer, 1st edition, 2005. [29] D. Mortari, M. Wilkins, and C. Bruccoleri. The Flower Constellations. John L. Junkins Astrodynamics Symposium, Texas A & M University, College Station, TX, AAS:03{274, May 23{24, 2003. [30] NASA. Communications and Navigation Systems Technology Area 05: TA05. NASA's DRAFT Space Technology Roadmaps, Nov 2010. [31] U.S. Government Accountability Oce. NASA's Deep Space Network: Current Management Structure Is Not Conducive to Eectively Matching Resources with Future Requirements. Apr 27, 2006. GAO{06{445. [32] U.S. Government Accountability Oce. Global Positioning System: Signicant Challenges in Sustaining and Upgrading Widely Used Capabilities. May 7, 2009. GAO{09{670T. [33] V. Romanov. Elemental Periodic Solutions of the Circular Restricted 3-Body Prob- lem. Master's thesis, Concordia Univesity, Montr eal, Queb ec, Canada, 2005. [34] R. Russell and N. Strange. Planetary Moon Cycler Trajectories. Journal of Guid- ance, Control, and Dynamics, 32(1):143{157, 2009. [35] H. Schaub, S.R. Vadali, J.L. Junkins, and K.T. Alfriend. Spacecraft Formation Flying Control Using Mean Orbital Elements. Journal of Astronautical Sciences, 48(1):69{87, 2000. [36] M. Spivak. Calculus on Manifolds. Westview Press, 1998. [37] M. Standish. JPL Planetary and Lunar Ephemerides, DE405/LE405. JPL IOM 312.F{98{048, Aug 1998. [38] V. Szebehely. Theory of Orbits. Academic Press, Inc., 1967. [39] B.D. Tapley, B.E. Schutz, and G.H. Born. Statistical Orbit Determination. Elsevier Inc., 2004. [40] B.F. Villac, C.C. Chow, M.W. Lo, and G.R. Hintz. A Dynamical System Approach to Orbit Down-Selection of Earth-Moon Autonomous Navigation Constellations. Jour- nal of Astronautical Sciences (accepted), 2011. [41] W. Wiesel. Earth Satellite Perturbation Theories as Approximate KAM Tori. Jour- nal of Astronautical Sciences, 58(2):153{165, 2011. [42] S. Wiggins. Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer{Verlag, 1990. 93 Appendix A Coordinate Transformations Often in trajectory design it is useful to view orbits from dierent coordinate systems. One benet is to highlight key features and structures that are apparent in one frame, that may otherwise be absent in another. As such, transformations between coordinate systems can be constructed as continuous mappings between basis vectors to allow smooth transitions from one frame to another. A distinction should be made clear that this section will describe the change of basis mappings and not the rotations of vectors within the same basis. Consider a linear map described by : R n ! R n (e.g. n = 3 for the standard three-dimensional Euclidean space). The relationship between the position, velocity, and acceleration vectors may be expressed as, r 0 = r _ r 0 = _ r + _ r r 0 = r + 2 _ _ r + r (A.1) 94 where the ( 0 ) indicates a dierent basis and the ( _ ) represents the total time derivative. The position,r, velocity, _ r, and acceleration, r, vectors can be expressed with standard Cartesian coordinates in their respective bases as, r = [x; y; z] T ; _ r = [ _ x; _ y; _ z] T ; r = [ x; y; z] T Conveniently, a matrix representation of Eq. (A.1) allows the formulation of one block matrix equation with which to transform a state vectorX = [r; _ r; r] T expressed in one basis into another. X 0 = 2 6 6 6 6 6 6 4 0 33 0 33 _ 0 33 2 _ 3 7 7 7 7 7 7 5 X (A.2) a) Synodic frame b) Sidereal frame Figure A.1 A sample Earth-Moon L2 vertical periodic orbit plotted in black for both the (a) synodic and (b) sidereal frames. Orthogonal projections are shown in light gray. 95 For example, a typical transformation of interest is going from a rotating frame to an inertial frame. 1 Figure A.1 illustrates an example of the dierent features uncovered from viewing an orbit in alternate reference frames. The corresponding mapping, , is the rotation matrix about the common angular momentum axis of the two frames. This mapping and its subsequent time derivatives are given as, = 2 6 6 6 6 6 6 4 cos!t sin!t 0 sin!t cos!t 0 0 0 1 3 7 7 7 7 7 7 5 _ = 2 6 6 6 6 6 6 4 sin!t cos!t 0 cos!t sin!t 0 0 0 0 3 7 7 7 7 7 7 5 ! = 2 6 6 6 6 6 6 4 cos!t sin!t 0 sin!t cos!t 0 0 0 0 3 7 7 7 7 7 7 5 ! 2 + 2 6 6 6 6 6 6 4 sin!t cos!t 0 cos!t sin!t 0 0 0 0 3 7 7 7 7 7 7 5 _ ! Sometimes it is desired to specify the state of the spacecraft with respect to a dierent origin. In this case the new state is obtained by adding the state vectors of the origins of intermediate systems. This origin relocation scheme can be summarized as, X s=q =X s + q X p=1 X p (A.3) where the subscripts s, p, and q denote spacecraft, intermediate origins, and nal origin, respectively. For example, to transform a spacecraft state given in the Earth-Moon 96 barycentered synodic frame, to a geocentric inertial coordinate system, one would apply Eq. (A.2) to map the state from the synodic to the sidereal frame, and then use Eq. (A.3), withX q = [ 0 0 0 0] T , to shift the origin of the system to the center of the Earth. If it is further desired to align the state to a particular epoch (e.g. in the ephemeris model), a rotation is necessary to match the Earth-Moon line in the geocentric frame to the actual position of the Moon at a given epoch. This rotation matrix can be conveniently constructed from the Moon's position vector, the angular momentum vector of the Moon's instantaneous orbital plane, and their cross product. Dening each vector and assembling the terms gives, ^ h $ = r $ _ r $ = r $ _ r $ ; ^ r $ =r $ = r $ ; ^ c $ = ^ h $ ^ r $ $ = h ^ r $ j ^ c $ j ^ h $ i With the epoch xingr $ and _ r $ , the associated nal state vector can be described by, X epoch = 2 6 6 4 $ 0 33 0 33 $ 3 7 7 5 X s=q The linear matrix operator in Eq. (A.2) can be truncated to its leading principal minor if only the position and velocity are considered as part of the state. 97 Appendix B Performance Index Partials During the formulation of pseudo-arclength continuation, nite dierence approx- imations of the implicit function are implemented rather than the explicit analytical derivatives themselves. This choice is based on experimental observations that suggest a slight deciency in these particular derivatives in fully capturing the appropriate trends. Nonetheless, the derivation of this performance index is described here for completeness. The parametrized nonlinear functions of interest are given by the Jacobian of Eq. (3.4) with the cost function set by Eq. (4.2), @(;) @ = 0 where 2 R N1 and ;2 R. Note that the performance index, , is dierent than the one posed by Hill and Born, in that the square-root is omitted from ave and the trace of the full covariance matrix is used rather than eigenvalues obtained from only the position covariance sub-matrix. These modications do not alter the trend, however, they do make computing the derivatives considerably easier. Here, note the distinction of 98 using only pairs of spacecraft. Thus, N spacecraft can be organized as N separate pairs (except when N = 2, when the pairs are degenerate). The following derivation will rely on evaluating separate components of the derivative and then piecing them back together at the end. The expansion of the implicit function is, @(;) @ = 1 mN N X j " m X i Trace @P i @ # j since the trace is a linear operator, it commutes with the derivative. Focusing on the rst component, @P i @ = @ @ h (t i ;t 0 )P 0 T (t i ;t 0 ) i = @(t i ;t 0 ) @ P 0 T (t i ;t 0 ) + (t i ;t 0 ) @P 0 @ T (t i ;t 0 ) + (t i ;t 0 )P 0 @ T (t i ;t 0 ) @ the quantities of interest are the partials of the state transition matrix, (t i ;t 0 ), and the partials of the initial covariance matrix, P 0 . Addressing the latter rst, @P 0 @ = @ 1 0 @ = 1 0 @ 0 @ 1 0 where the partials of the information matrix, 0 , obtained from the least-squares process becomes, @ 0 @ = @ @ h H T WH +P 1 0 i = @H T @ WH + H T W @H @ The two terms that vanish while applying the chain rule are the partials of W and P 1 0 because they are constants for this problem. TheH matrix accumulates during the batch 99 lter in the form of row \stacking", whereby each subsequent row is appended to the end of the matrix thus increasing its row dimension to match the number of observations. Each row, in turn, corresponds to a mapping of a measurement's eect on the initial state deviations, and conveniently, may be treated independently such that, @H i @ = @ @ h ~ H i (t i ;t 0 ) i = @ ~ H i @ (t i ;t 0 ) + ~ H i @(t i ;t 0 ) @ Here, again, the partials of the state transition matrix appear. But for the moment, its analysis will be postponed in order to concentrate on the partials of the observation-state mapping matrix, ~ H, @ ~ H i @ = @ ~ H i @X i @X i @ Since the observation-mapping matrix is the Jacobian of the measurement model, the partial above is then the product of the Hessian of the measurement model and the vector eld at timet i . Note that the change in phase angles,, only aect the placement of spacecraft other than the \rst" one. Thus, the vector eld of the state becomes @X i =@ = [0 16 ; @X () i =@] T , for = 1;:::;N 1, (assuming that the full state is organized to have the rst spacecraft as the rst entry block). Attention is now devoted to the Hessian of the measurement model. For context the measurement model for LiAISON (i.e. Euclidean norm) and its Jacobian are also given. Measurement Model: G(X i ) = p (x 1 x 2 ) 2 + (y 1 y 2 ) 2 + (z 1 z 2 ) 2 100 Jacobian, ~ H i : @G(X i ) @X i = " (x 1 x 2 ) G ; (y 1 y 2 ) G ; (z 1 z 2 ) G ; 0 13 ; ::: (x 1 x 2 ) G ; (y 1 y 2 ) G ; (z 1 z 2 ) G ; 0 13 # Hessian, @ ~ H i =@X i : @ 2 G(X i ) @X 2 i = " @ ~ H i (1) @X i @ ~ H i (2) @X i @ ~ H i (3) @X i 0 123 @ ~ H i (7) @X i @ ~ H i (8) @X i @ ~ H i (9) @X i 0 123 # where each non-trivial column is given as, @ ~ H i (1) @X i = " 1 G (x 1 x 2 ) 2 G 3 ; (x 1 x 2 )(y 1 y 2 ) G 3 ; (x 1 x 2 )(z 1 z 2 ) G 3 ; 0 13 ; ::: 1 G + (x 1 x 2 ) 2 G 3 ; (x 1 x 2 )(y 1 y 2 ) G 3 ; (x 1 x 2 )(z 1 z 2 ) G 3 ; 0 13 # T @ ~ H i (2) @X i = " (y 1 y 2 )(x 1 x 2 ) G 3 ; 1 G (y 1 y 2 ) 2 G 3 ; (y 1 y 2 )(z 1 z 2 ) G 3 ; 0 13 ; ::: (y 1 y 2 )(x 1 x 2 ) G 3 ; 1 G + (y 1 y 2 ) 2 G 3 ; (y 1 y 2 )(z 1 z 2 ) G 3 ; 0 13 # T @ ~ H i (3) @X i = " (z 1 z 2 )(x 1 x 2 ) G 3 ; (z 1 z 2 )(y 1 y 2 ) G 3 ; 1 G (z 1 z 2 ) 2 G 3 ; 0 13 ; ::: (z 1 z 2 )(x 1 x 2 ) G 3 ; (z 1 z 2 )(y 1 y 2 ) G 3 ; 1 G + (z 1 z 2 ) 2 G 3 ; 0 13 # T @ ~ H i (7) @X i = @ ~ H i (1) @X i ; @ ~ H i (8) @X i = @ ~ H i (2) @X i ; @ ~ H i (9) @X i = @ ~ H i (3) @X i 101 This formulation for the Hessian of the measurement model is, indeed, square and symmetric. Since the Jacobian of the measurement model is already being computed during the batch process at each time t i , it is an easy modication to include this com- putation alongside it. At this stage the only quantity remaining is the partials of the state transition matrix. For this mapping again apply the reasoning that a shift in will not alter the state of the rst spacecraft. Additionally, since the full state transition matrix for N-spacecraft is block diagonal, it is sucient to derive the formulation for a single block; i.e. a 6 6 matrix. The following analysis starts with taking the partials of the variational equations for the state transition matrix. @ @ d dt (t;t 0 ) =A(t) (t;t 0 ) ! ) d dt @(t;t 0 ) @ = @A(t) @ (t;t 0 ) +A(t) @(t;t 0 ) @ Here, a set of ordinary dierential equations is formulated that, in turn, contains a quan- tity which represents a partial derivative. That is, this form is developed for use with numerical integrators such that _ = F (;t), where the solution, , yields the partials of the state transition matrix. However, the last analytical quantity that needs to be further developed is the partials of the Jacobian of the vector eld. @A(t) @ = @ @ @F (X(t)) @X(t) = @ 2 F (X(t)) @X 2 (t) @X(t) @ This quantity is a complicated tensor product between a third-order tensor and a vector. The latter is the same vector eld of the state discussed previously. Prior to computing 102 the third-order tensor, a review of the Jacobian of the vector eld is necessary. Using the model of the CR3BP the Jacobian can be organized as the sub-divided matrix, 21 @F (X(t)) @X(t) = 2 6 6 4 0 33 I 33 U XX 2 3 7 7 5 where, = 2 6 6 6 6 6 6 4 0 1 0 1 0 0 0 0 0 3 7 7 7 7 7 7 5 The sub-matrix of second partials of the potential, U XX , is square and symmetric. As such, only the upper triangular values represent unique quantities. These six values are, U xx = 1 + (1) 3(x +) 2 d 5 1 1 d 3 1 + 3(x 1 +) 2 d 5 2 1 d 3 2 U yy = 1 + (1) 3y 2 d 5 1 1 d 3 1 + 3y 2 d 5 2 1 d 3 2 U zz = (1) 3z 2 d 5 1 1 d 3 1 + 3z 2 d 5 2 1 d 3 2 U xy = 3(1)(x +)y d 5 1 + 3(x 1 +)y d 5 2 U xz = 3(1)(x +)z d 5 1 + 3(x 1 +)z d 5 2 U yz = 3(1)yz d 5 1 + 3yz d 5 2 where the distances from the spacecraft to the larger and smaller primaries, respectively, are, d 1 = p (x +) 2 +y 2 +z 2 and d 2 = p (x 1 +) 2 +y 2 +z 2 103 To continue with the third-order tensor, this object is handled in parts. After a block division into sub-tensors, the one non-trivial sub-tensor is further partitioned into an array of second-order tensors. @ 2 F (X(t)) @X 2 (t) = 2 6 6 4 0 333 0 333 U XXX 0 333 3 7 7 5 Once again, because of the conservative dynamics, the order of the partial derivatives is path independent, and thus, the dierential operators become commutative. This arrangement leads to a simplication in the third-order sub-tensor, U XXX , where a reduction to a subset of unique values is found. The reduced set of ten unique values are, U xxx = (1) 9(x +) d 5 1 15(x +) 3 d 7 1 + 9(x 1 +) d 5 2 15(x 1 +) 3 d 7 2 U yyy = (1) 9y d 5 1 15y 3 d 7 1 + 9y d 5 2 15y 3 d 7 2 U zzz = (1) 9z d 5 1 15z 3 d 7 1 + 9z d 5 2 15z 3 d 7 2 U xxy = (1) 3y d 5 1 15(x +) 2 y d 7 1 + 3y d 5 2 15(x 1 +) 2 y d 7 2 U xxz = (1) 3z d 5 1 15(x +) 2 z d 7 1 + 3z d 5 2 15(x 1 +) 2 z d 7 2 U yyx = (1) 3(x +) d 5 1 15y 2 (x +) d 7 1 + 3(x 1 +) d 5 2 15y 2 (x 1 +) d 7 2 U yyz = (1) 3z d 5 1 15y 2 z d 7 1 + 3z d 5 2 15y 2 z d 7 2 U zzx = (1) 3(x +) d 5 1 15z 2 (x +) d 7 1 + 3(x 1 +) d 5 2 15z 2 (x 1 +) d 7 2 U zzy = (1) 3y d 5 1 15z 2 y d 7 1 + 3y d 5 2 15z 2 y d 7 2 U xyz = (1) 15(x +)yz d 7 1 + 15(x 1 +)yz d 7 2 104 Finally, this lengthy analytical derivative can be implemented along with the standard computations performed during a LiAISON simulation. Noting the dicult components in particular, the Hessian of the measurement model can be computed alongside the Jacobian of the measurement model, and the third-order tensor product can be appended to the standard set of dierential equations used to solve for the state and transition matrix (i.e. increasing the set of coupled rst-order ordinary dierential equations from 42 equations to 78.). The remaining partials and chain rule products follow routinely. 105 Appendix C Initial Conditions for Periodic Orbits Of the four subsets of Earth-Moon periodic orbit families considered in this disserta- tion, the initial conditions of twelve representative orbits from each family are provided here. To obtain the full orbits from this data, the equations of motion, Eq. (2.2), are numerically integrated for each set of initial conditions over the associated periods. All coordinates are given in the non-dimensional units of the Earth-Moon CR3BP y . The columns correspond to the 7-tuple parametrization of a periodic orbit, = (X(t 0 );T ). Table C.1: L1 Halo Family x y z _ x _ y _ z T +8.29883e-1 -4.33772e-2 +2.22229e-2 -3.49612e-2 +1.00728e-1 +3.97724e-2 +2.74864 +8.31503e-1 -5.14032e-2 +4.12597e-2 -4.08076e-2 +1.20197e-1 +6.93870e-2 +2.76058 +8.34343e-1 -6.04698e-2 +6.04896e-2 -4.80967e-2 +1.42444e-1 +9.37526e-2 +2.77486 +8.38426e-1 -6.93133e-2 +8.06699e-2 -5.63839e-2 +1.64726e-1 +1.13426e-1 +2.78604 +8.43685e-1 -7.68927e-2 +1.01986e-1 -6.54274e-2 +1.85329e-1 +1.28189e-1 +2.78319 +8.50163e-1 -8.17242e-2 +1.23808e-1 -7.45754e-2 +2.02439e-1 +1.37312e-1 +2.74495 +8.58446e-1 -8.15148e-2 +1.44182e-1 -8.18237e-2 +2.13134e-1 +1.39753e-1 +2.63979 +8.69121e-1 -7.48696e-2 +1.60091e-1 -8.37717e-2 +2.13519e-1 +1.35724e-1 +2.46107 +8.80827e-1 -6.48438e-2 +1.70692e-1 -7.99503e-2 +2.04039e-1 +1.28834e-1 +2.26583 +8.91936e-1 -5.50760e-2 +1.78002e-1 -7.32873e-2 +1.89208e-1 +1.22240e-1 +2.10145 +9.01954e-1 -4.67177e-2 +1.83801e-1 -6.59274e-2 +1.72191e-1 +1.16853e-1 +1.97683 +9.10852e-1 -3.98024e-2 +1.89114e-1 -5.88040e-2 +1.54530e-1 +1.12882e-1 +1.88861 Note, the table is transcribed in a copyable syntax, for ease of direct import into source code. y Specied in Chapter 4. 106 Table C.2: L1 Axial Family x y z _ x _ y _ z T +9.09542e-1 +1.97933e-3 -1.82591e-1 -7.19838e-2 -4.85375e-2 +2.28816e-1 +4.06280 +9.17542e-1 +2.38876e-2 -1.81276e-1 -6.14178e-2 -7.82879e-2 +2.19305e-1 +4.05744 +9.25167e-1 +4.58935e-2 -1.77381e-1 -5.02673e-2 -1.06263e-1 +2.06685e-1 +4.04931 +9.32316e-1 +6.76741e-2 -1.70810e-1 -3.86795e-2 -1.31999e-1 +1.91362e-1 +4.03889 +9.38900e-1 +8.88669e-2 -1.61508e-1 -2.68703e-2 -1.55133e-1 +1.73811e-1 +4.02675 +9.44845e-1 +1.09076e-1 -1.49477e-1 -1.51215e-2 -1.75421e-1 +1.54539e-1 +4.01359 +9.50092e-1 +1.27882e-1 -1.34783e-1 -3.77033e-3 -1.92747e-1 +1.34039e-1 +4.00012 +9.54597e-1 +1.44856e-1 -1.17567e-1 +6.80929e-3 -2.07111e-1 +1.12755e-1 +3.98709 +9.58330e-1 +1.59576e-1 -9.80496e-2 +1.62304e-2 -2.18603e-1 +9.10441e-2 +3.97517 +9.61274e-1 +1.71649e-1 -7.65357e-2 +2.41198e-2 -2.27371e-1 +6.91656e-2 +3.96503 +9.63419e-1 +1.80727e-1 -5.34103e-2 +3.01453e-2 -2.33581e-1 +4.72762e-2 +3.95719 +9.64760e-1 +1.86536e-1 -2.91310e-2 +3.40411e-2 -2.37384e-1 +2.54414e-2 +3.95209 Table C.3: L2 Vertical Family x y z _ x _ y _ z T +1.14785e+0 +2.58749e-3 -4.90716e-2 +1.34538e-2 +1.15410e-5 +8.60874e-2 +3.54111 +1.13726e+0 +6.20638e-3 -7.47542e-2 +3.12476e-2 +4.39284e-4 +1.27441e-1 +3.57789 +1.12235e+0 +1.15953e-2 -9.97864e-2 +5.55984e-2 +1.55617e-3 +1.62914e-1 +3.64131 +1.10338e+0 +1.90512e-2 -1.24094e-1 +8.48474e-2 +3.97803e-3 +1.90512e-1 +3.74578 +1.07944e+0 +2.97464e-2 -1.49524e-1 +1.18207e-1 +8.63434e-3 +2.10416e-1 +3.92713 +1.04206e+0 +5.05144e-2 -1.87644e-1 +1.60501e-1 +1.84829e-2 +2.27954e-1 +4.36576 +9.96435e-1 +8.29502e-2 -2.53909e-1 +1.98090e-1 +2.84695e-2 +2.67046e-1 +5.15156 +9.57122e-1 +1.12755e-1 -3.32972e-1 +2.37029e-1 +2.97597e-2 +3.35890e-1 +5.67149 +9.19402e-1 +1.43131e-1 -4.04720e-1 +2.89495e-1 +2.80168e-2 +4.06012e-1 +5.91565 +8.59147e-1 +1.95706e-1 -4.96322e-1 +3.90136e-1 +2.45729e-2 +4.98872e-1 +6.08611 +7.92939e-1 +2.56825e-1 -5.70212e-1 +5.10974e-1 +2.11569e-2 +5.74744e-1 +6.16730 +7.22518e-1 +3.23804e-1 -6.26009e-1 +6.44591e-1 +1.79830e-2 +6.32247e-1 +6.21151 Table C.4: L5 Vertical Family x y z _ x _ y _ z T +4.89886e-1 -8.62129e-1 -6.88673e-2 +2.30915e-3 -4.12316e-3 +6.88640e-2 +6.28332 +5.37656e-1 -7.60561e-1 -3.47325e-1 +6.60011e-2 -1.10556e-1 +3.48450e-1 +6.28681 +6.05305e-1 -5.66307e-1 -5.49626e-1 +2.02547e-1 -3.20069e-1 +5.59589e-1 +6.29333 +6.40913e-1 -4.15043e-1 -6.38345e-1 +3.19054e-1 -4.99261e-1 +6.52842e-1 +6.29797 +6.59527e-1 -3.04973e-1 -6.80979e-1 +4.07701e-1 -6.41638e-1 +6.90678e-1 +6.30096 +6.73048e-1 -1.98138e-1 -7.07669e-1 +4.95218e-1 -7.89735e-1 +7.00972e-1 +6.30348 +6.82854e-1 -9.56136e-2 -7.20418e-1 +5.78706e-1 -9.39952e-1 +6.82421e-1 +6.30556 +6.90197e-1 +2.13891e-3 -7.20716e-1 +6.55623e-1 -1.08850e+0 +6.33941e-1 +6.30723 +6.96089e-1 +9.47667e-2 -7.09631e-1 +7.23386e-1 -1.23100e+0 +5.54586e-1 +6.30855 +7.01317e-1 +1.81720e-1 -6.88007e-1 +7.79075e-1 -1.36203e+0 +4.43787e-1 +6.30958 +7.06458e-1 +2.62045e-1 -6.56756e-1 +8.19273e-1 -1.47467e+0 +3.02143e-1 +6.31033 +7.11868e-1 +3.34293e-1 -6.17300e-1 +8.40278e-1 -1.56036e+0 +1.32972e-1 +6.31082 107 Appendix D GMAT Conguration Settings The following parameter set-up congures the GMAT software for simulating the Earth-Moon-Sun system, using the JPL DE405 ephemeris library. To propagate the trajectory of a specic spacecraft, the corresponding initial conditions to change are noted below with \<----". The header le to HORIZONS' output for a selection of lunar ephemerides is given at the end of this appendix (most of the actual numerical data is suppressed for the sake of brevity). Create ForceModel DefaultProp_ForceModel; GMAT DefaultProp_ForceModel.CentralBody = Earth; GMAT DefaultProp_ForceModel.PrimaryBodies = {Earth}; GMAT DefaultProp_ForceModel.PointMasses = {Luna, Sun}; GMAT DefaultProp_ForceModel.Drag = None; GMAT DefaultProp_ForceModel.SRP = Off; GMAT DefaultProp_ForceModel.ErrorControl = RSSStep; GMAT DefaultProp_ForceModel.GravityField.Earth.Degree = 4; GMAT DefaultProp_ForceModel.GravityField.Earth.Order = 4; http://ssd.jpl.nasa.gov/horizons.cgi 108 GMAT DefaultProp_ForceModel.GravityField.Earth.PotentialFile = 'JGM2.cof'; Create Propagator DefaultProp; GMAT DefaultProp.FM = DefaultProp_ForceModel; GMAT DefaultProp.Type = RungeKutta89; GMAT DefaultProp.InitialStepSize = 60; GMAT DefaultProp.Accuracy = 9.999999999999999e-12; GMAT DefaultProp.MinStep = 0.001; GMAT DefaultProp.MaxStep = 2700; GMAT DefaultProp.MaxStepAttempts = 50; GMAT DefaultProp.StopIfAccuracyIsViolated = true; Create Spacecraft DefaultSC; GMAT DefaultSC.DateFormat = UTCGregorian; GMAT DefaultSC.Epoch = '22 Feb 2011 00:00:00.000'; GMAT DefaultSC.CoordinateSystem = EarthMJ2000Eq; GMAT DefaultSC.DisplayStateType = Cartesian; GMAT DefaultSC.X = {input X coordinate in km}; <---- GMAT DefaultSC.Y = {input Y coordinate in km}; <---- GMAT DefaultSC.Z = {input Z coordinate in km}; <---- GMAT DefaultSC.VX = {input Xdot coordinate in km/s}; <---- GMAT DefaultSC.VY = {input Ydot coordinate in km/s}; <---- GMAT DefaultSC.VZ = {input Zdot coordinate in km/s}; <---- GMAT DefaultSC.DryMass = 850; GMAT DefaultSC.Cd = 2.2; GMAT DefaultSC.Cr = 1.8; 109 GMAT DefaultSC.DragArea = 15; GMAT DefaultSC.SRPArea = 1; GMAT DefaultSC.NAIFId = -123456789; GMAT DefaultSC.NAIFIdReferenceFrame = -123456789; GMAT DefaultSC.Id = 'SatId'; GMAT DefaultSC.Attitude = CoordinateSystemFixed; GMAT DefaultSC.ModelFile = '../data/vehicle/models/aura.3ds'; GMAT DefaultSC.ModelOffsetX = 0; GMAT DefaultSC.ModelOffsetY = 0; GMAT DefaultSC.ModelOffsetZ = 0; GMAT DefaultSC.ModelRotationX = 0; GMAT DefaultSC.ModelRotationY = 0; GMAT DefaultSC.ModelRotationZ = 0; GMAT DefaultSC.ModelScale = 3; GMAT DefaultSC.AttitudeDisplayStateType = 'Quaternion'; GMAT DefaultSC.AttitudeRateDisplayStateType = 'AngularVelocity'; GMAT DefaultSC.AttitudeCoordinateSystem = 'EarthMJ2000Eq'; GMAT DefaultSC.Q1 = 0; GMAT DefaultSC.Q2 = 0; GMAT DefaultSC.Q3 = 0; GMAT DefaultSC.Q4 = 1; GMAT DefaultSC.EulerAngleSequence = '321'; GMAT DefaultSC.AngularVelocityX = 0; GMAT DefaultSC.AngularVelocityY = 0; GMAT DefaultSC.AngularVelocityZ = 0; 110 Propagate Synchronized DefaultProp(DefaultSC) ... {DefaultSC.ElapsedSecs = {input T in s}}; <---- Create EphemerisFile EphemerisFile1; GMAT EphemerisFile1.Spacecraft = DefaultSC; GMAT EphemerisFile1.Filename = 'EphemerisFile1.eph'; GMAT EphemerisFile1.FileFormat = CCSDS-OEM; GMAT EphemerisFile1.EpochFormat = UTCGregorian; GMAT EphemerisFile1.InitialEpoch = InitialSpacecraftEpoch; GMAT EphemerisFile1.FinalEpoch = FinalSpacecraftEpoch; GMAT EphemerisFile1.StepSize = 86400; GMAT EphemerisFile1.InterpolationOrder = 7; GMAT EphemerisFile1.CoordinateSystem = EarthMJ2000Eq; GMAT EphemerisFile1.WriteEphemeris = true; ******************************************************************************* Revised: Mar 11, 1998 Moon / (Earth) 301 PHYSICAL PROPERTIES: Radius, km = 1737.53+-0.03 Mass, 10^20 kg = 734.9 Density, gm cm^-3 = 3.3437 Geometric albedo = 0.12 V(1,0) = +0.21 GM, km^3/s^2 = 4902.798+-.005 Earth/Moon mass ratio = 81.300587 Surface gravity = 1.62 m s^-2 Nearside crust. thick.= 58+-8 km Farside crust. thick. = ~80 - 90 km Heat flow, Apollo 15 = 3.1+-.6 mW/m^2 Heat flow, Apollo 17 = 2.2+-.5 mW/m^2 111 Mean crustal density = 2.97+-.07g/cm^3 k2 = 0.0302+-.0012 Induced magnetic mom. = 4.23x10^22Gcm^3 Magnetometer moment = 435+-15 DYNAMICAL CHARACTERISTICS: Mean angular diameter = 31'05.2" Orbit period = 27.321582 d Obliquity to orbit = 6.67 deg Eccentricity = 0.05490 Semi-major axis, a = 384400 km Inclination = 5.145 deg Mean motion, rad/s = 2.6616995x10^-6 Nodal period = 6798.38 d Apsidal period = 3231.50 d Mom. of inertia C/MR^2= 0.3935+-.0011 beta (C-A/B), x10^-4 = 6.31(72+-15) gamma (B-A/C), x10^-4 = 2.278(8+-2) ******************************************************************************* ******************************************************************************* Ephemeris / WWW_USER Wed Jan 25 14:52:58 2012 Pasadena, USA / Horizons ******************************************************************************* Target body name: Moon (301) {source: DE405} Center body name: Earth (399) {source: DE405} Center-site name: BODY CENTER ******************************************************************************* Start time : A.D. 2011-Feb-22 00:00:00.0000 CT Stop time : A.D. 2012-Feb-22 00:00:00.0000 CT Step-size : 1440 minutes ******************************************************************************* Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)} Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)} 112 Center radii : 6378.1 x 6378.1 x 6356.8 km {Equator, meridian, pole} Output units : KM-S Output format : 02 Reference frame : ICRF/J2000.0 Output type : GEOMETRIC cartesian states Coordinate systm: Earth Mean Equator and Equinox of Reference Epoch ******************************************************************************* JDCT X Y Z VX VY VZ ******************************************************************************* $$SOE 2455614.500000000 = A.D. 2011-Feb-22 00:00:00.0000 (CT) -3.316107438701022E+05 -1.241813610448045E+05 -8.598752928570415E+04 3.917114730478682E-01 -9.335564991250414E-01 -3.675006740226504E-01 ... 2455979.500000000 = A.D. 2012-Feb-22 00:00:00.0000 (CT) 3.486490133736398E+05 -1.745407557903522E+05 -3.871489546880240E+04 4.886983723361418E-01 8.039685584527730E-01 3.640527800231241E-01 $$EOE ******************************************************************************* Coordinate system description: Earth Mean Equator and Equinox of Reference Epoch Reference epoch: J2000.0 113 xy-plane: plane of the Earth's mean equator at the reference epoch x-axis : out along ascending node of instantaneous plane of the Earth's orbit and the Earth's mean equator at the reference epoch z-axis : along the Earth mean north pole at the reference epoch Symbol meaning JDCT Epoch Julian Date, Coordinate Time X x-component of position vector (km) Y y-component of position vector (km) Z z-component of position vector (km) VX x-component of velocity vector (km/sec) VY y-component of velocity vector (km/sec) VZ z-component of velocity vector (km/sec) Geometric states/elements have no aberration corrections applied. Computations by ... Solar System Dynamics Group, Horizons On-Line Ephemeris System 4800 Oak Grove Drive, Jet Propulsion Laboratory Pasadena, CA 91109 USA Information: http://ssd.jpl.nasa.gov/ Connect : telnet://ssd.jpl.nasa.gov:6775 (via browser) telnet ssd.jpl.nasa.gov 6775 (via command-line) Author : Jon.Giorgini@jpl.nasa.gov ******************************************************************************* 114
Abstract (if available)
Abstract
According to NASA’s integrated space technology roadmaps, space-based infrastructures are envisioned as necessary ingredients to a sustained effort in continuing space exploration. Whether it be for extra-terrestrial habitats, roving/cargo vehicles, or space tourism, autonomous space networks will provide a vital communications lifeline for both future robotic and human missions alike. Projecting that the Moon will be a bustling hub of activity within a few decades, a near-term opportunity for in-situ infrastructure development is within reach. ❧ This dissertation addresses the anticipated need for in-space infrastructure by investigating a general design methodology for autonomous interplanetary constellations
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Optimal guidance trajectories for proximity maneuvering and close approach with a tumbling resident space object under high fidelity J₂ and quadratic drag perturbation model
PDF
Relative-motion trajectory generation and maintenance for multi-spacecraft swarms
PDF
Periodic motion analysis in spacecraft attitude dynamics
PDF
Application of the fundamental equation to celestial mechanics and astrodynamics
PDF
Trajectory mission design and navigation for a space weather forecast
PDF
Techniques for analysis and design of temporary capture and resonant motion in astrodynamics
PDF
The development of an autonomous subsystem reconfiguration algorithm for the guidance, navigation, and control of aggregated multi-satellite systems
PDF
Designing an optimal software intensive system acquisition: a game theoretic approach
PDF
Incorporation of mission scenarios in deep space spacecraft design trades
PDF
Long term lunar radiation degradation of potential lunar habitat composite materials
PDF
Quantifying the effect of orbit altitude on mission cost for Earth observation satellites
PDF
Unreal engine testbed for computer vision of tall lunar tower assembly
PDF
Systems engineering and mission design of a lunar South Pole rover mission: a novel approach to the multidisciplinary design problem within a spacecraft systems engineering paradigm
PDF
Numerical and experimental investigations of dust-plasma-asteroid interactions
PDF
Three-dimensional exospheric hydrogen atom distributions obtained from observations of the geocorona in Lyman-alpha
PDF
Laboratory investigations of the near surface plasma field and charging at the lunar terminator
PDF
Mission design study of an RTG powered, ion engine equipped interstellar spacecraft
PDF
Lobster eye optics: a theoretical and computational model of wide field of view X-ray imaging optics
PDF
Hydrogen peroxide vapor for small satellite propulsion
PDF
Advanced nuclear technologies for deep space exploration
Asset Metadata
Creator
Chow, Cornelius Channing, II
(author)
Core Title
Autonomous interplanetary constellation design
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Astronautical Engineering
Publication Date
04/06/2012
Defense Date
02/21/2012
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
astrodynamics,astronautics,Autonomous,bifurcation,constellation,CR3BP,Design,dynamical system,interplanetary,LiAISON,navigation,numerical continuation,OAI-PMH Harvest,optimization,pseudo-arclength,satellite,spacecraft
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Erwin, Daniel A. (
committee chair
), Gruntman, Michael (
committee member
), Lo, Martin W. (
committee member
), Rhodes, Edward J., Jr. (
committee member
), Villac, Benjamin F. (
committee member
)
Creator Email
channinc@usc.edu,channing.chow@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-3341
Unique identifier
UC11288184
Identifier
usctheses-c3-3341 (legacy record id)
Legacy Identifier
etd-ChowCornel-575.pdf
Dmrecord
3341
Document Type
Dissertation
Rights
Chow, Cornelius Channing, II
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
astrodynamics
astronautics
bifurcation
constellation
CR3BP
dynamical system
interplanetary
LiAISON
navigation
numerical continuation
optimization
pseudo-arclength
satellite
spacecraft