Dynamics and Control of Satellite Relative Motion in Elliptic Orbits using Lyapunov-Floquet Theory by Ryan Edward Sherrill A dissertation submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama May 4, 2013 Copyright 2013 by Ryan Edward Sherrill Approved by Andrew J. Sinclair, Chair, Associate Professor of Aerospace Engineering John E. Cochran Jr., Professor and Head of Aerospace Engineering David A. Cicci, Professor of Aerospace Engineering Subhash C. Sinha, Professor of Mechanical Engineering T. Alan Lovell, Research Aerospace Engineer, Air Force Research Laboratory George Flowers, Dean of the Graduate School Abstract Rendezvous and proximity operations involving satellites operating near each other have been performed since the Gemini missions. A better understanding of the dynamics of satellites in close proximity could be helpful for both mission planners and operators as well as for algorithm development. With this increased understanding, future satellite missions may involve an impromptu fly-around of another satellite and a greater focus on satellite autonomy. In satellite proximity operations, it is desirable to express the relative motion of a deputy satellite with respect to the chief satellite. For satellites in elliptic orbits, this relative motion can be described by the time-varying Linearized Equations of Relative Motion. Making an assumption that the chief satellite is in a circular orbit results in the linear time-invariant Hill-Clohessy-Wiltshire equations. These equations allow the relative-motion solution to be intuitively visualized. Due to the abundance of algorithms in the literature based on the Hill-Clohessy-Wiltshire equations, the focus of this dissertation is to explore the relationship between relative motion in elliptic orbits and the Hill-Clohessy-Wiltshire equations. The major contributions of this work are twofold. First, three time-varying coordinate transformations are derived which relate the Hill-Clohessy-Wiltshire equations to the Lin- earized Equations of Relative Motion. These transformations show that the Hill-Clohessy- Wiltshire equations are able to exactly capture the time-varying dynamics. Second, the literature does not contain infinite-horizon continuous-thrust control for satellites in elliptic orbits. A time-varying gain matrix was constructed using the time-varying coordinate trans- formations, optimal-control theory, and a control method from the literature. This control is able to drive the deputy?s position toward rendezvous with an elliptic chief. ii Acknowledgments The author is forever indebted to Dr. Andrew J. Sinclair for the guidance, patience, and encouragement shown during his graduate education, and for the help given while completing this dissertation. Our many conversations have lead to a deeper understanding of dynamics and established the foundation for my future career. The author also thanks the following members of his committee for their time, suggestions, and criticisms of this dissertation: Dr. John E. Cochran Jr., Dr. David A. Cicci, Dr. Subhash C. Sinha, Dr. T. Alan Lovell, and Dr. John Y. Hung. In addition, the author thanks Dr. R. Steven Gross for a decade of helpful advisement and the opportunity to teach. The author also acknowledges the financial support offered by the Department of Aerospace Engineering. Three summers of research experience at the Air Force Research Laboratory, Space Vehicles Directorate were provided by the American Society for Engineering Education?s Summer Faculty Fellowship Program and the Air Force Research Laboratory?s Space Scholars Program. Finally, the author thanks his beautiful fiancee Judith Ann Bailey for her constant support, love, and encouragement. iii When I Heard the Learn?d Astronomer When I heard the learn?d astronomer, When the proofs, the figures, were ranged in columns before me, When I was shown the charts and the diagrams, to add, divide, and measure them, When I, sitting, heard the astronomer, where he lectured with much applause in the lecture-room, How soon, unaccountable, I became tired and sick, Till rising and gliding out, I wanderd off by myself, In the mystical moist night-air, and from time to time, Lookd up in perfect silence at the stars. -Walt Whitman iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Orbital Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Orbital Motion of a Spacecraft . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 Kepler?s Laws of Planetary Motion . . . . . . . . . . . . . . . . . . . 3 2.1.2 The Two-Body Problem . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.3 The Orbit in Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.4 Orbital Integrals of Motion . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.5 Kepler?s Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.6 Perturbations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Relative Motion of Two Spacecraft . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.1 Nonlinear Equations of Relative Motion . . . . . . . . . . . . . . . . 16 2.2.2 Linear Equations of Relative Motion . . . . . . . . . . . . . . . . . . 19 2.2.3 Tschauner-Hempel Equations . . . . . . . . . . . . . . . . . . . . . . 20 2.2.4 Hill-Clohessy-Wiltshire Equations . . . . . . . . . . . . . . . . . . . . 24 2.2.5 Geometric Parameterization of the HCW Equations . . . . . . . . . . 26 2.2.6 Summary of Cartesian Linearized Relative-Motion Representations . 28 2.2.7 Additional Relative-Motion Representations . . . . . . . . . . . . . . 29 2.2.8 Calculation of Initial Conditions . . . . . . . . . . . . . . . . . . . . . 30 3 Preliminary Investigations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 v 3.1 Virtual-Chief Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Virtual-Time Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4 Lyapunov-Floquet Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.1 Floquet Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.1.1 Discussion of Floquet Theorem . . . . . . . . . . . . . . . . . . . . . 45 4.1.2 System Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2 Lyapunov-Floquet Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.3 Determining a Lyapunov-Floquet Transformation . . . . . . . . . . . . . . . 46 4.3.1 Commutative System . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.3.2 Non-commutative System . . . . . . . . . . . . . . . . . . . . . . . . 49 4.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.4 Lyapunov-Floquet Transformation of the TH Equations . . . . . . . . . . . . 51 4.5 Lyapunov-Floquet Generalization of the HCW Equations in the f-Domain . 53 5 Time-Varying Coordinate Transformations of the LERM . . . . . . . . . . . . . 56 5.1 Lyapunov-Floquet Transformations . . . . . . . . . . . . . . . . . . . . . . . 56 5.1.1 Periapse-Matching Transformation . . . . . . . . . . . . . . . . . . . 56 5.1.2 Apoapse-Matching Transformation . . . . . . . . . . . . . . . . . . . 60 5.1.3 Relative Motion Examples . . . . . . . . . . . . . . . . . . . . . . . . 63 5.2 Integral-Preserving Transformation . . . . . . . . . . . . . . . . . . . . . . . 68 5.3 Calibration of Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.3.1 Calculation of Initial Conditions . . . . . . . . . . . . . . . . . . . . . 70 5.3.2 Comparison of Methods . . . . . . . . . . . . . . . . . . . . . . . . . 71 6 Control of Relative Orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.1 Impulsive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.2 Linear Quadratic Regulator . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 6.3 Control of Time-Periodic Systems using Lyapunov-Floquet Theory . . . . . . 88 vi 6.3.1 Control Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6.3.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.4 Control of Satellites in Elliptic Orbits . . . . . . . . . . . . . . . . . . . . . . 94 6.4.1 Near-circular Orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.4.2 Highly Elliptic Orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.4.3 Effects of Eccentricity on Control Effort . . . . . . . . . . . . . . . . 122 6.4.4 Choice of Periapse-Matching or Apoapse-Matching Transformation . 127 7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 A Elements of P(f) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 B Elements of ?P(f) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 C Elements of ?(f) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 vii List of Figures 2.1 Kepler?s First and Second Laws. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Inertial coordinate frame. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 The Earth centered inertial (ECI) coordinate frame. . . . . . . . . . . . . . . . 7 2.4 Orbital Elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.5 The circumscribed auxiliary circle and eccentric anomaly. . . . . . . . . . . . . . 13 2.6 Magnitude of orbital perturbations vs satellite radius. . . . . . . . . . . . . . . . 15 2.7 Local-vertical local-horizonal (LVLH) coordinate frame . . . . . . . . . . . . . . 17 2.8 In-plane relative-orbit elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.9 Different Representations of Relative Motion Equations . . . . . . . . . . . . . . 29 3.1 Virtual-Chief Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2 NERM trajectories for chief eccentricities of 0, 0.01, 0.1, 0.2, and 0.4. . . . . . . 36 3.3 Solutions for the virtual time for several chief eccentricities. . . . . . . . . . . . 38 4.1 Different Representations of Relative Motion Equations . . . . . . . . . . . . . . 55 5.1 Different Representations of Relative Motion Equations . . . . . . . . . . . . . . 60 5.2 LERM, HCW, and LF transformation for case 1. . . . . . . . . . . . . . . . . . 64 viii 5.3 LERM, HCW, and LF transformation for case 2. . . . . . . . . . . . . . . . . . 64 5.4 LERM, HCW, and LF transformation for case 3. . . . . . . . . . . . . . . . . . 65 5.5 LERM, HCW, and LF transformation for case 4. . . . . . . . . . . . . . . . . . 65 5.6 LERM, HCW, and LF transformation for case 5. . . . . . . . . . . . . . . . . . 66 5.7 LERM, HCW, and LF transformation for case 6. . . . . . . . . . . . . . . . . . 66 5.8 Three dimensional representation of Case 5. . . . . . . . . . . . . . . . . . . . . 67 5.9 Three dimensional representation of Case 6. . . . . . . . . . . . . . . . . . . . . 67 5.10 Relative-motion trajectories for LERM and HCW propagation. . . . . . . . . . 70 5.11 Relative-motion trajectory for the LERM and three approximate HCW solutions for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.12 Relative-motion trajectory for the LERM and three approximate HCW solutions for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.13 Relative-motion trajectory for the LERM and three approximate HCW solutions for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.14 Relative-motion trajectory for the LERM and three approximate HCW solutions for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.15 Out-of-plane trajectory for the LERM and three approximate HCW solutions for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.16 Out-of-plane trajectory for the LERM and three approximate HCW solutions for case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 ix 5.17 Three dimensional representation of the relative-motion trajectory for case 5. . . 75 5.18 Three dimensional representation of the relative-motion trajectory for case 6. . . 75 5.19 Error distance as a function of time for the three underlying HCW trajectories in case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.20 Error distance as a function of time for the three underlying HCW trajectories in case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.21 Error distance as a function of time for the three underlying HCW trajectories in case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.22 Error distance as a function of time for the three underlying HCW trajectories in case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.23 Error distance as a function of time for the three underlying HCW trajectories in case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.24 Error distance as a function of time for the three underlying HCW trajectories in case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 6.1 Two-burn maneuver over one period. . . . . . . . . . . . . . . . . . . . . . . . . 84 6.2 Two-burn maneuver over five periods. . . . . . . . . . . . . . . . . . . . . . . . 84 6.3 Two-burn maneuver over one period. . . . . . . . . . . . . . . . . . . . . . . . . 85 6.4 Two-burn maneuver over half a period. . . . . . . . . . . . . . . . . . . . . . . . 85 6.5 Uncontrolled states of the system. . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.6 State feedback control of time-invariant auxiliary system. . . . . . . . . . . . . . 91 x 6.7 State feedback control of the time-varying system. . . . . . . . . . . . . . . . . . 92 6.8 Pole locations resulting in unstable time-varying system. . . . . . . . . . . . . . 93 6.9 Detailed study of the stability boundary. . . . . . . . . . . . . . . . . . . . . . . 93 6.10 Deputy trajectory case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.11 Deputy trajectory case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.12 Deputy trajectory case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.13 Deputy trajectory case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.14 Deputy trajectory case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.15 Deputy trajectory case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.16 Position error for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.17 Position error for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.18 Position error for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.19 Position error for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.20 Position error for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.21 Position error for case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.22 Control effort for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.23 Control effort for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.24 Control effort for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 xi 6.25 Control effort for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.26 Control effort for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.27 Control effort for case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.28 Deputy trajectory for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.29 Deputy trajectory for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.30 Deputy trajectory for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.31 Deputy trajectory for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.32 Deputy trajectory for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.33 Deputy trajectory for case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.34 Position error for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.35 Position error for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.36 Position error for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.37 Position error for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.38 Position error for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.39 Position error for case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 6.40 Control effort for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.41 Control effort for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.42 Control effort for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 xii 6.43 Control effort for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.44 Control effort for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.45 Control effort for case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.46 Position error for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.47 Position error for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.48 Position error for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.49 Control effort for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.50 Control effort for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 6.51 Control effort for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 6.52 Position error for periapse-matching transformation and apoapse-matching trans- formation for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.53 Position error for periapse-matching transformation and apoapse-matching trans- formation for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.54 Position error for periapse-matching transformation and apoapse-matching trans- formation for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.55 Position error for periapse-matching transformation and apoapse-matching trans- formation for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.56 Position error for periapse-matching transformation and apoapse-matching trans- formation for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 6.57 Position error for periapse-matching transformation and apoapse-matching trans- formation for case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 xiii 6.58 Control effort for periapse-matching transformation and apoapse-matching trans- formation for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.59 Control effort for periapse-matching transformation and apoapse-matching trans- formation for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.60 Control effort for periapse-matching transformation and apoapse-matching trans- formation for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.61 Control effort for periapse-matching transformation and apoapse-matching trans- formation for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.62 Control effort for periapse-matching transformation and apoapse-matching trans- formation for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.63 Control effort for periapse-matching transformation and apoapse-matching trans- formation for case 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 A.1 Elements of P(f) (1 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 A.2 Elements of P(f) (2 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 A.3 Elements of P(f) (3 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 A.4 Elements of P(f) (4 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 A.5 Elements of P(f) (5 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 A.6 Elements of P(f) (6 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 B.1 Elements of ?P(f) (1 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 B.2 Elements of ?P(f) (2 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 xiv B.3 Elements of ?P(f) (3 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 B.4 Elements of ?P(f) (4 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 B.5 Elements of ?P(f) (5 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 B.6 Elements of ?P(f) (6 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 C.1 Elements of ?(f) (1 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 C.2 Elements of ?(f) (2 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 C.3 Elements of ?(f) (3 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 C.4 Elements of ?(f) (4 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 C.5 Elements of ?(f) (5 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 C.6 Elements of ?(f) (6 of 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 xv List of Tables 3.1 Initial conditions for inplane test cases. . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 RMS error distances for inplane test cases. . . . . . . . . . . . . . . . . . . . . . 42 5.1 Orbital elements for the chief and deputy for six different relative motion cases. 63 5.2 Mean error (km) for each case over the simulation time. . . . . . . . . . . . . . 80 5.3 RMS error (km) for each case over the simulation time. . . . . . . . . . . . . . . 80 6.1 Orbital elements for the chief and deputy for six near circular relative motion cases. 97 6.2 Eigenvalues of ?x(P,0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.3 Eigenvalues of ??x(P,0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.4 Final position and total control effort for the auxiliary system tildewidez(t). . . . . . . . 98 6.5 Final position and total control effort for the LF transformed HCW system z(t). 98 6.6 Final position and total control effort for the LF transformed LERM system x(t). 98 6.7 Final position and total control effort for the direct application of the gain matrix to the LERM ?x(t). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.8 Orbital elements for the chief and deputy for six different relative motion cases. 109 6.9 Eigenvalues of ?x(P,0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.10 Final position and total control effort for the auxiliary system tildewidez(t). . . . . . . . 110 6.11 Final position and total control effort for the LF transformed HCW system z(t). 110 6.12 Final position and total control effort for the LF transformed LERM system x(t).110 6.13 Final position and total control effort for the auxiliary system tildewidez(t). . . . . . . . 123 6.14 Final position and total control effort for the LF transformed HCW system z(t). 123 6.15 Final position and total control effort for the LF transformed LERM system x(t).123 6.16 Difference in final position and control effort between the periapse-matching transformation and the apoapse-matching transformation . . . . . . . . . . . . . 127 xvi Chapter 1 Introduction Missions involving satellites operating in close proximity to each other are commonly referred to as rendezvous and proximity operations. Rendezvous and proximity operations are not a recent development; the Gemini and Apollo missions, as well as over forty percent of all Space Shuttle missions, involved on-orbit docking. However, missions currently being planned may involve an impromptu fly-around of a satellite and a greater focus on satellite autonomy. For these types of missions, thoroughly understanding the dynamics is helpful for both mission planners and operators, as well as for engineers designing the computer algorithms used onboard the satellite. In satellite proximity operations, it is often desirable to describe the motion of a ?deputy? satellite relative to a ?chief? satellite. Several models exist to predict these mo- tions. One choice is the Hill-Clohessy-Wiltshire equations, which are a set of linearized time- invariant ordinary differential equations, derived by assuming a circular chief orbit. Deriving linearized equations of motion without assuming a circular chief orbit results in time-varying differential equations. The relative motion can also be predicted by propagating the orbit of each satellite relative to the Earth. Among these approaches, the Hill-Clohessy-Wiltshire equations offer an advantage in the ability to understand and intuitively visualize the ge- ometry of the relative-motion solutions. The solutions obtained via linearized time-varying equations or by differencing the nonlinear orbit propagations of the chief and deputy are more accurate in describing elliptic orbits, but can be difficult to visualize conceptually. The focus of this dissertation is to explore the relationship between relative motion in elliptic orbits and the Hill-Clohessy-Wiltshire equations. This exploration is driven by the abundance of algorithms in the literature based on the Hill-Clohessy-Wiltshire equations and 1 the intuitive understanding these equations provide. By using Lyapunov-Floquet theory, several time-varying transformations are presented that relate the linearized time-varying equations of motion to the Hill-Clohessy-Wiltshire equations. These transformations show that the Hill-Clohessy-Wiltshire equations can be generalized to describe relative motion in elliptic orbits. Since the Hill-Clohessy-Wiltshire equations are able to capture the exact time-varying dynamics, existing methods for describing circular relative motion could be utilized for elliptic orbits. One application shown in this dissertation is continuous-thrust maneuvering in elliptic orbits based on infinite-horizon optimal control. The following chapter presents a review of orbital mechanics including basic formula- tions and several relative-motion equations. Next, preliminary investigations into the gen- eralization of the Hill-Clohessy-Wiltshire equations are discussed. Lyapunov-Floquet theory is then introduced, including an existing transformation from the literature and the author?s Lyapunov-Floquet generalization of the Hill-Clohessy-Wiltshire equations. The subsequent chapter presents three time-varying coordinate transformations which are able to relate the linearized time-varying equations to the Hill-Clohessy-Wiltshire equations. This chapter also includes using the coordinate transformations to select initial conditions to approxi- mate the motion using the Hill-Clohessy-Wiltshire equations. Finally, control algorithms and continuous-thrust maneuvering are discussed. The chapters in this dissertation are organized with background material and novel con- tributions conjointly presented. The specific contributions of the author to the field of satel- lite relative motion can be found in the following sections: Chapter 3, the Virtual-Chief and Virtual-Time methods; Section 4.5, Lyapunov-Floquet Generalization of the Hill-Clohessy- Wiltshire Equations in thef-Domain; Chapter 5, Time-Varying Coordinate Transformations of the linear equations of relative motion; Section 6.3.2, Discussion of the Lyapunov-Floquet Based Control Method; and Section 6.4, Continuous-Thrust Control of Satellites in Elliptic Orbits. 2 Chapter 2 Orbital Mechanics Classically, what is now called orbital mechanics was once referred to as celestial me- chanics, or the study of the motion of heavenly bodies. A few of the well-known classical philosophers and astronomers included: Hipparchus who cataloged over 1000 stars by bright- ness and expanded the theory of epicycles; Claudius Ptolemy who?s Almagest was a complete exposition of astronomy known to the Greeks and became the definitive work on astronomy for over 1000 years; and Copernicus who proposed the heliocentric model of the universe.1,2,3 Since the beginning of the ?Space Age? following the launch of Sputnik in 1957, the methods originally used to describe the motion of the planets have been adapted to artificial satellites. The first half of this chapter focuses on the basic laws of orbital mechanics which describe the motion of a satellite relative to the Earth developed by two seventeenth century mathematicians: Johannes Kepler and Isaac Newton. After discussing the motion of one body orbiting another, the later parts of the chapter will apply these principles to the problem of relative motion between objects in neighboring orbits. 2.1 Orbital Motion of a Spacecraft 2.1.1 Kepler?s Laws of Planetary Motion The Danish astronomer Tycho Brahe collected accurate measurements of the movement of the planets, and willed his observations to Johannes Kepler upon his death in 1601. Kepler was able to formulate his three comprehensive laws of planetary motion through intensive study of Brahe?s empirical data. In 1609, Kepler published Astronomia nova containing his first two laws of planetary motion. His third law was published in 1619 in Harmonices Mundi. Kepler?s three laws of planetary motion are stated below.4 3 Figure 2.1: Kepler?s First and Second Laws. 1. The orbit of each planet is an ellipse with the sun at one focus. 2. The line from the sun to a planet sweeps out equal areas inside the ellipse in equal lengths of time. 3. The squares of the orbital periods of the planets are proportional to the cubes of their mean distance from the Sun. Kepler?s first law describes the shape of a planet?s orbit. As depicted in Figure 2.1, the planets?s orbit is an ellipse in the orbital plane. The Sun occupies one focus while the other focus of the ellipse remains empty. As a planet travels in its orbit, its velocity changes with its position. The shaded regions in Figure 2.1 were traced out in equal amounts of time, demonstrating Kepler?s second law. While Kepler?s first two laws describes a planet?s orbit and how it traverses that orbit, the third law is a mathematical relationship describing how the motion of different planets are related as they revolve around the Sun. It should be 4 Y Z Xr? r? m? r m? Figure 2.2: Inertial coordinate frame. noted that Kepler?s laws apply to any body in orbit of another. Therefore, the same forces apply to a comet orbiting the Sun, or the Space Shuttle orbiting the Earth. 2.1.2 The Two-Body Problem In the Philosophiae Naturalis Principia Mathematica, Isaac Newton introduced the three fundamental laws of classical mechanics. The development of the two-body problem, which can describe the motion of a satellite with respect to the Earth, begins with a mathematical representation of Newton?s Second Law. summationdisplay F = ma (2.1) Here, F represents the forces acting on an object, m is the mass of the object, and a represents the acceleration of the object. Consider the motion of two point masses in an inertial coordinate system shown in Figure 2.2. The point masses are attracted to each other with a magnitude given by Newton?s Law of Universal Gravitation F = Gm1m2r2 where G is the gravitational constant, m1 and m2 are the masses of the bodies and r is the distance between the center of the masses. Equations of motion for each mass can be written in vector 5 form. m1?r1 = Gm1m2r3 r (2.2) m2?r2 = ?Gm1m2r3 r (2.3) Here, r represents the relative motion of the bodies given by r = r2 ?r1, and an overdot (?) represents a derivative with respect to time. Differentiating this equation twice gives the relative acceleration of the bodies. ?r = ?r2 ? ?r1 = ?G(m1 +m2)r3 r (2.4) Define the gravitational parameter ? such that ? = G(m1 +m2). Often when calculating the gravitational parameter, the mass of one of the bodies will be negligible such as the case of a satellite compared to the Earth. The relative equation of motion can be written in the familiar form. ?r = ??r3r (2.5) 2.1.3 The Orbit in Space In order to describe the position of a satellite in orbit, a coordinate system must be chosen. One choice is a Cartesian system shown in Figure 2.3. In this system, the origin is located at the center of the Earth where the X axis points toward the first point of Aries, symbolized by aries, the Z axis goes through the geographic North Pole, and the Y axis lies in the equatorial plane and completes the orthogonal triad. Using a Cartesian coordinate sys- tem, a satellite?s state can be completely described by six independent parameters: position (x,y,z) and velocity (?x, ?y, ?z). A disadvantage of selecting a Cartesian coordinate system is that, with the exception of a few select cases, visualizing a satellite?s orbit is difficult. Therefore, it is often advantageous to select another set of coordinates called orbital elements. When using orbital elements, five of the six parameters needed to describe a satellite?s state are constant and represent the 6 Y Z X Figure 2.3: The Earth centered inertial (ECI) coordinate frame. size, shape, and orientation of the orbit. The sixth parameter represents the motion along that orbit. The six orbital elements, shown Figure 2.4, are stated below. a: The semi-major axis describes the size of the ellipse. e: The eccentricity of the orbit which describes the shape of the ellipse. i: The inclination, or angle, of the orbital plane from the equator. ?: The right ascension of the ascending node, is the angle from the vernal equinox to the point on the equator where the satellite makes its south to north crossing. ?: The argument of perigee is the angle from the ascending node to perigee. f: The true anomaly represents the angle between the satellite?s current position in the orbit and perigee. 7 Y Z X ? i ? f Nodal Line periapse satellite geometric center a focus periapseapoapse satellite f ae line of apsides Figure 2.4: Orbital Elements. 8 As can be seen in Figure 2.4, the closest point on the ellipse to the focus is called periapse while apoapse is defined as the point on the ellipse furthest from the focus. The line of apsides passes through both periapse and apoapse. The focus is located a distance ae along the apse line from the geometric center of the ellipse. Note that the true anomaly is measured from periapse, i.e. the closest point on the ellipse to the focus occurs at f = 0. 2.1.4 Orbital Integrals of Motion In a polar coordinate frame (?ir ?i? ?iz) centered at m1, the position of the satellite is given by r = r?ir, and the velocity is given v = ?r = ?r?ir + r ?f?i?. The massless angular momentum can be determined by taking the cross product of these vectors. h=r?v = r2 ?f?iz (2.6) It is easy to see thathis a constant, and that the orbital motion lies in a plane perpendicular to h. Multiply both sides of Equation (2.5) by h (note that ?r is written as dvdt). dv dt ?h= ? ? r3r?h d dt (v?h) = ? df dt ?i? (2.7) Integrating both sides of Equation (2.7) and simplifying gives the following. v?h= ?rr+?e (2.8) Here, ?e is the constant of integration and e is the eccentricity vector. The eccentricity vector points toward periapse. The magnitude of the eccentricity vector is given by the following. e2 =e?e= 1?2 (v?h)?(v?h)? 2?rr?v?h+ 1 (2.9) 9 Using the properties of cross products, (v?h) = vh andr?v?h=r?v?h, Equation (2.9) can be written as shown below. e2 = 1?2h2v2 ? 2?rh2 + 1 e2 ?1 = h 2 ? parenleftbiggv2 ? ? 2 r parenrightbigg (2.10) 1?e2 = h 2 ? parenleftbigg2 r ? v2 ? parenrightbigg Define p = h 2 ? and a = parenleftbigg2 r ? v2 ? parenrightbigg?1 where p is called the parameter. From Equation (2.10), the parameter is related to a and e by the following. p = h 2 ? = a parenleftbig1?e2parenrightbig (2.11) The total energy of the system is related to a by the following. v2 2 ? ? r = ? u 2a (2.12) Where the kinetic energy is given by v 2 2 , potential energy by ? ? r, and ? ? 2a is the energy constant. The above expression can be rewritten as the energy integral or the vis-viva equation (latin for ?live force?). v2 = ? parenleftbigg2 r ? 1 a parenrightbigg (2.13) Solving Equation (2.8) for e and taking the dot product with position vector gives the following. e?r = 1? parenleftBig v?h??rr parenrightBig ?r ercosf = 1? (r?v?h??r) (2.14) ercosf = h 2 ? ?r 10 Solving Equation (2.14) for r and substituting for p gives the equation of an orbit. r = p1 +ecosf (2.15) 2.1.5 Kepler?s Equation Substituting the square of Equation (2.15) into Equation (2.6) gives the following ex- pression. h = p 2 (1 +ecosf)2 df dt ??p = p2 (1 +ecosf)2 df dt (2.16)radicalbigg ? p3dt = df (1 +ecosf)2 Integrating Equation (2.16) provides the relationship between time and true anomaly. radicalbigg? p3 (t?tp) = integraldisplay f 0 df (1+ecosf)2 (2.17) Here, the constant of integration is given by tp, and is the time when f = 0. For e = 0, Equation (2.17) can be easily integrated. radicalbigg? p3 (t?tp) = f (2.18) For e = 0, p = a, which allows Equation (2.18) to be written as the following. f = radicalbigg? a3 (t?tp) (2.19) Define the mean motion n = radicalbigg? a3. The mean motion is an average value of how fast a satellite progresses in its orbit. The period of an orbit P, is related to the mean motion 11 by P = 2?n . For circular orbits, the mean motion equals the satellite?s angular velocity. Equation (2.19) can be written in the following form which gives the relationship between time and true anomaly for circular orbits. f = n(t?tp) (2.20) Analogous to Equation (2.20), a relationship between time and true anomaly for elliptic orbits will now be developed. Evaluating the right hand side of Equation (2.17) for e negationslash= 0 results in the following.5 integraldisplay f 0 df (1 +ecosf)2 = 1 (1?e2)32 parenleftBigg 2tan?1 radicalbigg1?e 1 +e tan f 2 ? e?1?e2 sinf 1 +ecosf parenrightBigg (2.21) Define the mean anomaly, M, as the following. M = 2tan?1 radicalbigg1?e 1+e tan f 2 ? e?1?e2 sinf 1 +ecosf (2.22) Substituting the result of Equation (2.21) and (2.22) into Equation (2.17) gives the following relationship. radicalbigg ? p3 (t?tp) = 1 (1?e2)32 M (2.23) Solving Equation (2.23) for M, and substituting p = a(1 ? e2) and n = radicalbigg? a3 gives the relationship between time and mean anomaly. M = n(t?tp) (2.24) Note that for elliptic orbits, the mean motion does not equal the angular velocity, but does provide an average value of the satellite?s progression. The expression for mean anomaly given by Equation (2.22) can be simplified by introducing an auxiliary angle, E, called the 12 X Y E f ae r a Auxiliary Circle Ellipse satellite Figure 2.5: The circumscribed auxiliary circle and eccentric anomaly. eccentric anomaly shown in Figure 2.5. Since the eccentric anomaly is measured from the geometric center, the satellite?s position on the ellipse can expressed in parametric form by x = acosE, and y = bsinE = a?1?e2 sinE. In terms of the true anomaly, these expressions are given by x = ae+rcosf and y = rsinf. By equating the x components, the following expression is developed. acosE = ae+rcosf (2.25) Using Equation (2.15) to substitute for r in Equation (2.25) and simplifying gives the fol- lowing. cosE = e+ cosf1 +ecosf (2.26) The above expression for cosE is not ideal because for a given value of cosE, there are two corresponding values of E. To eliminate this ambiguity, the trigonometric identity 13 tan2 E2 = sin2 E2 cos2 E2 = 1?cosE1 + cosE is used. Here, 1?cosE = (1?e)(1?cosf)1 +ecosf and 1+cosE = (1 +e)(1+ cosf) 1 +ecosf . Upon simplification, the following relationship between true anomaly and eccentric anomaly can be determined. tanE2 = radicalbigg1?e 1 +e tan f 2 (2.27) Solving this expression for E gives the following. E = 2tan?1 parenleftBiggradicalbigg 1?e 1 +e tan f 2 parenrightBigg (2.28) Recall that in termsof themean anomalythedeputy?syposition isgiven byy = a?1?e2 sinE, while in terms of the true anomaly the position is given by y = rsinf. Equating these ex- pressions, solving for sinE, and using Equation (2.15) gives the following. sinE = ?1?e2 sinf 1 +ecosf (2.29) Using Equations (2.28) and (2.29), it is clear that Equation (2.22) can be simplified. M = n(t?tp) = E?esinE (2.30) Equation (2.30) is known as Kepler?s Equation which relates time to eccentric anomaly, which can be related to true anomaly through Equation (2.27). This relationship is considerably more complicated than the expression for circular orbits, and requires an iterative technique to transform time to true anomaly. 14 Figure 2.6: Magnitude of orbital perturbations vs satellite radius. (Reproduced from Refer- ence 6) 2.1.6 Perturbations The discussion presented in this chapter has only considered gravitational attraction of a point masses. While this basic model has allowed the development of Keplerian orbits, it is only an approximation of true satellite motion. The Earth, in fact, is not a point mass. 15 Due to the Earth?s rotation, the equatorial radius is approximately twenty kilometers greater than the polar radius. In addition, the mass of the Earth is not centered at a point, but distributed unequally under the surface. Models with increasing complexity up to seventy terms exist to model the zonal, tesseral, and sectorial harmonics. The mass distribution is also time-varying due to tidal effects. In addition to gravity, atmospheric effects such as drag perturb satellites in low Earth orbit. The motion of the solar system also effects a satellite due to the presence and movement of the Moon, Venus, and Jupiter. The Sun?s gravity effects a satellite as well as causing solar radiation pressure and solar wind. Figure 2.6, reproduced from Montenbruck and Gill, shows the relative magnitude of the acceleration caused by these forces as a function of orbital radius.6 While these perturbations will disturb a satellite from its orbit, these forces will not be considered in this dissertation. The idealization of point masses allows analytical solutions of the motion (and relative motion) to be constructed. These solutions provide geometrical interpretations of the orbit, and allow for intuitive understanding of the orbit geometry. In addition, the idealized problem provides a foundation for the future study of more compli- cated and realistic motion. 2.2 Relative Motion of Two Spacecraft In the previous section, equations of motion for an object in an inverse-squared gravity field were developed. While this expression is useful for describing the motion of a satellite around a central body, it can not directly describe the relative motion of two satellites in neighboring orbits. This section will develop three sets equations which directly describe the relative motion between two satellites. 2.2.1 Nonlinear Equations of Relative Motion Recall the nomenclature used to distinguish the satellites: one is often called the chief, and the other is referred to as the deputy. Note that the chief satellite is not necessarily 16 i jk Chief ^ ^ ^ r ? Deputy Figure 2.7: Local-vertical local-horizonal (LVLH) coordinate frame a physical object; in the case of a satellite formation it could be a useful reference point to describe the relative motion. The relative motion equations developed in this section utilize aCartesian local-vertical, local-horizontal(LVLH) frameattached to thechief satellite, as shown in Figure 2.7. This coordinate frame rotates with the chief?s radius vector and is a convenient reference frame to describe the relative motion.4 This reference frame is also sometimes referred to as the Hill frame or the CW frame. In this coordinate frame, x lies in the chief?s radial direction, z lies in the direction of the chief?s orbital angular momentum, and y completes the right-handed orthogonal triad. Note that the x and y directions correspond to the in-plane motion, and z corresponds to the out-of-plane motion. In the chief?s LVLH frame, the position of the deputy satellite is given by the following. rd =rc +?= (rc +x)?i+y?j+z?k (2.31) The angular velocity and acceleration of the LVLH frame are given by ? = ?fc?k = hr2 c ?k and ? = ?fc?k = ?2 ?rc ?fc rc ?k respectively. In the subsequent discussion, the subscript c will be omitted from r and f. From kinematics, the equation of motion for the deputy in the chief?s 17 frame is given by the following. ?rd = ?r+???+??(???) + 2?? ??+ ?? (2.32) Substitution and collecting terms in Equation (2.32) gives the following. ?rd = ?r?i+ ?f?k? parenleftBig x?i+y?j+z?k parenrightBig + ?f?k? parenleftBig ? f?k? parenleftBig x?i+y?j+z?k parenrightBigparenrightBig + 2 ?f?k? parenleftBig ?x?i+ ?y?j+ ?z?k parenrightBig + parenleftBig ?x?i+ ?y?j+ ?z?k parenrightBig (2.33) = ?r?i+ ?fx?j? ?fy?i? ?f2x?i? ?f2y?j+ 2?x ?f?j?2?y ?f?i+ ?x?i+ ?y?j+ ?z?k = parenleftBig ?r? ?fy? ?f2x?2?y ?f + ?x parenrightBig? i+ parenleftBig? fx? ?f2y+ 2?x ?f + ?y parenrightBig? j+ ?z?k Recall that the chief?s position in the LVLH frame is given by r = r?i. Substitution into Equation (2.5) allows the chief?s acceleration to be written as the following. ?r = ??r2 (2.34) Substituting Equation (2.34) into Equation (2.33) gives the following. ?rd = parenleftBig ??r2 ? ?fy? ?f2x?2?y ?f + ?x parenrightBig? i+ parenleftBig? fx? ?f2y+ 2?x ?f + ?y parenrightBig? j+ ?z?k (2.35) From Equation (2.5), the acceleration of the deputy is given by ?rd = ??r3 d rd where rd = [(r+x) y z]T. Equating coefficients in Equation (2.35) gives the full nonlinear equations of relative motion (NERM). ?x?2 ?f ?y? ?fy? ?f2x? ?r2 = ??r3 d (r+x) ?y+ 2 ?f ?x+ ?fx? ?f2y = ??r3 d y (2.36) ?z = ??r3 d z 18 The only assumption made in deriving Equation (2.36) was that both satellites obeyed Keplerian motion, i.e. the only force acting on each satellite was gravity. 2.2.2 Linear Equations of Relative Motion The nonlinear equations of motion presented in the previous section can be linearized for small separations between the chief and deputy. The deputy?s radius can be linearized as follows by neglecting higher order terms. rd = radicalBig (r+x)2 +y2 +z2 = r radicalbigg 1 + 2xr + x 2 +y2 +z2 r2 (2.37) ?r radicalbigg 1+ 2xr Using the binomial theorem, the term ?r3 d can be expanded in terms of r. Retaining only first order terms gives the following. ? r3d ? ?parenleftBig r radicalBig 1 + 2xr parenrightBig3 ? ?r3 parenleftbigg 1+ 2xr parenrightbigg?3 2 (2.38) ? ?r3 parenleftbigg 1? 3xr parenrightbigg Equation (2.38) is used to approximate the right hand side of the deputy?s acceleration given by Equation (2.5). ? ?r3 d rd = ??r3 d ? ?? ?? ? r+x y z ? ?? ?? ? ? ??r3 parenleftbigg 1? 3xr parenrightbigg ? ?? ?? ? r+x y z ? ?? ?? ? ? ??r3 ? ?? ?? ? r?2x y z ? ?? ?? ? (2.39) 19 Using Equation (2.39), Equation (2.36) can be linearized for small separations between the chief and deputy. These equations are here called the linearized equations of relative motion (LERM). ?x?2 ?f ?y? parenleftBig ? f2 + 2?r3 parenrightBig x? ?fy = 0 ?y+ 2 ?f ?x+ ?fx? parenleftBig ? f2 ? ?r3 parenrightBig y = 0 (2.40) ?z + ?r3z = 0 The system can be described by the state vector x = [x y z ?x ?y ?z]T. In state-space form, the LERM are given by the following. ?x= A(t)x= ? ?? ?? ?? ?? ?? ?? ?? ? 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1parenleftBig ?f2 + 2 ? r3 parenrightBig ? f 0 0 2 ?f 0 ??f parenleftBig ? f2 ? ?r3 parenrightBig 0 ?2 ?f 0 0 0 0 ? ?r3 0 0 0 ? ?? ?? ?? ?? ?? ?? ?? ? x (2.41) The state matrix A(t) is time varying, as r, ?f, and ?f vary with time. Significantly, for elliptic orbits these quantities are all time-periodic with the orbital period, T = P = 2?n . 2.2.3 Tschauner-Hempel Equations Solutions to Equation (2.40) can be found by first applying a coordinate scaling and a change of independent variable from time to true anomaly. The new state vector is ?x = [?x ?y ?z ?x? ?y? ?z?]T, where primes indicate derivatives with respect to true anomaly. The coordinate scaling is given by ?x = (1 +ecosf)x. Taking one derivative with respect to true anomaly gives ?x? = ?esinfx+ (1 +ecosf)x?, where x? = ?x?f = ?xr 2 h = p2 h(1 +ecosf). Generalizing this transformation for ?y, ?z, ?y?, and ?z? allows a transformation matrix between 20 x and ?x to be constructed as shown below. ?x=T (f)x (2.42) T (f) = ? ?? (1 +ecosf)I 0 ?esinfI p2h(1+ecosf)I ? ?? T?1 (f) = ? ?? 11+ecosfI 0 he p2 sinfI h p2 (1 +ecosf)I ? ?? (2.43) Here, I is a three-by-three identity matrix, and 0 is a three-by-three null matrix. Applying the transformation given by Equation (2.42) to Equation (2.41) results in the Tschauner- Hempel (TH) equations.7 The TH equations were first presented for relative motion by DeVries, who published an approximate solution in 1963.8 In fact, the inplane components of these equations were known and solved by Lawden in 1954.9 Lawden was not attempting to formulate equations for relative satellite motion, instead he was attempting to describe the primer vector for optimal spacecraft trajectory design. ?x?? ?2?y? ? 31 +ecosf ?x = 0 ?y?? + 2?x? = 0 (2.44) ?z?? + ?z = 0 In state-space form, the TH equations are given by the following. ?x? = B(f)?x= ? ?? ?? ?? ?? ?? ?? ?? ? 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 3 1 +ecosf 0 0 0 2 0 0 0 0 ?2 0 0 0 0 ?1 0 0 0 ? ?? ?? ?? ?? ?? ?? ?? ? x (2.45) 21 Because the TH equations use true anomaly as an independent variable, it is important to note that the true anomaly is not modulo 2? (in elliptic cases) but increases monotonically with time. The solution for the inplane states in TH equations reduces to solving a second-order, linear, nonhomogeneous, non-autonomous differential equation. The solution of this type of equation requires two homogeneous solutions and a particular solution, which were first presented by Tschauner and Hempel in 1965.7 This solution loses its independence at e = 0 and is undefined for e = 1. Yamanaka and Ankersen presented an alternative solution and showed how a recombination of terms eliminates the singularity at e = 0.10 Carter used the method of order reduction to eliminate both the singularities at e = 0 and e = 1, but this formulation is less compact than other solutions.11 A detailed review of the solutions to the TH equations can be found in Reference [12]. The specific form below consists of six fundamental solutions. For elliptic orbits, the fundamental solutions can be arranged in a matrix ?, and the vector of constants k = [k1 k2 k3 k4 k5 k6]T is defined. ?x= ?(f)k ?(f) = bracketleftbigg ?1 ?2 ?3 ?4 ?5 ?6 bracketrightbigg (2.46) ?1 = ? ?? ?? ?? ?? ?? ?? ?? ? sinf (1 +ecosf) 2cosf ?esin2f 0 cosf +ecos2f ?2sinf (1 +ecosf) 0 ? ?? ?? ?? ?? ?? ?? ?? ? ?2 = ? ?? ?? ?? ?? ?? ?? ?? ? cosf (1 +ecosf) ?2sinf ?esinf cosf 0 ?sinf ?esin2f e?2cosf (1 +ecosf) 0 ? ?? ?? ?? ?? ?? ?? ?? ? (2.47) 22 ?3 = ? ?? ?? ?? ?? ?? ?? ?? ? 1? 3e2 Ksinf (1+ecosf) ?32K(1 +ecosf)2 0 ?3e2 K(cosf +ecos2f)? 3esinf2(1+ecosf) 3eKsinf (1 +ecosf)? 32 0 ? ?? ?? ?? ?? ?? ?? ?? ? ?4 = ? ?? ?? ?? ?? ?? ?? ?? ? 0 1 0 0 0 0 ? ?? ?? ?? ?? ?? ?? ?? ? (2.48) ?5 = ? ?? ?? ?? ?? ?? ?? ?? ? 0 0 sinf 0 0 cosf ? ?? ?? ?? ?? ?? ?? ?? ? ?6 = ? ?? ?? ?? ?? ?? ?? ?? ? 0 0 cosf 0 0 ?sinf ? ?? ?? ?? ?? ?? ?? ?? ? (2.49) For elliptic orbits, K is related to the mean anomaly through the relationship K = M/(1? e2)3/2. The fundamental solutions correspond to four solutions for the in-plane motion and two solutions for the out-of-plane motion. General relative motions are composed of linear combinations of these six fundamental solutions in proportion to the constants k. These constants have units of length and are geometrically related to the size of the relative-motion contributions. The constants can also be considered integrals of Equation (2.40) and can be used to describe relative motion and plan maneuvers.13 By using Equation (2.43) and definingx(t0) =x0, the state-transition matrix for the original states, x, can be constructed as follows. k = ??1(f)T(f)x(t) (2.50) x(t) =T?1 (f)?(f)??1 (f0)T (f0)x0 = ?LERM(f,f0)x0 (2.51) 23 Note that time is the independent variable of the LERM, but the Tschauner-Hempel solution was used to construct the state transition matrix. Therefore, it is convenient to write the state transition matrix of the LERM, ?LERM(f,f0), in terms of true anomaly. 2.2.4 Hill-Clohessy-Wiltshire Equations Clohessy and Wiltshire published one of the most frequently used equations for relative satellite motion in 1960 when studying satellite rendezvous.14 These equations are frequently called the Hill-Clohessy-Wiltshire (HCW) equations, as Hill in 1878 was the first to linearize a set of equations to describe the motion of the moon relative to the Earth.15 Following the traditional derivation, assuming chief eccentricity equal to zero in Equation (2.40) results in the HCW equations. ?x?2n?y?3n2x = 0 ?y+ 2n?x = 0 ?z +n2z = 0 (2.52) Recall, n represents the mean motion. The HCW equations in state-space form are given by the following. ?x=Cx= ? ?? ?? ?? ?? ?? ?? ?? ? 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 3n2 0 0 0 2n 0 0 0 0 ?2n 0 0 0 0 ?n2 0 0 0 ? ?? ?? ?? ?? ?? ?? ?? ? x (2.53) The state matrix in this case is clearly time-independent. The state matrix has two eigen- values equal to zero that share a single independent eigenvector, and two pairs of imaginary eigenvalues both equal to ?in which have four independent eigenvectors. The state matrix 24 can be decomposed as shown below. C =MJM?1 (2.54) J = ? ?? ?? ?? ?? ?? ?? ?? ? 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ?n 0 0 0 0 n 0 0 0 0 0 0 0 0 ?n 0 0 0 0 n 0 ? ?? ?? ?? ?? ?? ?? ?? ? ; M = ? ?? ?? ?? ?? ?? ?? ?? ? 0 ?2n 0 2 0 0 3 3 4 0 0 0 0 0 0 0 ?1 0 0 0 2n 0 0 0 0 3 0 ?4n 0 0 0 0 0 0 0 n ? ?? ?? ?? ?? ?? ?? ?? ? (2.55) This decomposition is similar to the Jordan canonical form but instead of having the imag- inary eigenvalues on the diagonal in J, the imaginary magnitudes are placed in blocks on the diagonal. Associated with this, the imaginary and real parts of the eigenvectors appear separately in M. Note that the magnitude of each of the eigenvectors that make up M is arbitrary. However, the particular values chosen for M in Equation (2.55) will become important in later sections. Using this decomposition, a canonical form of the HCW equations is given by the following. w =M?1x ; ?w =Jw (2.56) This canonical form decomposes the system into three independent second-order subsystems. These subsystems are clearly related to a drifting in-plane (x and y directions) motion, a periodic in-plane motion, and a periodic out-of-plane (z direction) motion. Given an initial condition x0 at time t0, the HCW solution at time t can be expressed as x(t) = eC(t?t0)x0 = ?HCW(t,t0)x0. The solution of Equation (2.52) is given by the following.16 25 ?HCW(t,t0) = ? ?? ?? ?? ?? ?? ?? ?? ? 4?3cos(n(t?t0)) 0 0 6(sin(n(t?t0))?n(t?t0)) 1 0 0 0 cos(n(t?t0)) 3nsin(n(t?t0)) 0 0 6n(cos(n(t?t0))?1) 0 0 0 0 ?nsin(n(t?t0)) 1 n sin(n(t?t0)) 2 n (1?cos(n(t?t0))) 0 2 n (cos(n(t?t0))?1) 1 n (4sin(n(t?t0))?3n(t?t0)) 0 0 0 1n sin(n(t?t0)) cos(n(t?t0)) 2sin(n(t?t0)) 0 ?2sin(n(t?t0)) 4cos(n(t?t0))?3 0 0 0 cos(n(t?t0)) ? ?? ?? ?? ?? ?? ?? ?? ? (2.57) Approaches exist in the literature using the HCW equations as the basis for maneuver planning, see e.g. References [4, 17, 18, 19, 20, 21]. London extended the HCW equations for larger separation between the chief and deputy satellites in 1963.22 Like the HCW equations, London?s equations still assume a circular orbit, but the relative distance terms are retained to the second order. The solution of these equations is found through the method of successive approximations and includes secular and mixed-secular terms. 2.2.5 Geometric Parameterization of the HCW Equations An advantage of the HCW equations is the geometric insight into the solutions. Lovell and Tragesser found a geometric parametrization of the HCW equations called relative-orbit elements (ROEs).18 The ROEs allow relative motion described by the HCW equations to be described by geometric parameters similar to the concept of orbital elements for Keplerian 26 motion, which were discussed in Section 2.2.5. Lovell et al. have investigated using ROEs to extend the HCW equations to elliptic orbits and also for continuous-thrust maneuver- ing.23,24,25 The transformations between ROEs and Cartesian coordinates are given by the following. ae = 2 radicalBiggparenleftbigg ?x n parenrightbigg2 + parenleftbigg 3x+ 2 ?yn parenrightbigg2 x = ?ae2 cos? +xd xd = 4x+ 2 ?yn y = ae sin? +yd yd = y?2 ?xn z = zmax sin? ? = tan?1 parenleftbigg ?x 3nx+ 2?y parenrightbigg ?x = ae2 nsin? zmax = radicalBiggparenleftbigg ?z n parenrightbigg2 +z2 ?y = aencos?? 32nxd ? = tan?1 parenleftBignz ?z parenrightBig ?z = zmaxncos? (2.58) Here, the inverse tangent function is assumed to return a value between zero and 2?. An advantage of the ROEs is their simple behavior with respect to time. This behavior can be seen by substituting Equation (2.57) into Equation (2.58). ae = ae0 xd = xd0 yd = yd0 ? 32nxdt (2.59) ? = ?0 +nt (mod 2?) zmax = zmax0 ? = ?0 +nt (mod 2?) Combined, Equations (2.58) and (2.59) describe both the inplane and out-of-plane relative motion. As shown in Figure 2.8, the inplane motion traces an ellipse that has a semimajor axis ae in the y direction, a semiminor axis ae/2 in the x direction, centered at (xd, yd), and 27 x y Chief Satellite (xd, yd) Deputy Satellite ae ae Perigee Apogee Orbital Velocity Direction ?? Figure 2.8: In-plane relative-orbit elements. possibly drifts in the y direction. The out-of-plane motion oscillates about z = 0 with an amplitudezmax. Note that?? shown in Figure 2.8 is related to? through tan(?) = 2tan(??), and both angles are always in the same quadrant. In fact, ? is equal to the deputy?s mean anomaly. 2.2.6 Summary of Cartesian Linearized Relative-Motion Representations The previous sections presented three linearized descriptions of relative motion: the LERM, the TH equations, and the HCW equations. To summarize, the LERM are linear time-varying equations of motion which can describe the motion of a deputy satellite relative to a chief of any eccentricity. The assumptions used to derive the LERM was a chief and deputy satellite in close proximity, and both satellites obeying Keplerian motion. If an assumption is made that the chief is circular, the LERM reduce to the linear time-invariant HCW equations. Due to their simplicity, many solutions in the literature utilize the HCW equations. Applying a scaling and change of independent variable to the LERM results in the TH equations. The TH equations allow a state transition matrix of the LERM to be constructed. 28 LERM dx dt =A(t)xd79d79 ?x=T(f)x d15d15 eC=0 d47d47a95a95a95a95a95a95a95a95a95a95a95a95a95a95 HCWdx dt =Cx TH d?x df =B(f)?x Figure 2.9: Different Representations of Relative Motion Equations The relationships between the LERM, TH equations, and HCW equations are shown in Figure 2.9. Starting with the LERM, the assumption that the chief satellite is in a circular orbit is shown by a dashed line. One goal of this dissertation is to replace this assumption with a coordinate transformation. Note that the LERM and HCW equations share the same state vector. From the LERM, the transformation leading to the TH equations is shown by a solid line. Although an exact transformation exists, the state vector is now ?x. The following chapter discusses two approximate methods for applying the HCW equations to elliptic orbits. These methods result in time-varying equations of motion, but as will be seen, are lower fidelity than the LERM. 2.2.7 Additional Relative-Motion Representations The relative motion models presented above are not exclusive. A linearized set of equa- tions similar to the TH equations was developed by Brumberg and later expanded on by Kelly.26,27 The Brumberg-Kelly equations differ from the TH equations in that the inde- pendent variable is not transformed from time to true anomaly and that non-orthogonal coordinate axes are used. Kelly provides a state transition matrix in terms of time and eccentric anomaly, which reduces to a form similar to the HCW equations for circular orbits. However, these equations are not as compact as the TH equations and suffer from singularity issues. Representations for relative motion can also be performed with respect to various 29 parameter sets which are associated with an orbit. An element set defined using inclination and eccentricity vectors are attractive for defining the orbit of geostationary satellites and insuring safe collocation.28 Jones used this element set to describe rendezvous in circular orbits, while D?Amico and Montenbruck also used these vectors to describe relative motion and presented a case study for the TerraSAR-X/TanDEM-X formation.29,30 The relative motion can also be described through orbit element differences. In this method, the dif- ference between any two complete sets of orbit elements can be determined. Using orbit element differences, the parameterization of the relative motion at any instant of time can be determined by solving Kepler?s equation instead of a differential equation.31,32 2.2.8 Calculation of Initial Conditions Propagation of any of the equations of motion discussed in this section requires a set of initial conditions for the relative motion. The position and velocity of the deputy relative to the chief are related to the inertial position and velocity of each satellite through the following kinematic relationships. [?]C = [R]([rD]I ?[rC]I) [ ??]C = [R]([vD]I ?[vC]I)?[?]C [?]C (2.60) Here, inertial position vectors of the chief and deputy satellites are given by rC and rD respectively. The inertial velocities of the chief and deputy with respect an inertial frame are vC and vD. The skew-symmetric matrix of the angular velocity of the chief?s LVLH frame is given by [?], and the rotation matrix from inertial into the chief?s LVLH frame is given by [R]. The matrix operator [ ]I indicates coordinatization in an inertial frame, and [ ]C indicates coordinatization in the chief?s LVLH frame. In this dissertation, initial conditions for the chief and deputy will be specified using orbital elements. Orbital elements can be transformed into inertial position and velocity vectors (see, e.g. Reference [31] ?9.6.2), and using Equation (2.60), the relative position and velocity are calculated. As an example, 30 consider a chief satellite orbiting the Earth with the following orbital elements: aC = 8000 km, eC = 0.1, and iC = ?C = ?C = fC = 0. Additionally, a deputy satellite is specified in close proximity with the orbital elements: aD = 8000 km, eD = 0.10001, and iD = ?D = ?D = fD = 0. For this case, ?= [?0.080000 0 0]T km and ??= [0 0.0001655329 0]T km/s. The concept of calculating initial conditions will be revisited later in this dissertation. Specifically, Section 3.1 develops a virtual-chief method which propagates the motion in a coordinate frame which is more consistent with the assumptions of the HCW equations. Later, Section 5.1.3 compares the LERM to the HCW equations showing how the computed Lyapunov-Floquet transformations completely capture the dynamics, while in Section 5.3.2, the Lyapunov-Floquet transformations are used to select initial conditions which allows the underlying HCW equations to approximate the LERM solution. 31 Chapter 3 Preliminary Investigations? The assumptions made in deriving the HCW equations were a chief satellite in a circular orbit, a deputy satellite in close proximity to the chief, and both satellites obeying Keplerian motion. The assumption of a chief satellite with zero eccentricity is rarely achieved in practice. As the chief?s eccentricity increases, the HCW equations become less accurate. Alternative models exist in the literature, such as the LERM and TH equations discussed in the previous chapter, which are valid for arbitrary values of chief eccentricity. However, the analytical solutions of these alternative equations can be less intuitive to understand than the solutions of the HCW equations. This chapter presents two investigations into applying the HCW equations to elliptic chiefs. The Virtual-Chief (VC) model propagates the relative motion in a frame more consis- tent with the assumptions of the HCW equations. The Virtual-Time (VT) method evaluates the HCW equations at a virtual time to account for the time-varying effects of eccentric or- bital motion. Both of these methods define (somewhat ad hoc) transformations to the HCW solutions, and essentially construct approximate linear time-varying models for the relative motion. A comparison of both methods is presented in Section 3.3. ?Material in this chapter taken from the following papers: [33] R. E. Sherrill, A. J. Sinclair, T. A. Lovell, K. W. Johnson, and D. D. Decker, ?The Virtual-Chief Method for Modeling Relative Motion of Noncircular Satellites,? AAS 11-214, presented at the 21st AAS/AIAA Space Flight Mechanics Meeting, New Orleans, Louisiana, 2011. [34] R. E. Sherrill, A. J. Sinclair, and T. A. Lovell, ?A Virtual-Time Method for Modeling Rela- tive Motion of Noncircular Satellites,? AAS 11-208, presented at the 21st AAS/AIAA Space Flight Mechanics Meeting, New Orleans, Louisiana, 2011. 32 3.1 Virtual-Chief Method In the virtual-chief (VC) method, a fictional satellite with zero eccentricity is used as the chief satellite for the HCW equations, with both the actual chief and deputy satellites treated as deputies. This virtual chief satellite is a circularized version of the actual chief satellite, and is defined by setting its eccentricity to zero, and all other orbital element values equal to those of the chief.23 aVC = aC eVC = 0 iVC = iC ?V C = ?C ?VC = ?C fVC = MVC = MC (3.1) Here, a subscript VC refers to the virtual chief satellite. The above definitions for ?VC and fVC are somewhat arbitrary, because the virtual chief?s orbit is circular. Note that the virtual chief?s LVLH frame has constant angular velocity, unlike the chief?s LVLH frame. The three satellites are illustrated in Figure 3.1a. The VC method can also be called a dual-deputy approach, as both the real chief and deputy satellites are treated as deputies to the virtual chief. Because the circular orbit of the virtual chief has been defined, the motion of both the chief and deputy relative to the virtual chief can be described using the HCW equations without violating the circularity assumption. The motion of the deputy relative to the real chief is the difference of these solutions. However, the HCW equations do not need to be propagated twice in order to use the VC method. Since the HCW equations are linear, the difference of two solutions is also a solution to the equations. The state vector for the VC method is defined by u = [[?]TVC [??]TVC ]T, where [ ]VC indicates coordinatization in the virtual chief?s frame, and (?) indicates vector differentiation with respect to the virtual chief?s frame. The VC solution is propagated using 33 Earth Virtual Chief Chief Deputy Hill LVLH ? fVC = MC fC (a) Coordinate transformation between virtual-chief and chief frames. Earth Chief Deputy Hill LVLH ? fC fC-MC (b) Virtual-chief, chief, and deputy satellites. Figure 3.1: Virtual-Chief Method the HCW equations. u(t) = eCtu0 (3.2) Since the relative motion is propagated in the virtual chief?s frame, a coordinate trans- formation is necessary to convert the solution into the actual chief?s frame. For the position components, this coordinate transformation is a rotation in the chief?s orbital plane by the difference between the chief?s mean anomaly and the chief?s true anomaly, fC ?MC. Notice that only the in-plane motion is transformed since the chief and virtual chief are coplanar. For the velocity components, the angular velocity of the virtual chief?s frame relative to the chief?s frame is given by ?1 = ?2 = 0, and ?3 = n? ?fC (components are identical in the chief?s and virtual chief?s frame). Given an initial condition, the VC solution xVC =PVCu 34 is given by the following. R(t) = ? ?? ?? ? cos(fC ?MC) sin(fC ?MC) 0 ?sin(fC ?MC) cos(fC ?MC) 0 0 0 1 ? ?? ?? ? ; ?VC/C(t) = ? ?? ?? ? 0 ?fC ?n 0 n? ?fC 0 0 0 0 0 ? ?? ?? ? (3.3) P?1VC(t) = ? ?? R(t)T 0 ??VC/C(t)R(t)T R(t)T ? ?? (3.4) xVC =PVC(t)eC(t?t0)P?1VC(t0)x0 (3.5) Note that evaluating the VC solution does not require propagating the virtual chief?s position, and only requires knowledge of the orientation of the virtual chief?s frame. This suggests an alternative perspective for the VC method. The VC method can be thought of as propagating the relative motion in a constant-angular-velocity frame attached to the chief, labeled in Figure 3.1b as the Hill frame, which is more consistent with the HCW assumptions. In Reference 23, a set of governing equations for the VC method were presented. These equations are linear time-periodic differential equations. Equations (3.2-3.5) are a solution of those equations by Lyapunov-Floquet transformation (see e.g. References 35 and 36). The VC method applies a time-varying coordinate transformation to the solution of time- invariant ordinary differential equations. As will be seen, however, the VC equations are lower fidelity than the LERM. The next section describes the virtual-time method, which is another method for extending the HCW equations to eccentric chiefs. 3.2 Virtual-Time Method Consider a chief and deputy satellite orbiting the Earth, sharing the following orbital elements of a = 11,000 km and i = ? = ? = 0. Both satellites are initially at perigee. Figure 3.2 shows the NERM trajectory for a range of chief eccentricities, but in each case, 35 ?0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 y (km) x (km) Figure 3.2: NERM trajectories for chief eccentricities of 0, 0.01, 0.1, 0.2, and 0.4. the chief and deputy eccentricities are related by eD = eC + 10?5. In these cases, the trajectory shapes display only minor displacement from a 2?1 ellipse. In the above example it was seen that even for moderate chief eccentricity, the shape of the relative-motion trajectory is similar to the trajectory shape predicted by the HCW equations. This suggests that in some cases the error in applying the HCW equations to noncircular chiefs has largely to do with the motion along the trajectory, instead of the shape of the trajectory. This motivates an approach to describe the relative-motion state at time t, by evaluating the HCW solution at some virtual time ?. Define a state vector w = [[?]TC ?T ]T consisting of the virtual-time (VT) solution for the relative position, ?, and a set of velocity-like variables, ?. w(t) = eC(?(t)?t0)x0 (3.6) This solution can be evaluated at t0, and x0 can be eliminated. w(t) = eC(?(t)??(t0))w0 (3.7) 36 This VT solution forwpasses through the same points in state space as the HCW solution for x, and the trajectory shape for the VT solution is identical to the HCW solution. However, in the VT approach, each point along the trajectory is reached at a different time than in the HCW solution. Solving for the virtual time ?(t) as a function of the true time can be posed as an optimization problem over an interval 0 < t < T. The error in the virtual-time solution is defined as the difference from the NERM solution. J = integraldisplay T 0 (?V (t)??N(t))T (?V (t)??N(t))dt (3.8) = integraldisplay T 0 (?H(?(t))??N(t))T (?H(?(t))??N(t))dt Here, the chosen objective function focuses on the relative positions from the virtual-time solution, ?V , and the NERM solution, ?N. The value of ?(t) that minimizes Equation (3.8) is the optimal virtual time and will be denoted ??. This optimization problem can be solved numerically by an exhaustive search comparing the NERM and HCW solution at discrete time steps. Returning to the example motions shown in Figure 3.2, the optimal solution for the virtual time in each case is shown in Figure 3.3, with a time step of 15 seconds. For small chief eccentricities, ?? is approximately equal to t; but as the chief eccentricity increases, ?? begins to depart significantly from t. In each plot of Figure 3.3, two heuristic values of the virtual time are also shown. In these examples, the optimal virtual time is larger than the true time as the chief travels from periapse to apoapse, and the optimal virtual time is smaller than the true time as the chief travels from apoapse to periapse. This is the same behavior that eccentric anomaly and true anomaly exhibit relative to the mean anomaly. This motivates two heuristic values for the virtual time, ?E and ?f, that are related to the eccentric anomaly and true anomaly, respectively. n(t+tp) = n?E ?eC sinn?E (3.9) 37 0 2000 4000 6000 8000 100000 2000 4000 6000 8000 10000 t (s) ? (s) ?* ?E ?f 0 2000 4000 6000 8000 100000 2000 4000 6000 8000 10000 t (s) ? (s) ?* ?E ?f (a) eC = 0.01 (b) eC = 0.1 0 2000 4000 6000 8000 100000 2000 4000 6000 8000 10000 t (s) ? (s) ?* ?E ?f 0 2000 4000 6000 8000 100000 2000 4000 6000 8000 10000 t (s) ? (s) ?* ?E ?f (c) eC = 0.2 (d) eC = 0.4 Figure 3.3: Solutions for the virtual time for several chief eccentricities. tann?f2 = radicalbigg1 +e C 1?eC tan n?E 2 (3.10) Here, tp is the amount of time at t = 0 since the chief?s last periapse passage. If t = 0 occurs at the chief?s periapse, then at t = 0 all three solutions for the virtual time equal zero: ?(0) = 0. But this is not true for general initial locations, and in general w0 negationslash= x0. Solving for ?E (or ?f) at a particular instant t requires numerical solution of Equation (3.9). For completeness, the derivatives of ?E and ?f can be also be calculated. ??E = 11?e C cosn?E (3.11) The derivative of ?f is calculated by analogy with the true anomaly.37 ??f = (1 +eC cosf) 2 (1?eC)3/2 (3.12) 38 For the examples shown in Figure 3.3, ?f is a good approximation of the optimal virtual time, particularly at low chief eccentricities, and will be used to approximate the virtual time in the following subsection. It is possible to investigate a set of governing equations that admit the VT solution, by differentiating Equation (3.7) with respect to time. ?w =C??eC(?(t)??(0))w0 = ??Cw (3.13) The first three of these equations relate the velocity-like variables ? to the VT solution for the relative velocity. [ ??]C = ??? (3.14) These equations are a linear time-varying system where the state matrix varies with respect to time in a very specific manner: the time-varying state matrix is equal to a time-invariant matrix multiplied by a time-varying scalar. These equations can be compared to higher- fidelity linearized time-varying models where the state matrix varies with respect to time in more complicated manners, see References 10 and 38. The kinematics in Equation (3.14) can be further analyzed. ? = 1??[ ??]C ?? = 1??[??]C ? ????2[ ??]C (3.15) Substituting these relations into Equation (3.13) gives the virtual-time equations in second order form. ?x? ???? ?x?2n?? ?y?3n2 ??2x = 0 ?y+ 2n?? ?x? ???? ?y = 0 (3.16) ?z? ???? ?z +n2 ??2z = 0 39 Finally, the scalar form of the VT solution for the relative position can be written by eval- uating Equation (2.57) at the virtual time and taking into account Equation (3.15). For compactness, ?t?0 ??(t)??(0) and ??0 ? ??(0) are defined. x(t) = (4?3cosn?t?0)x0 + sinn?t?0n?? 0 ?x0 + 2n?? 0 (1?cosn?t?0) ?y0 y(t) = 6(sinn?t?0 ?n?t?0)x0 +y0 + 2n?? 0 (cosn?t?0 ?1) ?x0 + 1n?? 0 (4sinn?t?0 ?3n?t?0) ?y0 (3.17) z(t) = cosn?t?0 z0 + sinn?t?0n?? 0 ?z0 Equation (3.17) shows that by construction, the solutions of the linear time-varying system in Equation (3.16) have the same shapes as the linear time-invariant HCW equations but pass through these points at different times. 3.3 Results and Discussion To evaluate the performance of the virtual-chief and the virtual-time approaches relative to the HCW equations and the LERM, the motion of each method can be compared to the NERM. Here, initial conditions for the chief and deputy will be specified in terms of orbital elements. These orbital elements are transformed into the true initial conditions, x0, using Equation (2.60). These initial conditions were used for the NERM, the LERM, and the HCW equations. Recall that the VC method applies a coordinate transformation to the initial conditions before propagating the motion with the HCW equations: u0 = P?1VC(t0)x0. Here, the VT method will be evaluated using the same initial conditions: w0 = eC(?(t0)?t0)P?1VC(t0)x0 For each solution, the motion is propagated for one complete revolution of the chief?s orbit and at each time step, the position-error magnitude is computed. The root-mean-square 40 Case aC (km) eC aD (km) eD ?D (rad) 1 11,000 0.1 11,000 0.10001 0 2 11,000 0.4 11,000 0.40001 0 3 11,000 0.1 11,000.2 0.10001 0 4 11,000 0.4 11,000.2 0.40001 0 5 11,000 0.1 11,000 0.10001 2?10?5 6 11,000 0.4 11,000 0.40001 2?10?5 Table 3.1: Initial conditions for inplane test cases. (RMS) of these errors, ?, is given by the following. ?RMS = radicaltpradicalvertex radicalvertexradicalbt1 m msummationdisplay l=1 (?l ??N,l)T(?l ??N,l) (3.18) Here, ?l indicates the position vector at the lth time step obtained by propagating the LERM, HCW, VC, or VT method, ?N,l indicates the position vector at the lth time step obtained by propagating the NERM, and m represents the total number of time-steps in each solution. Six cases in particular were chosen to compare the performance of the LERM, HCW, VC, and VT methods. Table 3.1 gives the chief and deputy orbital elements for each case. Orbital elements not listed are set to zero. The RMS error was calculated for each case and is presented in Table 3.2. The RMS error indicates that both the VC and VT methods had considerably less error that the HCW solution for every case. For every case except cases 3 and 4, the VT method has less error than the VC method, although the improvement was not as significant as VC compared to HCW. For all cases, the LERM offered four orders of magnitude or more greater accuracy. The methods in this chapter selected transformations the the HCW equations which defined new approximate linear time-varying models. Whereas these models provide some accuracy improvement to the HCW equations, they are still lower fidelity than the LERM. 41 Case HCW Error (km) VC Error (km) VT Error (km) LERM Error (km) 1 0.4714 0.1625 0.1581 1.0460?10?5 2 3.2406 1.1377 1.0906 4.2539?10?5 3 0.4409 0.2559 0.2588 8.5585?10?5 4 0.8417 0.6887 0.6959 1.2905?10?4 5 0.4893 0.1294 0.1270 5.8095?10?5 6 3.3216 0.9411 0.8679 7.7002?10?5 Table 3.2: RMS error distances for inplane test cases. Therefore, Chapter 5 will investigate the opposite approach. Instead of defining a linear time- varying model by selecting a transformation, the transformation will be directly calculated by choosing the desired linear time-varying model. 42 Chapter 4 Lyapunov-Floquet Theory In 1883, Floquet published a theory on the solutions of linear systems with periodic coefficients.39 Floquet theory showed that once a fundamental solution to a linear periodic system over the principal period T has been found, it can be extended to all time. Addi- tionally, the solution at some time T + t is given by the product of the solution at time T and the solution at time t. In 1949, Lyapunov extended Floquet?s work by developing a transformation that converts a linear periodic system to a time-invariant form. 4.1 Floquet Theorem Consider a linear time-periodic system of the form, ?x=A(t)x, whereA(t) isT periodic such thatA(t+T) =A(t). Let?(t,t0) bethestate-transition matrixforthesystem. Because of the time periodicity of the system, clearly ?(t,t0) = ?(t+T,t0 +T). Additionally, if one considers the mapping from t0 to some time t+T, the state-transition matrix can be broken up arbitrarily. ?(t+T,t0) = ?(t+T,t0 +T)?(t0 +T,t0) = ?(t,t0)?(t0 +T,t0) (4.1) The constant matrix ?(t0 + T,t0) is often referred to as the discrete transition matrix, monodromy matrix, or the Floquet transition matrix. It is convenient to write this matrix as follows. ?(t0 +T,t0) =P0e?TP?10 (4.2) 43 Here, ? is a constant matrix, and P0 is also a constant, nonsingular matrix. Equation (4.1) can now be rewritten using Equation (4.2). ?(t+T,t0) = ?(t,t0)P0e?TP?10 (4.3) Next, define P(t) as follows. P(t) = ?(t,t0)P0e??(t?t0) (4.4) Evaluating Equation (4.4) at t0 gives P(t0) = P0. Using Eqs. (4.1) and (4.2), P(t) is seen to be T periodic. P(t+T) = ?(t+T,t0)P0e??(t+T?t0) = ?(t,t0)P0e?TP?10 P0e??Te??(t?t0) = ?(t,t0)P0e??(t?t0) (4.5) =P(t) Since ?(t,t0) is nonsingular by definition, P(t) is also nonsingular.40 Therefore, inverting Equation (4.4) leads to a formal definition of Floquet Theory. IfA(t) is a continuous, T-periodic matrix, then for all t ? R the state-transition matrix for ?x=A(t)x can be written in the form ?(t,t0) =P(t)e?(t?t0)P?1(t0) (4.6) where P(t) is a nonsingular, differentiable, T-periodic matrix and ? is a constant matrix.41 44 4.1.1 Discussion of Floquet Theorem Note that Floquet Theory only guarantees the existence of P(t) and ?. In general, there is no easy way to calculate ?(t,t0) for an arbitrary time-periodic system. Generally, a numeric algorithm is used to integrate the state equation over the principal period. See the work of Sinha et al. (e.g. References 42,43,44,45,46,47,48 and the references contained within) for more information on numerical methods. The previous paragraphs focused on evaluating the state transition matrix at the end of one period, T, and in general, both P(t) and ? will be complex. It is possible to instead evaluate the state transition matrix at the end of two periods, 2T. The resulting matrices, L(t) and ? (analogous to P(t) and ? respectively) will now be real. 4.1.2 System Stability The characteristic exponents, ?i, are defined as the eigenvalues of ?. The characteristic exponents can also be referred to as Floquet exponents. The characteristic exponents are related to the characteristic multipliers, ?i, by ?i = 1T ln|?i|. Note that the characteristic multipliers are the eigenvalues of ?(t0 + T,t0).36 System stability can be determined by examining the characteristic exponents or the multipliers. A necessary condition for the existence ofT-periodic solutions is that one or more of the characteristic exponents are purely imaginary (multiplier has modulus 1). A necessary and sufficient condition for asymptotic stability is that all of the exponents have negative real parts (multipliers have modulus less than one).49 4.2 Lyapunov-Floquet Theory Lyapunov extended Floquet Theory by introducing a change of coordinates that reduces the linear time-periodic system to a system of constant coefficients. Under the following periodic change of coordinates. x=P(t)z (4.7) 45 The system is reduced to the linear time-invariant form. ?z = ?z (4.8) This decomposition implies that the solution forxis a time-periodic coordinate transforma- tion of the solution to an underlying linear time-invariant system. The implication of this transformation is that the study of a time-periodic system can be reduced to the study of a much simpler time-invariant system. 4.3 Determining a Lyapunov-Floquet Transformation While Lyapunov-Floquet theory is a powerful tool in studying linear time-periodic equa- tions, it does not provide a framework for calculating P(t) and ?. This section will discuss two methods for calculating a Lyapunov-Floquet transformation for commutative and non- commutative systems. 4.3.1 Commutative System For a certain class of systems called commutative systems, the state-transition matrix and the Lyapunov-Floquet matrices can easily be found in closed form.50 A periodic system matrix A(t) is called commutative if there exists a matrix B(t), where A(t) and B(t) have the following relationship. dB(t) dt =A(t) (4.9) The matrix B(t) must satisfy the following properties. A(t)B(t) =B(t)A(t) B(T)B(0) =B(0)B(T) (4.10) 46 Here, B(t) is called the commutating antiderivative of A(t). The state-transition matrix is found through the matrix exponential. ?(t,t0) = eB(t)e?B(t0) (4.11) The term eB(t) from Equation (4.11) can be factored by defining another matrix BP(t). eB(t) = eBP(t)et? (4.12) The constant matrix ? is found by evaluating B(t) at two points. ? = 1T integraldisplay t0+T t0 A(?)d? = 1T (B(T +t0)?B(t0)) (4.13) From Equation (4.12), the matrix BP(t) can be determined by the following. BP(t) =B(t)?t? (4.14) TheT-periodic transformation matrixP(t) is found through the matrix exponential of Equa- tion (4.14). P(t) = eBP(t) (4.15) Equations (4.13) and (4.15) allow the state-transition matrix given by Equation (4.11) to be written as ?(t,t0) =P(t)e?(t?t0)P?1(t0). As an example, consider the periodic state matrix A(t) shown below. A(t) = ? ?? ?? ? 0 sint 1 sint 1 sint 1 sint 0 ? ?? ?? ? (4.16) 47 Here, T = 2?. The commutating antiderivative matrix B(t) can be determined which satisfies the conditions given by Equation (4.10). B(t) = ? ?? ?? ? 0 ?cost t ?cost t ?cost t ?cost 0 ? ?? ?? ? (4.17) Using Equation (4.17), the state-transition matrix for Equation (4.16) can be calculated from Equation (4.11). ?(t,t0) = ? ?? ?? ? 1 2 cosh(t?t0)(1 +a1) + 1 2 sinh(t?t0)(a1 ?1) ? a2? 2e t?t0 ?a2?2et?t0 et?t0a1 1 2 cosh(t?t0)(a1 ?1) + 1 2 sinh(t?t0)(1 +a1) ? a2? 2e t?t0 1 2 cosh(t?t0)(a1 ?1) + 1 2 sinh(t?t0)(1 +a1) ?a2?2et?t0 1 2 cosh(t?t0)(1 +a1) + 1 2 sinh(t?t0)(a1 ?1) ? ?? ?? ? (4.18) Here, a1 = coshparenleftbig?2(cost?cost0)parenrightbigand a2 = sinhparenleftbig?2(cost?cost0)parenrightbig. The matrices ? and BP(t) can be determined from Equations (4.13) and (4.14) respectively. ? = 12? (B(2?)?B(0)) = ? ?? ?? ? 0 0 1 0 1 0 1 0 0 ? ?? ?? ? (4.19) BP(t) =B?t? = ? ?? ?? ? 0 ?cost 0 ?cost 0 ?cost 0 ?cost 0 ? ?? ?? ? (4.20) 48 Using Equation (4.20), the periodic transformation matrixP(t) can be found through Equa- tion (4.15). P(t) = eBP(t) = ? ?? ?? ? cosh parenleftBig cost? 2 parenrightBig2 ?sinh( ?2cost) ?2 sinh parenleftBig cost? 2 parenrightBig2 ?sinh( ?2cost) ?2 cosh parenleftbig?2costparenrightbig ?sinh(?2cost) ?2 sinh parenleftBig cost? 2 parenrightBig2 ?sinh( ?2cost) ?2 cosh parenleftBig cost? 2 parenrightBig2 ? ?? ?? ? (4.21) 4.3.2 Non-commutative System For non-commutative systems, a Lyapunov-Floquet transformation can be determined from the state-transition matrix of a time-periodic system as discussed in Section 4.1. It is commonly assumed that P(t0) = I. Using this assumption, the constant matrix ? can be found from the monodromy matrix given by Equation (4.2). ? = 12? ln(?(t0 +T,t0)) (4.22) Once ? has been determined, P(t) is determined from Equation (4.4). Consider the previous example where the state matrix was given by Equation (4.16) and the state-transition matrix for this system was given by Equation (4.18). Evaluating the state-transition matrix over one period gives the following. ?(2?,0) = ? ?? ?? ? cosh2? 0 sinh2? 0 e2pi 0 sinh2? 0 cosh2? ? ?? ?? ? = ? ?? ?? ? 1 2 (e ?2pi +e2pi) 0 1 2 (?e ?2pi +e2pi) 0 e2pi 0 1 2 (?e ?2pi +e2pi) 0 1 2 (e ?2pi +e2pi) ? ?? ?? ? (4.23) The matrix ?(2?,0) can be written in canonical form as shown below. ?(2?,0) =MJM?1 49 J = ? ?? ?? ? e2pi 0 0 0 e2pi 0 0 0 e?2pi ? ?? ?? ? ; M = ? ?? ?? ? 1 0 ?1 0 1 0 1 0 1 ? ?? ?? ? (4.24) From Equation (4.22), the matrix ? can be determined. ? = 12? ln(?(t0 +T,t0)) = 12? lnparenleftbigMJM?1parenrightbig= 12?Mln(J)M?1 = ? ?? ?? ? 0 0 1 0 1 0 1 0 0 ? ?? ?? ? (4.25) From Equation (4.4), P(t) is determined and is shown below. P(t) = ? ?? ?? ? cosh parenleftBig cost? 2 parenrightBig2 ?sinh( ?2cost) ?2 sinh parenleftBig cost? 2 parenrightBig2 ?sinh( ?2cost) ?2 cosh parenleftbig?2costparenrightbig ?sinh(?2cost) ?2 sinh parenleftBig cost? 2 parenrightBig2 ?sinh( ?2cost) ?2 cosh parenleftBig cost? 2 parenrightBig2 ? ?? ?? ? (4.26) Note that Equations (4.25) and (4.26) give the same result as Equations (4.19) and (4.21) obtained in the previous subsection. 4.3.3 Discussion In Chapter 5, two Lyapunov-Floquet transformationswill bedetermined which relate the LERM to the HCW equations. Recall that the time-invariant HCW equations are a special case of the time-varying LERM when the chief satellite is in a circular orbit. In general, choosing a time-invariant special case of a time-varying system as ? does not produce a valid Lyapunov-Floquet transformation. Consider the Mathieu Equation without damping given by ?x+(a+?cos2t)x = 0. In state-space form, the Mathieu Equation is given by the following. ? ?? ?x ?x ? ??= ? ?? 0 1 ?(a+?cos2t) 0 ? ?? ? ?? x ?x ? ?? (4.27) 50 The following time-invariant form of Equation (4.27) is found by setting ? = 0. ? ?? ??x ??x ? ??= ? ?? 0 1 ?a 0 ? ?? ? ?? ?x ??x ? ?? (4.28) For constructing a Lyapunov-Floquet transformation of the Mathieu equation, the time- invariant system calculated in Equation (4.28) does not correspond with the time-invariant matrix ? calculated by Sinha et al. through perturbation methods or through expansion by Chebyshev polynomials.51 This indicates that in general, ? must be calculated using the time-periodic state matrix. Further, since the Lyapunov-Floquet transformation can- not change stability characteristics (being periodic), Equations (4.27) and (4.28) must have identical stability properties; which is obvious not the case. 4.4 Lyapunov-Floquet Transformation of the TH Equations A Lyapunov-Floquet (LF) transformation of the TH equations has previously been determined.7,52,53 The review of LF transformations in the previous sections generically used time for the independent variable. For the TH equations, of course, the independent variable is true anomaly. ?x=P(f)z ; z? = ?z (4.29) ? = ? ?? ?? ?? ?? ?? ?? ?? ? 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ?1 0 0 0 0 0 0 0 0 1 0 0 0 0 ?1 0 ? ?? ?? ?? ?? ?? ?? ?? ? (4.30) 51 P?1(f) = ? ?? ?? ?? ?? ?? ?? ?? ? c p1 13 ?q2 0 0 ?2q1 +e?? ?e? 0 ?q1 0 0 ?12esinf ?12(1 +ecosf) 12ecosf 0 0 0 ?12(3+ecosf) 0 ?12esinf ?12(2 +ecosf) 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? (4.31) where c = 1e bracketleftBig 1?parenleftbig1 + 2e2parenrightbig?1?e2 bracketrightBig sinf ?parenleftbig2 + 3ecosf +e2parenrightbigsin?1? p1 = ?16 bracketleftBig 1 + 3?1?e2 bracketrightBig ? 13e bracketleftBig 1?parenleftbig1?e2parenrightbig3/2 bracketrightBig cosf + 16 bracketleftBigparenleftbig 1+ 2e2parenrightbig?1?e2 ?1 bracketrightBig cos2f ?e?sin?1? q1 = (1 +ecosf)2 q2 = 13e bracketleftBigparenleftbig 2 +e2parenrightbig?1?e2 ?2 bracketrightBig sinf + 16 bracketleftBigparenleftbig 1 + 2e2parenrightbig?1?e2 ?1 bracketrightBig sin2f + (1 +ecosf)2 sin?1? ? = sinf (1+ecosf) ? = sinf1 +ecosf bracketleftBig e+ parenleftBig 1??1?e2 parenrightBig cosf bracketrightBig For e = 0, the elements c, p1, and q2 in Equation (4.31) must be calculated using L?Hospital?s Rule. For this case, P?1(f) is given by the following. P?1|e=0 = ? ?? ?? ?? ?? ?? ?? ?? ? 0 ?23 13 0 0 0 ?2 0 0 ?1 0 0 0 ?12 0 0 0 0 ?32 0 0 ?1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? (4.32) 52 Note that when deriving Equations (4.30) and (4.31), the state vector x = [?x ?x? ?y ?y? ?z ?z?]T was used. To remain consistent with the state vector used in this dissertation, a transfor- mation matrix F is required. F = ? ?? ?? ?? ?? ?? ?? ?? ? 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? (4.33) Using this LF transformation, the LERM solution at time t can be written as the following. x(t) =T(f)?1F?1P(f)e?(f?f0)P(f0)?1FT(f0)x(0) (4.34) 4.5 Lyapunov-Floquet Generalization of the HCW Equations in the f-Domain? The constant matrix?in Equation (4.30) strongly resemblesJ given in Equation (2.55). Due to this similarity, a transformation between ? and J should exist which would relate the HCW equations to the TH equations (and also the LERM). In order to relate ? and J, first define a virtual-time ? such that ? = 1nf and ddf = 1n dd?. Therefore, the derivative z? = ?z in Equation (4.29) can be rewritten in terms of ? and is given by dd?(z) = n?z. ?Material in this section taken from the following paper: [54] R. E. Sherrill, A. J. Sinclair, and T. A. Lovell, ?A Lyapunov-Floquet Generalization of the Hill- Clohessy-Wiltshire Equations,? AAS 12-103, presented at the 22nd AAS/AIAA Space Flight Mechanics Meeting, Charleston, SC, 2012. 53 The matrix n? can be written in canonical form n? =SJS?1 where S is given by S = ? ?? ?? ?? ?? ?? ?? ?? ? 0 1 0 0 0 0 1 n 0 0 0 0 0 0 0 ?1 0 0 0 0 0 0 1 0 0 0 0 0 0 ?1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? (4.35) Note that the magnitude of the eigenvectors that compriseS is somewhat arbitrary. Solving n? =SJS?1 for J and substituting for J in C =MJM?1 motivates another coordinate transformation, R=MS?1. y =Rz = ? ?? ?? ?? ?? ?? ?? ?? ? 0 ?2 0 2 0 0 3 0 ?4 0 0 0 0 0 0 0 1 0 0 0 ?2n 0 0 0 0 3n 0 ?4n 0 0 0 0 0 0 0 n ? ?? ?? ?? ?? ?? ?? ?? ? z (4.36) d d?(y) =R d d?(z) = nR?z = nR?R ?1y =MJM?1y =Cy (4.37) The virtual-time can be calculated from the true time through Kepler?s Equation.34 In order to derive the HCW equations presented earlier, an assumption was made that the chief satellite was in a circular orbit. However, Equation (4.37) shows that the same equations apply to any elliptic orbit through the appropriate coordinate and independent variable transformations. Figure (4.1) shows the different representations of the relative 54 LERM dx dt =A(t)xd79d79 ?x=T(f)x d15d15 HCW dy d? () =Cyd79d79 y=Rz d15d15 TH d?x df =B(f)?x d111d111 ?x=P(f)z d47d47 TH invariant formdz df() = ?z Figure 4.1: Different Representations of Relative Motion Equations motion equations, and the corresponding state vectors are shown below. x=T?1(f)?x =T?1(f)P(f)z (4.38) =T?1(f)P(f)R?1y The chief eccentricity equal to zero is the special case where t = ? and T?1PR?1 = I, for the choices of M and S shown in Equation (2.55) and Equation (5.13), respectively. Using Equation (4.38), the relative motion of two elliptic satellites can be related to an underlying circular relative motion, which can be described by the HCW equations. Given an initial state x0 at t0, the state x at time t can be determined as follows. The inverse of Equation (4.38) converts the initial condition x0 to y0, and Kepler?s Equation is used to convert t and t0 to ? and ?0. The solution given by the HCW equations is y(?(t)) = eC(???0)y0. Equation (4.38) is then used to convert y(?(t)) to x. 55 Chapter 5 Time-Varying Coordinate Transformations of the LERM? The previous chapter reviewed Lyapunov-Floquet theory and presented a Lyapunov- Floquet generalization of the HCW equations. Whereas this generalization was able to relate the relative motion in elliptic orbits to an underlying circular system, the HCW equations had to be evaluated at a virtual time as opposed to the true time. Lyapunov-Floquet theory clearly states the existence of a periodic coordinate change that relates the LERM to an underlying linear time-invariant system at the true time. Two Lyapunov-Floquet transformations are presented in this chapter that directly re- late the linearized equations of relative motion to the HCW equations with time as the independent variable. One of these transformations approximately matches satellite posi- tion at periapse, whereas the other approximately matches satellite position at apoapse. An integral-preserving transformation is also presented, which matches the relative motion in elliptic orbits to the HCW solution that shares certain integral values. 5.1 Lyapunov-Floquet Transformations 5.1.1 Periapse-Matching Transformation Whereas Lyapunov-Floquet theory guarantees the existence of P(t) and ?, it does not provide a mechanism for calculating them. Here, it will be assumed that for the LERM ?Material in this chapter taken from the following papers: [55] R. E. Sherrill, A. J. Sinclair, S. C. Sinha, and T. A. Lovell, ?Lyapunov-Floquet Transformation of Satellite Relative Motion in Elliptic Orbits,? AAS 13-466, presented at the 23nd AAS/AIAA Space Flight Mechanics Meeting, Kauai, Hawaii, 2013. [56] R. E. Sherrill, A. J. Sinclair, S. C. Sinha, and T. A. Lovell, ?Calibration of Hill-Clohessy- Wiltshire Initial Conditions for Elliptic Relative Motion,? submitted to the AAS/AIAA Astrodynamics Specialist Conference, Hilton Head, South Carolina, 2013. 56 it is possible to choose ? to equal the state matrix of the HCW equations, C, given in Equations. (2.53). Based on this assumption it must be determined if a corresponding P(t) can be computed which satisfies Equations. (4.6). This solution for P(t) will provide a coordinate transformation between the LERM solution x and an underlying HCW solution z. Recall that the state-transition matrix of the LERM is written as a function of f while the independent variable of the HCW equations is t. Defining f = 0 at t = 0 and arbitrarily choosing t0 = 0, Equation (4.6) can be rewritten as shown below. ?LERM(f,0) =P(f)eC(t?0)P?1(0) =P(f)?HCW(t,0)P?1(0) (5.1) As discussed in the previous section, P(f) must be periodic with the period of the system. Therefore, P(0) = P(2?). To determine P(0), the state-transition matrices of both the LERM and the HCW equations, given by Equations. (2.51) and (2.57) respectively, are evaluated at the end of one period. ?LERM(2?,0) = ? ?? ?? ?? ?? ?? ?? ?? ?? ? 1 0 0 0 0 0 ?6?(e+ 1)3(e+ 2) (1?e2)52 1 0 0 ?6?p2(e+ 1)2 h(1?e2)52 0 0 0 1 0 0 0 ?6?eh(e+ 1)4(e+ 2) p2(1?e2)52 0 0 1 ?6?e(e+ 1)3 (1?e2)52 0 0 0 0 0 1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ?? ? (5.2) 57 ?HCW parenleftbigg2? n ,0 parenrightbigg = ? ?? ?? ?? ?? ?? ?? ?? ? 1 0 0 0 0 0 ?12? 1 0 0 ?6?n 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? (5.3) Note that Equation (5.3) equals Equation (5.2) evaluated at e = 0. Equations (5.2) and (5.3) can be written in Jordan canonical form by defining ?LERM(2?,0) = MJM?1 and ?HCW parenleftbig2pin ,0parenrightbig=NJN?1. M = ? ?? ?? ?? ?? ?? ?? ?? ?? ? 0 2(1?e 2)52 (e+ 1)3(e+ 2) 0 ?1 0 0 ?12? 0 2 0 0 0 0 0 0 0 1 0 ?12?eh(e+ 1) p2 0 M43 0 0 0 0 0 0 h(e+ 1)(e+ 2)p2 0 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ?? ? (5.4) M43 = 2eh(e+ 1)p2 + n 2p2(e+ 1)2 2h(1?e2)52 J = ? ?? ?? ?? ?? ?? ?? ?? ? 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? ; N = ? ?? ?? ?? ?? ?? ?? ?? ? 0 1 0 ?1 0 0 ?12? 0 2 0 0 0 0 0 0 0 1 0 0 0 n2 0 0 0 0 0 0 2n 0 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? (5.5) 58 Note that the values ofM andN are not unique as the magnitude of each of the eigenvectors that comprise M and N are arbitrary. These eigenvectors correspond to three different second-order subsystems contained in J, related to a drifting in-plane motion, a periodic in-plane motion (here chosen to represent two different phasings of the motion on a 2 ? 1 ellipse), and a periodic out-of-plane motion. Since ?LERM(2?,0) and ?HCW parenleftbig2pin ,0parenrightbighave the same canonical form, J, a similarity transformation P0 can be determined. P0 =MN?1 MJM?1 =P0NJN?1P?10 (5.6) ?LERM(2?,0) =P0?HCW parenleftbigg2? n ,0 parenrightbigg P?10 P0 = ? ?? ?? ?? ?? ?? ?? ?? ?? ? 2(1?e2)52 (e+ 1)3(e+ 2) 0 0 0 (1?e2)52 n(e+ 1)3(e+ 2) ? 1 2n 0 0 1 0 0 0 0 0 0 1 0 0 0 0 eh(e+ 1)p2 0 np 2(e+ 1)2 h(1?e2)52 0 0 0 0 0 0 h(e+ 1)(e+ 2)2np2 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ?? ? (5.7) Due to the nonuniqueness of M and N, the value of P0 is also not unique. The scaling of M and N chosen here results in the following properties: det(P0) = 1, P0 =I for e = 0, and N equals M evaluated at e = 0. Note that the first three rows of M differ from the first three rows ofN by a single element. If the first three rows ofM equaled the first three rows of N, then the first three rows ofP0 would equal an identity matrix and a null matrix. This would provide exact matching between the position components ofxandz at periapse. The first row of Equation (5.7) shows that at periapse, the presented solution provides only approximate position matching in the radial direction. 59 LERM dx dt =A(t)xd79d79 ?x=T(f)x d15d15 d111d111 x=P(f)z d47d47 HCWdz dt =Cz TH d?x df =B(f)?x Figure 5.1: Different Representations of Relative Motion Equations It is now possible to determine P(f) by rewriting Equation (5.1) as shown below. P(f) = ?LERM(f,0)P0e?C(t?0) = ?LERM(f,0)P0??1HCW(t,0) (5.8) Note that P0 equals P evaluated at periapse. The elements of P(f) are shown in Ap- pendix A. As expected, P(f) = I when e = 0. Significantly, this solution for P(f) is T-periodic, thus forming a valid LF transformation, and ? = C is indeed a valid choice for the underlying linear time-invariant system. Therefore, Figure 2.9 which depicted the different representations of relative motion equations can be updated as shown in Figure 5.1 with a periodic transformation between the LERM and the HCW equations. 5.1.2 Apoapse-Matching Transformation Following a similar procedure, an alternative LF transformation ?P(f) can be deter- mined to provide approximate position matching at apoapse. To determine ?P(f), the state- transition matrices of both the LERM and the HCW equations are evaluated over one period, beginning and ending at apoapse (i.e. choosing t0 = pin). 60 ?LERM(3?,?) = ? ?? ?? ?? ?? ?? ?? ?? ?? ? 1 0 0 0 0 0 ?6?(e?1)3(e?2) (1?e2)52 1 0 0 ?6?p2(e?1)2 h(1?e2)52 0 0 0 1 0 0 0 ?6?eh(e?1)4(e?2) p2(1?e2)52 0 0 1 ?6?e(e?1)3 (1?e2)52 0 0 0 0 0 1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ?? ? (5.9) ?HCW parenleftbigg3? n , ? n parenrightbigg = ? ?? ?? ?? ?? ?? ?? ?? ? 1 0 0 0 0 0 ?12? 1 0 0 ?6?n 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? (5.10) Equations (5.9)and (5.10)can bewritten inJordan canonical form by defining?LERM(3?,?) = MJM?1 and ?HCW parenleftbig3pin , pinparenrightbig=NJN?1. M = ? ?? ?? ?? ?? ?? ?? ?? ?? ? 0 2(1?e 2)52 (e?1)3(e?2) 0 ?1 0 0 ?12? 0 2 0 0 0 0 0 0 0 1 0 ?12?eh(e?1) p2 0 2eh(e?1) p2 + n(4e+ 1) 2 0 0 0 0 0 0 h(e?1)(e?2)p2 0 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ?? ? (5.11) 61 J = ? ?? ?? ?? ?? ?? ?? ?? ? 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? ; N = ? ?? ?? ?? ?? ?? ?? ?? ? 0 1 0 ?1 0 0 ?12? 0 2 0 0 0 0 0 0 0 1 0 0 0 n2 0 0 0 0 0 0 2n 0 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ? (5.12) These matrices are used to determine ?P0. ?P0 =MN?1 = ? ?? ?? ?? ?? ?? ?? ?? ?? ? 2(1?e2)52 (e?1)3(e?2) 0 0 0 (1?e2)52 n(e?1)3(e?2) ? 1 2n 0 0 1 0 0 0 0 0 0 1 0 0 0 0 eh(e?1)p2 0 4e+ 1 0 0 0 0 0 0 h(e?1)(e?2)2np2 0 0 0 0 0 0 1 ? ?? ?? ?? ?? ?? ?? ?? ?? ? (5.13) Recall that the scaling of M and N given in Eqs. (5.11) and (5.12) is arbitrary, and note that for the scaling chosen here ?P0 = I for e = 0, and N equals M evaluated at e = 0. Similar to the periapse-matching transformation, ?P0 only provides approximate position matching in the radial direction at apoapse. It is now possible to determine ?P(f) by rewriting Equation (5.1) as shown below. ?P(f) = ?LERM(f,?) ?P0e?C(t?pin) = ?LERM(f,?) ?P0??1HCW parenleftBigt,? n parenrightBig (5.14) Note that ?P0 equals ?P evaluated at apoapse. The nonzero elements of ?P(f) are shown in Appendix B. Again, ?P(f) is T-periodic, providing a valid LF transformation. 62 Case aC (km) eC aD (km) eD iD (rad) ?D (rad) ?D (rad) 1 11,000 0.3 11,000 0.30001 0 0 0 2 11,000 0.3 11,000.2 0.30001 0 0 0 3 11,000 0.3 11,000 0.30001 0 0 4?10?5 4 11,000 0.6 11,000 0.60001 0 0 4?10?5 5 11,000 0.3 11,000 0.30001 4?10?5 0 0 6 11,000 0.3 11,000 0.30001 4?10?5 pi2 ?pi2 Table 5.1: Orbital elements for the chief and deputy for six different relative motion cases. 5.1.3 Relative Motion Examples Traditionally the HCW equations are thought of as a special case of the LERM obtained by setting the chief eccentricity equal to zero. However, by applying a LF transformation, motion described by the HCW equations can be transformed at any point to motion obeying the LERM and vice versa. This section will present several cases where this transformation is demonstrated. Six cases are presented here, with the orbital elements of the chief and deputy listed in Table 5.1. Orbital elements not listed are set to zero, i.e., iC = ?C = ?C = fC = fD = 0. Initial conditions for the LERM were calculated using Equation (2.60), while initial conditions for the HCW equations will use the method described in Section 5.3. Figures 5.2 through 5.6 show both a LERM trajectory and an HCW trajectory. Along the HCW trajectory, nine equally spaced points in time were selected and are designated by a circle. At each point, the periapse-matching transformation was applied to transform the motion to a corresponding point along the LERM trajectory given by a diamond. Since Cases 5 and 6 include an out-of-plane component, Figures 5.8 and 5.9 show a three-dimensional representation of the trajectory geometry. The figures clearly show that the HCW equations exactly capture the LERM dynamics and that the LERM and the HCW equations can be related by periodic coordinate transformation. In addition, the figures also show the effect of eccentricity on the relative motion. As can be seen, the equally-spaced points along the HCW trajectory are clustered near apoapse 63 ?0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 y (km) x (km) LERM HCW LF Figure 5.2: LERM, HCW, and LF transformation for case 1. ?7 ?6 ?5 ?4 ?3 ?2 ?1 0 ?2 ?1 0 1 2 3 y (km) x (km) LERM HCW LF Figure 5.3: LERM, HCW, and LF transformation for case 2. 64 0.2 0.3 0.4 0.5 0.6 0.7 ?0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 y (km) x (km) LERM HCW LF Figure 5.4: LERM, HCW, and LF transformation for case 3. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 0.25 y (km) x (km) LERM HCW LF Figure 5.5: LERM, HCW, and LF transformation for case 4. 65 0 0.2 0.4 0.6 0.8 1?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 time (P) z (km) LERM HCW LF Figure 5.6: LERM, HCW, and LF transformation for case 5. 0 0.2 0.4 0.6 0.8 1?0.4 ?0.3 ?0.2 ?0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 time (P) z (km) LERM HCW LF Figure 5.7: LERM, HCW, and LF transformation for case 6. 66 ?0.4 ?0.2 0 0.2 0.4 ?0.2 ?0.1 0 0.1 0.2 ?1 ?0.5 0 0.5 1 y (km)x (km) z (km) LERM HCW LF Figure 5.8: Three dimensional representation of Case 5. ?0.4 ?0.2 0 0.2 0.4 ?0.2 ?0.1 0 0.1 0.2 ?0.4 ?0.2 0 0.2 0.4 0.6 y (km)x (km) z (km) LERM HCW LF Figure 5.9: Three dimensional representation of Case 6. 67 on the LERM trajectory. The unequal distribution of points alludes to the motion progress- ing faster near periapse and slower near apoapse. The Virtual Time method described in Section 3.2 tried to emulate this natural progression of the motion by evaluating the HCW equations at a virtual time. While the VT method did significantly improve the error in directly applying the HCW equations to elliptic orbits, it was not able to capture the full time-varying dynamics. 5.2 Integral-Preserving Transformation Besides providing approximate position matching at periapse or apoapse, another trans- formation choice is to match the integral values of the x(t) and z(t) solutions. Specifically, the HCW solution that shares the same values for k is chosen. Recall from Equation (2.50) that the constant k for the LERM was given by k = ??1(f)T(f)x(t). Similarly for the HCW solution, the constant vector is given by k = ??1HCW(M)THCWz(t). THCW ? T(f)|e=0 = ? ?? ?? ?? ?? ?? ?? ?? ? 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1n 0 0 0 0 0 0 1n 0 0 0 0 0 0 1n ? ?? ?? ?? ?? ?? ?? ?? ? (5.15) ?HCW(M) ? ?(M)|e=0 = ? ?? ?? ?? ?? ?? ?? ?? ? sinM cosM 1 0 0 0 2cosM ?2sinM ?3M2 1 0 0 0 0 0 0 sinM cosM cosM ?sinM 0 0 0 0 ?2sinM ?2cosM ?32 0 0 0 0 0 0 0 cosM sinM ? ?? ?? ?? ?? ?? ?? ?? ? (5.16) 68 Note that the ?HCW is evaluated at the mean anomaly of the chief in the underlying circular orbit, which is equal to its true anomaly, and conveniently, also equal to the mean anomaly of the chief in the actual elliptic orbit. By equating the constants between the LERM and HCW solutions, the following transformation can be established. ??1(f)T(f)x(t) = ??1HCW(M)THCWz(t) (5.17) x(t) =T?1(f)?(f)??1HCW(M)THCWz(t) = ?(f)z(t) (5.18) By representing the LERM solution as a time-varying transformation of an HCW solution, the following integral-preserving (IP) transformation is proposed. x(t) = ?(f)eC(t?0)??1(f0)x0 = ?(f)?HCW(t,0)??1(f0)x0 (5.19) The nonzero elements of ?(f) are shown in Appendix C. Note that ?(f) =I when ec = 0. Recall that M increases monotonically with time. Since M appears outside of the trigono- metric functions in the elements of the first and fifth columns ?(f), this transformation is in fact aperiodic. This reveals that although ?(f) has a similar structure to an LF trans- formation, it is not an LF transformation. However, the simplicity and sparseness of ?(f) could make it attractive for certain relative-motion applications. 5.3 Calibration of Initial Conditions The previous subsections introduced several transformations between the LERM solu- tion and an HCW solution at any point in time. This subsection has a slightly different focus: using the underlying HCW solution as an approximation of the LERM solution. In this sense, the transformations are a means of calibrating initial conditions for the HCW solution. 69 ?1 ?0.8 ?0.6 ?0.4 ?0.2 0 ?0.4 ?0.3 ?0.2 ?0.1 0 0.1 0.2 0.3 0.4 0.5 y (km) x (km) LERM HCW Figure 5.10: Relative-motion trajectories for LERM and HCW propagation. 5.3.1 Calculation of Initial Conditions Section 2.2.8 gave an example of a chief and deputy satellite in close proximity. Using the initial conditions from this example the relative motion can be propagated using the LERM and HCW equations, and is presented in Fig. 5.10. Note that the instantaneous angular velocity is used in Equation (2.60), but the HCW equations assume constant angular velocity. Figure 5.10 shows that, in this case of an elliptic chief, significant error is introduced by propagating the initial conditions using the HCW equations. The shape of the LERM trajectory, however, is similar to a nondrifting HCW trajectory. This suggests that it should be possible to select an alternative set of initial conditions, that when propagated under the HCW equations, will better match the shape of the LERM trajectory. A novel strategy for selecting initial conditions would be to transform the linear time-varying system to a linear time-invariant system that can be described by the HCW equations. 70 The transformations presented in this chapter relate the LERM to the HCW equations. These transformations imply that there is an underlying circular system, z, which can be related to x through x=P(f)z, x= ?P(f)z, or x= ?(f)z. Since this underlying system can be described by the HCW equations, the previously-derived time-varying coordinate transformations can be used to determine a set of calibrated initial conditions. It is hoped that these calibrated initial conditions outperform the true initial conditions, given by Equa- tion (2.60), when applied to elliptic orbits. This method can be compared to other methods in the literature which attempt to match LERM and HCW trajectories.19 The first HCW trajectory uses the LF transformation which approximately matches satellite position at perigee and is obtained by post-multiplying the LF transformation evalu- ated at the initial time, i.e.,xHCW(t) = eCtP?1(f0)x0. The second HCW trajectory uses the apogee-matching LF transformation and is similarly obtained by xHCW(t) = eCt ?P?1(f0)x0. Finally, the IP transformation is given byxHCW(t) = eCt??1(f0)x0. The subsequent section will investigate the validity of using these transformation for several relative-motion cases. 5.3.2 Comparison of Methods To illustrate the effectiveness of the LF and IP transformations as initial condition calibrations, consider the case where a chief and deputy satellite are in close proximity around theEarth. Sixcases areconsidered, with thenon-zero orbital elements of thechief and deputy given in Table 5.1. The initial conditions were calculated using Equation (2.60). Figures 5.11 through 5.16 show the trajectory predicted by the LERM and three different approximate HCW trajectories. Figures 5.17 and 5.18 show a three-dimensional representation of the relative motion for cases 5 and 6. These figures show that each transformation approximates elliptic motion with some HCW solution. At each time step, k, the position error in the HCW approximations compared to the LERM solution can be computed. Let ??k represent the position of the HCW approximation obtained by propagating the perigee-matching LF transformation, the apogee-matching LF 71 ?0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 y (km) x (km) LERM Perigee LF Trans. Apogee LF Trans. IP Transformation Perigee LF Trans. Apogee LF Trans. IP Transformation LERM Figure5.11: Relative-motion trajectoryfor theLERM and three approximate HCW solutions for case 1. ?7 ?6 ?5 ?4 ?3 ?2 ?1 0 ?2 ?1 0 1 2 3 y (km) x (km) LERM Perigee LF Trans. Apogee LF Trans. IP Transformation Figure5.12: Relative-motion trajectoryfor theLERM and three approximate HCW solutions for case 2. 72 0.2 0.3 0.4 0.5 0.6 0.7 ?0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 y (km) x (km) LERM Perigee LF Trans. Apogee LF Trans. IP Transformation Figure5.13: Relative-motion trajectoryfor theLERM and three approximate HCW solutions for case 3. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 0.25 y (km) x (km) LERM Perigee LF Trans. Apogee LF Trans. IP Transformation Figure5.14: Relative-motion trajectoryfor theLERM and three approximate HCW solutions for case 4. 73 0 0.2 0.4 0.6 0.8 1 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 t (period) z (km) LERM Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.15: Out-of-plane trajectory for the LERM and three approximate HCW solutions for case 5. 0 0.2 0.4 0.6 0.8 1 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 t (period) z (km) LERM Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.16: Out-of-plane trajectory for the LERM and three approximate HCW solutions for case 6. 74 ?0.1 ?0.05 0 0.05 0.1 ?0.4 ?0.2 0 0.2 0.4 ?1 ?0.5 0 0.5 1 x (km)y (km) z (km) LERM Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.17: Three dimensional representation of the relative-motion trajectory for case 5. ?0.2 ?0.1 0 0.1 0.2 ?0.1 ?0.05 0 0.05 0.1 ?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4 x (km)y (km) z (km) LERM Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.18: Three dimensional representation of the relative-motion trajectory for case 6. 75 transformation, or the IP transformation at time k. Recall that the position of the LERM solution at timekis given by?k. The position error at each time step is given byek = ??k??k. The magnitude of the error, or the error distance, is given by dk = radicalbigeTkek. Figures 5.19 through 5.24 show the error magnitude, dk, for each case as a function of time for the perigee-matching LF transformation, the apogee-matching LF transformation, and the IP transformation. To further investigate the error of each method, the mean error, ?M, and the root-mean- square (RMS) of the error, ?RMS, can be calculated and are given by the following. ?M = 1m msummationdisplay k=1 dk; ?RMS = radicaltpradicalvertex radicalvertexradicalbt1 m msummationdisplay k=1 d2k (5.20) Here, m represents the total number of time-steps in each solution. The performance of the HCW approximation given by the perigee-matching LF transformation, the apogee-matching LF transformation, and the IP transformation can be compared to the HCW solution whose initial conditions equal the initial conditions of the LERM, i.e., x0. Table 5.2 gives the mean error for each case, while Table 5.3 gives the RMS error. As can be seen, the calibrated HCW solutions clearly outperform the nominal HCW trajectory for every case. The results also show that each transformation performed better for certain cases, but there was not a transformation which clearly outperformed the others for the cases considered. This perhaps indicates that a LF transformation can be derived for mission-specific goals instead of a general solution for all relative-motion cases. In addition, cases where the underlying relative-motion trajectory appeared to closely match the LERM trajectorydid not necessarily translate into low error. This can be best seen by examining case 2. Inspection of the relative-motion trajectory might initially indicate that the perigee-matching LF transformation might be an ideal choice to approximate the HCW trajectory. However, based on the error shown in Tables 5.2 and 5.3, the IP transformation is shown to be a better choice. This disparity is due to a satellite spending more time near 76 0 0.2 0.4 0.6 0.8 10 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 t (period) error (km) Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.19: Error distance as a function of time for the three underlying HCW trajectories in case 1. 0 0.5 1 1.5 2 2.5 30 0.5 1 1.5 2 2.5 3 3.5 4 t (period) error (km) Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.20: Error distance as a function of time for the three underlying HCW trajectories in case 2. 77 0 0.2 0.4 0.6 0.8 10 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 t (period) error (km) Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.21: Error distance as a function of time for the three underlying HCW trajectories in case 3. 0 0.2 0.4 0.6 0.8 10 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t (period) error (km) Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.22: Error distance as a function of time for the three underlying HCW trajectories in case 4. 78 0 0.2 0.4 0.6 0.8 10 0.05 0.1 0.15 0.2 0.25 t (period) error (km) Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.23: Error distance as a function of time for the three underlying HCW trajectories in case 5. 0 0.2 0.4 0.6 0.8 10 0.05 0.1 0.15 0.2 0.25 0.3 0.35 t (period) error (km) Perigee LF Trans. Apogee LF Trans. IP Transformation Figure 5.24: Error distance as a function of time for the three underlying HCW trajectories in case 6. 79 Perigee LF Trans. Apogee LF Trans. IP Transformation True IC 1 0.0489 0.0489 0.0489 1.5758 2 1.3537 0.9167 0.8183 1.9717 3 0.0644 0.0607 0.1199 1.7034 4 0.2106 0.1200 0.2869 6.9997 5 0.1300 0.0819 0.0642 1.5917 6 0.2077 0.2057 0.2046 1.6118 Table 5.2: Mean error (km) for each case over the simulation time. Perigee LF Trans. Apogee LF Trans. IP Transformation True IC 1 0.0530 0.0530 0.0530 1.9777 2 1.6619 1.2426 1.0098 2.4002 3 0.0775 0.0718 0.1210 2.0971 4 0.2237 0.1358 0.3093 8.6714 5 0.1476 0.0949 0.0698 1.9825 6 0.2255 0.2341 0.2105 1.9900 Table 5.3: RMS error (km) for each case over the simulation time. apogee. Since the perigee-matching LF transformation correctly matches position over a small portion of an orbit, it does not approximate the HCW trajectory over the entirety of the orbit. The apogee-matching transformation is a better choice, but trajectory lags behind the LERM causing significant error when the satellite is near perigee. The drift rate of the IP transformation is between the perigee-matching and apogee-matching LF transformations, and can better capture the motion on average over an orbit. This perhaps indicates the perigee-matching and apogee-matching transformations are upper and lower boundaries of reasonable LF transformations. 80 Chapter 6 Control of Relative Orbits This chapter will develop methods to control the relative motion of satellites in elliptic orbits. The literature contains impulsive methods and continuous-thrust control of satellite formations with fixed final time.16,57,58 This chapter will develop a method for continuous low-thrust maneuvering based on infinite-horizon optimal control. Since continuous-thrust control exists for circular orbits, it is proposed that the Lyapunov-Floquet transformations presented in the previous section be used to transform the control based on time-invariant equations to time-varying systems. This chapter will review an impulsive control method for two-burn solutions and linear quadratic regulators. Then, control of time-varying systems using Lyapunov-Floquet theory is introduced. Finally, this chapter will present four low- thrust scenarios comparing the results of different control aspects. 6.1 Impulsive Control Because the state-transition matrix for the LERM exists in closed-form, a simple two- burn solution to the two-point boundary-value problem can be constructed. This two-burn maneuver can also referred to as pseudo-Lambert method due to the similarity to Lambert?s problem for determining Keplerian transfer orbits. The first burn places the satellite on a trajectory intersecting the origin of the LVLH frame, while the second burn neutralizes any arrival velocity to keep the satellite at the origin. The state-transition matrix of the LERM, given by Equation (2.51), can be partitioned in the following manner. ? ?? r(t) v(t) ? ??= ? ?? ?rr(t,t0) ?rv(t,t0) ?vr(t,t0) ?vv(t,t0) ? ?? ? ?? r(t0) v(t0) ? ?? (6.1) 81 Each partition in Equation (6.1) is a 3 ? 3 matrix. Using the partitioned state-transition matrix, separate equations for the position and velocity at time t are given below. r(t) = ?rr(t,t0)r(t0) +?rv(t,t0)v(t0) (6.2) v(t) = ?vr(t,t0)r(t0) +?vv(t,t0)v(t0) (6.3) Note that in this case, t is a specified time in which to complete the maneuver. Since the first burn places the satellite on a trajectory toward the origin and the goal is rendezvous at time t, r(t) = 0. The velocity of this rendezvous trajectory can be determined from Equation (6.2). A superscript (+) indicates velocity after the impulsive burn. v+(t0) = ???1rv (t,t0)?rr(t,t0)r(t0) (6.4) The change in velocity necessary to put the satellite on the rendezvous trajectory is given by the difference between the pre-burn and post-burn velocities. ?v1 =v+(t0)?v(t0) (6.5) The final velocity can be found from Equation (6.3) where a superscript (?) indicates velocity before the impulsive burn. v?(t) = ?vr(t,t0)r(t0) +?vv(t,t0)v+(t0) = ?vr(t,t0)r(t0)??vv(t,t0)??1rv (t,t0)?rr(t,t0)r(t0) (6.6) =bracketleftbig?vr(t,t0)??vv(t,t0)??1rv (t,t0)?rr(t,t0)bracketrightbigr(t0) (6.7) The second velocity change brings the final velocity to zero. ?v2 = 0?v?(t) = ?v?(t) (6.8) 82 The total velocity magnitude is found from Equations (6.5) and (6.8). vtotal = bardbl?v1bardbl+bardbl?v2bardbl (6.9) A unique feature of the two-burn method is that the time to complete the maneuver needs to be specified. For certain initial conditions, changing the maneuver time has little effect on the deputy?s trajectory. Consider the case where the chief?s orbital elements are given by aC = 8000 km, eC = 0.1, and iC = ?C = ?C = fC = 0. The relative position and velocity of the deputy at time t0 is given by x(t0) = [0 ?2 0 0 0 0]T. Figure 6.1 shows the deputy?s trajectory for t = P with a total velocity change of 2.8872?10?4 km/s, while Figure 6.2 shows the deputy?s trajectory for t = 5P with a total velocity change of 2.1313?10?4 km/s. For this example of a chief and deputy satellite in a ?leader-follower? configuration, the trajectory is similar for both time spans and the difference in required velocity is on the order of centimeters per second. Consider a second case where the chief?s orbital elements are the same but the relative position and velocity of the deputy at time t0 is given by x(t0) = [0.1 0 0 0 0 0]T. Figure 6.3 shows the deputy?s trajectory for t = P with a total velocity change of 1.8051?10?2 km/s, while Figure 6.4 shows the deputy?s trajectory for t = 0.5P with a total velocity change of 2.5145 ? 10?4 km/s. For this case, changing the maneuver time to a fraction of a period dramatically changes the size of the relative trajectory and lowered the required velocity change by two orders of magnitude. As was seen with the above examples, the maneuver time needs to be chosen with care when employing the two-burn method. Okasha and Newman developed a method to specify intermediate points along the relative trajectory when using the two-burn method, which forces the relative motion to be bounded.59 83 ?2 ?1.5 ?1 ?0.5 0 ?1 ?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 x (km) y (km) Figure 6.1: Two-burn maneuver over one period. ?2 ?1.5 ?1 ?0.5 0 ?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 y (km) x (km) Figure 6.2: Two-burn maneuver over five periods. 84 ?35 ?30 ?25 ?20 ?15 ?10 ?5 0 ?10 ?5 0 5 10 y (km) x (km) Figure 6.3: Two-burn maneuver over one period. ?0.14 ?0.12 ?0.1 ?0.08 ?0.06 ?0.04 ?0.02 0 0.02 ?0.02 0 0.02 0.04 0.06 0.08 0.1 y (km) x (km) Figure 6.4: Two-burn maneuver over half a period. 85 6.2 Linear Quadratic Regulator Consider the linear time-invariant system ?x=Ax+Bu where x represents the states and u represents the control. It is desired to control the system while achieving some sort of system performance. This performance could include factors such as rise-time, overshoot, or to use a minimum control effort. To give a qualitative measure of system performance, a cost function J is introduced. By minimizing this cost function, the system can brought from an initial state at t0 to a final state at tf with acceptable levels of state error and using acceptable levels of control and along this trajectory. A common choice for the performance index is to use a quadratic form of the terminal state and an integral of the quadratic form of the current state and control. J = 12xT(tf)Sfx(tf) + 12 integraldisplay tf t0 parenleftbigxTQx+uTRuparenrightbigdt (6.10) Here, Sf and Qare symmetric positive semidefinite matrices and R is a symmetric positive definite matrix. The matrixSf represents the weighting on the final state, Q represents the weighting on the state, and R represents the weighting on the control. The Hamiltonian, H, can be formed with the addition of Lagrange multipliers or co-states, ?. The co-states represent sensitivity of the cost to the states. The form of the Hamiltonian is shown below. H = 12parenleftbigxTQx+uTRuparenrightbig+?T (Ax+Bu) (6.11) The Hamiltonian gives the three first-order necessary conditions for optimality. ??= ??H ?x = ?Qx?A T? (6.12) 0 = ?H?u =Ru+BT? (6.13) ?x= ?H?? =Ax+Bu (6.14) 86 From Equation (6.13), u= ?R?1BT?. This allows the closed loop system to be written in the following form. ?x=Ax?BR?1BT? (6.15) From Equations (6.12) and (6.14), the following system can be constructed. ? ?? ?x ?? ? ??= ? ?? A ?BR?1BT ?Q ?AT ? ?? ? ?? x ? ? ?? (6.16) Note that each entry in the state matrix in Equation (6.16) is n?n. Since each system in Equation (6.16) is linear, a transformation ?(t) = S(t)x(t) is proposed. Taking one derivative and using Equation (6.12) gives the following expression. ??= ?S(t)x(t) +S(t) ?x(t) = ?Qx(t)?AT?(t) ?S(t)x(t) +S(t)(Ax(t)?BR?1BTS(t)x(t)) = ?Qx(t)?ATS(t)x(t) parenleftBig ? S(t) +S(t)A+ATS(t)?S(t)BR?1BTS(t) +Q parenrightBig x(t) = 0 (6.17) Since x(t) negationslash= 0, Equation (6.17) gives the following. ?S(t) +S(t)A+ATS(t)?S(t)BR?1BTS(t) +Q= 0 (6.18) Equation (6.18) is known as the Ricatti Equation. As tf ? ?, stead-state solutions of the Ricatti Equation give ?S = 0. This steady-state value is a solution to the Algebraic Ricatti Equation. SA+ATS?SBR?1BTS+Q= 0 (6.19) 87 Defining K = ?R?1BTS and using ?(t) =Sx(t), the control law can be determined from Equation (6.15) and the system can be propagated forward in time. ?x=Ax?BR?1BTSx ?x= (A?BK)x (6.20) 6.3 Control of Time-Periodic Systems using Lyapunov-Floquet Theory 6.3.1 Control Method A control method described by Sinha and Joseph is here presented.60,61 Consider the state equation of a controllable linear time-periodic system. ?x(t) =A(t)x(t) +B(t)u(t) (6.21) Applying a Lyapunov-Floquet transformation of the control-free system to Equation (6.21) gives the following. ?z(t) = ?z(t) +P?1(t)B(t)u(t) (6.22) If the product P?1(t)B(t) is considered to be the control matrix, then Equation (6.22) has a constant system matrix and a time-varying control matrix. Sinha and Joseph proposed the following auxiliary system which has a more desirable form. ?tildewidez(t) = ?tildewidez(t) + tildewideBv(t) (6.23) Here, tildewideB is chosen such that ? and tildewideB are controllable, and a control law v(t) = ?tildewiderKtildewidez(t) is designed. The time-invariant gain matrix tildewiderK can be chosen through pole-placement or optimal control techniques such that the auxiliary system is asymptotically stable. Using this constant gain matrix, the closed-loop system is given by ?tildewidez = parenleftBig ?? tildewideBtildewiderK parenrightBig tildewidez, and the closed- loop state-transition matrix ?tildewidez(t2,t1) = e(??tildewideBtildewiderK)(t2?t1) is defined. The dynamics of the 88 state-transition matrix are given by ??tildewidez = parenleftBig ?? tildewideBtildewiderK parenrightBig ?tildewidez. Since tildewidez is chosen to approximate z, Equations (6.22) and (6.23) would be equivalent ifP?1(t)B(t)u(t) = ?tildewideBtildewiderKz(t). Instead, the equalityB(t)u(t) = ?P(t)tildewideBtildewiderKz(t) is chosen because this form better approximates the true time-varying system given by Equations (6.21). In general, no exact solution for u(t) exists. The approximate solution that minimizes the error norm is given by the following expression. u(t) = ?parenleftbigBT(t)B(t)parenrightbig?1BT(t)P(t)tildewideBtildewiderKz(t) (6.24) Using Equations (6.24) and (6.22), the closed-loop Lyapunov-Floquet system is given by the following. ?z(t) = ?z(t)?P?1(t)B(t)parenleftbigBT(t)B(t)parenrightbig?1BT(t)P(t)tildewideBtildewiderKz(t) = parenleftBig ??P?1(t)B(t)parenleftbigBT(t)B(t)parenrightbig?1BT(t)P(t)tildewideBtildewiderK parenrightBig z(t) (6.25) Applying the Lyapunov-Floquet transformation to Equation (6.25) results in the following. ?x(t) =A(t)x(t)?B(t)parenleftbigBT(t)B(t)parenrightbig?1BT(t)P(t)tildewideBtildewiderKP?1(t)x(t) = parenleftBig A(t)?B(t)parenleftbigBT(t)B(t)parenrightbig?1BT(t)P(t)tildewideBtildewiderKP?1(t) parenrightBig x(t) (6.26) Therefore, the control law for the time-varying system is given by u(t) = ?K(t)x(t), with the following gain matrix. K(t) =parenleftbigBT(t)B(t)parenrightbig?1BT(t)P(t)tildewideBtildewiderKP?1(t) (6.27) This time-varying gain matrix applies the calculated time-invariant gain to the time-varying system. Equation (6.27) allows the closed-loop system given by Equation (6.26) to be written in the following form. ?x(t) = (A(t)?B(t)K(t))x(t) (6.28) 89 The dynamics of the closed-loop state-transition matrix for Equation (6.28) is given by ??x = (A(t)?B(t)K(t))?x. Because an approximate solution for u(t) is used, stability of x is not guaranteed. Therefore, the eigenvalues of ?x(P,0) must be checked to verify stability. To demonstrate this control method, Sinha and Joseph presented the following example of a commutative system where ? is a parameter. ?x(t) = 2? ? ?? ?1 +?cos2 (2?t) 1??sin(2?t)cos(2?t) ?1??sin(2?t)cos(2?t) ?1 +?sin2 (2?t) ? ??x(t)+ ? ?? 0 t ? ??u(t) (6.29) The state transition matrix of this system is given by the following. ?(t,t0) = ? ?? e2pi(??1)(t?t0) cos(2?(t?t0)) e?2pi(t?t0) sin(2?(t?t0)) ?e2pi(??1)(t?t0) sin(2?(t?t0)) e?2pi(t?t0) cos(2?(t?t0)) ? ?? (6.30) Equation (6.30) has a known Lyapunov-Floquet transformation where P(t) is a rotation matrix and C is a diagonal matrix. P(t) = ? ?? cos(2?t) sin(2?t) ?sin(2?t) cos(2?t) ? ?? ; C = ? ?? 2?(??1) 0 0 ?2? ? ?? (6.31) Clearly, the system is unstable for ?> 1. Consider the case where ? = 1.2 and x(0) = [0.6 0.3]T. Figure 6.5 shows the uncon- trolled response of the system. The method presented in this section will be used to stabilize the system. Here, tildewideB was chosen to equal a 2?2 identity matrix, and tildewiderK was chosen to be a diagonal matrix. The poles, p1 and p2, of the closed-loop time-invariant system were selected to equal ?3.5 and ?5 respectively, with the gains determined through the pole-placement technique. Figure 6.6 shows the response of the closed-loop time-invariant auxiliary system 90 0 1 2 3 4 5?200 ?100 0 100 200 300 400 time (period) States State 1 State 2 Figure 6.5: Uncontrolled states of the system. 0 1 2 3 4 50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 time (period) States State 1 State 2 Figure 6.6: State feedback control of time-invariant auxiliary system. given by Equation (6.23), and Figure 6.7 shows the closed-loop response to the time-varying system where the time-varying gain matrix was calculated using Equation (6.27). 91 0 1 2 3 4 5?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 0.8 time (period) States State 1 State 2 Figure 6.7: State feedback control of the time-varying system. 6.3.2 Discussion An investigation was conducted to determine the relationship between the stability of the time-invariant auxiliary system and the time-varying system, i.e. are there pole locations that stabilize the time-invariant system but do not stabilize the time-varying system, and does the time-invariant system need to be stable for the time-varying system to be stable? Stability was determined by integrating the closed-loop state transition matrix of the time- varying system over one period. If the magnitude of all the eigenvalues are less than one, then the time-varying system is asymptotically stable. This investigation considered both positive and negative values of p1 and p2. These time-invariant pole locations were varied over the range ?10