Analysis of an Application where the Unscented Kalman Filter is Not Appropriate Except where reference is made to the work of others, the work described in this dissertation is my own or was done in collaboration with my advisory committee. This dissertation does not include proprietary or classi?ed information. Abby Anderson Certi?cate of Approval: Robert Dean, Co-Chair Assistant Professor Electrical and Computer Engineering A. Scottedward Hodel, Co-Chair Associate Professor Electrical and Computer Engineering Stanley Reeves Professor Electrical and Computer Engineering George Flowers Professor Mechanical Engineering George Flowers Dean Graduate School Analysis of an Application where the Unscented Kalman Filter is Not Appropriate Abby Anderson A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Ful?llment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama December 18, 2009 Analysis of an Application where the Unscented Kalman Filter is Not Appropriate Abby Anderson Permission is granted to Auburn University to make copies of this dissertation at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iii Dissertation Abstract Analysis of an Application where the Unscented Kalman Filter is Not Appropriate Abby Anderson Doctor of Philosophy, December 18, 2009 (M.S., Auburn University, 2006) (B.E.E., Auburn University, 2003) 241 Typed Pages Directed by A. Scottedward Hodel and Robert Dean In this work a spin stabilized rocket with a ring of lateral pulse jets for attitude correction and ?ns that open early in ?ight is simulated. The rocket is simulated with ?ve di?erent sensor packages: rate gyros only, rate gyros and an ideal magnetometer, rate gyros and a magnetometer, rate gyros and angle gyros, and rate gyros and angle gyros and a magnetometer. The gyros are microelectromechanical systems (MEMS) devices. All control e?ort must be applied in the ?rst seconds of ?ight because of the properties of the rocket. A comparison of the Extended Kalman Filter (EKF) andtheUnscented Kalman Filter (UKF) foreach sensorsuiteis presentedtodetermine the best approach to improve the circular probableerror (CEP) of the rocket. A solution to Wahba?s problem, the EStimator of the Optimal Quaternion (ESOQ) algorithm, is used for state estimation for certain sensor con?gurations. Wahba?s problem has traditionally been used for state estimation of orbiting satellites. iv Acknowledgments First, I would like to thank my family. Particularly my parents Je? and Rhonda Anderson. Without their support none of this would have been possible. I would also like to thank my committee: Dr. Hodel, Dr. Dean, Dr. Reeves, and Dr. Flowers. I would especially like to thank Dr. Dean for stepping into a di?cult situation to allow me to complete my work. But the person who is most directly responsible for the completion of this Ph.D. is Dr. Hodel. (Yes, Dr. Hodel, I realize I started the previous sentence with a conjunction, but note that there are no dangling modi?ers.) Over the years that I have been his student, we have had conversations ranging from work to our sometimes unorthodox religious views. No matter the topic, Babylon 5 always had a way of working itself into the discussion. As Delenn aptly said, ?The universe puts us in places where we can learn. They are never easy places, but they are right. Wherever we are, it?s the right place ... and the right time. Pain that sometimes comes is part of the process of constantly begin born.? I miss our conversations. I have been blessed not only to have Dr. Hodel as an adviser for both my graduate degrees but also to consider him my friend. One day in one of my undergraduate classes that Dr. Hodel was teaching, he said, ?My students? job is to make me look good.? I never quite let him forget that statement, but I hope this work (for which he suggested the title ?The Unscented Kalman Filter Stinks?) makes you look good, Dr. Hodel. It is dedicated to you. ?A student is not above his teacher, nor a servant above his master. It is enough for the student to be like his teacher and the servant like his master ...? (Matthew 10:24-25, NIV). ?And so it begins.? - Kosh, Babylon 5 v Style manual or journal used IEEE Transactions on Aerospace and Electronic Systems (together with the style known as ?auphd?). Bibliography follows IEEE Transactions. Computer software used The document preparation package TEX (speci?cally LATEX) together with the departmental style-?le auphd.sty. vi Table of Contents List of Figures x 1 Introduction 1 2 Wahba?s Problem 4 2.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Problem Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3 Davenport?s q-Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3.1 Rotation Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.2 q-Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 QUEST Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5 Kalman Filter Type Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5.1 Wahba?s Problem as Maximum Likelihood Estimation . . . . . . . . . . . . . 22 2.5.2 Filter QUEST Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.6 REQUEST Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.6.1 Time-Invariant REQUEST . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.6.2 Time-Varying REQUEST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.7 Extended QUEST Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.8 Energy Approach Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.9 Singular Value Decomposition Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 43 2.10 Fast Optimal Attitude Matrix (FOAM) Algorithm . . . . . . . . . . . . . . . . . . . 47 2.11 Alternative Quaternion Attitude Estimation Algorithm . . . . . . . . . . . . . . . . . 49 2.12 Euler-q Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.13 ESOQ Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.14 ESOQ2 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.15 A Slightly Sub-Optimal Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 2.16 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3 Magnetometer Navigation 65 3.1 Modeling Earth?s Magnetic Field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.1.1 Equation Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.1.2 Model Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.2 Magnetometer Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.3 Magnetometer Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.3.1 Attitude Independent Bias Estimation . . . . . . . . . . . . . . . . . . . . . . 72 3.3.2 TWOSTEP Algorithm Calibration . . . . . . . . . . . . . . . . . . . . . . . . 77 3.3.3 Recursive Least Squares Method . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.4 Angular Rate Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 vii 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4 Estimation And Control 88 4.1 Current Missile Control Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.2 Reaction Jet Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.2.1 Variable-Structure Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.2.2 Time-Optimal Control with A Single Reaction Jet . . . . . . . . . . . . . . . 92 4.2.3 Multiple-Reaction Jet Control . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.2.4 Projectile Linear Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5 Rocket Dynamics and State Estimation 105 5.1 Rotational Rate Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2 Pitch Rate and Yaw Rate Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.2.1 Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.3 State Estimation with a Direction Vector . . . . . . . . . . . . . . . . . . . . . . . . 111 5.3.1 Ideal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.3.2 Magnetometer Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.4 Estimation with Rate and Angle Gyros . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.4.1 EKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 5.5 Estimation with Magnetometer, Rate Gyros, and Angle Gyros . . . . . . . . . . . . 135 5.5.1 EKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 5.6 Wahba?s Problem Applied to a Rocket . . . . . . . . . . . . . . . . . . . . . . . . . . 140 5.6.1 Wahba?s Problem Solution as Input to an EKF and UKF . . . . . . . . . . . 141 5.6.2 Wahba?s Problem Solution Combined with an EKF and UKF . . . . . . . . . 143 6 Results 145 6.1 Rate Gyros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 6.2 Rate Gyros and an Ideal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 6.3 Rate Gyros and a Magnetometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 6.4 Rate Gyros and Angle Gyros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 6.5 Rate Gyros, Angle Gyros, and Magnetometer . . . . . . . . . . . . . . . . . . . . . . 167 6.5.1 ESOQ Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 6.5.2 EKF and UKF with All Sensors . . . . . . . . . . . . . . . . . . . . . . . . . 175 6.5.3 ESOQ Algorithm Combined with Kalman Filters . . . . . . . . . . . . . . . . 175 6.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 7 Conclusions and Future Work 197 7.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198 Bibliography 199 Appendices 205 viii A Rotation Sequences 206 A.1 Euler Angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 A.2 Quaternions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207 B Spherical Mathematics 210 B.1 Spherical Coordinate System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210 B.2 Legendre Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212 B.3 Laplace?s Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 B.3.1 Solution Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 B.3.2 Relationship with Legendre Functions and Spherical Harmonics . . . . . . . . 215 C Kalman Filtering 218 C.1 Linear Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 218 C.2 Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219 C.3 Unscented Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220 D Latitude and Longitude 224 D.1 Cartesian Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224 D.2 Geodetic Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 226 ix List of Figures 2.1 Inertial and Body Reference Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 5.1 Relationship Among Rotation Quaternion Components . . . . . . . . . . . . . . . . . 116 6.1 Scaled MSE of Roll Rates of KF and UKF Estimates from Rate Gyros . . . . . . . . 147 6.2 Scaled MSE of Pitch Rates of KF and UKF Estimates from Rate Gyros . . . . . . . 148 6.3 Scaled MSE of Yaw Rates of KF and UKF Estimates from Rate Gyros . . . . . . . . 148 6.4 Scaled MSE of Roll Angles of KF and UKF Estimates from Rate Gyros . . . . . . . 149 6.5 Scaled MSE of Pitch Angles of KF and UKF Estimates from Rate Gyros . . . . . . 149 6.6 Scaled MSE of Yaw Angles of KF and UKF Estimates from Rate Gyros . . . . . . . 150 6.7 Normalized CEP for KF Estimates from Rate Gyros . . . . . . . . . . . . . . . . . . 150 6.8 Normalized CEP for UKF Estimates from Rate Gyros . . . . . . . . . . . . . . . . . 151 6.9 Scaled MSE of Roll Rates of EKF and UKF Estimates from Rate Gyros and an Ideal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.10 Scaled MSE of Pitch Rates of EKF and UKF Estimates from Rate Gyros and an Ideal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.11 Scaled MSE of Yaw Rates of EKF and UKF Estimates from Rate Gyros and an Ideal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 6.12 Scaled MSE of Roll Angles of EKF and UKF Estimates from Rate Gyros and an Ideal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 6.13 Scaled MSE of Pitch Angles of EKF and UKF Estimates from Rate Gyros and an Ideal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 6.14 Scaled MES of Yaw Angles of EKF and UKF Estimates from Rate Gyros and an Ideal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 6.15 Normalized CEP for EKF Estimates from Rate Gyros and an Ideal Vector . . . . . . 156 x 6.16 Normalized CEP for UKF Estimates from Rate Gyros and an Ideal Vector . . . . . . 156 6.17 Scaled MSE of Roll Rates of EKF and UKF Estimates from Rate Gyros and a Magnetometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 6.18 Scaled MSE of Pitch Rates of EKF and UKF Estimates from Rate Gyros and a Magnetometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 6.19 Scaled MSE of Yaw Rates of EKF and UKF Estimates from Rate Gyros and a Magnetometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 6.20 Scaled MSE of Roll Angles of EKF and UKF Estimates from Rate Gyros and a Magnetometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 6.21 Scaled MSE of Pitch Angles of EKF and UKF Estimates from Rate Gyros and a Magnetometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 6.22 Scaled MSE of Yaw Angles of EKF and UKF Estimates from Rate Gyros and a Magnetometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 6.23 Normalized CEP for EKF Estimates from Rate Gyros and a Magnetometer . . . . . 161 6.24 Normalized CEP for UKF Estimates from Rate Gyros and a Magnetometer . . . . . 161 6.25 Scaled MSE of Roll Rates of EKF and UKF Estimates from Rate Gyros and Angle Gyros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 6.26 Scaled MSE of Pitch Rates of EKF and UKF Estimates from Rate Gyros and Angle Gyros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 6.27 Scaled MSE of Yaw Rates of EKF and UKF Estimates from Rate Gyros and Angle Gyros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 6.28 Scaled MSE of Roll Angles of EKF and UKF Estimates from Rate Gyros and Angle Gyros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 6.29 Scaled MSE of Pitch Angles of EKF and UKF Estimates from Rate Gyros and Angle Gyros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 6.30 Scaled MSE of Yaw Angles of EKF and UKF Estimates from Rate Gyros and Angle Gyros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 6.31 Normalized CEP for EKF Estimates from Rate Gyros and Angle Gyros . . . . . . . 166 6.32 Normalized CEP for UKF Estimates from Rate Gyros and Angle Gyros . . . . . . . 166 6.33 CEP for ESOQ Controlled Rocket with Ideal Angle Gyros and Magnetometer . . . . 168 xi 6.34 CEP for ESOQ Controlled Rocket with Nonideal Angle Gyros and Magnetometer . . 168 6.35 Actual Roll Rates (Normalized) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 6.36 ESOQ Estimated Roll Rates with Nonideal Sensors (Normalized) . . . . . . . . . . . 169 6.37 Actual Roll Angles (Normalized) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 6.38 ESOQ Estimated Roll Angles with Nonideal Sensors (Normalized) . . . . . . . . . . 170 6.39 Actual Pitch Rates (Normalized) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 6.40 ESOQ Estimated Pitch Rates with Nonideal Sensors (Normalized) . . . . . . . . . . 171 6.41 Actual Pitch Angles (Normalized) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 6.42 ESOQ Estimated Pitch Angles with Nonideal Sensors (Normalized) . . . . . . . . . . 172 6.43 Actual Yaw Rates (Normalized) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 6.44 ESOQ Estimated Yaw Rates with Nonideal Sensors (Normalized) . . . . . . . . . . . 173 6.45 Actual Yaw Angles (Normalized) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 6.46 ESOQ Estimated Yaw Angles with Nonideal Sensors (Normalized) . . . . . . . . . . 174 6.47 CEP for EKF Controlled Rocket with Rate Gyros, Angle Gyros, and Magnetometer 175 6.48 CEP for UKF Controlled Rocket with Rate Gyros, Angle Gyros, and Magnetometer 176 6.49 CEP for Rocket Controlled by Estimates from an EKF with Inputs from the ESOQ Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 6.50 CEP for EKF/ESOQ Hybrid Estimator Controlled Rocket . . . . . . . . . . . . . . . 178 6.51 CEP for UKF/ESOQ Hybrid Estimator Controlled Rocket . . . . . . . . . . . . . . 179 6.52 CEP for KF/ESOQ Hybrid Estimator Controlled Rocket . . . . . . . . . . . . . . . 179 6.53 Roll Rate MSE for ESOQ Estimates as Input to an EKF . . . . . . . . . . . . . . . 180 6.54 Pitch Rate MSE for ESOQ Estimates as Input to an EKF . . . . . . . . . . . . . . . 180 6.55 Yaw Rate MSE for ESOQ Estimates as Input to an EKF . . . . . . . . . . . . . . . 181 6.56 Roll Angle MSE for ESOQ Estimates as Input to an EKF . . . . . . . . . . . . . . . 181 6.57 Pitch Angle MES for ESOQ Estimates as Input to an EKF . . . . . . . . . . . . . . 182 xii 6.58 Yaw Angle MSE for ESOQ Estimates as Input to and EKF . . . . . . . . . . . . . . 182 6.59 Roll Rate MSE for EKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 184 6.60 Pitch Rate MSE for EKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 184 6.61 Yaw Rate MSE for EKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 185 6.62 Roll Angle MSE for EKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 185 6.63 Pitch Angle MSE for EKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . 186 6.64 Yaw Angle MSE for EKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 186 6.65 Roll Rate MSE for UKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 187 6.66 Pitch Rate MSE for UKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 188 6.67 Yaw Rate MSE for UKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 188 6.68 Roll Angle MSE for UKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 189 6.69 Pitch Angle MSE for UKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . 189 6.70 Yaw Angle MSE for UKF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 190 6.71 Roll Rate MSE for KF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . . 191 6.72 Pitch Rate MSE for KF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . . 192 6.73 Yaw Rate MSE for KF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . . 192 6.74 Roll Angle MSE for KF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . . 193 6.75 Pitch Angle MSE for KF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . 193 6.76 Yaw Angle MSE for KF/ESOQ Hybrid . . . . . . . . . . . . . . . . . . . . . . . . . . 194 B.1 The Spherical Coordinate System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 D.1 Derivation of Latitudes and Longitudes . . . . . . . . . . . . . . . . . . . . . . . . . 225 D.2 Geodetic Coordinate System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 227 xiii Chapter 1 Introduction Increasing the accuracy of munitions while reducing cost is of key interest to the military. One method is to use onboard guidance systems that can steer a missile to within meters or less of its target. However, this approach is expensive. A less costly but less accurate alternative is the spin stabilization approach. These systems lack closed-loop ?ight control and many of the sophisticated sensors found in onboard guidance systems. The accuracy of spin stabilized rockets is a?ected by many factors including motor misalignments, tip-o? error, and wind. One method to improve the accuracy of spin stabilized rockets is to include microelectrome- chanical system (MEMS) sensors. MEMS sensors have the advantages of low cost since they are batch producible using conventional IC technology and light weight. However, they are subject to various errors including constant biases, walking biases, and additive noise. In this work we compare state estimation techniques to improve the accuracy of a spin sta- bilized rocket equipped with various sensor suites which include MEMS devices. We compare the improvement in rocket accuracy when states are estimated by a Kalman Filter (KF), an Extended Kalman Filter (EKF), and an Unscented Kalman Filter (UKF) for various sensor con?gurations. We also use the EStimator of the Optimal Quaternion (ESOQ) algorithm, a solution to Wahba?s problem, alone and in conjunction with the EKF and UKF to estimate states with certain sensor suites. The ?rst sensor con?guration consists only of MEMS gyros to measure the rocket?s rota- tional rates. The second sensor suite consists of MEMS gyros and an ideal sensor that measures a known inertial frame vector in the body frame. The third sensor suite is comprised of MEMS gyros and a tri-axial magnetometer. The next sensor suite consists of MEMS rate gyros and MEMS gyros that directly measure rotational angles. The ?nal sensor suite combines the rate gyros, angle gyros, and magnetometer into one package. 1 For this work we make several assumptions about the rocket. The rocket is launched from a stationary ground launcher. A main thruster propels the rocket, and a ring of lateral thrusters located at the rocket?s rear supply torque for control. Fins are located at the back of the rocket that open shortly after launch. Due to aerodynamic forces, control of such a rocket is only feasible during the ?rst seconds of ?ight so state estimation is con?ned to that time interval. In this dissertation we present previous work pertaining to improving munition accuracy, de- velop system models for the various sensor suites, and present the results of state estimation. In Chapter 2 we present several solutions to Wahba?s problem, a method of determining a vehicle?s attitude using pairs of sensors such as accelerometers and magnetometers. Next, in Chapter 3 we review concepts associated with navigation via magnetometers. In Chapter 4 we present methods for controlling rockets with a focus on rockets with lateral jets. We develop the system equations for our various rocket models in Chapter 5. We present the results of simulating the systems from Chapter 5 in Chapter 6 and draw some conclusions. In this work we make the following contributions: ? We estimate the states of a rocket that is controlled solely with a ring of lateral thrusters. ? We concentrate on estimation of a rocket?s states during the ?rst seconds of ?ight. Other work focuses on controlling munitions later in ?ight. ? We compare the e?ects that various sensor suites have on munition accuracy. ? We use a magnetometer to estimate the states of a rocket. Magnetometers are conventionally used on satellites rather than rockets. Usually a satellite has a model of the magnetic ?eld of the body it is orbiting on board. A rocket, however, does not have a magnetic ?eld model available. ? We apply a solution of Wahba?s problem to estimate the states of a rocket. Wahba?s prob- lem is traditionally applied to orbiting satellites since magnetic ?eld and gravitational ?eld measurements can be found with magnetometers and accelerometers. 2 ? We compare the performances of various state estimators including the KF, EKF, UKF, and ESOQ algorithm to estimate the states of a rotating rocket and present conclusions about which estimator performs best. 3 Chapter 2 Wahba?s Problem Wahba?s problem is a minimization problem most often applied to ?nd the attitude of space- craft with a ?xed attitude. Measurements are taken in the spacecraft?s body frame that give the angles between known objects and the spacecraft. Often the known objects are stars and the sen- sors are star trackers. The di?erence between the body frame measurements and known inertial frame values are used to derive an attitude solution for the spacecraft. The relationship between the body frame and the inertial frame is illustrated in Figure 2.1. The origin of the inertial frame (with axes denoted by a subscript u1D456) is located at the center of the earth while the origin of the body frame (with axes denoted by a subscript u1D44F) is located at the center of gravity of the rocket. Various algorithms have been developed based on Wahba?s cost function. While the cost function as originally posed is for static systems, researchers have extended it for use with dynamic systems. 2.1 Problem Statement Grace Wahba ?rst presented what has become known as Wahba?s Problem in 1965 in SIAM Re- view [1]. The problem is as follows: Given two vector sets {v1, v2, ... ,vu1D45B} and {v?1, v?2, ... ,v?u1D45B} with u1D45B ? 2 entries, ?nd the rotation matrix u1D440 (an orthogonal matrix with determinant +1) that minimizes u1D45Buni2211.alt02 u1D457=1 ?v?u1D457 ?u1D440vu1D457?2. In other words, ?nd the rotation matrix u1D440 that minimizes the mean square error between the two vector sets. A solution to this problem means that given at least two noncolinear, nonzero measurements (such as magnetic ?eld and gravitational acceleration) in the body frame of a moving vehicle and known values of the measured quantities in a stationary reference frame (i.e. an inertial reference frame), the attitude of the vehicle can be found via the rotation matrix u1D440. 4 Figure 2.1: Inertial and Body Reference Frames 5 2.2 Problem Solutions Various solutions to Wahba?s Problem were published in 1966 [2]. The solution given by Farrell and Stuelpnagel, based on the polar decomposition, proceeds as follows. Let the column vectors v1, ... ,vu1D45B, v?1, ... ,v?u1D45B have dimension u1D458, and let u1D449 and u1D449? denote the two u1D458?u1D45B matrices formed from the two vector sets. Then the problem can be rewritten as u1D444(u1D440) = u1D45Buni2211.alt02 u1D457=1 ?v?u1D457 ?u1D440vu1D457?2 = tr((u1D449? ?u1D440u1D449)u1D447(u1D449? ?u1D440u1D449)), (2.1) where u1D444(u1D440) is de?ned as the sum of squares to be minimized. Proof. Let a u1D458?u1D45B matrix u1D434 be de?ned as [a1 a2 ... au1D45B] where au1D457 are column vectors of length u1D458. Then u1D434u1D447u1D434= ? ?? ?? ?? ?? au1D4471 au1D4472 ... au1D447u1D45B ? ?? ?? ?? ?? uni005B.alt03 a1 a2 ... au1D45B uni005D.alt03 = ? ?? ?? ?? ?? au1D4471 a1 au1D4471 a2 ... au1D4471au1D45B au1D4472 a1 au1D4472 a2 ... au1D4472au1D45B ... ... ... ... au1D447u1D45Ba1 au1D447u1D45Ba2 ... au1D447u1D45Bau1D45B ? ?? ?? ?? ?? tr{u1D434u1D447u1D434} = au1D4471 a1 +au1D4472 a2 +...+au1D447u1D45Bau1D45B = u1D45Buni2211.alt02 u1D457=1 ?au1D457?2 Expand expression (2.1) so that u1D444(u1D440) becomes u1D444(u1D440) = tr((u1D449?u1D447 ?u1D449u1D447u1D440u1D447)(u1D449? ?u1D440u1D449)) = tr(u1D449?u1D447u1D449?)+tr(u1D449?u1D447u1D440u1D449)?tr(u1D449u1D447u1D440u1D447u1D449?) +tr(u1D449u1D447u1D440u1D447u1D440u1D449). Since u1D440 is orthogonal (i.e. u1D440u1D447u1D440 = u1D440u1D440u1D447 = u1D43C) and tr(u1D434u1D435u1D436) = tr(u1D436u1D447u1D435u1D447u1D434u1D447), Q(M) can be written as u1D444(u1D440) = tr(u1D449?u1D447u1D449?) +tr(u1D449u1D447u1D449)?2tr(u1D449u1D447u1D440u1D447u1D449?). (2.2) 6 Since only the last term of Equation (2.2) is dependent on u1D440, u1D444(u1D440) is minimized when u1D439(u1D440) = tr(u1D449u1D447u1D440u1D447u1D449?) is maximized. u1D439(u1D440) may be written as u1D439(u1D440) = tr(u1D440u1D447u1D449?u1D449u1D447) because tr(u1D434u1D435u1D436) = tr(u1D435u1D436u1D434) = tr(u1D436u1D434u1D435) when the matrices? dimensions are conformable. By the polar decomposition, u1D449?u1D449u1D447 can be written as u1D448u1D443 where u1D448 is orthogonal (and unique if u1D449?u1D449u1D447 is nonsingular) and u1D443 is symmetric and positive semide?nite. So u1D439(u1D440) = tr(u1D440u1D447u1D448u1D443). Because u1D443 is a real symmetric matrix, it has a spectral decomposition u1D443 = u1D441u1D437u1D441u1D447 where u1D441 is an orthogonal matrix and u1D437 = diag(u1D4511,...,u1D451u1D45B) with u1D451u1D456 ?u1D451u1D456+1. De?ning u1D44B = u1D441u1D440u1D447u1D448u1D441u1D447 results in u1D439(u1D440) = tr(u1D440u1D447u1D448u1D441u1D447u1D437u1D441) = tr(u1D441u1D440u1D447u1D448u1D441u1D447u1D437) = tr(u1D44Bu1D437) = u1D458uni2211.alt02 u1D456=1 u1D451u1D456u1D465u1D456u1D456. u1D439(u1D440) attains its maximum value when the elements u1D465u1D456u1D456 are all at their maximum value. Since u1D44B is an orthogonal matrix, all of its elements lie between the values 1 and -1. This means the diagonal elements are maximized when they are equal to 1, which means the u1D439(u1D440) is the identity matrix. Since u1D440 is a rotation matrix, its determinant is required to be 1. Thus the determinant of u1D44B is ?u1D44B? = ?u1D441u1D440u1D447u1D448u1D441u1D447? = ?u1D441??u1D440u1D447??u1D448??u1D441u1D447? = ?u1D441?2?u1D440??u1D448? = ?u1D448?. If det(u1D448) = 1, then u1D44B = u1D43C maximizes u1D439(u1D440). If det(u1D448) = ?1, then det(u1D44B) = ?1, and a solution is u1D44B = ? ?? u1D43Cu1D458?1 0 0 ?1 ? ??. 7 De?ne u1D44B0 as the matrix which maximizes u1D439(u1D440) according to the determinant of u1D448 so that u1D44B0 = u1D441u1D440u1D4470 u1D448u1D441u1D447. Then the rotation matrix that minimizes the sum of squares of u1D444(u1D440) is u1D4400 = u1D448u1D441u1D447u1D44Bu1D4470 u1D441. The matrix is unique if u1D449?u1D449u1D447 is nonsingular. Wessner o?ers another solution to Wahba?s problem [2]. Like Farrell and Stuelpnagel, Wessner recasts the problem to maximize u1D439(u1D440) = tr(u1D440u1D447u1D449?u1D449u1D447). If u1D449?u1D449u1D447 is nonsingular, then u1D449?u1D449u1D447 has the polar decomposition u1D449?u1D449u1D447 = u1D434 = u1D448u1D443. De?ne u1D448 = (u1D434u1D447)?1(u1D434u1D447u1D434)1/2 and u1D443 = (u1D434u1D447u1D434)1/2 where (u1D434u1D447u1D434)1/2 is the symmetric square root of u1D434u1D447u1D434 with positive eigenvalues. From Farrell and Stuelpnagel, the optimal solution is u1D440u1D45C = u1D448u1D441u1D447u1D44Bu1D45Cu1D441 where u1D441 is orthogonal. Wessner assumes that the determinant of u1D434 is positive; thus, u1D44Bu1D45C = u1D43C. Then u1D440u1D45C = u1D448u1D441u1D447u1D441 = u1D448, which results in u1D440u1D45C = (u1D434u1D447)?1(u1D434u1D447u1D434)1/2 = (u1D449?u1D449u1D447)?1(u1D449?u1D449u1D447u1D449?u1D449u1D447)1/2 2.3 Davenport?s q-Method Another method for the solution of Wahba?s problem is the u1D45E-method developed by Paul Davenport in 1968 [3]. Davenport de?nes a rotation vector to solve the problem. We discuss the q-method in two parts. We ?rst present an overview of rotation vectors, followed by a discussion of Wahba?s problem. Davenport often uses the notation u1D44B2 to represent u1D44Bu1D447u1D44B where u1D44B is a vector. Thus, in Davenport?s notation, u1D44B2 is a scalar value not to be confused with the vector value u1D44B. For simplicity, we use Davenport?s notation in this discussion. 8 2.3.1 Rotation Vectors A rotation of an angle 0 ? u1D703 ? u1D70B about a unit vector u1D44B = uni005B.alt03 u1D4651 u1D4652 u1D4653 uni005D.alt03 ? ?3 can be represented as the 3?3 matrix operator u1D445u1D465(u1D703) = cos(u1D703) ? ?? ?? ? 1 0 0 0 1 0 0 0 1 ? ?? ?? ? +(1?cos(u1D703)) ? ?? ?? ? u1D46521 u1D4651u1D4652 u1D4651u1D4653 u1D4651u1D4652 u1D46522 u1D4652u1D4653 u1D4651u1D4653 u1D4652u1D4653 u1D46523 ? ?? ?? ? + sin(u1D703) ? ?? ?? ? 0 u1D4653 ?u1D4652 ?u1D4653 0 u1D4651 u1D4652 ?u1D4651 0 ? ?? ?? ? (2.3) and u1D445 = (u1D435u1D447)?1u1D435 (2.4) where u1D435 = ? ?? ?? ? 1 u1D4653 tanuni0028.alt01u1D7032uni0029.alt01 ?u1D4652 tanuni0028.alt01u1D7032uni0029.alt01 ?u1D4653 tanuni0028.alt01u1D7032uni0029.alt01 1 u1D4651 tanuni0028.alt01u1D7032uni0029.alt01 u1D4652 tanuni0028.alt01u1D7032uni0029.alt01 ?u1D4651 tanuni0028.alt01u1D7032uni0029.alt01 1 ? ?? ?? ? . (2.5) In order to simplify notation, de?ne u1D44C = tan uni0028.alt03u1D703 2 uni0029.alt03 u1D44B and u1D44D = sin uni0028.alt03u1D703 2 uni0029.alt03 u1D44B from which we obtain sinuni0028.alt01u1D7032uni0029.alt01 = ?u1D44C2 ?1+u1D44C2 = ?u1D44D2, cosuni0028.alt01u1D703 2 uni0029.alt01 = 1 ?1+u1D44C2 = ?1?u1D44D2, sin(u1D703) = 2?u1D44C2 1+u1D44C2 = 2uni221A.alt01u1D44D2(1?u1D44D2), and cos(u1D703) = 1?u1D44C21+u1D44C2 = 1?2u1D44D2. Hence Equation (2.3) can be rewritten as u1D445 = 11+u1D44C2 ? ?? ?? ? (1?u1D44C2) ? ?? ?? ? 1 0 0 0 1 0 0 0 1 ? ?? ?? ? +2 ? ?? ?? ? u1D46621 u1D4661u1D4662 u1D4661u1D4663 u1D4661u1D4662 u1D46622 u1D4662u1D4663 u1D4661u1D4663 u1D4662u1D4663 u1D46622 ? ?? ?? ? +2 ? ?? ?? ? 0 u1D4663 ?u1D4662 ?u1D4663 0 u1D4661 u1D4662 ?u1D4661 0 ? ?? ?? ? ? ?? ?? ? (2.6) 9 or u1D445 = (1?2u1D44D2) ? ?? ?? ? 1 0 0 0 1 0 0 0 1 ? ?? ?? ? +2 ? ?? ?? ? u1D46721 u1D4671u1D4672 u1D4671u1D4673 u1D4671u1D4672 u1D46722 u1D4672u1D4673 u1D4671u1D4673 u1D4672u1D4673 u1D46723 ? ?? ?? ? +2 uni221A.alt01 1?u1D44D2 ? ?? ?? ? 0 u1D4673 ?u1D4672 ?u1D4673 0 u1D4671 u1D4672 ?u1D4671 0 ? ?? ?? ? . (2.7) Alternatively, u1D435 may be written in terms of the elements of the vector u1D44C as u1D435 = ? ?? ?? ? 1 u1D4663 ?u1D4662 ?u1D4663 1 u1D4661 u1D4662 ?u1D4661 1 ? ?? ?? ? . (2.8) Let u1D445 be a rotation matrix, and separate u1D445 into symmetric and skew symmetric parts u1D445 = 12(u1D445+u1D445u1D447)+ 12(u1D445?u1D445u1D447) (2.9) where u1D445+u1D445u1D447 is symmetric and u1D445?u1D445u1D447 is skew symmetric. Then a comparison of Equation (2.9) to Equation (2.7) yields nine conditions for the vector u1D44D. Taking the trace u1D70E of Equation (2.7) yields u1D70E = 3?4u1D44D2 = 1+2cos(u1D703) (2.10) and inspection of the skew symmetric portion of u1D445 yields 2 uni221A.alt01 1?u1D44D2u1D44D = 12 ? ?? ?? ? u1D45F23 ?u1D45F32 u1D45F31 ?u1D45F13 u1D45F12 ?u1D45F21 ? ?? ?? ? . (2.11) Rearrange Equations (2.10) and (2.11) to obtain u1D44D2 = 3?u1D70E4 and u1D44D = 12?1+u1D70E ? ?? ?? ? u1D45F23 ?u1D45F32 u1D45F31 ?u1D45F13 u1D45F12 ?u1D45F21 ? ?? ?? ? . 10 We write u1D44D2 in terms of u1D703 which results in u1D44D2 = 12 ? 12 cos(u1D703), which means that u1D44D2 ? 1. Thus, the mapping de?ned by Equation (2.7) is a mapping of three- dimensional vectors over the ?eld of real numbers whose Euclidean length is less than or equal to one (denote this set of vectors byu1D701) onto the group of rotation matrices. The mapping is one-to-one except when u1D44D2 = 1 (u1D70E = ?1). Since u1D44C = u1D44Dcos(u1D703 2) , u1D44C = 1?1?u1D44D2u1D44D. Then using the same methods as above u1D44C = 11+u1D70E ? ?? ?? ? u1D45F23 ?u1D45F32 u1D45F31 ?u1D45F13 u1D45F12 ?u1D45F21 ? ?? ?? ? . u1D44C is unde?ned when u1D70E = ?1. To avoid this singularity, allow vectors of in?nite magnitude whose direction is given by a unit vector u1D44B. When u1D70E = ?1, Equation (2.6) reduces to u1D445 = ?u1D43C +2u1D44Bu1D44Bu1D447. (2.12) Let u1D702 denote the set of all real three-dimensional vectors augmented by the vectors of in?nite magnitude just discussed. Then Equations (2.6) and (2.12) de?ne a mapping from u1D702 onto the group of rotation matrices. Thus, either u1D701 or u1D702 may be used to parametrize the group of rotations. To distinguish between the new type of vector and ordinary vectors, Davenport de?nes the rotation vector: Given a set u1D6FF of real three-dimensional vectors and a mappingu1D70F that maps u1D6FF onto the group of rotation matrices, then two elements of u1D6FF are said to be equivalent if they map into the same rotation matrix. Normal vector operations are applied to rotation vectors. Rotation vectors can be combined to yield new rotation vectors. Given two rotation vectors u1D44C1 and u1D44C2, there are two associated rotation matrices u1D4451 and u1D4452. Then u1D445 = u1D4452u1D4451 is also a rotation 11 matrix, and there exists an associated rotation vector u1D44C. This vector u1D44C is given by u1D44C = 11?u1D44C 1u1D44C2 (u1D44C1 +u1D44C2 +u1D44C1 ?u1D44C2). (2.13) Similarly, for the two rotation vectors u1D44D1 and u1D44D2 which de?ne rotations u1D4451 and u1D4452, there exists a u1D44D that gives u1D445 = u1D4452u1D4452. u1D44D is de?ned by u1D44D = sgn uni0028.alt03uni221A.alt02 1?u1D44D21 uni221A.alt02 1?u1D44D22 ?u1D44D1 ?u1D44D2 uni0029.alt03 u1D44D0 (2.14a) u1D44D0 = uni221A.alt02 1?u1D44D22u1D44D1 + uni221A.alt02 1?u1D44D21u1D44D2 +u1D44D1 ?u1D44D2. (2.14b) This leads to the de?nition of Davenport?s rotation product: Let u1D6FF be a set of vectors, let ? be a binary operation on u1D6FF, and let u1D70F be a mapping of u1D6FF onto the group of rotation matrices. Then ? is said to be a rotation product if it is preserved by u1D70F, i.e. if u1D70F(u1D449 ?u1D44A) = u1D70F(u1D449)u1D70F(u1D44A) for every u1D449 and u1D44A in u1D6FF. Davenport also develops other useful relationships between rotation matrices, vectors, and rotation vectors. Assume a vector u1D449 is rotated by u1D445 to yield u1D449? = u1D445u1D449. If u1D44C and u1D44D are the rotation vectors de?ning u1D445, then u1D449? is u1D449? = u1D44Cu1D449 = 11+u1D44C2[(1?u1D44C2)u1D449 +2(u1D449 ?u1D44C)u1D44C +2u1D449 ?u1D44C], u1D44C2 2. When det(u1D435) = 0, u1D7061 = 0, and the eigenvectors associated with the eigenvalues ?u1D7061 = 0 cannot be discriminated. This can be overcome due to the orthogonality of u1D448 and u1D437. Once (u2,d2) and (u3,d3) are found u1 = ?2u2 ?u3 d1 = ?2d2 ?d3. 40 An alternate solution form that does not have singularities can be found by performing an eigenanalysis on Equations (2.63a)-(2.63b). The resulting u1D448 and u1D437 matrices are orthonormal, and the conditions of Equation (2.64) do not hold. The ?rst alternative solution form is u1D434u1D43Cu1D43C = u1D448u1D437?1 = u1D448u1D437u1D447. The eigenvectors u1D448 = uni005B.alt03 u1 u2 u3 uni005D.alt03 andu1D437 = uni005B.alt03 d1 d2 d3 uni005D.alt03 must refer to the same eigenvalues sequence 0 ?u1D7061 ?u1D7062 ?u1D7063, and the condition det(u1D448) = det(u1D437) must be satis?ed. Since u1D435u1D435u1D447 and u1D435u1D447u1D435 have the same eigenvalues, let u1D440 = u1D435u1D435u1D447 (or u1D440 = u1D435u1D447u1D435). Then the characteristic equation of u1D440 is u1D7063 +u1D4502u1D7062 +u1D4501u1D706+u1D4500 where u1D4502 = ?tr(u1D440) = ?u1D45A11 ?u1D45A22 ?u1D45A33 u1D4501 = tr(adj(u1D440)) = u1D45A11u1D45A22 +u1D45A11u1D45A33 +u1D45A22u1D45A33 ?u1D45A212 ?u1D45A213 ?u1D45A223 u1D4500 = ?det(u1D440) = u1D45A11u1D45A22u1D45A13 +2u1D45A12u1D45A13u1D45A23 ?u1D45A222u1D45A213 ?u1D45A11u1D45A223 ?u1D45A33u1D45A212. u1D440 is symmetric so all of its eigenvalues are real. Setting u1D45D = uni221A.alt01 (u1D4502/3)2 ?u1D4501/3 u1D45E = [u1D4501/2?(u1D4502/3)2]u1D4502/3?u1D4500/2 u1D467 = 13 cos?1(u1D45E/u1D45D3), then ? ?? ?? ? u1D7061 u1D7062 u1D7063 ? ?? ?? ? = ?u1D45D ? ?? ?? ? cos(u1D467) +?3sin(u1D467) cos(u1D467)??3sin(u1D467) ?2cos(u1D467) ? ?? ?? ? ? u1D45023 . 41 When u1D45B = 2, u1D4500 = ?det(u1D440) = 0. Then (u1D440 ?u1D706u1D456u1D43C)fu1D456 = 0 where fu1D456 is uu1D456 or du1D456 depending on the choice of u1D440. The row vectors of (u1D440 ?u1D706u1D456u1D43C) must lie on the same plane. The solution fu1D456 can be computed by normalizing the cross product between two row vectors of the matrix (u1D440 ?u1D706u1D456u1D43C). It is desired to compute fu1D456 from the cross product having the highest modulus. The 3 ? 3 symmetric matrix (u1D440 ?u1D706u1D456u1D43C) can be written as u1D440 ?u1D706u1D456u1D43C = ? ?? ?? ? mu1D4471 mu1D4472 mu1D4473 ? ?? ?? ? = ? ?? ?? ? u1D45Au1D44E u1D45Au1D465 u1D45Au1D466 u1D45Au1D465 u1D45Au1D44F u1D45Au1D467 u1D45Au1D466 u1D45Au1D467 u1D45Au1D450 ? ?? ?? ? . Then the eigenvector u1D453u1D456 is chosen from among e1 = m2 ?m1 = uni005B.alt03 u1D45Au1D44Fu1D45Au1D450 ?u1D45A2u1D467 u1D45Au1D466u1D45Au1D467 ?u1D45Au1D465u1D45Au1D450 u1D45Au1D465u1D45Au1D467 ?u1D45Au1D466u1D45Au1D44F uni005D.alt03u1D447 e2 = m3 ?m1 = uni005B.alt03 u1D45Au1D466u1D45Au1D467 ?u1D45Au1D465u1D45Au1D450 u1D45Au1D44Eu1D45Au1D450 ?u1D45A3u1D466 u1D45Au1D465u1D45Au1D466 ?u1D45Au1D467u1D45Au1D44E uni005D.alt03u1D447 e3 = m1 ?m2 = uni005B.alt03 u1D45Au1D465u1D45Au1D467 ?u1D45Au1D466u1D45Au1D44F u1D45Au1D465u1D45Au1D466 ?u1D45Au1D467u1D45Au1D44E u1D45Au1D44Eu1D45Au1D44F ?u1D45A2u1D465 uni005D.alt03u1D447 which are parallel. The eu1D456 with highest modulus has the maximum element u1D45Du1D458 u1D45D1 = (u1D45Au1D44Fu1D45Au1D450 ?u1D45A2u1D467)2 u1D45D2 = (u1D45Au1D44Eu1D45Au1D450 ?u1D45A2u1D466)2 u1D45D3 = (u1D45Au1D44Eu1D45Au1D44F ?u1D45A2u1D465)2. Then fu1D456 = eu1D458/ uni221A.alt02 eu1D447u1D458eu1D458. By this method the evaluation of the eigenvectors f2 and f3 is su?cient. Then f1 = f2 ?f3. Mortari gives several other possible solution forms. Other previously known solutions can be derived from the R-equation as well. Appropriate manipulation of the R-equations gives the ?direct 42 solution? developed by Stuelpnagel as well as a solution based on the singular value decomposition of u1D435. Some algorithms become singular when the u1D435 matrix is singular (as is always the case when u1D45B = 2). Mortari o?ers a general method to eliminate the singularity for the case u1D45B = 2. When u1D45B = 2, it is possible to add the unit vector s3 = (s1 ? s2)/sin(u1D717u1D460) and the associated vector v3 = (v1 ?v2)/sin(u1D717u1D463) to the data without a?ecting the computed optimal attitude. So s3 and v3 can be added to the data to avoid a singularity. New weights u1D6FD1 and u1D6FD2 that replace u1D6FC1 and u1D6FC2 are such that u1D6FD1/u1D6FD2 = u1D6FC1/u1D6FC2 and u1D6FD1 +u1D6FD2 +u1D6FD3 = 1. Mortari suggests the values u1D6FD1 = 2u1D6FC1/3, u1D6FD2 = 2u1D6FC2/3, andu1D6FD3 = 1/3 because these values maximize the distance from the singularity provided by the value of ?det(u1D435)?. 2.9 Singular Value Decomposition Algorithm In [10] a singular value decomposition approach is used to solve Wahba?s problem by ?nding the matrix u1D434u1D45Cu1D45Du1D461 that minimizes u1D43F(u1D434) = 12 u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456?bu1D456 ?u1D434ru1D456?2 (2.66) where ru1D456 are inertial frame unit vectors, bu1D456 are the corresponding unit vectors in the body frame, u1D44Eu1D456 are weights, and u1D45B is the number of observations. The weights are normalized so that u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456 = 1. (2.67) 43 Then u1D43F(u1D434) = 12 u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456(bu1D456 ?u1D434ru1D456)u1D447(bu1D456 ?u1D434ru1D456) = 12 u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456(bu1D447u1D456 ?ru1D447u1D456 u1D434u1D447)(bu1D456 ?u1D434ru1D456) = 12 u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456[bu1D447u1D456 ?2bu1D447u1D456 u1D434ru1D456 +ru1D447u1D456 u1D434u1D447u1D434ru1D456] = 12 u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456[2?2bu1D447u1D456 u1D434ru1D456] = u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456 ? u1D45Buni2211.alt02 u1D456=1 bu1D447u1D456 u1D434ru1D456 = 1? u1D45Buni2211.alt02 u1D456=1 bu1D447u1D456 u1D434ru1D456. Then the cost function can then be written as u1D43F(u1D434) = 1?tr(u1D434u1D435u1D447) (2.69) where u1D435 = u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456bu1D456ru1D447u1D456 . (2.70) The singular value decomposition of the matrix u1D435 is given by u1D435 = u1D448u1D446u1D449u1D447 where u1D448 and u1D449 are orthogonal matrices and u1D446 = diag(u1D4601,u1D4602,u1D4603) u1D4601 ?u1D4602 ?u1D4603 ? 0. 44 De?ne u1D448+ = u1D448[diag(1,1,det(u1D448))] u1D449+ = u1D449[diag(1,1,det(u1D449))] and u1D44A = u1D448u1D447+u1D434u1D449+ = cos(?)u1D43C +(1?cos(?))eeu1D447 ?sin(?)[u1D452?] where [u1D452?] = ? ?? ?? ? 0 ?u1D4523 u1D4522 u1D4523 0 ?u1D4521 ?u1D4522 u1D4521 0 ? ?? ?? ? . Thus u1D44A can be represented by a Euler axis/angle rotation with unit vector e and rotation angle ?. De?ne u1D446? = diag(u1D4601,u1D4602,u1D451u1D4603) (2.71) u1D451 = det(u1D448)det(u1D449) = ?1. (2.72) Then u1D435 can be written as u1D435 = u1D448+u1D446?u1D449u1D447+. (2.73) Substitute Equation (2.73) into Equation (2.69) to get u1D43F(u1D434) = 1?tr(u1D446?u1D44A) = 1?cos(?)tr(u1D446?)?(1?cos(?))[u1D4601u1D45221 +u1D4602u1D45222 +u1D451u1D4603u1D45223]+tr(u1D446?)?tr(u1D446?) = 1?tr(u1D446?) +(1?cos(?))[u1D4601 +u1D4602 +u1D451u1D4603 ?u1D4601u1D45221 ?u1D4602u1D45222 ?u1D451u1D4603u1D45223] = 1?tr(u1D446?) +(1?cos(?))[u1D4601(u1D45221 +u1D45222 +u1D45233) +u1D4602 +u1D451u1D4603 ?u1D4601u1D45221 ?u1D4602u1D45222 ?u1D451u1D4603u1D45223] = 1?tr(u1D446?) +(1+cos(?))[u1D4602 +u1D451u1D4603 +(u1D4601 ?u1D4602)u1D45222 +(u1D4601 ?u1D451u1D4603)u1D45223]. 45 u1D43F(u1D434) is minimized when ? = 0, which gives u1D44A = u1D43C so u1D43F(u1D434u1D45Cu1D45Du1D461) = 1?tr(u1D446?) = 1?u1D4601 ?u1D4602 ?u1D451u1D4603 and u1D434u1D45Cu1D45Du1D461 = u1D448+u1D449u1D447+ = u1D448[diag(1,1,u1D451)]u1D449u1D447. (2.74) The minimum is unique unless u1D4602 +u1D451u1D4603 = 0, in which case there is a family of minimizing u1D44A matrices given by setting e2 = e3 = 0 u1D44A = ? ?? ?? ? 1 0 0 0 cos(?) sin(?) 0 ?sin(?) cos(?) ? ?? ?? ? . Equation (2.74) is a transformation from the inertial frame to the body frame via two transfor- mations. The matrix u1D449u1D447+ transforms from the inertial frame to an intermediate from (S-frame), and u1D448+ transforms from the S-frame to the body frame. The rank of u1D435 determines the solution?s uniqueness. If u1D435 has rank less than two, the solution is not unique. The sensitivity in the attitude as a function of the variations u1D6FFru1D456 and u1D6FFbu1D456 is given by u1D467 = u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456[(u1D448u1D447+u1D6FFbu1D456)?(u1D449u1D447+ ru1D456)+(u1D448u1D447+bu1D456)?(u1D449u1D447+u1D6FFru1D456)]. In summary the algorithm is as follows: 1. Compute u1D435 from Equation (2.70). 2. Find the singular value decomposition of u1D435. 3. Compute u1D451 from Equation (2.72). 4. Compute u1D434u1D45Cu1D45Du1D461 from Equation (2.74). 5. Compute u1D43F(u1D434u1D45Cu1D45Du1D461) and any desired statistics. 46 2.10 Fast Optimal Attitude Matrix (FOAM) Algorithm In [11] Markley develops the FOAM algorithm to solve Wahba?s problem. Markley presents Wahba?s loss function as detailed in ?2.9. The cost function can be rewritten as u1D43F(u1D434) = u1D7060 ?tr(u1D434u1D435u1D447) u1D7060 = u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456 u1D435 = u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456bu1D456ru1D447u1D456 . The orthogonal matrix u1D434 that maximizes tr(u1D434u1D435u1D447) minimizes the expression ?u1D434?u1D435?2 = tr[(u1D434?u1D435)(u1D434?u1D435)u1D447] = tr(u1D43C)?2tr(u1D434u1D435u1D447) +?u1D435?2 soWahba?s problemis equivalent to ?ndingtheorthogonal matrixu1D434that isclosest tou1D435in Euclidean norm. The matrix u1D435 has the decomposition u1D435 = u1D448+ diag[u1D4461,u1D4462,u1D4463]u1D449u1D447+ where u1D448+ and u1D449+ are orthogonal matrices required to have positive det u1D4461 ?u1D4462 ? ?u1D4463?. The optimal attitude estimate is then u1D434u1D45Cu1D45Du1D461 = u1D448+u1D449u1D447+. 47 Matrix u1D435 has the properties ?u1D435?2u1D439 = u1D44621 +u1D44622 +u1D44623 det(u1D435) = det(u1D448+ diag[u1D4461,u1D4462,u1D4463]u1D449u1D447+ ) = det(u1D448+)det(diag[u1D4461,u1D4462,u1D4463])det(u1D449u1D447+ ) = (1)(u1D4461u1D4462u1D4463)(1) = u1D4461u1D4462u1D4463 adj(u1D435u1D447) = u1D448+ diag[u1D4462u1D4463,u1D4463u1D4461,u1D4461u1D4462]u1D449u1D447+ u1D435u1D435u1D447u1D435 = u1D448+ diag[u1D44631,u1D44632,u1D44633]u1D449u1D447+ all of which can be evaluated without computing the SVD. They are used to compute u1D434u1D45Cu1D45Du1D461 u1D434u1D45Cu1D45Du1D461 = [(u1D705+?u1D435?2)u1D435+u1D706adj(u1D435u1D447)?u1D435u1D435u1D447u1D435]/u1D701 (2.75) u1D705 = u1D4462u1D4463 +u1D4463u1D4461 +u1D4461u1D4462 (2.76) u1D706 = u1D4461 +u1D4462 +u1D4463 (2.77) u1D701 = (u1D4462 +u1D4463)(u1D4463 +u1D4461)(u1D4461 +u1D4462). (2.78) The SVD is not necessary to calculate the coe?cients u1D705, u1D706, and u1D701. The coe?cients u1D705 and u1D701 can be written in terms of u1D706 u1D705 = 12(u1D7062 ??u1D435?2) (2.79) u1D701 = u1D705u1D706?det(u1D435) (2.80) so u1D434(u1D706) = u1D434u1D45Cu1D45Du1D461 when u1D706= u1D4461 +u1D4462 +u1D4463. u1D706 can be found with the quartic polynomial 0 = u1D45D(u1D706) = (u1D7062 ??u1D435?2)2 ?8u1D706det(u1D435)?4?adj(u1D435)?2. 48 The roots are all real, and the maximum root is required (it is unique unless u1D4462 +u1D4463 = 0 in which case the attitude solution is not unique). u1D706 can be found by applying Newton?s method. The covariance matrix u1D443 for u1D434u1D45Cu1D45Du1D461 can be given as u1D443 = u1D7060u1D70E2u1D461u1D45Cu1D461(u1D705u1D43C +u1D435u1D435u1D447)/u1D701 (2.81) where u1D70E2u1D461u1D45Cu1D461 = uni0028.alt04 u1D45Buni2211.alt02 u1D456=1 u1D70E?2u1D456 uni0029.alt04?1 . (2.82) The u1D456u1D461? measurement error vector components are assumed to have a Gaussian distribution with respect to the actual vector, and phase is assumed to have a uniform distribution with variance u1D70E2u1D456 per axis. The two choices for u1D7060 are 1 and u1D70E?2u1D461u1D45Cu1D461. The normalized (u1D7060 = 1) form is useful with ?xed-point arithmetic or if measurement weights are arbitrarily assigned. The unnormalized form (u1D7060 = u1D70E?2u1D461u1D45Cu1D461) is good if weights are computed with measurement variances. 2.11 Alternative Quaternion Attitude Estimation Algorithm In [12] Markley presents another method for attitude estimation based on Wahba?s problem. This method also uses SVD the formulation of Wahba?s problem presented in?2.9. The loss function can be rewritten as u1D43F(u1D434) = u1D7060 ?tr(u1D434u1D435u1D447) (2.83) u1D7060 = u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456 (2.84) u1D435 = u1D45Buni2211.alt02 u1D456=1 u1D44Eu1D456bu1D456ru1D447u1D456 . (2.85) The assumed orthogonality of A gives ?u1D434?u1D435?2 = tr[(u1D434?u1D435)(u1D434?u1D435)u1D447] = 3?2tr(u1D434u1D435u1D447)+?u1D435?2. This norm is minimized by the same matrix that maximizes tr(u1D434u1D435u1D447). 49 Let u1D701(u1D706,u1D435) = 12u1D706(u1D7062 ??u1D435?2)?det(u1D435). Then the optimal matrix u1D434u1D45Cu1D45Du1D461 is given by u1D434(u1D706) = uni005B.alt031 2(u1D706 2 +?u1D435?2)u1D435+u1D706adj(u1D435u1D447)?u1D435u1D435u1D447u1D435 uni005D.alt03 /u1D701(u1D706,u1D435) (2.86) where u1D706 is the largest root of the quartic equation resulting from u1D706 = tr[u1D434(u1D706)u1D435u1D447]. Several methods exist to solve for u1D434u1D45Cu1D45Du1D461 that involve the computation of u1D706. u1D706 should be close to the value of u1D7060 from Equation (2.84) since u1D43F(u1D434u1D45Cu1D45Du1D461) = u1D7060 ?u1D706? 0, and with small measurement errors the loss function should be close to zero. The resulting attitude estimate is u1D4340 = u1D440/u1D701(u1D7060,u1D435) (2.87) u1D440 = 12(u1D70620 +?u1D435?2)u1D435+u1D7060 adj(u1D435u1D447)?u1D435u1D435u1D447u1D435. (2.88) The estimate u1D4340 is only approximately orthogonal. A variant of Shepperds? algorithm can be used to obtain the normalized attitude quaternion u1D45E from u1D4340 to construct an orthogonal attitude matrix u1D434u1D452u1D460u1D461 u1D434u1D452u1D460u1D461 = ? ?? ?? ? u1D45E24 +u1D45E21 ?u1D45E22 ?u1D45E23 2(u1D45E1u1D45E2 +u1D45E3u1D45E4) 2(u1D45E1u1D45E3 ?u1D45E2u1D45E4) 2(u1D45E2u1D45E1 ?u1D45E3u1D45E4) u1D45E24 ?u1D45E21 +u1D45E22 ?u1D45E23 2(u1D45E2u1D45E3 +u1D45E1u1D45E4) 2(u1D45E3u1D45E1 +u1D45E2u1D45E4) 2(u1D45E3u1D45E2 ?u1D45E1u1D45E4) u1D45E24 ?u1D45E21 ?u1D45E22 +u1D45E23 ? ?? ?? ? . (2.89) 50 If u1D456, u1D457, u1D458 is a cyclic permutation of 1, 2, 3, the quaternion components obey the relations 4u1D701(u1D7060,u1D435)u1D45E2u1D456 ? u1D701(u1D7060,u1D435) +u1D440u1D456u1D456 ?u1D440u1D457u1D457 ?u1D440u1D458u1D458 = u1D463u1D456 4u1D701(u1D7060,u1D435)u1D45Eu1D456u1D45Eu1D457 ? u1D440u1D456u1D457 +u1D440u1D457u1D456 = u1D463u1D457 4u1D701(u1D7060,u1D435)u1D45Eu1D456u1D45Eu1D458 ? u1D440u1D456u1D458 +u1D440u1D458u1D456 = u1D463u1D458 4u1D701(u1D7060,u1D435)u1D45Eu1D456u1D45E4 ? u1D440u1D457u1D458 ?u1D440u1D458u1D457 = u1D4634 = u1D464u1D456 4u1D701(u1D7060,u1D435)u1D45Eu1D457u1D45E4 ? u1D440u1D458u1D456 ?u1D440u1D456u1D458 = u1D464u1D457 4u1D701(u1D7060,u1D435)u1D45Eu1D458u1D45E4 ? u1D440u1D456u1D457 ?u1D440u1D457u1D456 = u1D464u1D458 4u1D701(u1D7060,u1D435)u1D45E24 ? u1D701(u1D7060,u1D435) +u1D440u1D456u1D456 +u1D440u1D457u1D457 +u1D440u1D458u1D458 = u1D4644. Let u1D456 correspond to the index of the largest diagonal element of u1D440, and de?ne the quaternion components for u1D459 = 1,...,4 by u1D45Eu1D459 = u1D463u1D459/?u1D463? for u1D440u1D457u1D457 +u1D440u1D458u1D458 < 0 or u1D45Eu1D459 = u1D464u1D459/?u1D464? for u1D440u1D457u1D457 +u1D440u1D458u1D458 ? 0 where ?u1D463? and ?u1D464? are the Euclidean norms of the u1D463 and u1D464. 2.12 Euler-q Algorithm In [13] Mortari develops the Euler-q algorithm. It is desired to ?nd the attitude matrix u1D434 that rotates unit vectors vu1D456 from the inertial frame to corresponding unit vectors su1D456 in the body frame. u1D434 can be expressed in terms of the Euler axis e, a unit 3-vector about which a rotation takes place, and the the Euler angle u1D719, the angle of rotation about the Euler axis, as follows u1D434 = u1D43Ccos(u1D719) +(1?cos(u1D719))eeu1D447 ? ?u1D452sin(u1D719) 51 where ?u1D452 is the 3?3 skew-symmetric, cross-product matrix ?u1D452 = ? ?? ?? ? 0 ?u1D4523 u1D4522 u1D4523 0 ?u1D4521 ?u1D4522 u1D4521 0 ? ?? ?? ? . Let u1D6FDu1D456 represent the accuracy of the u1D456u1D461? sensor. Then the sensor relative precision can be de?ned as u1D6FCu1D456 = 1/ uni0028.alt04 u1D6FDu1D456 u1D45Buni2211.alt02 u1D458=1 (1/u1D6FDu1D458) uni0029.alt04 . Wahba?s problem is to ?nd the 3?3 matrix u1D434 that minimizes the loss function u1D43Fu1D44A(u1D434) = 12 u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456?su1D456 ?u1D434vu1D456?2 = 1? u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456su1D447u1D456 u1D434vu1D456 or equivalently maximizes the gain function u1D43Au1D44A(u1D434) = 1?u1D43Fu1D44A(u1D434) = u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456su1D447u1D456 u1D434vu1D456 = tr(u1D434u1D435u1D447) u1D435 = u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456su1D456vu1D447u1D456 . The gain function can be rewritten as u1D43Au1D44A(e,u1D719) = cos(u1D719)tr(u1D435) +(1?cos(u1D719))eu1D447u1D435e+sin(u1D719)zu1D447e z = u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456su1D456 ?vu1D456. 52 De?ne a loss function u1D43Fu1D440(e) = u1D45Buni2211.alt02 u1D456=1 u1D709u1D456u1D6FF2u1D456 = u1D45Buni2211.alt02 u1D456=1 u1D709u1D456eu1D447du1D456du1D447u1D456 e = eu1D447u1D43Be du1D456 = vu1D456 ?su1D456/?vu1D456 ?su1D456? u1D6FFu1D456 = eu1D447du1D456 = du1D447u1D456 e u1D43B = u1D43Bu1D447 = u1D45Buni2211.alt02 u1D456=1 u1D709u1D456du1D456du1D447u1D456 where u1D709u1D456 are relative weights. The worst case for the du1D456 direction deviation vector occurs at the angle u1D6FD?u1D456 when Su1D456 is displaced from the measured su1D456 by the angle u1D6FDu1D456 and the spherical triangle is right. The value of u1D6FD?u1D456 is obtained from sin(u1D714u1D456)sin(u1D6FD?u1D456) = sin(u1D6FDu1D456) where u1D714u1D456 is the angle between the su1D456 and vu1D456 directions. The relative weights u1D709u1D456 are then derived from u1D6FD?u1D456 in the same way that u1D6FCu1D456 are obtained from u1D6FDu1D456 u1D709u1D456 = 1/ uni0028.alt04 u1D6FD?u1D456 u1D45Buni2211.alt02 u1D458=1 (1/u1D6FD?u1D458) uni0029.alt04 . The augmented cost function is de?ned as u1D43F?u1D440(e) = eu1D447u1D43Be?u1D706(eu1D447e?1) which leads to u1D43Be = u1D706e. Thus, the Euler axis is the eigenvector of the u1D43B matrix associated with the eigenvalue u1D706 so u1D43Fu1D440(e) = u1D706= eu1D447u1D43Be. 53 Since u1D43Fu1D440(e) has to be a minimum, the unknown eigenvalue has to be the smallest u1D706 u1D43Beu1D45Cu1D45Du1D461 = u1D706u1D45Au1D456u1D45Beu1D45Cu1D45Du1D461. u1D706u1D45Au1D456u1D45B can be found using the characteristic equation of u1D43B 0 = u1D7063 +u1D44Eu1D7062 +u1D44Fu1D706+u1D450 u1D44E= ?tr(u1D43B) u1D44F= tr(adj(u1D43B)) u1D450 = ?det(u1D43B). Since u1D43B is symmetric positive semide?nite, its eigenvalues are real and nonnegative. Let 0 ?u1D7061 ? u1D7062 ?u1D7063, and de?ne the variables u1D45D2 = (u1D44E/3)2 ?(u1D44F/3) u1D45E = [(u1D44F/2)?(u1D44E/3)2](u1D44E/3)?(u1D450/2) u1D464 = 13 cos?1(u1D45E/u1D45D3) which leads to u1D7061 = ?u1D45D(?3sin(u1D464) +cos(u1D464))?u1D44E/3 u1D7062 = u1D45D(?3sin(u1D464)?cos(u1D464))?u1D44E/3 u1D7063 = 2u1D45Dcos(u1D464)?u1D44E/3. The calculation of eu1D45Cu1D45Du1D461 proceeds as follows (u1D43B ?u1D706u1D45Au1D456u1D45Bu1D43C)eu1D45Cu1D45Du1D461 = u1D440eu1D45Cu1D45Du1D461 = ? ?? ?? ? u1D45Au1D4471 u1D45Au1D4472 u1D45Au1D4473 ? ?? ?? ? eu1D45Cu1D45Du1D461 = ? ?? ?? ? u1D45Au1D44E u1D45Au1D465 u1D45Au1D466 u1D45Au1D465 u1D45Au1D44F u1D45Au1D467 u1D45Au1D466 u1D45Au1D467 u1D45Au1D450 ? ?? ?? ? eu1D45Cu1D45Du1D461 = 0. 54 The direction of the optimal Euler axis can be found from the cross product of any two row vectors of u1D440 e1 = m2 ?m3 e2 = m3 ?m1 e3 = m1 ?m2. The most accurate eu1D456 is the one with the greatest modulus, which is determined by ?nding the greatest u1D45Du1D456 u1D45D1 = ?u1D45Au1D44Fu1D45Au1D450?u1D45A2u1D467? u1D45D2 = ?u1D45Au1D44Eu1D45Au1D450 ?u1D45A2u1D466? u1D45D3 = ?u1D45Au1D44Eu1D45Au1D44F ?u1D45A2u1D465?. Let u1D45Du1D458 = max(u1D45D1,u1D45D2,u1D45D3). Then Euler axis with the greatest modulus is eu1D45Cu1D45Du1D461 = eu1D458/?eu1D458?. The optimal Euler angle can be derived from sin(u1D719u1D45Cu1D45Du1D461) = (1/u1D437)zu1D447eu1D45Cu1D45Du1D461 cos(u1D719u1D45Cu1D45Du1D461) = (1/u1D437)(tr(u1D435)?eu1D447u1D45Cu1D45Du1D461u1D435eu1D45Cu1D45Du1D461) u1D4372 = (zu1D447eu1D45Cu1D45Du1D461)2 +(tr(u1D435)?eu1D447u1D45Cu1D45Du1D461u1D435eu1D45Cu1D45Du1D461)2. The matrix u1D440 is singular or ill-conditioned under one of three conditions: a) if the row vectors ofu1D440 are parallel (colinear), b) the Euler axis e and all observed vectors are approximately coplanar (one vector is nearly linearly dependent on the other two), or c) the row vectors of u1D440 become 0. For case (a) there is no solution, but for cases (b) and (c) the method of sequential rotations (MSR) can be applied. MSR states that if an optimal attitude matrix exists for the u1D45B unit vector pairs 55 (su1D456,vu1D456), then there exists u1D45B unit vector pairs (su1D456,wu1D456) that imply an optimal attitude matrix u1D439. The vectors su1D456 and wu1D456 are related by any rotation matrix u1D445: wu1D456 = u1D445su1D456. u1D434 is related to u1D439 by u1D439 = uni005B.alt03 f1 f2 f3 uni005D.alt03 = u1D434u1D445u1D447. If the (su1D456,vu1D456) data set implies a singularity, the set (su1D456,wu1D456) would not necessarily imply a singularity. So MSR evaluates the attitude of u1D439 by using the rotated unit vector wu1D456 and then computes u1D434 as u1D434 = u1D439u1D445. By using one of the following rotation matrices u1D4451 = ? ?? ?? ? 1 0 0 0 ?1 0 0 0 ?1 ? ?? ?? ? u1D4452 = ? ?? ?? ? ?1 0 0 0 1 0 0 0 ?1 ? ?? ?? ? u1D4453 = ? ?? ?? ? ?1 0 0 0 ?1 0 0 0 1 ? ?? ?? ? , no new computation is required, only sign changes. 2.13 ESOQ Algorithm In [14] Mortari presents the commonly used EStimator of the Optimal Quaternion (ESOQ) algorithm. The ESOQ algorithm is a singularity-free algorithm that ?nds a closed form expression for u1D706u1D45Au1D44Eu1D465 and u1D45Eu1D45Cu1D45Du1D461 of the q-method equation (Equation (2.94) below). Let su1D456 be unit vectors measured in the body frame, vu1D456 be the actual vectors in the inertial frame, and u1D434 be the matrix that rotates vu1D456 to the body frame. Let u1D6FDu1D456 be the precision of the u1D456u1D461? sensor such that the angle between the true and observed direction is smaller than u1D6FDu1D456. The sensor relative precision is de?ned as u1D6FCu1D456 = 1/ uni0028.alt04 u1D6FDu1D456 u1D45Buni2211.alt02 u1D458=1 (1/u1D6FDu1D458) uni0029.alt04 . (2.90) 56 Wahba?s problem can be solved by computation of the 3?3 rotation matrix u1D434 that minimizes u1D43Fu1D44A(u1D434) = 12 u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456?su1D456 ?u1D434vu1D456?2 = 1? u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456su1D447u1D456 u1D434vu1D456 (2.91) or equivalently maximizes the gain function u1D43Au1D44A(u1D434) = 1?u1D43Fu1D44A(u1D434) = u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456su1D447u1D456 u1D434vu1D456 = tr(u1D434u1D435u1D447) (2.92) u1D435 = u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456su1D456vu1D447u1D456 . (2.93) The q-method solution equation is u1D43Eu1D45Eu1D45Cu1D45Du1D461 = u1D706u1D45Au1D44Eu1D465u1D45Eu1D45Cu1D45Du1D461 (2.94) meaning u1D45Eu1D45Cu1D45Du1D461 is the eigenvector of the symmetric matrix u1D43E u1D43E = ? ?? u1D446 z zu1D447 u1D461 ? ??= ? ?? u1D435+u1D435u1D447 ?u1D43Ctr(u1D435) z zu1D447 tr(u1D435) ? ?? (2.95) z = u1D45Buni2211.alt02 u1D456=1 u1D6FCu1D456su1D456 ?vu1D456. (2.96) associated with its maximum eigenvalue assuming that the vectors su1D456 and vu1D456 are normalized. The characteristic equation of the matrix u1D43E can be written as 0 = u1D7064 +u1D44Eu1D7063 +u1D44Fu1D7062 +u1D450u1D706+u1D451 (2.97) u1D44E = tr(u1D43E) = 0 (2.98) u1D44F = ?2(tr(u1D435))2 +tr(adj(u1D435+u1D435u1D447))?zu1D447z (2.99) u1D450 = ?tr(adj(u1D43E)) (2.100) u1D451 = det(u1D43E). (2.101) 57 The third order auxiliary equation of the characteristic equation is u1D4623 ?u1D44Fu1D4622 ?4u1D451u1D462+4u1D44Fu1D451?u1D4502 = 0 (2.102) which has the solution u1D4621 = 2?u1D45Dcos uni0028.alt031 3 cos ?1(u1D45E/u1D45D3/2) uni0029.alt03 +u1D44F/3 (2.103) u1D45D = (u1D44F/3)2 +4u1D451/3 (2.104) u1D45E = (u1D44F/3)3 ?4u1D451u1D44F/3 +u1D4502/2. (2.105) Then u1D7061 = 12 uni0028.alt02 ?u1D4541 ? uni221A.alt01 ?u1D4621 ?u1D44F+u1D4542 uni0029.alt02 (2.106) u1D7062 = 12 uni0028.alt02 ?u1D4541 + uni221A.alt01 ?u1D4621 ?u1D44F+u1D4542 uni0029.alt02 (2.107) u1D7063 = 12 uni0028.alt02 u1D4541 ? uni221A.alt01 ?u1D4621 ?u1D44F?u1D4542 uni0029.alt02 (2.108) u1D7064 = 12 uni0028.alt02 u1D4541 + uni221A.alt01 ?u1D4621 ?u1D44F?u1D4542 uni0029.alt02 (2.109) u1D4541 = uni221A.alt01 u1D4621 ?u1D44F (2.110) u1D4542 = 2 uni221A.alt02 u1D46221 ?4u1D451 (2.111) where ?1 ?u1D7061 ?u1D7062 ?u1D7063 ?u1D7064 = u1D706u1D45Au1D44Eu1D465 ? 1. (2.112) If u1D45B = 2 u1D7064 = ?u1D7061 = (u1D4543 +u1D4544)/2 (2.113) u1D7063 = ?u1D7062 = (u1D4543 ?u1D4544)/2 (2.114) u1D4543 = uni221A.alt02 2 ? u1D451?u1D44F (2.115) u1D4544 = uni221A.alt02 ?2 ? u1D451?u1D44F. (2.116) 58 The q-method solution equation is equivalent to (u1D43E?u1D706u1D45Au1D44Eu1D465u1D43C)u1D45Eu1D45Cu1D45Du1D461 = u1D43Bu1D45Eu1D45Cu1D45Du1D461 = 0 (2.117) which meansu1D45Eu1D45Cu1D45Du1D461 is perpendicular to each row vector hu1D456 of symmetricu1D43B. Four di?erent and parallel cross-products u1D45Eu1D458 can be evaluated using the 4-D cross-product u1D45Eu1D458(u1D456) = (?1)u1D458+u1D456det(u1D43Bu1D458u1D456) where u1D456 = 1,...,4 and the 3?3 submatrices u1D43Bu1D458u1D456 are obtained from matrix u1D43B by deleting the u1D458u1D461? row and the u1D456u1D461? column. The optimal u1D45E is then u1D45Eu1D45Cu1D45Du1D461 = u1D45Eu1D43A/?u1D45Eu1D43A? where u1D43A is the index associated with the element u1D45Eu1D458(u1D458) with the largest magnitude. 2.14 ESOQ2 Algorithm In[15] Mortari presents theESOQ2algorithm, an Euler angle variation of the ESOQalgorithm. The same procedure as outlined in ?2.13 in Equations (2.90)-(2.117) is used to develop the ESOQ2 algorithm. The quaternion solution to the q-method equation (Equation (2.94)) can be expressed in terms of an Euler axis and angle u1D45E = {eu1D447 sin(u1D719/2) cos(u1D719/2)}u1D447. It can be shown that the equation for determining the optimal Euler axis is [(u1D461?u1D706u1D45Au1D44Eu1D465)(u1D446?u1D706u1D45Au1D44Eu1D465u1D43C)?zzu1D447]e = u1D440e = 0 59 where u1D440 = ? ?? ?? ? mu1D4471 mu1D4472 mu1D4473 ? ?? ?? ? = ? ?? ?? ? u1D45Au1D44E u1D45Au1D465 u1D45Au1D466 u1D45Au1D465 u1D45Au1D44F u1D45Au1D467 u1D45Au1D466 u1D45Au1D467 u1D45Au1D450 ? ?? ?? ? . All row vectors of u1D440 are perpendicular to e. The optimal Euler axis e can be found by taking the cross-product of two rows of u1D440. This fails in two cases: 1. When the rows of u1D440 become zero which occurs if u1D719? 0 (the Euler angle singularity). 2. When the rows of u1D440 become parallel (an unresolvable case). The three choices for e are e1 = m2 ?m3 e2 = m3 ?m1 e3 = m1 ?m2 where the eu1D456?s di?er only in modulus. The most accurate eu1D456 is the one with the greatest modulus, which is determined by ?nding the greatest u1D45Du1D456 u1D45D1 = ?u1D45Au1D44Fu1D45Au1D450?u1D45A2u1D467? u1D45D2 = ?u1D45Au1D44Eu1D45Au1D450 ?u1D45A2u1D466? u1D45D3 = ?u1D45Au1D44Eu1D45Au1D44F ?u1D45A2u1D465?. Let u1D45Du1D458 have the largest value. Then the most reliable Euler axis is eu1D45Cu1D45Du1D461 = eu1D458/?eu1D458?. 60 The optimal Euler angle is then derived from u1D44B(u1D458) = ?sin(u1D719/2) (2.118) u1D44C(u1D458) = ?cos(u1D719/2) (2.119) where ? is a proportional constant and u1D458 identi?es the element of vector u1D44B with the greatest modulus. The optimal quaternion can then be computed from e and u1D719. To avoid the singularity that occurs when u1D719? 0, the MSR can be applied. 2.15 A Slightly Sub-Optimal Algorithm In [16] an algorithm for the solution of Wahba?s problem that requires no iterations and little computation but results in a slightly nonorthogonal matrix solution is presented. Wahba?s loss function can be written as u1D43F(u1D434) = 12 u1D45Buni2211.alt02 u1D457=1 ?vu1D457 ?u1D434uu1D457?2 where uu1D457 are unit vectors in the inertial frame of u1D45B observations and vu1D457 are the corresponding unit vectors in the body frame. The loss function can be generalized by using an u1D45B ?u1D45B symmetric positive-de?nite weight matrix u1D44A u1D43F(u1D434) = 12 tr(u1D44A(u1D434u1D448 ?u1D449)u1D447(u1D434u1D448 ?u1D449)) u1D448 = [u1,u2,...,uu1D45B] u1D449 = [v1,u1D4632,...,vu1D45B]. The loss function can be written as u1D43F(u1D434) = 12 tr(u1D44A(u1D448u1D447u1D448 +u1D449u1D447u1D449))?tr(u1D434u1D435u1D447) u1D435 = u1D449u1D44Au1D448u1D447. 61 The loss function is minimized when tr(u1D434u1D435u1D447) is maximized. This occurs when u1D434 is the orthogonal matrix that is closest to u1D435 in the Euclidean (or Frobenius) norm. The loss function can be further generalized to u1D43F(u1D434) = 12 tr(u1D44A(u1D434u1D448 ?u1D449)u1D447u1D44D(u1D434u1D448 ?u1D449)) where u1D44D is a 3?3 symmetric positive-de?nite weight matrix. Let u1D4340 be an extremum of the loss function and u1D716 an arbitrary variation in the direction of an arbitrary non-zero matrix u1D43B. Then u1D43F(u1D4340 +u1D716u1D43B) = 12?u1D44D1/2(u1D4340u1D448 ?u1D449)u1D44A1/2?2u1D439 +u1D716tr(u1D43Bu1D447u1D44D(u1D4340u1D448 ?u1D449)u1D44Au1D448u1D447] (2.120) + 12u1D7162?u1D44D1/2u1D43Bu1D448u1D44A1/2?2u1D439. (2.121) In order for Equation (2.121) have a solution, the following condition must hold (u1D4340u1D448 ?u1D449)u1D44Au1D448u1D447 = 0. If u1D448 is full rank u1D4340 = u1D449u1D44Au1D448u1D447(u1D448u1D44Au1D448u1D447)?1 = u1D435(u1D448u1D44Au1D448u1D447)?1. The solution provides an unbiased estimate of u1D434 under some assumptions: 1. Error-free reference vectors 2. Additive random measurement errors vu1D457 = u1D434u1D461u1D45Fu1D462u1D452uu1D457 +nu1D457 ?u1D457 u1D449 = u1D434u1D461u1D45Fu1D462u1D452u1D448 +u1D441 3. E[N] = 0. 62 The attitude estimate is given by u1D4340 = (u1D434u1D461u1D45Fu1D462u1D452u1D448 +u1D441)u1D44Au1D448u1D447(u1D448u1D44Au1D448u1D447)?1 = u1D434u1D461u1D45Fu1D462u1D452 +u1D6FFu1D434 u1D6FFu1D434 = u1D441u1D44Au1D448u1D447(u1D448u1D44Au1D448u1D447)?1. An estimate of the deviation of u1D4340 from u1D434u1D461u1D45Fu1D462u1D452 is u1D443 = u1D438[(u1D6FFu1D434)u1D447u1D6FFu1D434] = (u1D448u1D44Au1D448u1D447)?1u1D448u1D44Au1D445u1D44Au1D448u1D447(u1D448u1D44Au1D448u1D447)?1 u1D445 = u1D438[u1D441u1D447u1D441]. If we choose u1D44A = u1D445?1, then u1D443 = (u1D448u1D44Au1D448u1D447)?1 = (u1D448u1D445?1u1D448u1D447)?1. If {nu1D457} is a white sequence, then u1D445 and u1D44A are diagonal matrices. This method requires a minimum of three observations (u1D45B = 3). In the unconstrained problem, a third observation can be added (the cross-product v1?v2 as a measurement of u1?uu1D460) to make u1D448 rank 3. 2.16 Conclusions Many algorithms have been developed to ?nd solutions to Wahba?s problem. While many mathematical solutions are o?ered to solve the problem, Davenport is the ?rst to present a solution in the form of a practical algorithm. Davenport?s algorithm develops what has become known as Davenport?s equation but does not provide a quick way to solve the equation. However, Davenport?s equation is the foundation of many other algorithms. The QUEST algorithm is one such algorithm. Faster than Davenport?s method, the QUEST algorithm has the major disadvantage of using only single measurements. Past information is not taken into account so anomalous data can cause 63 bad estimates. To overcome these limitations, variations of the QUEST algorithm have been developed to make the algorithm more like Kalman ?ltering. These algorithms cannot compute sensor biases and accounting for process noise proves to be too computationally expensive. The REQUEST algorithm makes the QUEST algorithm recursive, but REQUEST only works if high quality sensors are used to get attitude measurements and does not estimate biases. The extended QUEST algorithm is developed to estimate biases but is slow. The energy approach algorithm o?ers a new formulation of Wahba?s problem that is analogous to a physical system; therefore, it is easy to understand. The algorithm has the drawbacks of assuming small errors in measurements and of often encountering singularities when ?nding a solution. The SVD algorithm, which requires taking a computationally expensive SVD, provides the basis for the FOAM algorithm, which does not require computing an SVD. The FOAM algorithm and most algorithms developed after the FOAM algorithm provide the same solution, but the speeds at which the solution is derived di?er. The ESOQ algorithm provides the same optimal solution as the FOAM algorithm but is faster and guaranteed singularity free. The ESOQ2 algorithm di?ers from the ESOQ algorithm in that it provides a solution in terms of Euler angles rather than a quaternion. The sub-optimal algorithm is computationally inexpensive and the fastest of all the algorithms presented; however, it has the drawback of computing a rotation matrix that is nonorthogonal. The ESOQ algorithm is best suited for estimating the states of a ballistic rocket in the early stages of ?ight. This algorithm is fast while providing an orthogonal rotation matrix estimate. As much accuracy in an attitude estimate as possible as quickly as possible is necessary for the system presented in this work. The ESOQ algorithm also accounts for past measurements rather than being a single-point algorithm like other Wahba?s problem algorithms. No assumptions are made about the quality of sensors used, which is essential since MEMS sensors are not as accurate as more conventional sensors. 64 Chapter 3 Magnetometer Navigation Magnetometers are a key component of solving Wahba?s problem as discussed in Chapter 2. In this chapter we present information relevant to the use of a magnetometer for state estimation of a vehicle. In ?3.1 we develop equations to model the earth?s magnetic ?eld and present a model that implements those equations. Next, in ?3.2 we present a model for a magnetometer, a device that measures magnetic ?elds. We also present several possible error sources that corrupt magnetometer measurements. In ?3.3 we present methods to calibrate a magnetometer in order to minimize measurement errors. In ?3.4 we present an algorithm that estimates angular rates based on magnetometer measurements. Finally, in ?3.5 we evaluate the presented algorithms as they pertain to the problem in this work. 3.1 Modeling Earth?s Magnetic Field To simulate a magnetometer?s behavior, a model of the earth?s magnetic ?eld is necessary. In [17] Roithmayr develops equations that are necessary to construct such a model. 3.1.1 Equation Development An in?nite series of spherical harmonics can be used to describe magnetic ?elds in general. Then at a point u1D444 above the earth?s surface, the magnetic ?eld vector B is given by Bu1D45B,u1D45A = u1D43Eu1D45B,u1D45Au1D44E u1D45B+2 u1D445u1D45B+u1D45A+1 uni007B.alt02u1D454u1D45B,u1D45Au1D436u1D45A +?u1D45B,u1D45Au1D446u1D45A u1D445 [(u1D460u1D706u1D434u1D45B,u1D45A+1 +(u1D45B+u1D45A+1)u1D434u1D45B,u1D45A)?r?u1D434u1D45B,u1D45A+1?e3] ? u1D45Au1D434u1D45B,u1D45A[(u1D454u1D45B,u1D45Au1D436u1D45A?1 +?u1D45B,u1D45Au1D446u1D45A?1)?e1 +(?u1D45B,u1D45Au1D436u1D45A?1 ?u1D454u1D45B,u1D45Au1D446u1D45A?1)?e2] uni007D.alt02 (3.1) where B = ?uni2211.alt02 u1D45B=1 u1D45Buni2211.alt02 u1D45A=0 Bu1D45B,u1D45A. (3.2) 65 In Equation (3.1) u1D44E is the earth?s average radius (6371 km), u1D445 is the magnitude of the position vector R from the earth?s center to point u1D444, ?r is a unit vector in the direction of R, and u1D454u1D45B,u1D45A and ?u1D45B,u1D45A are degree u1D45B, order u1D45A Gauss coe?cients. ?e1, ?e2, and ?e3 are unit vectors of an Earth-?xed coordinate system. The coe?cients u1D43Eu1D45B,u1D45A are determined recursively from either u1D43Eu1D45B,u1D45A = uni0028.alt03u1D45B?u1D45A u1D45B+u1D45A uni0029.alt031/2 u1D43Eu1D45B?1,u1D45A u1D45A = 1,...,?, u1D45B?u1D45A+1 or u1D43Eu1D45B,u1D45A = [(u1D45B+u1D45A)(u1D45B?u1D45A+1)]?1/2u1D43Eu1D45B,u1D45A?1 u1D45A = 2,...,?, u1D45B?u1D45A whereu1D43Eu1D45B,0 = 1 wheneveru1D45A = 0 andu1D43E1,1 = 1. u1D434u1D45B,u1D45A andu1D434u1D45B,u1D45A+1 are degreeu1D45BLegendre polynomials (see ?B.2) with orders u1D45A and u1D45A+1, respectively. When u1D45B = u1D45A, the polynomial is given by u1D434u1D45B,u1D45B = (1)(3)(5)...(2u1D45B?1) u1D45B = 1,...,? = (2u1D45B?1)u1D434u1D45B?1,u1D45B?1 u1D45B = 2,...,? and u1D4340,0 = 1. In general, u1D434u1D45B,u1D45A(u1D462) = 1u1D45B?u1D45A[(2u1D45B?1)u1D462u1D434u1D45B?1,u1D45A(u1D462)?(u1D45B+u1D45A?1)u1D434u1D45B?2,u1D45A(u1D462)] u1D45A = 0,...,?, u1D45B? (u1D45A+1). Arguments of the Legendre polynomials are u1D460u1D706 = sin(u1D706) = ?r??e3 where u1D706 is the geographic latitude. The variables u1D446u1D45A and u1D436u1D45A are de?ned recursively as u1D446u1D45A = u1D4461u1D436u1D45A?1 +u1D4361u1D446u1D45A?1 u1D436u1D45A = u1D4361u1D436u1D45A?1 ?u1D4461u1D446u1D45A?1 66 where u1D4460 = 0 u1D4461 = R??e2 u1D4360 = 1 u1D4361 = R??e1. the earth?s magnetic ?eld can be modeled as a tilted dipole as long as the point u1D444 is not near the magnetic poles. For the dipole model u1D4341,0(u1D460u1D706) = u1D460u1D706 u1D4341,1(u1D460u1D706) = 1 u1D4341,2(u1D460u1D706) = 0 u1D43E1,1 = u1D43E1,0 = 1. Then B1,0 = uni0028.alt02u1D44E u1D445 uni0029.alt023 u1D4541,0(3u1D460u1D706?r??e3) B1,1 = uni0028.alt02u1D44E u1D445 uni0029.alt023 [3(u1D4541,1?r??e1 +?1,1?r??e2)?r?(u1D4541,1?e1 +?1,1?e2)]. So B1,0 +B1,1 = uni0028.alt02u1D44E u1D445 uni0029.alt023 [3(?u1D45F?M)?u1D45F?M] where the terrestrial dipole moment M is de?ned as M = u1D4541,1?e1 +?1,1?e2 +u1D4541,0?e3. 67 3.1.2 Model Implementation Two implementations of models of the earth?s magnetic ?eld are commonly used: the Inter- national Geomagnetic Reference Field (IGRF) [18] and the World Magnetic Model (WMM) (an eighth-order model) [19] which is used by the U.S. DoD and NATO. Due to the changes in the be- havior of the earth?s magnetic ?eld such as from secular variation (slow changes in time of the main magnetic ?eld), these models are updated every ?ve years (last updated in 2005). Our development is based on the WMM, and further discussion is concentrated solely on the WMM. the earth?s magnetic ?eld is generated mainly by three sources: 1. The main ?eld generated by the earth?s outer core (Bu1D45A). 2. The crustal ?eld generated by the earth?s crust and upper mantle (Bu1D450). 3. Electrical currents in the upper magnetosphere and atmosphere which induce electrical cur- rents in the ground and oceans (Bu1D451). The total magnetic ?eld can then be written as B(r,u1D461) = Bu1D45A(r,u1D461) +Bu1D450(r) +Bu1D451(r,u1D461) where r is a position vector and u1D461 is time. To model the earth?s magnetic ?eld the WMM uses data gathered from the Danish ?rsted and German CHAMP satellites, which have good global coverage and low noise, as well as data from ground observatories. The WMM only takes into account the contributions of Bu1D45A. This introduces errors into the model since Bu1D450 and Bu1D451 are ignored. The WMM also has other error sources. A magnetic sensor will not match the WMM in all locations on the earth; it may observe spatial and temporal anomalies. Spatial anomalies on land are caused by mountain ranges, ore deposits, geological faults, trains, railroad tracks, power lines, and other such conditions. Disturbances in the atmosphere from ionic activity from space will also cause variations in the magnetic ?eld. Seven quantities describe the geomagnetic ?eld vector B. These quantities are northerly inten- sity u1D44B, easterly intensity u1D44C, vertical intensity (positive downwards)u1D44D, total intensity u1D439, horizontal 68 intensity u1D43B, inclination (or dip) u1D43C which is the angle between the horizontal plane and the ?eld vector (measured positive downwards), and declination (or magnetic variation) u1D437 which is the hor- izontal angle between true north and the ?eld vector (measured positive eastwards). These values are calculated as follows u1D43B = uni221A.alt01 u1D44B2 +u1D44C2 u1D439 = uni221A.alt01 u1D43B2 +u1D44D2 u1D437 = arctan(u1D44C,u1D44B) u1D43C = arctan(u1D44D,u1D43B) u1D43Au1D449 = u1D437?u1D706 for u1D719> 55? u1D43Au1D449 = u1D437+u1D706 for u1D719 0 0 if u1D45E?u1D45E < 0 (4.32) where u1D43B(y) = ??u1D453/?y. The boundary conditions of the Hamiltonian and costate vector are u1D43B(u1D447u1D453) = 0 (4.33) q(u1D4470) = q0 = uni005B.alt03 free free free free free uni005D.alt03u1D447 (4.34) q(u1D447u1D453) = qu1D453 = uni005B.alt03 free free free free 0 uni005D.alt03u1D447 . (4.35) A control history can be generated with the following steps: 1. Initialize yu1D453 and qu1D453. 2. Integrate the ?y and ?q equations backward in time, and at each u1D447?u1D45B = u1D45B?u1D447 obtain y(u1D447?u1D45B) where ?u1D447 is the time interval chosen to give desired data point spacing between y(u1D447?u1D45B) and y(u1D447?u1D45B+1) and u1D45B is a positive integer. 3. Transform y(u1D447?u1D45B) to obtain the boundary conditions in the original form u1D4651,0, u1D4652,0, u1D4653,u1D451, and u1D4654,u1D451. Store the thruster switching times as functions of these four variables. 4.2.3 Multiple-Reaction Jet Control Costello et al. [58], [59], [60] are the only authors that investigate control of a projectile solely through multiple reaction jets. In [58] a rocket is simulated that has a main thruster that has a limited burn time and is stabilized by three pop-out ?ns. Lateral pulse jets, each of which can only be ?red once, are located toward the front of the rocket. Lateral pulse jet control is investigated to improve the dispersion pattern of direct-?re atmospheric rockets. This method of control falls under 96 the category of impulse control, unlike other common rocket control methods such as proportional navigation guidance (PNG), which is continuous. In [58] Jitpraphai et al. investigate PNG, parabolic and proportional navigation guidance (PAPNG), and trajectory tracking (TT) applied to a direct-?re atmospheric rocket equipped with a lateral pulse jet control mechanism. All three control laws use the same control logic to ?re the lateral pulse jets. Proportional Navigation Guidance (PNG) The guidance law for PNG is u1D434u1D436 = u1D441u1D449u1D450?u1D706 (4.36) where u1D441 is the proportional navigation constant (typically valued between 3 and 5), u1D449u1D436 is rocket closing velocity, and ?u1D706 is line-of-sight (LOS) angular rate. PNG has horizontal plane and vertical plane control law components and three reference frames are used: the LOS frame, the target frame, and the inertial frame. For the following equations, the target is assumed to be stationary. All axes in [58] are labeled (u1D43C,u1D43D,u1D43E) with a subscript to indicate the reference frame. The u1D43C-axis of the LOS frame points from the rocket to the target. The horizontal component of u1D706 is then u1D706u1D43B = tan?1 uni0028.alt03u1D466 u1D447 ?u1D466 u1D465u1D447 ?u1D465 uni0029.alt03 (4.37) and ?u1D706u1D43B = ??u1D466(u1D465u1D447 ?u1D465) + ?u1D465(u1D466u1D447 ?u1D466) (u1D465u1D447 ?u1D465)2 +(u1D466u1D447 ?u1D466)2 . (4.38) Then the horizontal acceleration command is ?u1D434u1D44Cu1D436 = u1D441u1D43B ?u1D706u1D43Bu1D462u1D43B (4.39) u1D462u1D43B = cos(u1D706u1D43B)?u1D465+sin(u1D706u1D43B)?u1D466 (4.40) 97 where the ? indicates the LOS frame. The vertical plane component is given as u1D706u1D449 = tan?1 uni0028.alt03u1D467 u1D447 ?u1D467 u1D465u1D447 ?u1D465 uni0029.alt03 (4.41) ?u1D706u1D449 = ??u1D467(u1D465u1D447 ?u1D465)+ ?u1D465(u1D467u1D447 ?u1D467) (u1D465u1D447 ?u1D465)2 +(u1D467u1D447 ?u1D467)2 (4.42) ?u1D434u1D44Du1D436 = u1D441u1D449 ?u1D706u1D449u1D462u1D449 (4.43) u1D462u1D449 = cos(u1D706u1D449)?u1D465+sin(u1D706u1D44C)?u1D467 (4.44) The total acceleration command is then ?u1D434u1D436 = ?u1D434u1D44Cu1D436u1D43Du1D43F + ?u1D434u1D44Du1D436u1D43Eu1D43F = u1D434u1D44Bu1D436u1D43Cu1D435 +u1D434u1D44Cu1D436u1D43Du1D435 +u1D434u1D44Du1D436u1D43Eu1D435 (4.45) where (u1D43C,u1D43D,u1D43E) are unit vectors. The input to the lateral pulse ring ?ring logic is the magnitude of the o?-axis command acceleration. This magnitude is given by ? = uni221A.alt02 u1D4342u1D44Cu1D436 +u1D4342u1D44Du1D436 (4.46) and, the error signal phase is given by u1D6FE = tan?1(u1D434u1D44Du1D436/u1D434u1D44Cu1D436). (4.47) Parabolic and Proportional Navigation Guidance (PAPNG) PAPNG uses the same horizontal control law as PNG. In the vertical plane, a desired parabolic trajectory is described as ?u1D467u1D443 = ?u1D467u1D447 +u1D43E1?u1D465u1D45D +u1D43E2?u1D4652u1D443 (4.48) where the ?u1D465u1D443 and ?u1D467u1D443 are components of the rocket position in the target frame, ?u1D465u1D447 and ?u1D467u1D447 are components of the target position in the target frame, and u1D43E1 and u1D43E2 are de?ned by Equations (4.49)-(4.50). The desired angle at which the rocket passes through the target is de?ned as u1D6FDu1D453. 98 Then u1D43E1 = tan(u1D6FDu1D439) (4.49) u1D43E2 = ?(u1D43E1 uni221A.alt01(u1D465 u1D447 ?u1D465)2 +(u1D466u1D447 ?u1D466)2 + ?u1D467u1D447 ? ?u1D467u1D443) (u1D465u1D447 ?u1D465)2 +(u1D466u1D447 ?u1D466)2 . (4.50) The command acceleration in the vertical LOS plane is ?u1D434u1D44Du1D436 = u1D449u1D43F(u1D6FDu1D437 ?u1D6FD) u1D70F (4.51) where u1D449u1D43F is the rocket velocity magnitude in the LOS frame, u1D6FD is the ?ight path angle in the LOS frame, u1D70F is the acceleration command time constant, and u1D6FDu1D437 is the rocket?s desired ?ight path angle. The acceleration command is then given by ?u1D434u1D436 = ? ?u1D434u1D44Du1D436 sin(u1D6FD)u1D43Cu1D43F + ?u1D434u1D44Cu1D436u1D43Du1D43F ? ?u1D434u1D44Du1D436 cos(u1D6FD)u1D43Eu1D43F (4.52) = u1D434u1D44Bu1D436u1D43Cu1D435 +u1D434u1D44Cu1D436u1D43Du1D435 +u1D434u1D44Du1D436u1D43Eu1D435. (4.53) Trajectory Tracking (TT) In this method a command trajectory is assumed to be known prior to launch. From this infor- mation the error between the desired path and the actual path can be computed. The magnitude of this error is used as control input to the thruster ?ring logic. The magnitude of the error and the error phase angle are ? = uni221A.alt02 u1D4522u1D44C +u1D4522u1D44D (4.54) u1D6FE = tan?1(u1D452u1D467/u1D452u1D44C). (4.55) Lateral Pulse Jet Firing Logic In order for a jet to ?re, four requirements must be met. The following two requirements apply to all lateral jets. 99 ? The error magnitude must be greater than some tolerance. ? The time between thruster ?rings must be greater than some threshold. The next two requirements are for a speci?c lateral jet. ? The di?erence between the error phase angle and the individual pulse jet force must be less than some threshold angle. ? The pulse jet has not been previously ?red. The last condition enforces the requirement that each jet can only be ?red once. In [59] the path control system is designed to track a speci?c trajectory. To do this the rocket?s measured position is compared to the commanded position to create an error signal, which in the rocket body frame is given by ? ?? ?? ? u1D452u1D44B u1D452u1D44C u1D452u1D44D ? ?? ?? ? = ? ?? ?? ? u1D450u1D703u1D450u1D713 u1D450u1D703u1D460u1D713 ?u1D460u1D703 u1D460u1D719u1D460u1D703u1D450u1D713 ?u1D450u1D719u1D460u1D713 u1D460u1D719u1D460u1D703u1D460u1D719 +u1D450u1D719u1D450u1D713 u1D460u1D719u1D450u1D703 u1D450u1D719u1D460u1D703u1D450u1D719 +u1D460u1D719u1D460u1D713 u1D450u1D719u1D460u1D703u1D460u1D713 ?u1D460u1D719u1D450u1D713 u1D450u1D719u1D450u1D703 ? ?? ?? ? ? ?? ?? ? u1D465u1D436 ?u1D465 u1D466u1D436 ?u1D466 u1D467u1D436 ?u1D467 ? ?? ?? ? (4.56) where u1D719, u1D703, and u1D713 are the roll, pitch, and yaw Euler angles. Magnitude and phase error are de?ned as ? = uni221A.alt02 u1D4522u1D44C +u1D4522u1D44D (4.57) u1D6FE = tan?1(u1D452u1D467/u1D452u1D44C). (4.58) With this information available at each computation cycle, four criteria are checked to see if an individual thruster should be ?red (lateral jet thrust is modeled as a constant of a known duration). The four conditions are 1. The tracking error magnitude must be greater than some threshold ? >u1D452u1D461?u1D45Fu1D452u1D460. (4.59) 100 2. A certain amount of time must have passed since the last jet ?ring u1D461?u1D461? > ?u1D461u1D461?u1D45Fu1D452u1D460 (4.60) where u1D461? is the time of the most recent pulse jet ?ring. 3. The angle between the tracking error and the individual pulse jet under consideration must be less than a threshold ?u1D703u1D456 ?u1D45Du1D456?u1D6FE? ?u1D6FE(?u1D443u1D43D/2)?