Course Hero - We put you ahead of the curve!

You have requested the below document.

Coursehero >> Texas >> Texas >> BAMFORDW 82161

by Copyright William Alfred Bamford Jr. 2004 The Dissertation Committee for William Alfred Bamford Jr. certifies that this is the approved version of the following dissertation: Navigation and Control of Large Satellite Formations Committee: E. Glenn Lightsey, Supervisor Robert H. Bishop Wallace T. Fowler Cesar Ocampo Russell Carpenter Navigation and Control of Large Satellite Formations by William Alfred Bamford Jr., B.S., M.S. DISSERTATION Presented to the Faculty of the Graduate School of The University of Texas at Austin in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY THE UNIVERSITY OF TEXAS AT AUSTIN December 2004 Acknowledgments First and foremost, this research was partially funded by the NASA Cooperative Agreement No: NCC5-732. To my advisor: Dr. Glenn Lightsey: Thank you for sticking with me, even though it seemed as if I would never go away. You kept pulling for me long after most would have stopped, and I appreciate that. To the members of my committee: Thank you for the time you spent reading, re-reading, and in some cases, re-re-reading my dissertation. To Dr. Takuji Ebinuma: After all of your help, this dissertation is almost as much yours as it is mine. Thanks for always being around to answer my questions. To Larry Jackson, Rich Burns, Dr. Mike Moreau, Dr. Dave Gaylor, and the guys at Emergent Space Technologies: Thanks for the remote support, and the hours on the telephone trying to make the Spirent simulator work correctly. To the research group: Greg Holt, Jacob Williams, and Shaun Stewart: Thanks for the tips and the white board art (though I almost never understood it). Also, thanks for always being there when I needed a favor. To my parents and in-laws: Thanks for the constant support, even from half a continent away. iv To my children: Ashley, Tyler and Zachary: You always understood when I had to miss your soccer practices or baseball games. You tolerated seeing me for one hour a day, and that was usually on the drive to Mommy's work. The nights that I was home, you usually saw more of the back of my head then the front (as I worked on the computer). Yet, through it all, you were always there with a smile, a hug, and a kiss on the cheek. You three provided the motivation to keep going when things got hard. I thank you for being understanding, and still loving your Daddy through it all. Finally, to my wife Katie: It took seven long years, but it is finally over. You bore each sacrifice alongside of me, and you never complained, no matter how tough things got. (And things couldn't get much worse than University Housing). You took care of the kids when I was gone, and took care of me when I was home. You spent seven years pulling the weight for the both of us at home, while I chased my dreams in a classroom. I couldn't have accomplished any of this without you. Thank you, and I love you. v Navigation and Control of Large Satellite Formations Publication No. William Alfred Bamford Jr., Ph.D. The University of Texas at Austin, 2004 Supervisor: E. Glenn Lightsey In recent years, there has been substantial interest in autonomous satellite formations, driven by the new technologies that enable smaller and cheaper spacecraft. Formation flying allows for mission designs, such as stereoscopic imaging, that are impractical or impossible for a single satellite. Much of the current work focuses upon small formations, which can be defined as four or less satellites in a relatively tight grouping. Next generation formations may be composed of more satellites spanning greater spatial distances. The large formation problem becomes more difficult for several reasons, including an increased amount of communication required between the satellites, and orbit perturbations, which become more important as the formation size grows. The purpose of this dissertation is to examine formation flying for large formations, and determine whether or not generalizations can be made linking the large and small formation regimes. vi In order to model formations with many satellites, a simulation environment was constructed in which different observers, controllers, and formation architectures can be modelled. This dissertation focuses on a decentralized control scheme, but the software is general enough to accommodate a variety of control architectures. Validation of the large formation models is accomplished by initially modelling only a pair of satellites and comparing the results against those found in the literature. As a demonstration of the theoretical results, a real-time, closed-loop, hardware-in-the-loop simulation was constructed using GPS receivers as the measurement source. A large constellation, real-time simulation system was developed that utilized the Internet to connect simulation equipment from research centers in different locations. vii Table of Contents Acknowledgments Abstract List of Tables List of Figures Chapter 1. Introduction 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 TECHSAT-21 . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 A-TRAIN . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Research Contributions . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Small Formation to Large Formation Generalization . . 1.3.2 Extended Kalman Filter Design . . . . . . . . . . . . . . 1.3.3 Demonstration of a Generalized Kalman Filter . . . . . 1.3.4 Construction of the Extended Formation Flying Testbed 1.3.5 GPS Specific Instance of the EFFTB . . . . . . . . . . . 1.4 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 2. Reference Systems 2.1 Time Systems . . . . . . . . . . . . 2.1.1 Solar and Sidereal Times . . 2.1.2 Atomic time . . . . . . . . . 2.1.3 Coordinated Universal Time 2.1.4 GPS Time . . . . . . . . . . 2.1.5 GPS Roll Over . . . . . . . . 2.2 Reference Frames . . . . . . . . . . viii iv vi xiii xiv 1 1 2 3 3 5 5 6 6 7 7 7 9 9 9 10 10 10 11 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Earth Centered Inertial Reference Frame (ECI) . . . . . 2.2.2 Earth Centered Earth Fixed Reference Frame (ECEF) . 2.2.3 Spacecraft Centered Frame . . . . . . . . . . . . . . . . Chapter 3. GPS Models and Measurements 3.1 GPS Observables . . . . . . . . . . . . . . 3.1.1 Pseudorange Measurement . . . . . 3.1.1.1 User Range Error . . . . . . 3.1.1.2 Ionospheric Delay . . . . . . 3.1.2 Carrier Phase Measurement . . . . . 3.1.3 Doppler Measurement . . . . . . . . 3.2 Double-Differenced Observables . . . . . . 3.2.1 Double-Differenced Carrier Phase . . 3.2.2 Double-Differenced Doppler . . . . . 3.2.3 GPS Satellite Propagation . . . . . . 11 12 12 14 14 14 16 16 17 19 19 21 22 22 25 25 26 27 28 28 29 29 30 30 32 32 33 33 35 35 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 4. Kalman Filter Design 4.1 General Formation Description . . . . . . . . . . . . . . . . . . 4.2 Filter State . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Absolute And Relative States . . . . . . . . . . . . . . . 4.2.2 Receiver Clock and Frequency Offset State . . . . . . . 4.2.3 Double-Differenced Ambiguity State . . . . . . . . . . . 4.2.4 GPS Clock Offset . . . . . . . . . . . . . . . . . . . . . 4.3 Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 State Vector Construction . . . . . . . . . . . . . . . . . 4.3.2 Propagation . . . . . . . . . . . . . . . . . . . . . . . . 4.3.3 Initialization . . . . . . . . . . . . . . . . . . . . . . . . 4.3.4 Measurement Update . . . . . . . . . . . . . . . . . . . 4.4 Dynamic Model and State Transition Matrix . . . . . . . . . . 4.4.1 Dynamic Model . . . . . . . . . . . . . . . . . . . . . . 4.4.2 State Transition Matrix . . . . . . . . . . . . . . . . . . 4.4.2.1 Absolute State Acceleration Partial Derivatives Abs . . . . . . . . . . . . . . . . . . . . . . . . ix 4.4.2.2 Relative State Acceleration Partial Derivatives Rel . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2.3 Receiver Clock States Partial Derivatives b . . 4.4.2.4 Double-Differenced Ambiguity Partial Derivatives N . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2.5 GPS Clock Error Model B . . . . . . . . . . . 4.5 Process Noise Matrix . . . . . . . . . . . . . . . . . . . . . . . 4.5.1 Absolute and Relative State Process Noise . . . . . . . . 4.5.2 GPS Receiver Clock Process Noise . . . . . . . . . . . . 4.5.3 Double-Differenced Ambiguities Process Noise . . . . . . 4.5.4 GPS Satellite Clock Bias . . . . . . . . . . . . . . . . . 4.6 GPS Measurement Equations . . . . . . . . . . . . . . . . . . . 4.6.1 Undifferenced Measurement Equations . . . . . . . . . . 4.6.2 Double-Differenced Measurement Equations . . . . . . . 4.7 Measurement Noise Covariance . . . . . . . . . . . . . . . . . . 4.7.1 Pseudorange Measurement . . . . . . . . . . . . . . . . 4.7.2 Doppler Measurement . . . . . . . . . . . . . . . . . . . 4.7.3 Double-Differenced Carrier Phase Measurement . . . . . Chapter 5. Control System 5.1 Controller Development . . . . . 5.2 Controller Implementation . . . 5.3 Offline Controller Tuning . . . . 5.4 Formation Control Architecture 36 36 36 37 37 38 39 39 39 39 40 40 41 42 42 42 45 45 47 48 51 57 57 57 59 60 61 62 62 63 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 6. Extended Formation Flying 6.1 Hardware Description . . . . . . . . . 6.1.1 GPS Signal Simulator . . . . . 6.1.2 Orion GPS Receivers . . . . . 6.1.3 Master Control Computer . . . 6.1.4 Remote Control Computer . . 6.1.5 Hub Computer . . . . . . . . . 6.1.6 Node Computer . . . . . . . . 6.2 Software Description . . . . . . . . . x Testbed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 7. Results 7.1 Offline Filter Validation . . . . . . . . . . . . 7.1.1 Filter Performance at 1 Hz . . . . . . 7.1.2 Filter Performance at 5 Seconds . . . 7.1.3 Cycle Slips . . . . . . . . . . . . . . . 7.1.4 Best Case Filter Performance . . . . . 7.1.5 Receiver Clock Model Validation . . . 7.2 Extended Formation Flying Testbed Results 7.2.1 Two Satellite Local Simulation . . . . 7.2.2 Multi-Satellite Formation . . . . . . . Chapter 8. Conclusion 8.1 Summary of Results . . . 8.2 Future Work . . . . . . . 8.2.1 Ionosphere Study . 8.2.2 Multi-Path Study 8.2.3 Attitude Modelling 8.2.4 Formation Control 8.2.5 EFFTB Extension Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 65 67 74 79 80 86 88 88 92 99 99 102 102 103 103 103 104 105 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . and Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix A. Dynamic Equation Partial Derivatives 106 A.1 Dynamic Equation Partial Derivatives . . . . . . . . . . . . . . 106 Appendix B. Measurement Equation Partial Derivatives B.1 Pseudorange Partial Derivatives . . . . . . . . . . . . . . B.2 Doppler Partial Derivatives . . . . . . . . . . . . . . . . . B.3 Double-Differenced Carrier Phase Partial Derivatives . . B.4 Double-Differenced Doppler Partial Derivatives . . . . . . 110 110 111 112 114 . . . . . . . . . . . . xi Appendix C. Filter Verification Procedures C.1 Verification of the Dynamic Partial Derivative Matrix . . C.2 Verification of the Measurement Partial Derivative Matrix C.3 Verification of Measurement Noise Settings . . . . . . . . C.4 Covariance Update Verification . . . . . . . . . . . . . . . Bibliography Vita . . . . . . . . . . . . 118 118 120 122 123 125 129 xii List of Tables 1.1 2.1 2.2 2.3 3.1 5.1 5.2 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.10 Actual and Proposed Formation Flying Missions.[22] . . . . . J2000 Reference System Definition . . . . . . . . . . . . . . . ECEF Reference System Definition . . . . . . . . . . . . . . . UVW Reference System Definition . . . . . . . . . . . . . . . Almanac Information . . . . . . . . . . . . . . . . . . . . . . . 2 12 12 13 23 Controller Design Parameters . . . . . . . . . . . . . . . . . . 49 Initial Node Absolute State Offset From The Desired Trajectory 54 66 66 67 68 73 73 75 78 86 89 95 98 119 120 121 123 Simulator Environment Settings. . . . . . . . . . . . . . . . . True Initial Conditions: November 6 2001 00:00:00 GPS. . . . Error Covariance Initialization . . . . . . . . . . . . . . . . . . Absolute State Error in u, v , w frame (mean std) . . . . . . ^ ^ ^ Process Noise Settings: 1Hz Update . . . . . . . . . . . . . . . Relative State Error: 1 Sec. Update Rate (mean std) . . . . Process Noise Settings: 0.2 Hz Update . . . . . . . . . . . . . Relative State Error: 5 Sec. Update Rate (mean std) . . . . Relative State Error: Best Case Filter Runs (mean std) . . Relative State Error: EFFTB 2 Satellite Simulation (mean std) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.11 Relative State Error: EFFTB 5 Satellite Simulation (mean std) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.12 Control Effort in the u, v , w frame . . . . . . . . . . . . . . . . ^ ^ ^ C.1 C.2 C.3 C.4 Dynamic Partial Differential Equation Validation . . . Dynamic Partial Differential Equation Validation Cont. Measurement Partial Differential Equation Validation . Measurement Noise Verification. . . . . . . . . . . . . . xiii . . . . . . . . . . . . . . . . List of Figures 3.1 4.1 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 6.1 6.2 6.3 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.10 7.11 7.12 Double-Difference. . . . . . . . . . . . . . . . . . . . . . . . . Problem Setup. . . . . . . . . . . . . . . . . . . . . . . . . . . Tracking Error With Too Much Control Penalty. . . . . Tracking Error With Too Little Control Penalty. . . . . Tracking Error With Appropriate Penalty. . . . . . . . Tracking With 2 m Of Measurement Noise at 300 Sec. . Tracking With .5 m Of Measurement Noise at 300 Sec. Tracking With 2 m Of Feedback Noise at 10 Sec. . . . Control Architectures . . . . . . . . . . . . . . . . . . . Implemented Decentralized Architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 26 50 51 52 52 53 53 55 56 58 59 64 69 69 70 70 71 71 72 72 74 75 76 76 Hardware Setup. . . . . . . . . . . . . . . . . . . . . . . . . . Orion Receivers. . . . . . . . . . . . . . . . . . . . . . . . . . . Software Flowchart. . . . . . . . . . . . . . . . . . . . . . . . . Ensemble Mean Absolute Position Actual and Formal Ensemble Mean Absolute Velocity Actual and Formal 1 km Ensemble Mean Relative Position . . . . . . . . 1 km Ensemble Mean Relative Velocity . . . . . . . 100 km Ensemble Mean Relative Position . . . . . . . 100 km Ensemble Mean Relative Velocity . . . . . . 500 km Ensemble Mean Relative Position . . . . . . 500 km Ensemble Mean Relative Velocity . . . . . . Common View GPS Satellites. . . . . . . . . . . . . . 1 km Ensemble Mean Relative Position . . . . . . . 1 km Ensemble MeanRelative Velocity . . . . . . . . 100 km Ensemble Mean Relative Position . . . . . . . xiv Errors. . Errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.13 7.14 7.15 7.16 7.17 7.18 7.19 7.20 7.21 7.22 7.23 7.24 7.25 7.26 7.27 7.28 7.29 7.30 7.31 7.32 7.33 7.34 7.35 7.36 7.37 100 km Ensemble Mean Relative Velocity . . . . . . . . . . . . 500 km Ensemble Mean Relative Position . . . . . . . . . . . . 500 km Ensemble Mean Relative Velocity . . . . . . . . . . . Best Case 1 km Relative Position . . . . . . . . . . . . . . . . Best Case 1 km Relative Velocity . . . . . . . . . . . . . . . . Best Case 100 km Relative Position . . . . . . . . . . . . . . . Best Case 100 km Relative Velocity . . . . . . . . . . . . . . . Best Case 1 km Relative Position . . . . . . . . . . . . . . . . Best Case 1 km Relative Velocity . . . . . . . . . . . . . . . . Best Case 100 km Relative Position . . . . . . . . . . . . . . . Best Case 100 km Relative Velocity . . . . . . . . . . . . . . . Cycle Slip Visualization. . . . . . . . . . . . . . . . . . . . . . Best Case Cycle Slip Profile. . . . . . . . . . . . . . . . . . . . GPS Clock Offset Estimation. . . . . . . . . . . . . . . . . . . GPS Frequency Offset Estimation. . . . . . . . . . . . . . . . Ensemble Mean Relative Position . . . . . . . . . . . . . . . . Ensemble Mean Relative Velocity . . . . . . . . . . . . . . . . Ensemble Mean Tracking Error with 300 Sec Window . . . . . Ensemble Mean Actuation Effort in the u, v , w frame . . . . . ^ ^ ^ Full Simulation Environment Ensemble Mean Relative Position Full Simulation Environment Ensemble Mean Relative Velocity Ensemble Mean Relative Position . . . . . . . . . . . . . . . . Ensemble Mean Relative Velocity . . . . . . . . . . . . . . . . Ensemble Mean Tracking Error . . . . . . . . . . . . . . . . . Ensemble Mean Actuation Effort . . . . . . . . . . . . . . . . 77 77 78 81 81 82 82 83 83 84 84 85 85 87 88 90 90 91 91 93 94 96 96 97 97 122 124 124 C.1 Measurement Residuals . . . . . . . . . . . . . . . . . . . . . . C.2 Absolute Position Covariances Updates . . . . . . . . . . . . . C.3 Relative State Covariance Updates . . . . . . . . . . . . . . . xv Chapter 1 Introduction 1.1 Background Formation flying can be defined as the use of two or more vehicles moving with specified dynamics for scientific, military, or commercial purposes. Specific examples could be the docking of the Space Shuttle to the Space Station, or the Space Shuttle capturing the Hubble Space Telescope. While there have been numerous examples of formations with humans in-theloop, the number of unmanned satellite missions have been relatively few. The disparity between these two mission types is due, in part, to the high cost of manufacturing and maintaining multiple satellites. In recent years there has been substantial interest in autonomous satellite formations, driven by the new technologies that enable smaller and cheaper spacecraft. The advantages of using formations are numerous, and include the ability to distribute instruments among several spacecraft, thus eliminating the possibility of a single point failure. Additionally, by placing sensors on multiple platforms, new classes of missions that require simultaneous sensing from different angles become possible. Table 1.1 lists several formation missions that have been considered by the National Aeronautics and Space Administration (NASA). Two of the missions listed, TECHSAT-21 and the AQUA mission, are designed to be Earth observing formations whose mission plans lie in tandem with the work documented in this paper. 1 Table 1.1: Actual and Proposed Formation Flying Missions.[22] Projected Launch Mission 2000 2004 2004 2006 2008 2011 15+ * Launched in March 2002 The sophistication and complexity of the missions grow as one descends the time line of Table 1.1. To ensure proper functionality of these formations, our simulation ability must keep pace with the satellite design technology. This work is designed to open the door to new formation testing strategies that will accommodate flight hardware in a real-time setting. 1.1.1 TECHSAT-21 The Technology Satellite of the 21st Century (TECHSAT-21) is a formation being driven by the Air Force Research Laboratory (AFRL). The Air Force hopes to fully take advantage of the benefits of formation flying by constructing an array of small, lightweight satellites that can be reconfigured to allow for multiple mission profiles. One of the keys to the mission is having several platforms flying in a precise formation providing simultaneous measurements from multiple positions [25]. While the proof-of-concept mission may only consist of three satellites, the proposed Techsat-21 mission will be 2 Gravity Recovery and Climate Recovery (GRACE)* Techsat-21/AFRL ESSP-3-Cena (w/ Aqua) Magnetosphere Imaging Constellation (MAGIC) Laser Interferometer Space Antenna (LISA) Terrestrial Planet Finder (TPF) Planet Imager (PI) comprised of 35 clusters of 8 satellites each [17]. The simulation system demonstrated in this work would be ideal for the modelling of a control system with real-time hardware integration. 1.1.2 A-TRAIN The "A-Train" is an Earth observing formation that will consist of five satellites: Aqua, CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations), CloudSat, Parasol (Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar), and Aura. Mission constraints require an in-line formation with intersatellite distances ranging from 100 km out to 6000 km [23]. The long baseline relative navigation filter documented in this work could provide a good starting point for developing a formation-wide autonomous control system. 1.2 Previous Work The number of formation flying papers has increased dramatically over the past five years. There are also numerous works available detailing the relative navigation problem, which has been studied for decades. This section focuses on some of the more recent works that are directly relevant to this thesis. Several papers have been written outlining the basic fundamentals of formation flying. Bauer defines many of the key developments required for formations, including sensor development and formation control strategies [2]. Sabol expands on this work by defining some basic satellite formation types and the corresponding vehicle dynamics [24]. 3 Ebinuma demonstrated the use of GPS for rendezvous by developing a real-time hardware-in-the-loop feedback control system at the University of Texas at Austin. The filter in this work estimated the absolute state of the each satellite, and then differenced these states to generate the relative separation. His development of the mathematical models is rigorous [9]. Ebinuma, Bishop, and Lightsey published a condensed account of these results[10]. In another work, Ebinuma, Montenbruck, and Lightsey demonstrated 20 cm relative navigation accuracies for a pair of satellites separated by 10 km [11]. They found that at this separation distance, only minimal gains were achieved when their Kalman filter employed a 10x10 gravity model as opposed to a simple model based upon Hill's equations. Binning used single-differenced carrier phase measurements and an integer ambiguity search algorithm to explore navigation accuracies for 5 and 50 km baselines. His hardware-in-the-loop simulations generated an accuracy of 3 cm over a 50 km baseline using a dual frequency widelane integer resolution method. [3]. Wolfe and Speyer examined the relative navigation problem with vehicle separations near 100 km. They used differential carrier phase measurements coupled with a Wald test to aid in the resolution of the integer ambiguities. Their procedure yielded accuracies near 5 cm for a GPS carrier phase signal with 1 cm of noise [29]. Gramling discusses the extension of formation flying algorithms from low Earth orbits to high Earth orbits using GPS signals when available. The algorithms detailed describe relative position estimation by differencing point solutions as well as by differencing filtered solutions [14]. Further work by Long compares accuracies in relative navigation for eccentric medium and 4 high Earth orbits. Their filter uses GPS measurements, cross link data, and celestial object measurements to compute the states [19]. How has published several papers tackling the fundamental issues behind formation design and control. Busse and How demonstrated sub-centimeter navigation accuracy for a 1 km baseline using an adaptive Kalman filter and differential carrier phase measurements [5]. Several papers were written in conjunction with the EO-1/Landsat-7 mission. Folta provides a good overview of the mission as well as the results from the Autocon flight control system [13]. 1.3 Research Contributions The primary objective of this research was to construct a simulation tool that could be used to evaluate estimation and control systems for large formations in real-time utilizing actual hardware. In addition to the simulation environment, a basic Kalman filter was developed for use in a wide variety of satellite formations. A simple linear-quadratic control system was also employed in order to close the loop. 1.3.1 Small Formation to Large Formation Generalization In order to understand the differences between the small and large formation cases, research was done to determine the strengths and weaknesses of both types. A literature review was conducted in order to determine the best way to construct a general Kalman filter that could leverage the simplicity of the small formation plant against the increasingly complex relative dynamics encountered with longer baselines. 5 1.3.2 Extended Kalman Filter Design An extended Kalman filter (EKF) was constructed for use in the Ex- tended Formation Flying Testbed (EFFTB) developed in this work. The filter uses the local Global Positioning System (GPS) measurements for its own absolute state estimate, and measurements from a distant receiver for relative navigation. The filter utilizes all three GPS observables: the pseudorange, the carrier phase, and the Doppler shift. To reduce the effect of common error sources, the relative measurements are all double-differenced prior to processing. Special care was taken to include error sources such as the differential ionospheric delay and the differential transmit time between receivers. These error sources have a greater impact on missions with longer baselines. 1.3.3 Demonstration of a Generalized Kalman Filter The Kalman filter design was tested over a variety of baselines ranging from 1 km to 500 km. The 1 km simulation yielded sub-centimeter accuracies, similar to works by Ebimuma and Binning [11] [3]. By slightly shifting the gains, the filter exhibited 3-4 centimeter steady-state error levels at 100 km spacing, comparable to Speyer's work [29]. It was determined that the 500 km baseline scenario exceeded the threshold of the filter's ability due to the breakdown of the models implemented. By implementing a simple gain-scheduling algorithm into the flight computer, this single filter could allow for a group of satellites launched together to separate from the launch vehicle and position themselves up to 400 km away with less than 50 cm of error. 6 1.3.4 Construction of the Extended Formation Flying Testbed A real-time, closed-loop, hardware in-the-loop simulation tool was de- veloped to model formations with several satellites. Initially, this tool consisted of a generic communications package that could accept any number of hardware or software measurement sources. The observer and controller portions of the code were added as subroutines, making it easy to swap in and out the different modules. This allows for rapid implementation of different estimator and controller combinations. The EFFTB could be used for Low Earth Orbits (LEO), High Earth Orbits (HEO), or even deep space missions, with the appropriate propagator and measurement source substitutions. 1.3.5 GPS Specific Instance of the EFFTB The general EFFTB architecture was adapted to use Spirent GPS simulators and Orion type GPS receivers in order to simulate a large LEO formation. The propagation machine was configured to feed the simulator reference data at 10Hz, and subroutines were developed to read in and process the GPS data correctly. The testbed was connected to a similar configuration of machines at Goddard Space Flight Center (GSFC) to generate a multiple satellite formation transferring real-time measurements over the Internet. 1.4 Overview Chapter 2 details the different time systems and reference frames used in this work, and includes the mathematical descriptions on how to convert back and forth between these systems. Chapter 3 describes the GPS measurement models and differentiates 7 between the measurements used for the absolute and relative navigation portions. The GPS satellite propagation models are also outlined. Chapter 4 presents the extended Kalman filter equations and the dynamic models used. Chapter 5 diagrams the equations and implementation of the linearquadratic controller. Chapter 6 outlines the hardware and software development of the Extended Formation Flying Testbed. Chapter 7 describes the formation considered for this study, and presents the results of the real-time simulation environment. Preliminary results validating the filter and controller operation, as well as the validation of the testbed architecture, are also presented. Chapter 8 summarizes the research and includes conclusions and topics for possible future work. 8 Chapter 2 Reference Systems 2.1 Time Systems There are three basic time systems utilized in this work. Coordinated Universal Time (UTC) is the basic time used for computations and integration. In order to convert from the Earth Centered Inertial (ECI) coordinate frame and the Earth Centered Earth Fixed (ECEF), Greenwich Sidereal Time (GST) is used. Finally, GPS time is used in conjunction with the corresponding ephemerides to propagate the GPS satellites forward in time. The following is a short list defining each system and the transformation equations. 2.1.1 Solar and Sidereal Times Universal time (UT) is a solar time system defined by the Greenwich Hour Angle (GHA) augmented by 12 hours of a fictitious sun uniformly orbiting in the equatorial plane [15]. Sidereal time is defined by the hour angle of the vernal equinox, and the sidereal time associated with the Greenwich Meridian is called the Greenwich Sidereal Time [27]. Since the Earth's rotation rate is not constant, neither of these time systems are uniform. In order to correct for some of this non-uniformity, the UT1 time system was introduced, which accounts for the polar motion of the Earth. 9 2.1.2 Atomic time In order to meet the demands of the scientific community, the atomic time standard was introduced. International Atomic Time (IAT) is based on the quantum transitions of the Cesium-133 atom. The atomic second is based upon a fixed number of these cycles. 2.1.3 Coordinated Universal Time In 1972, UTC was introduced to tie the irregularity of the UT1 system to the rigidity of the IAT system. In order to align the two systems, leap seconds are added to both the UTC and IAT. UTC always differs from IAT by an integer number of seconds, such that IAT = U T C + 1.000n seconds (2.1) The United States Naval Observatory (USNO) is tasked with the determination of this integer n, as well as the UT1-UTC corrections, and publishes them in the International Earth Rotation Service (IERS) Bulletin-A. 2.1.4 GPS Time GPS time is an atomic time system that is based on the number of weeks elapsed from the GPS epoch, as well as the number of seconds into the current GPS week. The GPS standard epoch date is January 6, 1980 at midnight. GPS time is related to the IAT time system by a nominal 19 second offset such that: IAT = GP S + 19.000 seconds Substituting equation 2.2 into 2.1 results in: U T C = GP S + 19.000 - 1.000n seconds 10 (2.3) (2.2) According to the USNO, in 2003 the integer value was n = 32, which yields a difference between UTC and GPS time of 13 seconds. 2.1.5 GPS Roll Over In order to accurately convert between GPS time and other time systems, the GPS week must be accounted for. The week number starts at the standard GPS epoch and is incremented until it reaches 1024, at which time it is reset to 0. The first rollover occurred at midnight on August 21, 1999. The GPS week number can be obtained from the first subframe of the GPS navigation message. 2.2 Reference Frames There are three different references systems used throughout this work. The following subsections define those system and provide the appropriate transformations. 2.2.1 Earth Centered Inertial Reference Frame (ECI) The J2000 reference system is used predominantly throughout this dissertation. This system is a geocentric equatorial frame defined at the epoch January 1, 2000. The reference plane is constructed from the mean equatorial plane of the Earth at the J2000 epoch. All of the orbit propagation and the generation of the state transition matrix is performed in this frame. Table 2.1 describes the specifics of the coordinate system. 11 Table 2.1: J2000 Reference System Definition X Axis Y Axis Z Axis Mean vernal equinox at epoch Normal to x and z to form a right-handed system Normal to the mean equatorial plane of J2000 in the direction of the Earth's mean spin 2.2.2 Earth Centered Earth Fixed Reference Frame (ECEF) The ECEF frame differs from the ECI frame in that it is fixed in the Earth, and the primary axis is always aligned with the Greenwich Meridian. Since it is rotating, it is not an inertial frame. The ECEF frame is used for communicating with the GPS signal generator. Table 2.2 describes the ECEF system unit vectors. Table 2.2: ECEF Reference System Definition X Axis Y Axis Z Axis Intersection of the Greenwich Meridian with the equator Normal to x and z to form a right-handed system Adopted geographic pole direction 2.2.3 Spacecraft Centered Frame The U,V,W frame is a spacecraft centered frame that is used only for presenting results. Table 2.3 describes the coordinate system. The unit vectors comprising the UVW frame, denoted u, v , and w can be generated from a ^ ^ ^ J2000 position using equations 2.4. r |r| (2.4) u= ^ v =w u ^ ^ ^ 12 w= ^ where: r r |r r| r = the position vector in the J2000 frame r = the velocity vector in the J2000 frame Table 2.3: UVW Reference System Definition U Axis V Axis W Axis Points along the radius vector to the satellite Normal to U and W to complete a right hand system Direction of the orbital angular momentum 13 Chapter 3 GPS Models and Measurements 3.1 GPS Observables There are three main observables provided by the GPS system, and all three are utilized within this work. This section details the pseudorange, carrier phase, and doppler measurements and the errors associated with them. 3.1.1 Pseudorange Measurement The pseudorange is simply the range between the GPS receiver and the GPS satellite with offsets due to clock uncertainties and other measurement error sources. For the measurement equations, superscripts will be used for the GPS satellites, and subscripts will indicate the GPS receiver. Since the true time is not known, the satellite transmit time (tS ) and the receiver acquire time (tR ) are given by: tS = T S + tS tR = TR + tR (3.1) (3.2) where the capital T's are indicative of the true times, and tS and tR denote the satellite and receiver clock error respectively. The pseudorange measurement is simply the speed of light C times the difference between the receiver clock and the satellite clock. P (tR ) = C(tR - tS ) 14 (3.3) Substituting equations 3.1 and 3.2 into 3.3 and rearranging yields: P (tR ) = C(TR - T S ) - C(tR - tS ) (3.4) The first term on the right hand side of the equality is simply the true range between the GPS satellite and receiver, which will be denoted as S . Equation R 3.4 simplifies to: P (tR ) = S - C(tR - tS ) R (3.5) The measurement generated in the receiver is at the acquired time tR not the true time TR . Since the true time is not known, the geometric distance is linearized about the acquired time by [9]: S = S (tR - tR ) R R = S (tR ) - S (tR )tR R R Inserting this result into the equation 3.5 gives: P (tR ) = S (tR ) + (C - S )tR - CtS R R (3.7) (3.6) Equation 3.7 is the mathematical equation for the pseudorange in a vacuum with no other errors. When these additional error sources are accounted for, the equation becomes: P (tR ) = S + (C - S )tR - CtS + E + Iono + R R where: E = User Range Error Iono = Ionospheric Delay = Pseudorange measurement noise 15 (3.8) Three possible error sources have been omitted from the equation, but are documented here for completeness. Since this work only examines spacecraft orbits, the tropospheric delay can be neglected. Selective ability, the slewing of the GPS clock and ephemeris data, was turned off on May 1, 2000, and is therefore not accounted for. Finally the multipath error is not accounted for, as this error is highly dependent on the application in which the receiver is used. 3.1.1.1 User Range Error The user range error is a combination of two error sources: the GPS satellite clock correction error, and GPS satellite ephemeris error. Each GPS satellite has an atomic clock on board for precise timekeeping. However, there are slight drifts and biases to these clocks which are tracked by the GPS ground segment. These offsets are broadcast in the form of curve fit coefficients in the GPS message. However, since only terms through the first order are broadcast, there is a slight error due to truncation. The GPS ground segment also provides the GPS satellite ephemerides which are broadcast in the almanac file. Slight errors in the satellite's orbital elements results in satellite position errors when propagated. 3.1.1.2 Ionospheric Delay The ionosphere, which extends from approximately 50 km to 1000 km above the Earth's surface, acts as a dispersive medium with respect to the GPS signals. The net affect on the measurement is that the pseudorange mea- 16 surements are delayed, and therefore measured too long. The opposite affect is observed the carrier phase measurements. Due to the variable nature of the ionosphere, precise modelling of these group delays and advances is difficult. 3.1.2 Carrier Phase Measurement The carrier phase measurement is a direct measurement of the phase of the received signal. If we denote the phase of the received carrier signal as S and the phase of the reference carrier signal generated by the receiver as R , we can write [15]: R (tR ) = R (TR ) + R tR S (tS ) = S (T S ) + S tS (3.9) (3.10) The time t is an epoch reckoned from an initial epoch t0 = 0. The carrier beat phase S (tR ) can be defined as: R S (tR ) = S (tS ) - R (tR ) R = (S (T S ) - R (TR )) + f tS - f tR (3.11) where the L1 carrier frequency has been replaced with f, and the appropriate terms from equations 3.9 and 3.10 have substituted in. The term f /C can be used to convert the geometric distance into cycles. Doing so allows the true 17 beat phase to be represented by: f S (T S ) - R (TR ) = - S (TR ) C R f R = - (S (tR ) - S (tR )tR ) C R (3.12) Substituting 3.12 into the first term in 3.11 gives the beat phase f S (tR ) S (tR ) = - S (tR ) - f 1 - R tR + f tS R R C C (3.13) When the receiver is switched on the initial epoch is set, and the beat phase is measured. The receiver is unable to determine the integer number of whole phase cycles between the receiver and GPS satellite, only the portion of the cycle it has measured. This cycle ambiguity N, remains a constant as long as signal lock is not lost to the GPS satellite. Therefore, at some epoch t, the beat phase can be represented by: S (t) = S |t0 + N R R t (3.14) The term S represents the phase measured between the times t and R t0 . The carrier phase measurement can be realized by substituting 3.13 into 3.14, denoting the observation as , and assigning it a negative sign: S (tR ) = R f S S (tR ) R (tR ) + f (1 - R )tR - f tS + N C C (3.15) The carrier phase measurement is often represented in terms of range, so we can scale it by = C f where is the carrier wavelength. (3.16) R S (tR ) = S (tR ) + (C - S (tR ))tR - CtS + N R R 18 As with the pseudorange equation, when the error sources are added in, the carrier phase becomes: S (tR ) = S (tR ) + (C - S (tR ))tR - CtS + N + E - Iono + (3.17) R R R where is the carrier phase measurement noise. It is worth noting that the only qualitative differences between the pseudorange measurement and the carrier phase measurement is the addition of the integer ambiguity term at the end, and the opposite sign of the ionosphere delay. 3.1.3 Doppler Measurement Since the transmitting satellite in the GPS constellation is always in motion relative to the receiving body, the received frequency will be Doppler shifted. In other words, the received frequency will differ from the transmitted frequency due to the relative motion. The equation for the Doppler measurement scaled to units of range can be obtained by differentiating equation 3.17 to obtain: D(t) = = + C(fR - f S ) - Iono + (3.18) Where term is the noise in the Doppler measurement. The rate of change of the ionospheric delay is often neglected in practical work. 3.2 Double-Differenced Observables Many of the error sources described above are common mode, and will cancel out when differenced, leaving a more accurate relative measurement. For this reason, double-differenced carrier phase and doppler measurements are utilized for the relative positioning portion of the filter. A schematic of the double-difference process is diagramed in Figure 3.1. In the figure, SAT 19 A receives measurements from GPS A and GPS B at the same epoch. SAT A subtracts these measurements, creating a single-difference. At the same epoch, SAT B performs the same single-difference. The double-difference is performed by subtracting the two single-differences. The double-differenced carrier phase and doppler measurements will be developed mathematically in the following subsections. Figure 3.1: Double-Difference. 20 3.2.1 Double-Differenced Carrier Phase If the carrier phase measurement between SAT A and GPS I is denoted as I , the measurement can be represented by: A I I (t) = I + (c - I )tR,A - CtS,I + E,I + Iono,IA + NA + I A A A A similar equation can be written for SAT A and GPS J : J J (t) = J + (c - J )tR,A - CtS,J + E,J + Iono,JA + NA + J A A A Subtracting these equations yields the single difference equation for SAT A: I J I J IJ = (I -J )-CtS,I +CtS,J +(Iono,IA -Iono,JA )+(NA -NA )+(A -A ) A A A The receiver clock error terms have dropped out due to the differencing. A similar equation can be written for the SAT B receiver: I J I J IJ = (I -J )-CtS,I +CtS,J +(Iono,IB -Iono,JB )+(NB -NB )+(B -B ) B B B performing the second differencing yields: IJ = ((I - J ) - (I - J )) + ((Iono,IA - Iono,JA ) - (Iono,IB Iono,JB )) AB A A B B I J I J +(NA - NA ) - ((NAB - NB )) + where is the double differenced random errors. The GPS satellite clock errors and the user range errors have dropped out, leaving the double differenced ionosphere as the sole non-random error. For short intersatellite distances, the path between the GPS satellites and the receivers are close enough that this error is negligible. 21 3.2.2 Double-Differenced Doppler Following the above derivation, the Doppler equation is written for the measurement between SAT A and GPS I. I I Iono,A + A DA (t) = I + C(fR,A - f S,I ) - I A (3.19) For the SAT A and GPS J combination: J J Iono,A + A DA (t) = J + C(fR,A - f S,J ) - J A (3.20) Taking the single-difference yields: IJ I J J DA (t) = (I - J ) - C(f S,I + f S,J ) - (I A A Iono,A - Iono,A ) + (A - A ) Writing the complementary equation for SAT B and subtracting gives the double-differenced Doppler equation: IJ Iono,A - J Iono,A ) DAB (t) = ((I - J ) - (I - J )) - ((I A A B B I J I J J -(I Iono,B - Iono,B )) + (A - A ) - (B - B ) 3.2.3 GPS Satellite Propagation In order to predict the location of the GPS satellites, the GPS receiver uses the almanac data which is superimposed upon the GPS signal. The almanac data is updated at least every six days and contains orbital ephemerides and satellite clock corrections. Table 3.1 provides an overview of the almanac data elements [15]. The propagation of the elements from the almanac is straightforward, and the ECEF satellite coordinates at time t can be readily obtained using the following algorithm [20]. 22 Table 3.1: Almanac Information Parameter Description ID Week ta a e M0 i a0 a1 Satellite PRN number Current GPS week Reference epoch in seconds within current GPS week Square root of the semi major axis Eccentricity Mean anomaly at reference epoch Argument of Perigee Inclination Longitude of the node at the weekly epoch Drift of the right ascention of the node Satellite clock offset Satellite clock drift Determine the instantaneous longitude of the ascending node: = 0 + (t - ta ) - E (t - t0 ) where t0 denotes the start of the current GPS week, and E is the rotation rate of the Earth. Solve for the eccentric anomaly E using Kepler's equation E - e sin E = M0 + GME (t - ta ) a3 (3.21) The ECEF position is obtained via the following rotations: a(cos(E) - e) = Rz (-)Rx (-i)Rz (-) a 1 - e2 sin(E) 0 rECEF where: (3.22) 23 cos (-) sin (-) 0 Rz (-) = - sin (-) cos (-) 0 0 0 1 1 0 0 Rx (-i) = 0 cos (-i) sin (-i) 0 - sin(-i) cos (-i) cos (-) sin (-) 0 Rz (-) = - sin (-) cos (-) 0 0 0 1 Special care must be taken to ensure that the geometric range used in the measurement equations is the range generated using the receiver position at tR and the GPS satellite position at time tS . This is accomplished by applying the following iterative scheme [9]. The signal transmission time is solved using: tS,n+1 = tR - |rS (tS,n ) - rR (tR )| C where rS , rR are the GPS and receiver states respectively and tS,n is the nth approximation for tS , the signal transmit time from the GPS satellite. The iteration continues until: |tS,n+1 - tS,n | < The default value for the tolerance required for convergence. is 1e-12 which yields less than 1 mm of accuracy for the range measurement. Usually no more than 3 iterations are 24 Chapter 4 Kalman Filter Design 4.1 General Formation Description The design of any particular formation is mission dependent, and sev- eral very interesting formation designs have been developed. In order to generalize this work, it is assumed that each satellite will generate a relative navigation solution with respect to one common satellite. This vehicle, which will be labelled the master satellite, shares its measurements with each of the remaining satellites, which are designated as nodes[5]. The nodes use these measurements to generate a relative state estimate at each epoch. While this assumption may be impractical in some missions, many formations, such as a multi-vehicle space telescope, emphasize the relative distances from some central craft. One of the advantages in this scenario is that, for estimation purposes, the formation flying problem breaks down into several relative navigation problems that can be examined independently. As described in Section 4.6, all of the measurement equations are rewritten to remove the necessity of estimating the absolute position of the hub, though its state is reconstructed from the node state and the relative position at every epoch. Figure 4.1 details the layout of the absolute and relative vectors. 25 r B GPS r r A GPS rN r H Figure 4.1: Problem Setup. 4.2 Filter State The estimation state vector X at each node can be written as: r r b X= b N B 26 where: r = Absolute state of the node r = Relative state between the node and hub b = Receiver clock offset b = Receiver clock frequency offset N = Double-Differenced integer ambiguities B = GPS clock offset Each of the states are developed more completely in the following subsections. 4.2.1 Absolute And Relative States Each node satellite estimates its own absolute state as well as the relative state to the hub: = x y z x y z x y z x y z r r where each of these elements are realized in the ECI reference frame. For simplicity, the GPS receiver is assumed to be located at the center of mass of each satellite. While estimation of the attitude is crucial to a formation mission, it will be left for future research. 27 4.2.2 Receiver Clock and Frequency Offset State The GPS receiver typically has a quartz crystal clock on board, making it much less accurate than the ones in the GPS satellites. In order to compensate for these clock errors in the pseudorange and doppler measurements, the receiver clock and frequency offsets must be estimated. This may be achieved by means of a second order Gauss-Markov process. The Orion GPS receivers used for this research are enabled with a clock-steering algorithm. This code slightly changes the TIC, the receiver's fundamental time unit, so that the GPS receiver is always aligned to the GPS time. This makes exchanging measurements much more feasible. Since the receiver clock is no longer drifting freely, the commonly used receiver clock model no longer applies. Instantaneous changes in the receiver clock are not reflected in the clock's frequency, creating a situation where the two states vary independently. It was empirically determined that the clock error could be successfully modelled using a constant offset b and a constant frequency offset b. These states are modelled as uncorrelated random walk processes. 4.2.3 Double-Differenced Ambiguity State The double-differenced ambiguity state 12 NAB 13 NAB N = . . . NAB 1(n-1) is defined by: where n is the number of GPS satellites tracked by both GPS receivers. 28 4.2.4 GPS Clock Offset As described in Section 3.1.1.1, there is a slight error in the transmitted GPS clock model due to the truncation of the model down to the first order. This error can be accurately modelled as a constant bias for each GPS satellite tracked. B= 1 BA 2 BA . . . n BA This error is common mode, and cancels out for the double-differenced measurements. Consequently, the B state is needed for the pseudorange measurement only. 4.3 Kalman Filter The feedback control system that makes up the software portion of the EFFTB requires the implementation of an estimator and a controller. The backbone of the software is general enough to allow for many combinations of those elements. Many different estimator architectures have been developed for the relative navigation problem, such as the centralized, decentralized, and partially centralized formation. Instead of developing a filter that provides the "best" results for a single baseline length, a different approach was taken in this dissertation. A Kalman filter has been developed which yields good results over a wide variety of baselines. Thus some accuracy has been traded for a more versatile filter design. Each satellite in the formation will run its own EKF for absolute and relative state estimation. The filter functions as a weighting scheme between the current measurement set and the dynamic model. The filtering can be 29 broken up into several distinct stages for clarity, and the mathematical descriptions will be developed accordingly. 4.3.1 State Vector Construction Due to the chosen estimation strategy, the number of estimated states changes with the addition or removal of each GPS satellite. Therefore, if the number of GPS satellites tracked changes, the state vector must be resized to account for the appropriate number of tracked satellites. The error-covariance matrix must be rescaled accordingly as well. The hub state is also reconstructed from the current state parameters. 4.3.2 Propagation The propagation step is used to bridge the gap in time between successive measurements. In an extended Kalman filter, the state is propagated by the direct integration of the non-linear dynamic equations. The state error covariance is propagated forward using a state transition matrix. If we denote X k (+) as the state vector at time k after the measurement, and X k+1 (-) as the state vector at time k+1 before the measurement, then the propagation algorithm is presented as: 1. Propagate the position and velocity vectors for the absolute states by integrating the non-linear dynamical model, which is provided in section 4.4.1. The states are integrated using a Runge-Kutta 4th order numerical integrator. = aE + ap r (4.1) where aE is the acceleration due to the Earth's gravity and ap is the 30 desired perturbative forces. 2. The relative state is not propagated directly. Instead, the reconstructed hub state is propagated forward to the current epoch via Equation 4.1. The relative state is then formed by subtracting the current node state from the hub state. 3. The remaining terms in the state are modelled as constants, so their values are changed only by the inclusion of process noise during the propagation step. 4. The state error covariance matrix contains the estimates for the closeness of the fit with the actual observations [27]. The diagonal terms are the variances of the estimate of the state parameters. The square roots of the variances are the standard deviations for each element of the state space. For a properly working filter, the errors in the problem must be accurately reflected by the entries in the covariance matrix. The error covariance must be propagated during this step. The state error covariance matrix, Pk at time tk is defined by the expectation operator as: Pk = E Xk - Xk Xk - Xk T (4.2) The update equation for the covariance is given by: Pk+1 (-) = (tk+1,k )Pk (+)T (tk+1,k ) + Qk+1 (4.3) The generation of the state transition matrix (tk+1,k ) and the covariance propagation strategy are discussed in Section 4.4.2. The process noise matrix Q is used as a tuning parameter. 31 4.3.3 Initialization Each state of the EKF must be initialized before it is processed the first time. The node state is initialized from the GPS position fix solution at the initial epoch. The initial hub position is also set to its position fix solution, and the initial relative state is set to the difference between the two. The double-differenced ambiguity term is initialized as the difference between the double-differenced pseudorange and the double-differenced carrier phase measurements. The remainder of the states are initialized to zero. 4.3.4 Measurement Update Once the state vector and error covariance matrix are propagated to the current time, the updated state vector X k (+) and the updated covariance matrix Pk (+) can be computed. The update requires the construction of the measurement noise covariance matrix Rk . The measurement equation matrix G is covered in section 5.5. 1. The predicted measurement is generated by evaluating the non-linear measurement models along the current state: (Y k ) = G(X, tk ) (4.4) The residual vector is computed as the difference between the measurement and the predicted measurement: z k = Yk - (Y k ) (4.5) where the hat denotes an estimated quantity. The Hk matrix is generated by taking the partial derivative of the measurement equations with 32 respect to the state. These partial derivatives are presented in Appendix B. Hk = G | X X=X k (-) (4.6) 2. The error covariance is updated using the numerically stable Joseph formula: Pk (+) = (I - Kk Hk )Pk (-)(I - Kk Hk )T + Kk Rk KT k where the Kalman gain is given by: Kk = Pk (-)HT (Hk Pk (-)HT + Rk )-1 k k The state update also uses the Kalman gain: X k (+) = X k (-) + Kk z (4.9) (4.8) (4.7) 4.4 Dynamic Model and State Transition Matrix In an EKF, the state is propagated forward using the non-linear dif- ferential equations. The covariance matrix, however, uses a state transition matrix for propagation. The following sections provides the full non-linear differential equations and the state transition matrix structure. 4.4.1 Dynamic Model One of the potential benefits for formation flying is to be able to create several smaller, cheaper, less-capable nodes to replace large monolithic satellites. With this philosophy in mind, the dynamic model was intentionally left relatively simple, which reduces the computational burden. The acceleration 33 equations contain the gravitational attraction of the Earth, the J2 zonal harmonic, and the atmospheric drag effect with an exponential air density model. These equations are detailed in Equation 4.10. 3 x x = - 3 1 - J2 r 2 3 y y = - 3 1 - J2 r 2 z 3 z = - 3 1 - J2 r 2 RE r RE r RE r 2 5 2 z2 -1 r2 z2 -1 r2 z2 -3 r2 1 A - CD V (x + E y) 2 M 1 A - CD V (y - E x)(4.10) 2 M 1 A - CD V z 2 M 5 2 5 where E is the rotation rate of the Earth, RE is the Earth's radius, Cd is the drag coefficient of the satellite, m is the satellite's mass, and A is the instantaneous cross-sectional area. The magnitude of the relative wind, V, can be computed from: V = 2 2 x2 + 2E xy + E y 2 + y 2 - 2E yx + E y 2 + z 2 (4.11) The value of , the atmospheric density, is generated using an exponentially varying model [27]: = 0 e-(r-r0 /H) Here: 0 = 3.725x10-12 (kg/m3 ) r = current altitude(m) r0 = reference altitude = 400000(m) H = scale height = 58515(m) (4.12) 34 4.4.2 State Transition Matrix While much of the state is propagated directly via the dynamic equa- tions, the covariance matrix must be propagated forward through the use of a state transition matrix. The state transition matrix used to propagate between the two times tk and tk+1 satisfies: (tk+1 , tk ) = A(tk )(tk+1 , tk ) (tk , tk ) = I F A(tk ) = | X X=X k (-) (4.13) (4.14) (4.15) where I is the appropriately sized identity matrix and F is the set of dynamical equations. The state transition matrix has the block diagonal form: = Abs 0 0 0 Rel 0 . . . . . b . 0 N 0 0 0 B (4.16) each submatrix is described in the following subsections. 4.4.2.1 Absolute State Acceleration Partial Derivatives Abs The acceleration partial derivatives are found by taking the partial derivatives of equation 4.10 with respect to the state. The equations for development of the absolute partial derivatives are provided in Appendix A. 35 4.4.2.2 Relative State Acceleration Partial Derivatives Rel The relative state is not propagated directly, and the following approximation is made in order to generate the state transition matrix. The relative state at time tk+1 can be written as: r (tk+1 ) = N (tk+1,k )xN - H (tk+1,k )xH (4.17) If we assume that, for separations under a certain threshold, the accelerations felt by the hub spacecraft are approximately equal to those of the node spacecraft, then we can say that the relative state transition matrix block can be generated from the partial differential equations of either the node or the hub. This approximation is validated by evaluating the partial differential equations numerically and comparing them to the non-linear models. The results of this analysis are presented in Appendix C. 4.4.2.3 Receiver Clock States Partial Derivatives b The receiver clock and frequency offsets are modelled as random-walk processes. The contribution to the state transition matrix is given by: b (tk+1 , tk ) = 1 0 0 1 (4.18) 4.4.2.4 Double-Differenced Ambiguity Partial Derivatives N The double-differenced ambiguity terms are modelled as random-walk processes. Though theory describes these ambiguities as integers, their values are not constrained as such. Their partial derivative submatrix has the form: N (tk+1 , tk ) = I 36 (4.19) For N GPS satellites, there are (N-1) double-differenced ambiguity terms to estimate at every epoch. 4.4.2.5 GPS Clock Error Model B The GPS satellite clock error is modelled as a random-walk process. B (tk+1 , tk ) = I Each GPS satellite tracked requires one offset parameter. (4.20) 4.5 Process Noise Matrix As the time progresses into the simulation, the covariance matrix will begin to converge. The filter will be increasingly less sensitive to the new measurements. To overcome this effect, process noise is added in the form of the process noise matrix Q. This matrix also accounts for the errors in the dynamical model. The Q matrix can be represented in the block diagonal form: Qk = QR 0 0 0 0 Qb 0 . . . . QB . . 0 QN (4.21) Each submatrix is described below. 37 4.5.1 Absolute and Relative State Process Noise The process noise covariance matrix for the absolute and relative state is defined by: QR = tk+1 tk (tk+1 , )S( )T (tk+1 , )d (4.22) where S(t) is the process noise spectral density matrix. If the state transition matrix can be assumed constant over the time span t, the process noise matrix can be written in discrete form: QR = S( ) (4.23) and the absolute and relative state process noise contributions are defined as: QRa 0 (4.24) 0 QRr correspond to the absolute state and relative state process QR = where QRa and QRr noise respectively. These matrices can be written as: = 2 xABS 0 0 0 0 0 2 0 yABS 0 0 0 2 0 0 zABS 0 0 0 2 0 0 0 xABS 0 0 2 0 0 0 0 yABS 0 2 0 0 0 0 0 zABS QRa (4.25) and: = 2 xREL 0 0 0 0 0 2 0 yREL 0 0 0 2 0 0 zREL 0 0 0 2 0 0 0 xREL 0 0 2 0 0 0 0 yREL 0 2 0 0 0 0 0 zREL QRr (4.26) 38 4.5.2 GPS Receiver Clock Process Noise The process noise covariance for the receiver clock model is: Qb = 2 b 0 2 0 b (4.27) 2 2 The b and b represent the variance for the clock and frequency offsets. This model is valid only for the receiver used in this experiment, due to the decoupling of the clock states by the clock steering algorithm. 4.5.3 Double-Differenced Ambiguities Process Noise The process noise covariance block for the double-differenced ambiguity terms have the form: 2 QN = IN (4.28) 2 where I is an (n - 1) (n - 1) identity matrix, and N is the variance for the double differenced ambiguities. 4.5.4 GPS Satellite Clock Bias The process noise covariance entries for the GPS satellite clock bias also have a diagonal form: 2 QB = IB (4.29) 2 where I is an (n) (n) identity matrix, and B is the variance of the GPS clock bias. 4.6 GPS Measurement Equations Sections 3.1 and 3.2 describe the basic formulation for the measure- ments utilized by the filter. The partial derivative of these equations, which 39 are required for the measurement update portion of the filter, are provided in Appendix B. 4.6.1 Undifferenced Measurement Equations The undifferenced measurement equations are very similar to the theoretical ones presented in Section 3.1. The pseudorange filter model is: P (tR ) = S + b + B + Iono R (4.30) where B is the scaled GPS clock state offset, b is the scaled receiver clock state offset, and PIono is the pseudorange ionospheric delay. The ionosphere calculation uses a simplified model which consists of a shell of constant thickness and uniform electron distribution. The delay is calculated as [26]: Iono = where: T EC = 2e17 Constant Total Electron Count FC = 1575.42e6 Frequency of the carrier Signal EL = Elevation of the GPS satellite from a local horizon 82.1T EC 2 FC sin (EL) + 0.076 + sin2 (EL) 2 (4.31) The equation for the Doppler model is given by: D(t) = S + b R where b is the receiver clock frequency state. 4.6.2 Double-Differenced Measurement Equations In order to remove the necessity of estimating the absolute state of the hub satellite, the measurement for the relative states are rewritten in terms of 40 (4.32) the node state and the relative state. If we define our system such that: r(t) = rH (t) - rN (t) then the absolute state of the hub can be represented by: rH (t) = r(t) + rN (t) (4.34) (4.33) Plugging this expression for the hub into the double-differenced measurement equations yields the following: IJ (t) = AB (|rA S (t - ) - r(t) - rN (t)|) - (|rA S (t - ) - rN (t)|) - GP GP (|rB S (t - ) - r(t) - rN (t)|) - (|rB S (t - ) - rN (t)|) - GP GP Iono + N where the notation signifies a double-difference such that: A A B B N = (NH - NN ) - (NH - NN ) (4.35) Similarly, the double-differenced Doppler equation has the form: IJ DAB (t) = (rA S (t - ) - r(t) - rN (t)) GP (rA S (t - ) - r(t) - rN (t)) GP - (|rA S (t - ) - r(t) - rN (t)|) GP (rB (t - ) - r(t) - rN (t)) (rB S (t - ) - r(t) - rN (t)) GP S GP (|rB S (t - ) - r(t) - rN (t)|) GP 4.7 Measurement Noise Covariance The measurement noise covariance matrix represents the uncertainty in the measurement model. The measurement noise covariance matrix, R has the block diagonal structure: RP 0 0 0 0 RD 0 0 Rk = 0 0 R 0 0 0 0 RD 41 where RP ,RD ,R ,RD , correspond to the pseudorange, doppler, doubledifferenced carrier phase, and double-differenced doppler measurements respectively. 4.7.1 Pseudorange Measurement The raw pseudorange measurements are inherently uncorrelated, yielding a diagonal block matrix: 2 0 . . RP = . . . . . . . 2 0 2 where is the standard deviation of the pseudorange error. 4.7.2 Doppler Measurement The doppler measurements are also uncorrelated, and have a similar diagonal structure: 2 0 . . RD = . . . . . . . 2 0 2 where is the standard deviation of the doppler error. 4.7.3 Double-Differenced Carrier Phase Measurement While the raw carrier phase measurements are uncorrelated, the double- differenced measurements are correlated with each other. If we write two 42 double-differenced carrier phase measurements in matrix representation: A N A H B AB 1 -1 1 -1 0 0 NH N = AC 1 -1 0 0 1 -1 B NH H C N C H 43 The covariance matrix would be generated by: 2 0 0 0 0 0 2 0 0 0 0 0 2 0 1 -1 1 -1 0 0 0 0 0 0 R = 2 0 0 0 1 -1 0 0 1 -1 0 0 2 0 0 0 0 0 2 0 0 0 0 0 1 -1 1 -1 1 0 -1 0 0 1 0 -1 2 where is the spectral density of the double-differenced carrier phase mea- surement. Performing the matrix multiplication yields: 2 R = 4 2 2 4 This pattern continues as more measurements are added. The double-differenced doppler measurements are likewise correlated, and their measurement error covariance block structures must be generated with 4 on the diagonals and 2s in all other elements. 44 Chapter 5 Control System The EFFTB uses a linear quadratic regulator (LQR) to maintain the desired separation between the formation satellites. The following section outlines the development and implementation of the LQR used in this work. 5.1 Controller Development While most natural processes can be represented with continuous func- tions, controlling these systems is usually a discrete time process. Since an optimal controller is desired, a cost function J(t) is constructed such that: tf J(t) = t [xT (t)W(t)x (t) + uT (t)V(t)u(t)]dt + xT (tf )Sf (tf )x (tf ) (5.1) The control vector u(t) is said to be optimal if it minimizes this cost function. In this equation, x is the difference between the current trajectory and some reference trajectory. The matrices W(t), V(t), and Sf (t) represent the state penalty matrix, the control penalty matrix, and the terminal penalty matrix respectively. These matrices are used as design parameters to produce the desired performance. The solution is given by [4]: u(t) = L(t)x(t) 45 where x is the state vector at time t. The matrix L(t) is determined by: L(t) = V-1 (t)BT S(t) (5.2) where B is a matrix mapping the control to the state. The matrix S(t) can be solved using a Ricatti equation: S(t) = S(t)BV-1 (t)BT S(t) - S(t)A(t) - AT (t)S(t) + W(t) S(tf ) = Sf It is worth noting that the Riccati equation can be solved by backwards propagation as a discrete process. Letting ti+1 correspond to the final epoch, and ti be the epoch just previous, the recursion relationship can be written as [8]: Li = - Vd (ti+1 , ti ) + [Bd (ti+1 , ti )]T Si+1 Bd (ti+1 , ti ) [Bd (ti+1 , ti )]T Si+1 (ti+1 , ti ) Si = [(ti+1 , ti )]T Si+1 (ti+1 , ti ) - (Li )T Vd (ti+1 , ti ) + [Bd (ti+1 , ti )]T Si+1 Bd (ti+1 , ti ) Li + Wd (ti+1 , ti ) If u(t) is held constant over the interval t - ti then: t -1 Bd (t, ti ) = Wd (ti+1 , ti ) = (t, )B( )d ti ti+1 ti ti+1 (t, ti )T W(t)(t, ti )d [Bd (t, ti )T W(t)Bd (t, ti ) + V(t)]d Vd (ti+1 , ti ) = ti The matrix (ti+1 , ti ) is the state transition matrix, and its generation is discussed in Section 4.4.2. 46 5.2 Controller Implementation Since the type of controller was not critical to the overall simulation, a simple tracking controller was implemented. In this scenario, the controller's job is to maintain a constant distance between the node and the hub satellite. The controller makes use of the absolute position and velocity estimates, and returns a maneuver to the node. The algorithm for the generation of the control vector is: 1. The absolute position of both the hub and the node spacecraft are propagated 30 seconds into the future. This is the epoch at which the control will be applied, and will be designated tc . The 30 second buffer window allows ample time for the controller to run and the control vector to be distributed to the Remote Control Machine (defined in section 6.1.4). 2. The node spacecraft constructs the absolute state of the hub by differencing its own absolute state estimate and the relative state estimate. The reconstructed hub's absolute state is then used to generate a UVW reference frame centered at the hub spacecraft. 3. The desired relative state is constructed as an offset in this UVW frame. In this work, a 100 km separation was added in the v direction to form ^ the desired state. 4. The node's estimate of the hub position is then used to rotate the desired state back into the inertial reference frame. It should be noted that each node will have a slightly different estimate of the hub's position, and therefore a different UVW frame. The result of this offset is that the desired trajectory for each node will vary slightly from the ideal 100 km separation. 47 5. The desired state is propagated forward to some epoch in the future, tf , at Sf which is defined. An arbitrary coasting period of five minutes was chosen for this controller. Accordingly, the desired state is integrated 5 minutes into the future, and the Ricatti equation is back propagated 5 minutes as well. At the end of the coasting window, a new control is calculated and the process continues. 6. The Ricatti Equation is solved by back propagating from the epoch tf to the control epoch tc in one second increments. The state transition matrix used in this procedure is generated using the same plant formulation that the filter used. 7. The back propagation yields the L(tc ) matrix, which is used to find the control vector at time tc . The algorithm presented was selected arbitrarily, and may not be efficient enough for actual flight hardware. One option to conserve computational time would be to generate an entire orbit worth of gains and store them in RAM. This would alleviate having to perform the back propagation every 5 minutes. 5.3 Offline Controller Tuning An empirical approach was selected to determine the control parame- ters that would optimize the linear controller. A "cheap" controller was used to simplify the search space for the desired values. In this situation, the terminal penalty matrix Sf was set to zero. The state penalty matrix was set to an identity matrix, and the control penalty matrix was constructed as the 48 appropriately sized identity matrix premultiplied by the scalar parameter . This scalar was adjusted until the desired results were obtained. The controller was tested offline by providing it with simulated absolute state estimates. In order to verify that the controller could track the desired position without error, the input states were noise free. The node satellite was offset from the desired initial state according to the values given in Table 5.2. The simulation spanned one orbit, and a control vector was output every 300 seconds. The maneuvers were assumed to be ideal, and always along the inertial axes. Figure 5.1 depicts a case where the control penalty matrix has a value of 5e6. This value is too high and does not allow enough actuation effort to drive the tracking error to zero. The tracking error is defined as the difference between the estimated state and the desired state. With an value of 5e5, there is too little penalty applied and the controller consistently overshoots the desired trajectory. This can be seen in Figure 5.2. With the desired control solution bracketed, the value was adjusted slightly until the results shown in Figure 5.3 were achieved. The design parameters are summarized in Table 5.1. Table 5.1: Controller Design Parameters Parameter Value Sf W V 0 0 1e6 49 Tracking Error 100 80 60 Tracking Error (m) 40 20 0 -20 -40 -60 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 5.1: Tracking Error With Too Much Control Penalty. The controller targets a desired absolute position based upon the propagated absolute state estimates. A typical pseudorange-based filtered GPS solution contains approximately two to three meters of error per axis. Therefore, to study the predicted controller performance, the simulated absolute states in this offline test were augmented with two meters of white noise before being sent to the controller. Figure 5.4 demonstrates that the controller has a 20 meter limit on its tracking accuracy. It should be emphasized that this is due to the relatively low accuracy of the absolute states sent to the controller. If more emphasis was placed upon improving this accuracy, the overall tracking error would decrease accordingly. This fact is verified in Figure 5.5, where the estimated state noise is approximately one-half meter and the tracking error drops to within a few meters. 50 Tracking Error 100 80 60 Tracking Error (m) 40 20 0 -20 -40 -60 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 5.2: Tracking Error With Too Little Control Penalty. As a final verification step, the maneuver window was decreased to ten seconds, and the noise on the absolute state was returned to two meters. With a ten second control window, the errors induced by lengthy propagation of inaccurate states are minimized, and the controller is able to track the desired trajectory to within the two meter input noise. Figure 5.6 demonstrates this tracking capability, and verifies that no stray tracking errors are leaking into the solution. 5.4 Formation Control Architecture Just as there are many different types of controllers, there are several different ways to implement a control architecture for a formation. Three such architectures are the centralized controller, decentralized controller, and a hy- 51 Tracking Error With No Feedback 100 80 60 Tracking Error (m) 40 20 0 -20 -40 -60 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 5.3: Tracking Error With Appropriate Penalty. Tracking Error With 2 Meters of Feedback Error 80 60 40 Tracking Error (m) 20 0 -20 -40 -60 -80 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 5.4: Tracking With 2 m Of Measurement Noise at 300 Sec. 52 Tracking Error With .5 Meters of Feedback Error 100 80 60 Tracking Error (m) 40 20 0 -20 -40 -60 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 5.5: Tracking With .5 m Of Measurement Noise at 300 Sec. Tracking Error With 2 Meters of Feedback Error 80 60 40 Tracking Error (m) 20 0 -20 -40 -60 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 5.6: Tracking With 2 m Of Feedback Noise at 10 Sec. 53 Table 5.2: Initial Node Absolute State Offset From The Desired Trajectory State Deviation u v w u v w 75 m -50 m 50 m .01 m/s -.05 m/s .2 m/s brid controller that leverages the strengths of the centralized and decentralized controllers[7]. These architectures are illustrated in Figure 5.7. A form of decentralized architecture was implemented for this work, and is depicted in Figure 5.8. The simulation environment was designed to allow for rapid adaptation to any of these architecture types, and would make an excellent tool to study the benefits and drawbacks of each technique. 54 Figure 5.7: Control Architectures 55 Figure 5.8: Implemented Decentralized Architecture 56 Chapter 6 Extended Formation Flying Testbed The EFFTB is a combination of hardware and software that allows for the simulation of multiple vehicles in a real-time closed-loop hardware-in-theloop situation. Each of the major components is described in the following subsections. 6.1 Hardware Description The EFFTB hardware was customized for a Low Earth Orbit (LEO) using GPS navigation signals. The design of the testbed was generalized to allow for multiple measurement sources and dynamic inputs. For example, to accommodate a deep space formation, the GPS simulator/receiver measurement source could be replaced with a simulated star tracker, and the dynamic modelling software could be replaced with the appropriate equations of motion. In this way, the EFFTB can be used to simulate nearly any proposed formation. The following descriptions detail the hardware as configured for this work. Figure 6.1 illustrates this hardware design. 6.1.1 GPS Signal Simulator The heart of the simulation setup is the GPS signal simulator. Each of the participating installations is required to have at least one such simulator 57 58 Figure 6.1: Hardware Setup. onsite. The GPS laboratory at the University of Texas at Austin/Center for Space Research (UTCSR) uses a Spirent STR4760 GPS signal simulator. This simulator is capable of simulating L1 signals for two vehicles on up to 16 channels each. All of the environmental variables affecting the GPS signals, such as ionospheric delay and GPS clock errors, are controlled by the simulator. The main goal of the EFFTB is to overcome the limitations of two RF outputs on this simulator by connecting several such simulators together. 6.1.2 Orion GPS Receivers The Orion type GPS receiver was chosen to provide measurements. The Orion receiver is based on the terrestrial GPS receiver built around the Zarlink GP2000 chipset [1]. The Orion receiver is capable of tracking 12 satellites on Figure 6.2: Orion Receivers. the L1 frequency. In order to optimize the receiver performance for space use, several modifications were made to the receiver [21]. These modifications include the fixes related to the implicit assumption of a low speed vehicle in the Doppler prediction and the correction of the time tagging error [11]. In order to facilitate measurement exchange between receivers, a clock steering algorithm has been implemented in the receiver which aligns the measurement epochs 59 and the navigation solutions to the integer GPS second. The Orion also utilizes a 3rd order phase lock loop (PLL) assisted by a 2nd order frequency lock loop (FLL) to improve the accuracy of the integrated carrier phase measurement. 6.1.3 Master Control Computer The Master Control Computer (MCC) functions as the nerve center of the EFFTB. As denoted in 6.3 this is the only machine that doesn't have to be collocated with the simulation equipment. The functions of the MCC are described below: 1. Formation Initialization The MCC possesses the initialization data for all of the GPS receivers and the Remote Control Machines. When the simulation is initialized, the MCC sends out the appropriate two-line element set to the each node. The ECEF position of each node as well as its mass and area characteristics are then sent to the Remote Control Computer for propagation. 2. Simulation Activation In order to keep the dynamic environments at all facilities within a second of each other, an electronic activation has been designed. Once all of the initialization is completed, the MCC sends out a pulse alerting each machine that the simulation has begun. Upon receiving this pulse, the Node Computers upload the two-line elements into their respective receivers, and the Remote Control Computers begin to propagate the equations of motion and feed the simulator the trajectory data. Accordingly, the difference between the start times at each facility is simply the 60 signal transmit time through the internet. 3. Data Relay In order to emulate a true broadcast environment via hardwired cables, a mock communication system has been developed. In this simulation, no node has the ability to communicate directly with any other node. All packets are instead sent to the MCC, which then rebroadcasts them to all of the nodes. The communication system on each node must be robust enough to recognize the information sent to it, identify the relevant data and reject the rest. The MCC also routes the control information directly to the Remote Control Computer. 4. Data Logging The MCC computer allows for a convenient central point for data logging and archiving. This requires extra effort on the part of the node to distribute the position and covariance information for the state, but it relieves the user of the inconvenience of collecting multiple files on several different systems. 6.1.4 Remote Control Computer The Remote Control Computer (RCC) is responsible for generating the dynamics for the simulation. While the Spirent simulator is capable of performing this task, the closed-loop environment requires an external data source. The RCC is initialized with each of the node's initial positions, and propagates them forward at 10 Hz according to the user defined environment variables. At the appropriate time, the RCC augments the node's velocity vector with the v provided by the Node Computer. The RCC smooths this 61 impulse by integrating it over one second. Ten times each second the RCC sends the trajectory information to the simulator, which then generates the appropriate GPS signals. The RCC at each location runs independently of the remainder of the RCC's. As long as the filter update intervals are sufficiently large, strict time synchronization of these machines is not necessary. However, if a rapid update rate is required, the RCC or the simulators may need to be driven by a precise clock in order to maintain system wide synchronization. 6.1.5 Hub Computer The satellite that is designated the hub is only required to broadcast its measurements at each epoch. The Hub Computer (HC) extracts the pseudorange, doppler and carrier phase measurements from the receiver, as well as the position fix solution, and forwards the data on to the MCC. The MCC then distributes the measurements to the formation 6.1.6 Node Computer The Node Computer (NC) performs all of the estimation and control computations. One NC can administer multiple GPS receivers, each utilizing its own instance of the node software. When the NC communication routine detects measurements with identical time-tags, it invokes the estimator to determine the current state. If it is time to perform a station keeping maneuver, the node calculates this maneuver and sends it to the RC. 62 6.2 Software Description The EFFTB code was developed in standard C. The primary utility is the generality of the code structure. The NC software was developed with the estimator and controller being subroutines called by the main shell. The subroutines can easily be swapped out with other, more sophisticated programs. The RC requires only a control vector, and is impartial to how it is generated. The RC software can easily be modified to write to a device other than the GPS signal generator. The overall flow of the software is illustrated in 6.3, where the dashed arrows represent data transfer. 63 Figure 6.3: Software Flowchart. 64 Chapter 7 Results This chapter details the results obtained from the simulation environment for three satellite separation distances. These tests include runs made in offline simulations as well as the filter and controller performances when inserted into the EFFTB. Tuning parameters for both the observer and controller are also documented. 7.1 Offline Filter Validation The Kalman filter was tuned using measurement data from the GPS receivers and processed in an offline environment. In these tests only two satellites were used, the node satellite and the hub satellite. The simulator settings for the "truth" environment are listed in Table 7.1. Since the goal was to focus on how the intersatellite distances affect the navigation error, the GPS clock and ephemerides errors in the GPS simulator were disabled. To evaluate the filter's performance over various distances, three candidate distances were chosen. The satellites were placed in the same orbit plane with an inclination of 87 degrees and altitude of approximately 600 km. Spacing of the nodes was accomplished by adjusting the the true anomaly. The true ECI initial conditions for each of these cases are provided in Table 7.2. The initial error covariance matrix is diagonal, and the initialization 65 Parameter Table 7.1: Simulator Environment Settings. Setting 10x10 2 m2 2 100 Kg Off Off Constant Thickness with Sinusoidal Variation Earth Gravity Model Cross Sectional Area CD Mass GPS Clock Error GPS Ephemerides Errors Ionosphere Table 7.2: True Initial Conditions: November 6 2001 00:00:00 GPS. Hub Node 1 km Node 100 km Node 500 km X Y Z Vx Vy Vz m -4.81976e6 4.81976e6 0.00000 -0.28313e3 -0.28313e3 7.64046e3 -4.81980e6 4.81972e6 0.00099e6 -0.28234e3 -0.28393e3 7.64046e3 -4.82294e6 4.81555e6 0.09975e6 -0.20390e3 -0.36231e3 7.63964e3 -4.82531e6 4.78837e6 0.49837e6 0.11332e3 -0.67808e3 7.61998e3 m/s parameters are listed in Table 7.3. The error covariance terms for the GPS ephemerides errors and the doubledifferenced integer ambiguities must be reset with the inclusion of each new GPS satellite. This is accomplished by resetting the corresponding diagonal element and zeroing out the associated row and column elements in the error covariance matrix. To examine the filter performance, and perform the necessary tuning, 66 Table 7.3: Error Covariance Initialization State Value u, v, w u, v, w u, v, w u, vw b b N B m m/s m m/s m/s m/s2 m m 400 1 400 1 10 1 25 100 GPS data was analyzed offline. The data were generated using the internal propagator in the GPS simulator. Two Orion receivers were utilized to process the measurements, which were then saved to text files. These text files were used as input to the filter to generate the results in this section. 7.1.1 Filter Performance at 1 Hz To determine the optimal performance, the measurement update rate was set at 1 Hz. Figures 7.1 and 7.2 show the absolute position and velocity errors for the node spacecraft. Since the absolute states are independent of the relative distances, these figures are not repeated for the longer baselines. It was observed that the accuracy of the absolute position is typical for a GPS pseudorange based navigation solution. The jumps in the error and error covariance are caused by the introduction and removal of satellites from the solution. The accuracies for the absolute states are listed in Table 7.4. The accuracy of the relative state is not strongly correlated with the accuracy of the absolute states. For short baselines, it was shown that the estimation of 67 Table 7.4: Absolute State Error in u, v , w frame (mean std) ^ ^ ^ State Accuracy u v w u v w m -0.0034 3.4934 -0.0450 3.0136 -0.0004 2.4191 0.0703 0.0335 -0.0185 0.0267 -0.0092 0.0415 m/s the absolute state could be replaced by simply propagating the absolute state forward with little lack of accuracy [28]. The relative accuracies were determined by generating an ensemble average over five simulation runs for each baseline length. The results are depicted in Figures 7.3 - 7.8. The filter parameters for each baseline are detailed in Table 7.5. It should be noted that the large process noise values for the clock bias and drift reflect the inability to accurately model these states. Table 7.6 details the steady-state accuracies for these simulations. The errors at the one kilometer baseline are comparable to those found in current literature [29] [11]. The convergence for the relative position is dependent on the correct resolution of the double differenced integer ambiguity states. For the one second update rate, this convergence takes approximately 1000 seconds. For the 1 km and 100 km case, once these integers are determined, the solution is only slightly disturbed by the inclusion of new satellite 68 UVW Absolute Position Error 20 0 -20 0 500 1000 1500 2000 2500 3000 3500 4000 4500 20 Error (m) 0 -20 0 500 1000 1500 2000 2500 3000 3500 4000 4500 20 0 -20 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.1: Ensemble Mean Absolute Position Actual and Formal Errors. UVW Absolute Velocity Error 0.5 0 -0.5 0.5 Error (m/s) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.2: Ensemble Mean Absolute Velocity Actual and Formal Errors. 69 UVW Relative Position Error Plots: 1 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.3: 1 km Ensemble Mean Relative Position Actual and Formal Errors/1 Second Update. UVW Relative Velocity Error Plots: 1 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.4: 1 km Ensemble Mean Relative Velocity Actual and Formal Errors/1 Second Update. 70 UVW Relative Position Error Plots: 100 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.5: 100 km Ensemble Mean Relative Position Actual and Formal Errors/1 Second Update. UVW relative velocity Error Plots: 100 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.6: 100 km Ensemble Mean Relative Velocity Actual and Formal Errors/1 Second Update. 71 UVW Relative Position Error Plots: 500 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.7: 500 km Ensemble Mean Relative Position Actual and Formal Errors/1 Second Update. UVW Relative Velocity Error Plots: 500 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.8: 500 km Ensemble Mean Relative Velocity Actual and Formal Errors/1 Second Update. 72 Table 7.5: Process Noise Settings: 1Hz Update Parameter 1 km Case 100 km Case 500 km Case 2 rABS 2 vABS 2 rREL 2 vREL 2 b 2 b 2 N 2 B m2 m m m2 s2 2 m2 s2 2 m2 s2 2 m m2 1e-2 5e-5 1e-6 1e-7 1e6 1e6 1e-6 1e-5 1e-2 5e-5 1e-6 1e-7 1e6 1e6 1e-5 1e-5 1e-1 1e-5 1e-7 1e-6 1e6 1e6 1e-5 1e-5 Table 7.6: Relative State Error: 1 Sec. Update Rate (mean std) State 1 km Case 100 km Case 500 km Case u v w u v w pairs. The 500 km ensemble does not converge as well as the other two scenarios. In two of the simulations, the filter struggled to accurately estimate the ambiguities in the face of the greater differential ionosphere delay error. The ensemble averaging reflects this difficulty in larger errors. Another factor limiting the accuracy at the long baseline is the number of common view GPS satellites available. The number of satellites utilized in the double difference is illustrated in Figure 7.9. While the 500 km baseline has five or more satellites m -0.0024 0.0044 0.0041 0.0105 0.012 0.0063 0.00001 0.0002 -0.00001 0.0002 0.00002 0.0001 -0.0478 0.0149 -0.0114 0.0324 0.0310 0.02376 -0.0000 0.0001 -0.0001 0.0002 0.0000 0.0002 -0.0307 0.0644 0.0989 0.0627 -0.0501 0.0281 0.0001 0.0005 0.0000 0.0005 -0.0000 0.0001 m/s 73 Number of Common GPS Satellites for Double Differencing 12 10 500 Km 100 Km 1 Km Available Common View Satellties 8 6 4 2 0 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.9: Common View GPS Satellites. in view for a majority of the time, there is a ten minute time period where there are four or fewer satellites available. It is apparent that the filter as constructed has a maximum relative baseline length near 500 km. 7.1.2 Filter Performance at 5 Seconds In an actual distributed system, a one second update rate may pose too heavy of a burden on the computer and communication systems. Accordingly, the measurement update rate of the filter was slowed to 5 seconds. The filter gains were adjusted to account for the larger emphasis placed upon the dynamical models. The process noise covariance values are listed in 7.7. The ensemble results for the five second update rate are displayed in 7.10-7.15. 74 Table 7.7: Process Noise Settings: 0.2 Hz Update Parameter 1 km Case 100 km Case 500 km Case 2 rABS 2 vABS 2 rREL 2 vREL 2 b 2 b 2 N 2 B m2 m m m2 s2 2 m2 s2 2 m2 s2 2 m m2 1e-2 5e-4 1e-5 1e-5 1e6 1e6 1e-6 1e-5 1e-1 5e-4 1e-4 1e-6 1e6 1e6 1e-5 1e-5 1e-1 1e-4 1e-4 1e-5 1e6 1e6 1e-4 1e-5 UVW Relative Position Error Plots: 1 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 5000 Figure 7.10: 1 km Ensemble Mean Relative Position Actual and Formal Errors/5 Second Update. 75 UVW relative velocity Error Plots: 1 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 5000 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Figure 7.11: 1 km Ensemble MeanRelative Velocity Actual and Formal Errors/5 Second Update. UVW Relative Position Error Plots: 100 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.12: 100 km Ensemble Mean Relative Position Actual and Formal Errors/5 Second Update. 76 UVW relative velocity Error Plots: 100 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.13: 100 km Ensemble Mean Relative Velocity Actual and Formal Errors/5 Second Update. UVW Relative Position Error Plots: 500 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.14: 500 km Ensemble Mean Relative Position Actual and Formal Errors/5 Second Update. 77 UVW relative velocity Error Plots: 500 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.15: 500 km Ensemble Mean Relative Velocity Actual and Formal Errors/5 Second Update. Table 7.8: Relative State Error: 5 Sec. Update Rate (mean std) State 1 km Case 100 km Case 500 km Case u v w u v w m 0.0295 0.0129 0.0035 0.0048 0.0122 0.0001 0.0000 0.0001 -0.00001 0.0002 0.00000 0.0001 -0.0116 0.0566 -0.0310 0.0317 0.0718 0.0252 -0.0000 0.0001 0.0001 0.0000 0.0001 0.0003 -0.0127 0.2237 0.2308 0.2141 0.2308 0.0676 0.0000 0.0016 0.0004 0.0005 -0.0001 0.0008 m/s 78 Decreasing the update rate had little impact on the overall accuracies of 1 km and 100 km simulations. The steady state position noise levels are still better than 2 cm and 6 cm respectively. The inability of the filter to correctly resolve the double differenced integer ambiguities became more pronounced for the 500 km simulation. The large jump around 4000 seconds occurred in one of the runs when the number of common view satellites dropped to three, causing the relative position estimate to jump dramatically. 7.1.3 Cycle Slips The integer ambiguity, as described in Chapter 3, is the number of whole carrier waves between the GPS satellite and the receiver. As long as the receiver keeps locked on the carrier, this integer is constant. If the receiver momentarily loses the carrier signal, the ambiguity may change by one or more whole cycles, resulting in a cycle slip. This slip would then have to be detected and accounted for in the filter in order to correctly process the carrier-phase measurement. The Orion receivers manufactured at UT/CSR are prone to incur cycle slips, possibly due to the hardware fabrication process, or perhaps the hardware or software itself. These jumps usually manifest themselves in the form of a half-cycle slip. Since the carrier signal has a 19 cm wavelength, a halfcycle slip causes a 10 cm jump in the observed measurement. These slips are shown in Figure 7.24. Each spike occurring at a multiple near 10 cm indicates a cycle slip. Longer spikes indicate a poorly initialized ambiguity, which was reset during the next measurement update. In order to minimize the impact of these slips, the filter checks the residuals for outliers at each measurement update. If a residual is greater than 79 9 cm, it is excluded from the solution. If this satellite combination generates two successive outliers, the filter assumes that a cycle slip has occurred. On the next measurement update, the ambiguity term and the corresponding error covariance element are reinitialized. The drawback to this method occurs when the one reference GPS satellite, the satellite against which all of the others are differenced, experiences a cycle slip. That error then saturates all of the double-differenced carrier phase measurements. When passing through the filter, all of the phase measurements are considered outliers, and the filter is left with no carrier-phase measurements for two time steps. A double-differencing algorithm which computed the differences successively would rectify this design weakness. 7.1.4 Best Case Filter Performance In order to establish the best case filter performance, the two most stable Orion receivers were used, minimizing the impact of cycle slips. Figures 7.16-7.23 show these single run simulations for both the one second and the five second update rate. The accuracy for each scenario is listed in Table 7.9. 80 UVW Relative Position Error Plots: 1km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.16: Best Case 1 km Relative Position Actual and Formal Errors/1 Second Update. UVW Relative Velocity Error Plots: 1 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.17: Best Case 1 km Relative Velocity Actual and Formal Errors/1 Second Update. 81 UVW Relative Position Error Plots: 100 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.18: Best Case 100 km Relative Position Actual and Formal Errors/1 Second Update. UVW Relative Velocity Error Plots: 100 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.19: Best Case 100 km Relative Velocity Actual and Formal Errors/1 Second Update. 82 UVW Relative Position Error Plots: 1 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.20: Best Case 1 km Relative Position Actual and Formal Errors/5 Second Update. UVW Relative Velocity Error Plots: 1 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.21: Best Case 1 km Relative Velocity Actual and Formal Errors/5 Second Update. 83 UVW Relative Position Error Plots: 100 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.22: Best Case 100 km Relative Position Actual and Formal Errors/5 Second Update. UVW relative velocity Error Plots: 100km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.23: Best Case 100 km Relative Velocity Actual and Formal Errors/5 Second Update. 84 Cycle Slips 0.4 0.3 0.2 Carrier-Phase Residuals (m) 0.1 0 -0.1 -0.2 -0.3 -0.4 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.24: Cycle Slip Visualization. Best-Case Cycle Slips 0.2 0.15 0.1 Carrier-Phase Residuals (m) 0.05 0 -0.05 -0.1 -0.15 -0.2 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.25: Best Case Cycle Slip Profile. 85 Table 7.9: Relative State Error: Best Case Filter Runs (mean std) State 1 km Case 100 km Case 1 Second u v w u v w 5 Second u v w u v w m -0.0025 0.0029 -0.0033 0.0024 -0.0097 0.0008 0.0000 0.0002 -0.0000 0.0005 0.00000 0.0002 0.0131 0.0122 -0.0042 0.0155 -0.0069 0.0124 0.0001 0.0003 -0.0006 0.0006 0.0003 0.0004 m -0.0018 0.0034 -0.0015 0.0026 -0.0053 0.0025 0.0000 0.0001 -0.000 0.0002 0.0000 0.0002 0.0062 0.0197 -0.0292 0.0157 -0.0030 0.0123 -0.0000 0.0002 0.0000 0.0003 0.0000 0.0001 m/s m/s Table 7.9 illustrates the filter performance when the cycle slips are kept to a minimum, as depicted in Figure 7.25. These accuracies range from two to five times lower than the ensemble runs, which combined data with varying amounts of cycle slips. Accordingly, in order to maximize the filter's efficiency, GPS receivers which have a minimal cycle slip rate should be chosen. 7.1.5 Receiver Clock Model Validation As described in Chapter 4, the GPS receiver clock and receiver offsets are modelled as random variables that are disconnected. Traditionally, the frequency offset is modelled as the derivative of the clock offset state. The clock steering performed by the Orion disconnects these states, which requires 86 them to be modelled independently. Figures 7.26 and 7.27 show the typical estimation values for the two clock parameters. It is apparent that there is little correlation between the figures, which validates the modelling approach. Scaled Clock Offset Estimation 0 -5 Scaled Clock Offset Estimation (m/s) -10 -15 -20 -25 -30 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.26: GPS Clock Offset Estimation. 87 Scaled Clock Frequency Estimation -230 -232 Scaled Clock Frequency Estimation (m/s ) 2 -234 -236 -238 -240 -242 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.27: GPS Frequency Offset Estimation. 7.2 Extended Formation Flying Testbed Results In order to test robustness of the navigation and control routines, they were inserted into the EFFTB communication package. The fidelity of the algorithms were tested by requiring them to process real time data at a 5 second update rate, and generate control vectors at 300 second intervals. A benchmark test was performed in the GPS Laboratory at UT/CSR. This result was then used to analyze the results generated from tests spanning the internet. 7.2.1 Two Satellite Local Simulation Since the GPS simulator at UT/CSR is limited to two spacecraft, initial EFFTB tests were conducted with one hub spacecraft and one node craft. This testing allowed for the verification of the communication links, and the validation of the filter and controller in the real-time environment. Figure 7.28 88 and Figure 7.29 demonstrate the relative position and velocity errors, which are summarized in Table 7.11. Table 7.10: Relative State Error: EFFTB 2 Satellite Simulation (mean std) State Steady State Error u v w u v w (m) 0.0435 0.0621 0.0300 0.0345 -0.0126 0.0604 -0.0000 0.0002 -0.0001 0.0003 0.0000 0.0004 (m/s) The filter performs very closely to the predicted laboratory results. Excluding the large jump at 3500 seconds, the steady state noise is around 3 cm. The spike in estimation error is caused by a large (35 cm/s) maneuver commanded by the controller. The filter treats the delta-v as an impulse, while the truth environment integrates the maneuver over one second. Additionally, the large maneuver causes the receiver to lose lock on all of the GPS satellites, so the filter processes no relative measurements for 10 seconds while it resets the ambiguity terms. During these 10 seconds, the filter's dynamic model and the truth model drift apart, yielding the error shown in the plot. It is also shown that within 400 seconds, the solution has once again converged. The tracking errors are shown in Figure 7.30. The controller performed as in the offline cases, exhibiting a 20 meter tracking error over the simulation period. Figure 7.31 illustrates the control effort generated by the controller. The maneuvers are typically on the order of 5 cm/s in any direction, but can grow as high as 30 cm/s. It is apparent that such a control scheme would be very costly in fuel, and may not be an optimal control strategy. 89 UVW Relative Position Error Plots: 100 Km Separation 0.5 0 -0.5 0.5 Error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.28: Ensemble Mean Relative Position Error Actual and Formal Errors at 100 km. UVW Relative Velocity Error Plots: 100 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 Error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.29: Ensemble Mean Relative Velocity Actual and Formal Errors at 100 km. 90 EFFTB Tracking Error: 100 Km 60 40 20 Error (m) 0 -20 -40 -60 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.30: Ensemble Mean Tracking Error with 300 Sec Window Actuation Effort 0.2 0.1 0 -0.1 -0.2 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Actuation (m/s) 0.2 0.1 0 -0.1 -0.2 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0.2 0.1 0 -0.1 -0.2 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.31: Ensemble Mean Actuation Effort in the u, v , w frame ^ ^ ^ 91 The controller passes the maneuver back to the filter to update its dynamical model. Since no error is superimposed upon this velocity vector, the covariance envelope is not widened during the burns. This can be visualized in Figure 7.28, especially around the 3600 second mark. A degree of uncertainty should be added to the solution to compensate for the the misrepresentation of the maneuver, requiring an increase of the covariance around the control time. 7.2.2 Multi-Satellite Formation The full EFFTB was tested by coupling the Spirent GPS signal generator located at UT/CSR with the 4 RF signal simulator hosted at Goddard Space Flight Center's Formation Flying Testbed. Several small modifications were performed on the communication environment in order to comply with network security policies at both sites. The formation scenario implemented for the final phase of testing was similar to the set of conditions presented earlier. The four node spacecraft were distributed approximately 100 km from the hub spacecraft in either a leading or trailing position. The results of four real-time data sets were averaged to generate Figures 7.32 and 7.33. Both the relative position and velocity errors were consistent with the previous results, implying that it is not necessary to maintain exact synchronization between the two simulation environments. The typical offset between the two environment computers was on the order of one second, well below the 5 second filter update rate. For applications requiring estimates closer to 1Hz, time synchronization devices would be required at both laboratories. A set of four simulations was run with the controller enabled. The 92 UVW Relative Position error Plots: 100 Km Separation 0.5 0 -0.5 0.5 error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.32: Full Simulation Environment Ensemble Mean Relative Position Actual and Formal Errors scenario contained five satellites, the hub and four node spacecraft. Such a situation would have been previously impossible without a 5 RF GPS simulator. The controller was again set to burn every 300 seconds, and the initial offsets from the desired trajectory were also held the same. The relative position and velocity errors are illustrated in Figures 7.34 and 7.35 respectively. It is apparent that the ensemble averaging of the 16 individual runs helped to smooth out many of the large jumps seen in Figure 7.28. The tracking error shown in 7.36, also demonstrates a somewhat smoother nature. The EFFTB communication system had no difficulties in routing messages to any of the nodes, and seamlessly linked the two laboratories. Figure 7.37 illustrates the average control effort exhausted over the four controlled simulations. The total control effort is detailed in Table 7.12. On average, it would require three meters/second of control per orbit to maintain 93 UVW relative velocity error Plots: 100 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.33: Full Simulation Environment Ensemble Mean Relative Velocity Actual and Formal Errors the desired spacing of this orbit. This expenditure would be prohibitively high for most formation missions. The point of this example was to demonstrate a proof of control capability using this navigation system rather than to recommend a specific guidance method, which would be application dependent. 94 Table 7.11: Relative State Error: EFFTB 5 Satellite Simulation (mean std) Non-Controlled State Steady State Error u v w u v w Controlled State u v w u v w (m) 0.0114 0.0699 -0.0559 0.0458 0.0323 0.0312 0.00000 0.0002 -0.00002 0.0001 0.00006 0.0003 Steady State Error -0.0612 0.0710 -0.0943 0.0556 0.0037 0.0962 0.00000 0.0005 0.00001 0.0004 0.00016 0.0007 (m/s) (m) (m/s) 95 UVW Relative Position error Plots: 100 Km Separation 0.5 0 -0.5 0.5 error (m) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0.5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 -0.5 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.34: Ensemble Mean Relative Position Actual and Formal Errors With Control Effort UVW relative velocity error Plots: 100 Km Separation 0.01 0.005 0 -0.005 -0.01 0.01 error (m/s) 0.005 0 -0.005 -0.01 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 7.35: Ensemble Mean Relative Velocity Actual and Formal Errors With Control Effort 96 Tracking Error 80 60 U error V error W error 40 Tracking Error (m) 20 0 -20 -40 -60 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure 7.36: Ensemble Mean Tracking Error in the u, v , w frame ^ ^ ^ Control Effort 0.2 0.1 0 -0.1 -0.2 0 Control Effort (m/s) 500 1000 1500 2000 2500 3000 3500 0.2 0.1 0 -0.1 -0.2 0 500 1000 1500 2000 2500 3000 3500 0.2 0.1 0 -0.1 -0.2 0 500 1000 1500 Time (s) 2000 2500 3000 3500 Figure 7.37: Ensemble Mean Actuation Effort in the u, v , w frame ^ ^ ^ 97 Table 7.12: Control Effort in the u, v , w frame ^ ^ ^ State Control Effort u v w T otal m/s 0.8952 0.777885 1.144059 2.817145 98 Chapter 8 Conclusion The purpose of this work was to develop a general formation navigation filter capability. The filter possesses high accuracy relative position estimation over a variety of inter-satellite distances and number of satellites included in the formation. This filter was tested in a real-time closed-loop hardwarein-the-loop tesbed which was developed to utilize the internet to exchange measurements. A summary of the results is presented in the following section. Finally, directions for future work are indicated. 8.1 Summary of Results The study of large satellite formations can be divided into two cat- egories: formations with large intersatellite spacings and formations with a large number of satellites. It was desired to look at both of these problems as generally as possible, and show that both are feasible for future formation missions using GPS sensor measurements. The problem of relative navigation over large baselines was examined first. In order to strike a balance between filter complexity and accuracy, a plant model with atmospheric drag and J2 gravity perturbation was implemented. All of the GPS observables were used in a filter state that estimated both the relative position and the absolute position of the host satellite. It was found that the filter performed quite well over a large range of 99 baselines. It was determined that for shorter baselines, the accuracy of the absolute position had little impact on the overall relative solution. For the relative solution at the 1Hz update rate, both the 1 km separation and the 100 km separation had steady state noise of less than a few centimeters over an ensemble of runs. The 1 km results were consistent with other works, validating the filter design, and displaying its robustness at multiple baselines. At 500 km, the models, primarily the ionospheric delay model, began to degrade in such a way that the integer ambiguity states became difficult to determine. This led to a drop in performance, signaling the limit of the design outlined here. When the measurement update rate was increased to 5 seconds, neither the 1 km nor the 100 km estimator accuracies decreased significantly. In general, the steady state errors doubled, but left them still below 10 centimeters. One of the problems that was faced in processing the actual GPS data was the handling of cycle slips. The slips occurred sporadically in some of the receivers that were used in the testing, and more frequently in others. Simulations with a large number of slips had a much slower relative solution convergence rate, and higher overall navigation error. It was shown that for a run with minimal cycle slips, the 1 km spacing had steady state error no larger than 1 centimeter, and the 100 km case maintained error less than 2 centimeters. These results validate the estimator design and emphasize the importance of the hardware performance on the final results. A simple linear controller was developed in order to close the loop in the EFFTB. The purpose of the controller was to demonstrate the closed-loop hardware-in-the-loop system performance for a representative formation flying mission. The designed controller was able to accurately control the system with no feedback error, and when the controller update interval was reduced 100 to 10 seconds, it was able to track small amounts of measurement error precisely. When the absolute state filter feedback was simulated by adding 3 meters of feedback noise, the accuracy of the controller dropped to 20 meters. This was determined to be the result of the less accurate absolute position solution being propagated forward in time to generate the desired trajectory. It was shown that if the absolute position error was reduced to .5 meters, the controller accuracy improved to within a couple of meters. The Extended Formation Flying Testbed (EFFTB) was developed to link hardware from remote sites in order to simulate formations with many satellites. This was done because it is unrealistic to locate all of the testing resources that are required for a large formation at one facility. The communication package was designed to allow for many combinations of hardware or software components required for testing. Once developed, it was specialized to accommodate Orion GPS receivers and the Spirent GPS signal simulators used in this test program. This instance of the EFFTB makes use of several different programs running simultaneously on different machines to emulate a satellite formation in orbit. Data is passed back and forth in a broadcast fashion to mimic space communications. Every effort was made to keep the simulation environment as authentic as possible. The long baseline filter and the linear controller were inserted into the testbed architecture in order to verify the offline results. Since the 500 km baseline demonstrated to be at the limit of applicability for this filter, the focus was placed on the 100 km filter. It was found that the steady state error was on the order of 5 centimeters, only 3 times as noisy as the offline case. The decrease in performance is due to the addition of the controller, which caused small jumps in the relative position error when the control effort was large. The tracking error had a 20 meter error band, as predicted in the offline 101 testing. Both the controller and the estimator performed comparably in the real-time testing and in the offline studies, verifying the design and utility of this navigation system. There was no appreciable decrease in accuracy when the simulation environment was extended to include other GPS simulators located at a different facility. The synchronization between the two simulation systems did not play a critical role in the relative navigation, assuming the update rate was greater than the anticipated time offset. Using the newly-developed tool, closed-loop simulations can be created which span multiple environments, allowing the study of navigation and control on truly large formations. 8.2 Future Work The results presented in the previous section suggest that the filter de- sign is robust up to 500 km, but a few simplicities were included that may increase the errors from the values stated here. 8.2.1 Ionosphere Study Future work should be done examining the effects of the ionosphere, es- pecially at longer baselines. The Spirent simulator can produce various smooth ionosphere models, which may not accurately represent the actual space environment. This filter design should be tested against a more realistic ionosphere model to further examine its effects on the ambiguity resolution. 102 8.2.2 Multi-Path Study GPS multi-path error was neglected during this study, but could serve as a major error source. This error is highly dependent on the satellite's shape, and cannot easily be modelled in a general sense. It may be helpful to take a single satellite, and analyze the impact that multi-path would have on it, generating some insight into the effect on the filter's performance. 8.2.3 Attitude Modelling and Control The main piece that is missing from this formation simulation is the at- titude determination and control capabilities. It was assumed that the control was always performed along the X,Y,Z inertial axis, regardless of the satellite's orientation. The inclusion of attitude modelling would open a new realm of reality to the testbed. This work is anticipated to be performed in the near future at UT/CSR. 8.2.4 Formation Control The EFFTB can be used as a tool to investigate the control of a for- mation as a single entity. The different control architectures can be studied to determine which is most appropriate for a given mission. With the addition of more accurate thruster models, formation fuel consumption and mission life can also be studied. 103 8.2.5 EFFTB Extension The EFFTB was designed to be a general purpose simulation tool. Future work could include the developing of different instances that use other simulated measurements instead of the GPS simulator/receiver combination. In doing so, missions such as High Earth Orbit (HEO) or deep space missions can be considered and modelled. The last step could be to connect flight computers to the EFFTB in order to encompass the limited computation power experienced in space. 104 Appendices 105 Appendix A Dynamic Equation Partial Derivatives A.1 x x y x z x y x x y Dynamic Equation Partial Derivatives =0 =0 =0 =0 =0 x x y y z y z x x z =0 =0 =0 =0 =0 x x y z z z z y y z =0 =0 =0 =0 =0 x = 1 x y = 1 y z = 1 z x 3 = - 3 1 - J2 x r 2 +3 RE r 2 5 2 z2 -1 r2 7 z2 -1 r2 x2 5 1 - J2 r5 2 RE r 1 A x(x + E y) + CD V 2 m rH A E (x + E y)(y - E x) 1 + CD 2 m V 106 x xy 5 = 3 5 1 - J2 y r 2 RE r 2 7 z2 -1 r2 A y(x + E y) 1 + CD V 2 m rH 1 A E (x + E y)2 - CD V E + 2 m V x xz 5 = 3 5 1 - J2 z r 2 RE r 2 z2 7 2 -3 r x x x y x z y x = = = = 1 A z(x + E y) + CD V 2 m rH 1 A (x + E y)2 - CD V + 2 m V 1 A (y - E x)(x + Ey) - CD 2 m V 1 A z(x + E y) - CD 2 M V 2 xy 5 RE z2 3 5 1 - J2 7 2 -1 r 2 r r 1 A x(y - E x) + CD V 2 m rH 1 A E (y - E x)2 - CD V E - 2 m V y 3 = - 3 1 - J2 y r 2 RE r 2 5 2 z2 -1 r2 7 z2 -1 r2 y 2 5 +3 5 1 - J2 r 2 RE r A y(y - E x) 1 + CD V 2 m rH 1 A E (x + E y)(y - E x) - CD 2 m V 107 y yz 5 = 3 5 1 - J2 z r 2 RE r 2 7 z2 -3 r2 y x y y y z z x = = = = A z(y - E y) 1 + CD V 2 m rH 1 A (y - E x)(x + Ey) - CD 2 m V 1 (y - E x)2 A - CD V + 2 m V 1 A z(y + E x) - CD 2 M V 2 xz 5 RE z2 3 5 1 - J2 7 2 -3 r 2 r r 1 A zx + CD V 2 m rH A E z(y - E x 1 + CD 2 m V 2 z yz 5 RE z2 = 3 5 1 - J2 7 2 -3 y r 2 r r 1 A zy + CD V 2 m rH A E z(x - E y 1 - CD 2 m V 2 z 3 RE z2 = - 3 1 - J2 5 2 -3 z r 2 r r z 2 5 +3 5 1 - J2 r 2 1 A zz + CD V 2 m rH RE r 2 7 z2 -5 r2 108 z 1 A z(x + E y) = - CD x 2 m V z 1 A z(y + E x) = - CD y 2 m V z 1 A z2 = - CD V + z 2 m V 109 Appendix B Measurement Equation Partial Derivatives B.1 Pseudorange Partial Derivatives P =0 y P =0 y P =0 y P =0 N P z P z P z P B P =0 x P =0 x P =0 x P = 0 b =0 =0 =0 =0 P x P y P z P b (rx,GP S - rx ) |rGP S - r| (ry,GP S - ry ) = - |rGP S - r| (rz,GP S - rz ) = - |rGP S - r| = - = 1 110 B.2 Doppler Partial Derivatives D y D y P N D =0 x D =0 x P =0 b =0 =0 =0 D z D z P B =0 =0 =0 D x D y D z D x D y D z D b = = = = = = (rx - rx,GP S ) |rGP S - r| (ry - ry,GP S ) |rGP S - r| (rz - rz,GP S ) |rGP S - r| (rx - rx,GP S ) |r - rGP S | (ry - ry,GP S ) |r - rGP S | (rz - rz,GP S ) |r - rGP S | = 1 111 B.3 x x b Double-Differenced Carrier Phase Partial Derivatives =0 =0 =0 y y b =0 =0 =0 z z B =0 =0 =0 = x A A -rGP S,x + rx + rx (rGP S,x - rx ) + |rA S - r - r| |rA S - r| GP GP - = y B B -rGP S,x + rx + rx (rGP S,x - rx ) + |rB S - r - r| |rB S - r| GP GP A A -rGP S,y + ry + ry (rGP S,y - ry ) + |rA S - r - r| |rA S - r| GP GP - = z B B -rGP S,y + ry + ry (rGP S,y - ry ) + |rB S - r - r| |rB S - r| GP GP A A -rGP S,z + rz + ry (rGP S,z - rz ) + |rA S - r - r| |rA S - r| GP GP - = x B B -rGP S,z + rz + rz (rGP S,z - rz ) + |rB S - r - r| |rB S - r| GP GP A -rGP S,x + rx + rx | - rGP S + r + r - B -rGP S,x + rx + rx | - rB S + rx + rx GP 112 = y A -rGP S,y + ry + ry | - rGP S + r + r - = z B -rGP S,y + ry + ry | - rB S + ry + ry GP A -rGP S,z + rz + rz | - rGP S + r + r - N = 1 B -rGP S,z + rz + rz | - rB S + rz + rz GP 113 B.4 Double-Differenced Doppler Partial Derivatives =0 AB DIJ b AB DIJ b =0 AB DIJ N =0 AB DIJ B =0 AB DIJ x = -rGP S,x + rx + rx A A (| - rGP S + r + r|) A -rGP S,x + rx + rx + (-rA S + r + r) (-rA S + r + r) GP GP A 3 (| - rGP S + r + r|) - A -rGP S,x + rx -|rA S GP + r| + A rGP S,x - rx -|rA S GP + r| (-rA S + r) (-rA S + r) GP GP -rGP S,x + rx + rx B B (| - rGP S + r + r|) B -rGP S,x + rx + rx (-rB S + r + r) (-rB S + r + r) GP + GP (| - rB S + r + r|)3 GP - AB DIJ y B -rGP S,x + rx -|rB S + r| GP + B rGP S,x - rx -|rB S + r| GP (-rB S + r) (-rB S + r) GP GP = -rGP S,y + ry + ry A A (| - rGP S + r + r|) A -rGP S,y + ry + ry + (-rA S + r + r) (-rA S + r + r) GP GP (| - rA S + r + r|)3 GP - A -rGP S,y + ry -|rA S + r| GP + A rGP S,y - ry -|rA S + r| GP (-rA S + r) (-rA S + r) GP GP -rGP S,y + ry + ry B (| - rB S + r + r|) GP B -rGP S,y + ry + ry + (-rB S + r + r) (-rB S + r + r) GP GP (| - rB S + r + r|)3 GP - B -rGP S,y + ry -|rB S GP + r| + B rGP S,y - ry -|rB S GP + r| (-rB S + r) (-rB S + r) GP GP 114 AB DIJ z -rGP S,z + rz + rz A = (| - rA S + r + r|) GP A -rGP S,z + rz + rz + (-rA S + r + r) (-rA S + r + r) GP GP (| - rA S + r + r|)3 GP - A -rGP S,z + rz -|rA S + r| GP + A rGP S,z - rz -|rA S + r| GP (-rA S + r) (-rA S + r) GP GP -rGP S,z + rz + rz B (| - rB S + r + r|) GP B -rGP S,z + rz + rz + (-rB S + r + r) (-rB S + r + r) GP GP (| - rB S + r + r|)3 GP - AB DIJ x B -rGP S,z + rz -|rB S GP + r| + B rGP S,z - rz -|rB S GP + r| (-rB S + r) (-rB S + r) GP GP = A A -rGP S,x + rx + rx -rGP S,x + rx - | - rA S + r + r| |rA S + r| GP GP - AB DIJ y B B -rGP S,x + rx + rx -rGP S,x + rx - | - rB S + r + r| |rB S + r| GP GP = A A -rGP S,y + ry + ry -rGP S,y + ry - | - rA S + r + r| |rA S + r| GP GP - AB DIJ z B B -rGP S,y + ry + ry -rGP S,y + ry - | - rB S + r + r| |rB S + r| GP GP = A A -rGP S,z + rz + rz -rGP S,z + rz - | - rA S + r + r| |rA S + r| GP GP - B B -rGP S,z + rz + rz -rGP S,z + rz - | - rB S + r + r| |rB S + r| GP GP 115 AB DIJ x -rGP S,x + rx + r A = (| - rA S + r + r|)3 GP A -rGP S,x - rx - rx + (rA S - r - r) (rA S - r - r) GP GP A | - rGP S + r + r| -rGP S,x + rx + r B B (| - rGP S + r + r|)3 + AB DIJ y B -rGP S,x - rx - rx (rB S - r - r) (rB S - r - r) GP GP | - rB S + r + r| GP -rGP S,y + ry + r A = (| - rA S + r + r|)3 GP A -rGP S,y - ry - ry (rA S - r - r) (rA S - r - r) GP + GP | - rA S + r + r| GP -rGP S,y + ry + r B B (| - rGP S + r + r|)3 + AB DIJ z B -rGP S,y - ry - ry (rB S - r - r) (rB S - r - r) GP GP | - rB S + r + r| GP -rGP S,z + rz + r A = A (| - rGP S + r + r|)3 + A -rGP S,z - rz - rz (rA S - r - r) (rA S - r - r) GP GP | - rA S + r + r| GP -rGP S,z + rz + r B B (| - rGP S + r + r|)3 B -rGP S,z - rz - rz + (rB S - r - r) (rB S - r - r) GP GP B | - rGP S + r + r| 116 AB DIJ x AB DIJ y AB DIJ z = = = A B -rGP S,x + rx + rx -rGP S,x + rx + rx - -|rA S + r + r| -|rB S + r + r| GP GP A B -rGP S,y + ry + ry -rGP S,y + ry + ry - -|rA S + r + r| -|rB S + r + r| GP GP B A -rGP S,z + rz + rz -rGP S,z + rz + rz - -|rA S + r + r| -|rB S + r + r| GP GP 117 Appendix C Filter Verification Procedures To ensure proper functionality of the EKF, several tests were performed. These tests and their results are provided for completeness. C.1 Verification of the Dynamic Partial Derivative Matrix In order to verify the proper construction of the dynamic partial deriva- tive matrix, its entries are compared to the complimentary quantities generated from taking the finite differences of the non-linear differential equations. From the discussion in Chapter 5, the dynamic partial derivative matrix is generated from: A(t) = F |x x 0 This can be compared to the approximate value computed from the forward difference formula: F fi (x1 , x2 , . . . , xi + h, . . . , xn ) - fi (x1 , x2 , . . . , xi , . . . , xn ) xi h where fi is the ith state element that is being examined, and h is a small number. The non-trivial results of this comparison, performed for a 100 Km satellite separation, are detailed in Table C.1 and Table C.2. 118 Table C.1: Dynamic Partial Differential Equation Validation Derivative A Matrix Entry Numerical Derivative Difference x x x y x z x x x y x z y x y y y z y x y y y z z x z y z z z x z y z z 4.43337762E-07 -1.63758144E-06 -9.26468001E-07 -3.18904346E-10 1.96743068E-11 -7.35739287E-11 -1.63757938E-06 3.15532103E-07 8.91018701E-07 1.96743068E-11 -3.16468679E-10 6.91435351E-11 -9.26494804E-07 8.91045651E-07 -7.58869552E-07 -7.35739287E-11 6.91435351E-11 -5.56539317E-10 4.43339843E-07 -1.63757952E-06 -9.26466901E-07 -3.18917300E-10 1.96751108E-11 -7.35734433E-11 -1.63757938E-06 3.15532141E-07 8.9101873243E-07 1.96751123E-11 -3.1647196713E-10 6.91432569544E-11 9.26494814E-07 8.91045641E-07 -7.58869605E-07 -7.35734404E-11 6.91432596E-11 -5.56536974E-10 -2.08117E-12 -1.92069E-12 -1.10075E-12 1.29539E-14 -8.04001E-16 -4.85365E-16 1.30301E-17 -3.82064E-14 -3.15294E-14 -8.05459E-16 3.28765E-15 2.78177E-16 9.75885E-15 1.02093E-14 5.29765E-14 -4.88296E-16 2.75518E-16 -2.34303E-15 119 Table C.2: Dynamic Partial Differential Equation Validation Cont. Derivative A Matrix Entry Numerical Derivative Difference x x x y x z x x x y x z y x y y y z y x y y y z z x z x z x z x z x z x 4.25051620E-07 -1.61680315E-06 -9.58182648E-07 -3.20395142E-10 2.12700884E-11 -7.59046096E-11 -1.61680111E-06 2.92706778E-07 9.19767672E-07 2.12700884E-11 -3.17854490E-10 7.15372075E-11 -9.58210700E-07 9.19795839E-07 -7.17760550E-07 -7.58941111E-11 7.15267089E-11 -5.53001301E-10 4.25051613E-07 -1.61680318E-06 -9.58182615E-07 -3.20378886E-10 2.12727721E-11 -7.58928249E-11 -1.61680339E-06 2.92708944E-07 9.19769109E-07 2.12727633E-11 -3.17858283E-10 7.15292774E-11 -9.58210693E-07 9.19795864E-07 -7.17760552E-07 -7.58930910E-11 7.15293743E-11 -5.52997605E-10 6.58033E-15 3.25156E-14 -3.27764E-14 -1.62555E-14 -2.68365E-15 -1.17847E-14 2.28546E-12 -2.16553E-12 -1.43745E-12 -2.67489E-15 3.79343E-15 7.93010E-15 -6.92870E-15 -2.43985E-14 2.18525E-15 -1.02002E-15 -2.66539E-15 -3.69589E-15 C.2 Verification of the Measurement Partial Derivative Matrix The measurement partial derivative matrix was verified as in Section C.1. The non-trivial results are listed in Table C.3. 120 Table C.3: Measurement Partial Differential Equation Validation Derivative H Matrix Entry Numerical Derivative Difference P x P y P z D x D y D z D x D y D z x y z x y z D x D y D z D x D y D z D x D y D z D x D y D z 7.84860376E-01 -3.44975293E-01 -5.14768327E-01 1.64229470E-04 -1.73000219E-04 3.66335628E-04 7.84860260E-01 -3.44975283E-01 -5.14768323E-01 -2.02226288E-03 6.38032661E-04 5.45557807E-04 -1.73277916E+00 1.76946906E-01 2.38128883E-01 6.16708027E-06 -3.09217499E-07 -1.81840626E-06 -2.02246147E-03 6.37983912E-04 5.45673581E-04 -2.44516126E-04 4.78947935E-06 2.90376513E-05 -1.73277918E+00 1.76946888E-01 2.38128946E-01 7.84860260E-01 -3.44975283E-01 -5.14768323E-01 1.64229470E-04 -1.73000200E-04 3.66335629E-04 7.84860260E-01 -3.44975283E-01 -5.14768323E-01 -2.02246146E-03 6.37983933E-04 5.45673565E-04 -1.73277916E+00 1.76946894E-01 2.38128962E-01 6.16705936E-06 -3.09234916E-07 -1.81841582E-06 -2.02246146E-03 6.37983933E-04 5.45673565E-04 -2.44516132E-04 4.78946488E-06 2.90376732E-05 -1.73277918E+00 1.76946888E-01 2.38128946E-01 1.15618E-07 -1.00706E-08 -4.35421E-09 1.65339E-13 -1.96636E-11 -1.16779E-12 -5.08993E-12 -1.92002E-12 1.13600E-11 1.98578E-07 4.87274E-08 -1.15758E-07 -4.49096E-09 1.24717E-08 -7.85587E-08 2.09169E-11 1.74173E-11 9.56467E-12 -7.15000E-12 -2.14815E-11 1.56880E-11 5.69595E-12 1.44737E-11 -2.18656E-11 6.41998E-12 -5.44001E-12 -5.42000E-12 121 C.3 Verification of Measurement Noise Settings Verification of the measurement noise settings in the filter was achieved by comparing the standard deviation on the measurement residuals to the onesigma setting used in the filter. Figure C.1 illustrates the residuals over one orbit, and Table C.4 provides the comparison of the filter settings. Pesudorange 15 10 5 Meters 0 -5 -10 -15 0 1000 2000 3000 Time (s) 4000 -2 0 1000 2000 3000 Time (s) 4000 Meters/Second 1 2 Doppler 0 -1 Double-Differenced Phase 1 0.015 0.01 Meters 0.005 0 -0.005 -0.01 -0.015 0 1000 2000 3000 Time (s) 4000 -1 0 Meters/Second 0.5 Double-Differenced Doppler 0 -0.5 1000 2000 3000 Time (s) 4000 Figure C.1: Measurement Residuals 122 Table C.4: Measurement Noise Measurement Residuals P(m) 0.8099 D(m/s) 0.1981 (m) 0.0012 D (m/s) 0.1055 Verification. Filter Setting 1 .25 .004 .15 C.4 Covariance Update Verification In order to verify the covariance update process was reducing the co- variance, it was plotted before and after the update. A typical result for the absolute case are depicted in Figure C.2. Figure C.3 illustrates a relative position covariance update. In both sets of figures, the blue dots are pre-update covariances and the red dots are after the update. 123 Absolute Position Covariance 40 30 20 10 Covariance 0 -10 -20 -30 -40 0 500 1000 1500 Time (s) 2000 2500 3000 3500 Figure C.2: Absolute Position Covariances Updates -4 8 x 10 Relative Position Covariance Example 7 6 5 Covariance (m) 4 3 2 1 0 0 500 1000 1500 2000 2500 Time (s) 3000 3500 4000 4500 Figure C.3: Relative State Covariance Updates 124 Bibliography [1] GPS Orion 12 Channel GPS Reference Design, 2.0 edition, October 2001. AN4808. [2] Bauer, F. et al. "Enabling Spacecraft Formation Flying Through Spaceborn GPS and Enhanced Automony Technologies". Navigation and Control, September 1999. [3] Binning, P. Absolute and Relative Satellite to Satellite Navigation Using GPS. PhD thesis, University of Colorado, 1997. [4] Brogran, W. Modern Control Theory, volume 3rd. Prentice Hall, 1991. [5] Busse, F. and J. How. "Real-Time Experimental Demonstration of Precise Decentralized Relative Navagation for Formation Flying Spacecraft". AIAA Guidance Navigation and Control, August 2002. [6] Carpenter, J. "Decentralized Control - Theoretical Development". Internal Memo. [7] Carpenter, J. "Partially Decentralized Control Architectures for Satellite Formations". Proc. of AIAA GNC Conference, Monterey, CA, August 2002. [8] Carpenter, J. Decentralized control of satellite formations. International Journal of Robust and Nonlinear Control, (12):141 161, 2002. AIAA Guidance 125 [9] Ebinuma, T. Precision Spacecraft Rendezvous Using Global Positioning System: An Integrated Hardware Approach. PhD thesis, The University of Texas at Austin, August 2001. [10] Ebinuma T., R. Bishop, and G. Lightsey. "Spacecraft Rendezvous Using GPS Relative Navigation". AAS/AIAA Spaceflight Mechanics Meeting, February 2001 Santa Barbara, CA. [11] Ebinuma, T., O. Montenbruck, and G. Lightsey. ION GPS-2002, September 2002. [12] Farrell, J. and M. Barth. The Global Positioning System and Inertial "Precise Spacecraft Relative Navigation Using Kinematic Inter-Spacecraft State Estimates". Navigation. McGraw-Hill, New York, 1999. [13] Folta, D. and A. Hawkins. hanced Formation Flying". "E0-1 Technology Validation Report Enhttp://eo1.gsfc.nasa.gov/miscPages/ TechForumReports/EFF-GSFC-2.pdf, 2001. [14] Gramling, C. et al. "Autonomous Relative Navigation For FormationFlying Sattelites Using GPS". http://geons.gsfc.nasa.gov/library_ docs/ISSFDrelnavfinal.PDF. [15] Hoffman-Wellenhof, H., H. Lichtenegger, and J. Collins. Global Positioning System: Theory and Practice. Springer-Verlag/Wien, New York, fifth edition, 2001. [16] Kaplan, E. Understanding GPS: Principles and Applications. House Publishers, Norwood, MA, 1996. Artech 126 [17] Kong, E. "Minimum Energy Trajectories for Techsat-21 Earth Orbiting Clusters". http://web.mit.edu/hilstad/www/school/16.851/FF_ 2-resizing.ppt, May 2001. [18] Lightsey, G. Development and Flight Demonstration of a GPS Receiver for Space. PhD thesis, Stanford University, February 1997. [19] Long, A. et al. "Relative Navigation Of Formation-Flying Satellites". http://geons.gsfc.nasa.gov/library_docs/ISFF02-relnav.pdf. [20] Montenbruck, O. and G. Eberhard. Satellite Orbits: Models, Mehtods, and Applications. Springer-Verlag, Berlin, Germany, 2000. [21] Montenbruck, O., M. Markgraf, S. Leung, and E. Gill. "A GPS Receiver for Space Application". ION GPS 2001, Salt Lake City, UT, Sep. 12-14 2001. [22] NASA. "Distributed Space Systems Mission List". December 2003. [23] NASA. "Formation Flying: The Afternoon `A-Train' Satellite Constellation". http://eospso.gsfc.nasa.gov/ftp_docs/A-Train_Fact_Sheet. pdf, December 2003. [24] Sabol C., R. Burns, and C. McLaughlin. "Satellite Formation Flying Design and Evolution". Spaceflight Mechanics, v 102 pp 265-384, February 1999. [25] AFRL Space Vehicles Directorate. html, june 1999. "Techsat-21 Space Missions Using Satellite Clusters". http://www.vs.afrl.af.mil/Factsheets/techsat21. 127 [26] Spirent. Spirent GSS Users Manual. [27] Vallado, D. Fundamentals of Astrodynamics and Applications. Microcosm Press, El Segundo, CA, second edition, 2001. [28] Bamford W., G. Lightsey, and T. Ebinuma. "Navigation of Large Autonomously Controlled Formations". TX, August 11, 2003. [29] Wolfe, J., M. Abdel-Hafez, and J. Speyer. "Effective Estimation of Relative Positions in Orbit Using Differential GPS". AIAA GNC Conference, Austin, 128 Vita William Alfred Bamford Jr. was born July 19, 1974 in San Lorenzo California to Margaret and William Bamford Sr. William graduated from San Lorenzo High School in 1992 and was awarded the position of Salutitorian. Following high school, he attended San Jose State University where he graduated Cum Laude with a B.S. in Aerospace Engineering . While at San Jose State, William worked on the Spartnik microsatellite, focusing on structure testing and assembly. Upon graduation, William enrolled at the University of Texas at Austin in the Aerospace Engineering Department. For his Master's Thesis, William worked on orbit determination of the GLONASS satellites using NASA's GIPSY/OASIS II software. William continued on for his PhD studying navigation of formations using GPS receivers as primary sensors. Permanent address: 440 Retford Dr. Severna Park, MD 21146 Email: bamford@csr.utexas.edu A This dissertation was typeset with L TEX by the author. A LT EX is a document preparation system developed by Leslie Lamport as a special version of Donald Knuth's TEX Program. 129

Find millions of documents here - Study Guides, Homework Solutions, Papers, Exam Answer Keys and more. Course Hero has millions of course related materials that will enable you to learn better, faster and get an A in all your courses.
Below is a small sample set of documents:

russellr74662.pdf
Path: Texas >> RUSSELLR >> 74662 Fall, 2009

lab9.pdf
Path: CSU San Bernardino >> CS >> 201 Fall, 2009
mukadama15106.pdf
Path: Texas >> MUKADAMA >> 15106 Fall, 2009
kellerkm71167.pdf
Path: Texas >> KELLERKM >> 71167 Fall, 2009
oxfordwt32223.pdf
Path: Texas >> OXFORDWT >> 32223 Fall, 2009
bennettl81291.pdf
Path: Texas >> BENNETTL >> 81291 Fall, 2009
engelas504835.pdf
Path: Texas >> ENGELAS >> 504835 Fall, 2009
curranma71134.pdf
Path: Texas >> CURRANMA >> 71134 Fall, 2009
stanleyk74304.pdf
Path: Texas >> STANLEYK >> 74304 Fall, 2009
protsenkode026.pdf
Path: Texas >> PROTSENKOD >> 026 Fall, 2009
Chapter07.outline.pdf
Path: Concordia NE >> PHYS >> 110 Fall, 2009
rutherfordg022.pdf
Path: Texas >> RUTHERFORD >> 022 Fall, 2009
auerbachs13838.pdf
Path: Texas >> AUERBACHS >> 13838 Fall, 2009
dechapanyaw029.pdf
Path: Texas >> DECHAPANYA >> 029 Fall, 2009
shoemakerdb042.pdf
Path: Texas >> SHOEMAKERD >> 042 Fall, 2009
johnsonam71217.pdf
Path: Texas >> JOHNSONAM >> 71217 Fall, 2009
sampselld77810.pdf
Path: Texas >> SAMPSELLD >> 77810 Fall, 2009
complex.txt
Path: CSU San Bernardino >> CS >> 330 Fall, 2009
cadenheadjk046.pdf
Path: Texas >> CADENHEADJ >> 046 Fall, 2009
benjaminsmr042.pdf
Path: Texas >> BENJAMINSM >> 042 Fall, 2009
simpsonal13317.pdf
Path: Texas >> SIMPSONAL >> 13317 Fall, 2009
hamiltont84490.pdf
Path: Texas >> HAMILTONT >> 84490 Fall, 2009
kotrlaka518287.pdf
Path: Texas >> KOTRLAKA >> 518287 Fall, 2009
harrisont86130.pdf
Path: Texas >> HARRISONT >> 86130 Fall, 2009
brandonjc99738.pdf
Path: Texas >> BRANDONJC >> 99738 Fall, 2009
MATH107A46024536.doc
Path: MD University College >> ASIA >> 2092 Fall, 2009
crawforda65881.pdf
Path: Texas >> CRAWFORDA >> 65881 Fall, 2009
achacosom07761.pdf
Path: Texas >> ACHACOSOM >> 07761 Fall, 2009
jarroldwl86380.pdf
Path: Texas >> JARROLDWL >> 86380 Fall, 2009
sharyginany026.pdf
Path: Texas >> SHARYGINAN >> 026 Fall, 2009
goncalvesac026.pdf
Path: Texas >> GONCALVESA >> 026 Fall, 2009
zieglerkj47418.pdf
Path: Texas >> ZIEGLERKJ >> 47418 Fall, 2009
burtnerjc90760.pdf
Path: Texas >> BURTNERJC >> 90760 Fall, 2009
alvarezla07232.pdf
Path: Texas >> ALVAREZLA >> 07232 Fall, 2009
MATH012A46124534.doc
Path: MD University College >> ASIA >> 2092 Fall, 2009
bonningew86532.pdf
Path: Texas >> BONNINGEW >> 86532 Fall, 2009
CMIS141AA44024445.doc
Path: MD University College >> ASIA >> 2092 Fall, 2009
CMIS102AA42086692.doc
Path: MD University College >> ASIA >> 2088 Fall, 2009
STAT200A42186896.doc
Path: MD University College >> ASIA >> 2088 Fall, 2009
kulkarnis86095.pdf
Path: Texas >> KULKARNIS >> 86095 Fall, 2009
chapmanbg60287.pdf
Path: Texas >> CHAPMANBG >> 60287 Fall, 2009
slattonkc78713.pdf
Path: Texas >> SLATTONKC >> 78713 Fall, 2009
michalskylo026.pdf
Path: Texas >> MICHALSKYL >> 026 Fall, 2009
batemanmt33508.pdf
Path: Texas >> BATEMANMT >> 33508 Fall, 2009
lodowskid97061.pdf
Path: Texas >> LODOWSKID >> 97061 Fall, 2009
raichlend29983.pdf
Path: Texas >> RAICHLEND >> 29983 Fall, 2009
perkinsjd44616.pdf
Path: Texas >> PERKINSJD >> 44616 Fall, 2009
mehdiabadinj026.pdf
Path: Texas >> MEHDIABADI >> 026 Fall, 2009
borisovasa86653.pdf
Path: Texas >> BORISOVASA >> 86653 Fall, 2009
Abuhakema504399.pdf
Path: Texas >> ABUHAKEMA >> 504399 Fall, 2009
hw03_solution.doc
Path: Penn State >> ME >> 581 Fall, 2009
oestreichj19588.pdf
Path: Texas >> OESTREICHJ >> 19588 Fall, 2009
evstatieve01477.pdf
Path: Texas >> EVSTATIEVE >> 01477 Fall, 2009
paschvaldesg042.pdf
Path: Texas >> PASCHVALDE >> 042 Fall, 2009
alvaradocg86236.pdf
Path: Texas >> ALVARADOCG >> 86236 Fall, 2009
martinssonpj026.pdf
Path: Texas >> MARTINSSON >> 026 Fall, 2009
makowitza504694.pdf
Path: Texas >> MAKOWITZA >> 504694 Fall, 2009
andersonmw81540.pdf
Path: Texas >> ANDERSONMW >> 81540 Fall, 2009
martinezrs39334.pdf
Path: Texas >> MARTINEZRS >> 39334 Fall, 2009
elshayebta87380.pdf
Path: Texas >> ELSHAYEBTA >> 87380 Fall, 2009
cowmeadowr17589.pdf
Path: Texas >> COWMEADOWR >> 17589 Fall, 2009
schougaardsb029.pdf
Path: Texas >> SCHOUGAARD >> 029 Fall, 2009
kordoskyma87090.pdf
Path: Texas >> KORDOSKYMA >> 87090 Fall, 2009
metcalfets016-x.pdf
Path: Texas >> METCALFETS >> 016 Fall, 2009
bocknackbm84986.pdf
Path: Texas >> BOCKNACKBM >> 84986 Fall, 2009
mahdjoubid26824.pdf
Path: Texas >> MAHDJOUBID >> 26824 Fall, 2009
vanderveenaa029.pdf
Path: Texas >> VANDERVEEN >> 029 Fall, 2009
crabtreejc17037.pdf
Path: Texas >> CRABTREEJC >> 17037 Fall, 2009
steubingdm73657.pdf
Path: Texas >> STEUBINGDM >> 73657 Fall, 2009
johnsonhl692102.pdf
Path: Texas >> JOHNSONHL >> 692102 Fall, 2009