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