01jun Palermo

Calhoun: The NPS Institutional Archive
DSpace Repository
Theses and Dissertations Thesis and Dissertation Collection
2001-06
Angular rate estimation for multi-body

spacecraft attitude control
Palermo, William J.
Monterey, California. Naval Postgraduate School
http://hdl.handle.net/10945/2389
Downloaded from NPS Archive: Calhoun

NAVAL POSTGRADUATE SCHOOL
Monterey, California
THESIS
ANGULAR RATE ESTIMATION FOR MULTI-BODY
SPACECRAFT ATTITUDE CONTROL
by
William J. Palermo
June 2001
Thesis Advisor: Brij N. Agrawal

Second Reader: Harold A Titus
Approved for public release; distribution is unlimited

REPORT DOCUMENTATION PAGE Form Approved OMB No. 0704-0188
Public reporting burden for this collection of information is estimated to average 1 hour per response, including
the time for reviewing instruction, searching existing data sources, gathering and maintaining the data needed, and
completing and reviewing the collection of information. Send comments regarding this burden estimate or any
other aspect of this collection of information, including suggestions for reducing this burden, to Washington
headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite
1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project
(0704-0188) Washington DC 20503.
1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED
June 2001 Engineers Thesis
4. TITLE AND SUBTITLE: Angular Rate Estimation for Multi-Body Spacecraft 5. FUNDING NUMBERS
Attitude Control
6. AUTHOR(S) William J. Palermo
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING
Naval Postgraduate School ORGANIZATION REPORT
Monterey, CA 93943-5000 NUMBER
9. SPONSORING / MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING / MONITORING
N/A AGENCY REPORT NUMBER
11. SUPPLEMENTARY NOTES The views expressed in this thesis are those of the author and do not reflect the official
policy or position of the Department of Defense or the U.S. Government.
12a. DISTRIBUTION / AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE
13. ABSTRACT (maximum 200 words)
Spacecraft with high performance attitude control systems requirements have traditionally relied on imperfect mechanical
gyroscopes for primary attitude determination. Gyro bias errors are corrected with a Kalman filter algorithm that uses updates
from precise attitude sensors like star trackers. Gyroscopes, however, have a tendency to degrade or fail on orbit, becoming a
life limiting factor for many satellites. When errors become erratic, pointing accuracy may be lost during short star gaps.
Unpredictable gyro degradations have impacted NASA spacecraft missions such as Skylab and Hubble Space Telescope as
well as several DoD and ESA satellites. An alternative source of angular rate information is a software implemented real time
dynamic model. Inputs to the model from internal sensors and known spacecraft parameters enable the tracking of total system
angular momentum from which body rates can be determined. With this technique, the Kalman filter algorithm provides error
corrections to the dynamic model. The accuracy of internal sensors and input parameters determine the effectiveness of this
angular rate estimation technique. This thesis presents the background for understanding and implementation of this technique
into a representative attitude determination system. The system is incorporated into an attitude simulation model developed in
SIMULINK to evaluate the effects of dynamic modeling errors and sensor inaccuracies. Results are presented that indicate that
real time dynamic modeling is an effective method of angular rate determination for maneuvering multi-body spacecraft
attitude control systems.
14. SUBJECT TERMS dynamic gyro, Kalman filter, attitude determination, rate estimation, star 15. NUMBER OF
trackers, attitude simulation, multi-body dynamics, quaternions, MATLAB, SIMULINK PAGES
16. PRICE CODE

17. SECURITY 18. SECURITY 19. SECURITY 20. LIMITATION
CLASSIFICATION OF CLASSIFICATION OF THIS CLASSIFICATION OF OF ABSTRACT
REPORT PAGE ABSTRACT
Unclassified Unclassified Unclassified UL
NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)
Prescribed by ANSI Std. 239-18
i
THIS PAGE INTENTIONALLY LEFT BLANK
ii
ANGULAR RATE ESTIMATION FOR MULTI-BODY SPACECRAFT

ATTITUDE CONTROL
William J. Palermo
Lieutenant, United States Navy
B.S., United States Naval Academy, 1992
Submitted in partial fulfillment of the

requirements for the degree of
MASTER OF SCIENCE IN ASTRONAUTICAL ENGINEERING

and
AERONAUTICAL AND ASTRONAUTICAL ENGINEER
from the
NAVAL POSTGRADUATE SCHOOL

June 2001
Author: ___________________________________________
William J. Palermo
Approved by: ___________________________________________

Brij N. Agrawal, Thesis Advisor
___________________________________________
Harold A. Titus, Second Reader
___________________________________________
Max F. Platzer, Chairman
Department of Aeronautics and Astronautics
iii
iv
ABSTRACT
Spacecraft with high performance attitude control systems requirements have

traditionally relied on imperfect mechanical gyroscopes for primary attitude
determination. Gyro bias errors are corrected with a Kalman filter algorithm that uses
updates from precise attitude sensors like star trackers. Gyroscopes, however, have a
tendency to degrade or fail on orbit, becoming a life limiting factor for many satellites.
When errors become erratic, pointing accuracy may be lost during short star gaps.
Unpredictable gyro degradations have impacted NASA spacecraft missions such as
Skylab and Hubble Space Telescope as well as several DoD and ESA satellites. An
alternative source of angular rate information is a software implemented real time
dynamic model. Inputs to the model from internal sensors and known spacecraft
parameters enable the tracking of total system angular momentum from which body rates
can be determined. With this technique, the Kalman filter algorithm provides error
corrections to the dynamic model. The accuracy of internal sensors and input parameters
determine the effectiveness of this angular rate estimation technique. This thesis presents
the background for understanding and implementation of this technique into a
representative attitude determination system. The system is incorporated into an attitude
simulation model developed in SIMULINK to evaluate the effects of dynamic modeling
errors and sensor inaccuracies. Results are presented that indicate that real time dynamic
modeling is an effective method of angular rate determination for maneuvering multi-
body spacecraft attitude control systems.
v
vi
TABLE OF CONTENTS
I. INTRODUCTION........................................................................................................1
A. ANGULAR RATE ESTIMATION ................................................................1
B. OPTIMAL ESTIMATION..............................................................................1
C. ANGULAR RATE ESTIMATION FROM VECTOR
MEASUREMENTS .........................................................................................2
D. FULL ORDER ESTIMATORS......................................................................4
E. REDUCED ORDER ESTIMATORS.............................................................5
F. STAR SENSOR BASED ESTIMATORS ......................................................7
G. ANALYTICAL RATE DETERMINATION.................................................8
H. ATTITUDE ESTIMATION FROM CALCULATED RATE .....................9
I. PURPOSE .........................................................................................................9
II. SIMULATION OVERVIEW....................................................................................11
A. CONCEPTUAL ATTITUDE SIMULATION MODEL ............................11
B. SPACECRAFT MOTION.............................................................................13
C. CONTROL SYSTEM....................................................................................14
D. ATTITUDE DETERMINATION SYSTEM ...............................................14
III. DYNAMICS AND CONTROL.................................................................................17
A. BIFOCAL RELAY MIRROR SATELLITE DYNAMICS........................17
1. Rigid Body Equations of Motion ......................................................19
2. Multi-Body Equations of Motion......................................................20
3. Moments of Inertia.............................................................................21
a. Rate of Change of Spacecraft Inertia Matrix ........................23
4. Angular Velocity ................................................................................24
5. Angular Momentum ..........................................................................25
6. Solving for Spacecraft Angular Rate ...............................................25
B. COMMAND ...................................................................................................26
C. CONTROL .....................................................................................................28
D. DISTURBANCE TORQUES ........................................................................29
1. Gravity Gradient Torque ..................................................................29
2. Other External Disturbance Torques ..............................................31
E. MAGNETIC MOMENTS .............................................................................31
1. Magnetic Field Model........................................................................31
2. Magnetic Control Torque..................................................................33
IV. ATTITUDE REPRESENTATION AND KINEMATICS......................................35
A. QUATERNION DEFINITION AND CHARACTERISTICS ...................35
1. Meaning of Quaterions ......................................................................36
2. Attitude Matrix ..................................................................................36
3. Quaternion Multiplication ................................................................37
4. Coordinate Rotations .........................................................................38
5. Quaternion Inverse and Identity ......................................................39
vii
6. Quaternion Error...............................................................................39
7. Vector Transformations ....................................................................40
B. QUATERNION KINEMATICS...................................................................40
1. Continuous Kinematics .....................................................................41
2. Discrete Kinematics ...........................................................................42
V. KALMAN FILTER....................................................................................................45
A. RECURSIVE DISCRETE KALMAN FILTER..........................................45
B. ERROR STATE EXTENDED KALMAN FILTER
IMPLEMENTATION ...................................................................................47
1. State Variables....................................................................................47
2. Attitude Propagation Error Correction Methods ...........................47
3. Plant Model.........................................................................................48
4. State Transition Matrix .....................................................................49
5. Kalman Filter Prediction Equations ................................................49
6. Plant Noise Covariance......................................................................50
7. Kalman Filter Initialization ..............................................................51
8. Sensor Measurement Update ............................................................51
a. Measurement Vector...............................................................51
b. Predicted Measurement ..........................................................52
c. Measurement Residual ...........................................................53
9. Observation Matrix............................................................................54
10. Kalman Gain ......................................................................................54
11. Kalman Filter Correction Equations ...............................................55
C. CONTROLLER DESIGN IMPLICATIONS..............................................56
VI. DYNAMIC RATE CALCULATION.......................................................................57
A. DISCRETE EQUATIONS OF MOTION....................................................57
B. MOMENTUM CORRECTION FROM KALMAN FILTER
UPDATES .......................................................................................................58
C. INPUTS AND ERROR SOURCES ..............................................................58
1. External Control and Disturbance Torques....................................59
2. Reaction Wheel Relative Momentum...............................................60
3. Moment of Inertia Calculations ........................................................61
4. Appendage Relative Momentum ......................................................61
D. CALCULATED ANGULAR RATE ERROR CORRECTION ................63
VII. RESULTS ...................................................................................................................65
A. BASELINE SIMULATION ..........................................................................65
1. Simulation Input Parameters............................................................65
2. Command Attitude Profile ................................................................67
3. Attitude Determination System Performance Results....................70
4. Attitude Control System Performance Results ...............................73
B. DYNAMIC GYRO VS GYRO PERFORMANCE .....................................75
1. Continuous Star Tracker Coverage .................................................76
2. Gapped Star Tracker Coverage........................................................78
C. DYNAMIC GYRO PLANT ERROR ANALYSIS......................................79
viii
1. Disturbance Torque Modeling ..........................................................80
2. Rotation Axis Alignment Error........................................................81
3. Potentiometer Quantization..............................................................82
4. Moment of Inertia Update Frequency .............................................82
D. STAR TRACKER MEASUREMENT ERROR ANALYSIS ....................83
1. Star Tracker Accuracy ......................................................................83
2. Update Interval ..................................................................................84
3. Star Tracker Selection.......................................................................84
VIII. CONCLUSIONS ........................................................................................................87
A. SUMMARY ....................................................................................................87
B. RECOMMENDATIONS...............................................................................88
APPENDIX A: KALMAN FILTER BACKGROUND .....................................................89
A. PROBABILITY..............................................................................................89
B. EXPECTATION ............................................................................................90
C. LEAST-SQUARES ESTIMATION CONCEPT.........................................91
D. STATE ESTIMATE AND COVARIANCE ................................................92
E. STOCHASTIC LINEAR DYNAMIC MODEL ..........................................93
F. PROPAGATION OF ERRORS ...................................................................94
G. MEASUREMENT UPDATES ......................................................................95
APPENDIX B: MATLAB CODE........................................................................................99
A. INPUT PARAMETERS AND SIMULATION CALL ...............................99
B. SIMULATION RESULTS PLOTS ............................................................104
1. Commanded Profile Plots................................................................104
2. Attitude Determination Performance Results...............................105
3. Attitude Control Performance Results ..........................................106
LIST OF REFERENCES ....................................................................................................109
INITIAL DISTRIBUTION LIST .......................................................................................111
ix
x
LIST OF FIGURES
Figure 1. Basic Kalman Filter Block Diagram..................................................................2

Figure 2. Conceptual Attitude Simulation Block Diagram.............................................12
Figure 3. Top Level SIMULINK Attitude Simulation Model ........................................13
Figure 4. SIMULINK Attitude Determination Subsystem..............................................15
Figure 6. SIMULINK Subsystem Diagram: X-Axis Rotation Matrix ...........................22
Figure 7. SIMULINK Subsystem Diagram: Spacecraft Moment of Inertia Matrix.......23
Figure 8. SIMULINK Subsystem Diagram: Rate of Change of Spacecraft Inertia
Matrix...............................................................................................................24
Figure 9. SIMULINK Subsystem Diagram: Spacecraft Dynamics................................26
Figure 10. SIMULINK Subsystem Diagram: Feed Forward Torque Command .............27
Figure 11. SIMULINK Subsystem Diagram: Reaction Wheel Control...........................29
Figure 12. SIMULINK Subsystem Diagram: Orbital Reference Propagation.................30
Figure 13. SIMULINK Subsystem Diagram: Gravity Gradient Torque Model ..............30
Figure 14. SIMULINK Subsystem Diagram: Earth Magnetic Field Vector in Orbital
Coordinates ......................................................................................................32
Figure 15. SIMULINK Subsystem Diagram: Magnetic Field Conversion to Body
Coordinates ......................................................................................................33
Figure 16. SIMULINK Subsystem Diagram: Magnetic Torquers for Momentum
Dumping ..........................................................................................................33
Figure 17. SIMULINK Subsystem Diagram: Attitude Matrix from Quaternions ...........37
Figure 18. SIMULINK Subsystem Diagram: Quaternion Multiplication........................38
Figure 19. SIMULINK Subsystem Diagram: Error Quaternion ......................................40
Figure 20. SIMULINK Subsystem Diagram: Continuous Kinematics............................42
Figure 21. SIMULINK Subsystem Diagram: Discrete Attitude Propagator....................43
Figure 22. Discrete Kalman Filter Loop............................................................................46
Figure 23. SIMULINK Subsystem Diagram: State Transition Matrix ............................49
Figure 24. SIMULINK Subsystem Diagram: Kalman Filter Prediction Step..................50
Figure 25. SIMULINK Subsystem Diagram: Star Tracker Measurement.......................52
Figure 26. SIMULINK Subsystem Diagram: Inertial Star Vector...................................53
Figure 27. SIMULINK Subsystem Diagram: Measurement Residual .............................54
Figure 28. SIMULINK Subsystem Diagram: Observation Matrix ..................................54
Figure 29. SIMULINK Subsystem Diagram: Kalman Gain ............................................55
Figure 30. SIMULINK Subsystem Diagram: Covariance Updeate .................................56
Figure 31. SIMULINK Subsystem Diagram: Dynamic Gyro..........................................58
Figure 32. SIMULINK Subsystem Diagram: Error Corrupted Reaction Wheel
Momentum Measurement ................................................................................61
Figure 33. SIMULINK Subsystem Diagram: Appendage Relative Angle and Rate
Measurement....................................................................................................62
Figure 34. SIMULINK Subsystem Diagram: Dynamic Gyro Inputs ...............................63
Figure 35. Baseline Commanded Attitude Profile ............................................................68
Figure 36. Baseline Commanded Angular Rate Profile ....................................................68
xi
Figure 37. Baseline Commanded Relative Angle to Receive Telescope Profile ..............69
Figure 38. Total Spacecraft Angular Momentum Profile..................................................69
Figure 39. Baseline Dynamic Gyro Angular Momentum Error........................................70
Figure 40. Baseline Kalman Filter Attitude and Rate Bias Errors ....................................71
Figure 41. Baseline Estimated Attitude Quaternion Error ................................................72
Figure 42. Baseline Estimated Angular Rate Error...........................................................73
Figure 43. Reaction Wheel Control Torques.....................................................................74
Figure 44. Reaction Wheel Angular Momentum ..............................................................74
Figure 45. Baseline Control Attitude Quaternion Error....................................................75
Figure 46. Baseline Control Angular Rate Error...............................................................75
Figure 47. SIMULINK Subsystem Diagram: Gyro Measurement and Bias Error
Tracker. ............................................................................................................76
Figure 48. Results: Dynamic Gyro vs. Gyro with Continuous Star Coverage ..................77
Figure 49. Results: Dynamic Gyro vs. Gyro with Gapped Star Coverage ........................79
Figure 50. Baseline Attitude Determination Performance Results ...................................80
Figure 51. Gravity Gradient Disturbance Modeled in the Dynamic Gyro ........................81
Figure 52. Reaction Wheel Alignment Error Effects on the Dynamic Gyro ....................81
Figure 53. Reduced Potentiometer Quantization Effect on Dynamic Gyro
Performance .....................................................................................................82
Figure 54. 10 Hz Moment of Inertia Calculation..............................................................83
Figure 55. Star Tracker 1 Sigma Variance: 1x10 -4 [left] vs 2.5x10-5 [right] ....................84
Figure 56. Star Tracker Update Interval: 2 Seconds [left] vs 6 Seconds [right] ..............84
Figure 57. Star Tracker Selection Comparison: Even Distribution [left] vs One
Tracker Favored [right] ....................................................................................85
xii
LIST OF TABLES
Table 1. List of Symbols ...............................................................................................xvi

Table 2. Discrete Kalman Filter Equations ....................................................................46
Table 3. Simulation Input Parameters............................................................................67
Table 4. Plant Error Covariance Matricies.....................................................................78
xiii
xiv
LIST OF SYMBOLS
MATLAB/
Symbol SIMULINK Description Units
Variable
t dt Computer Step Size s
mT m1 Mass of Primary Body (Transmit Telescope) kg
mR m2 Mass of Secondary Body (Receive Telescope) kg
IT I1 Primary Body Inertia Matrix kgm2
IR I2 Secondary Body Inertia Matrix kgm2
I I Spacecraft Inertia Matrix kgm2
rT d1 Primary Body Center of Mass Offset m
rR d2 Secondary Body Center of Mass Offset m
mu Earth Gravitational Constant km3/s2
Re Re Earth Radius km
h h Orbit Altitude km
r r Orbit Radius km
e we Earth Rotation Rate rad/s
o w0 Orbital Rate rad/s
nu True Anomaly rad
i inc Orbital Inclination rad
e1 Earth Magnetic Diploe Tilt Angle rad
Angle of Magnetic Dipole Normal to Line of
u u rad
Ascending Nodes
v
B B Earths Magnetic Field Vector N/Am
Km Magnetic Field Constant Nm2/A
v
M cmd Mcmd Magnetic Dumping Control Torque Nm
v
MG Mdg Gravity Gradient Disturbance Torque Nm
Mdp Periodic Disturbance Torque Nm
Mds Secular Disturbance Torque Nm
v
H H Total System Angular Momentum kgm/s
v System Angular Momentum Neglecting Relative
HS- H kgm/s
Internal Momentum
v
Hw hw Reaction Wheel Relative Angular Momentum kgm/s
v
Hrel I2bd Secondary Body Relative Angular Momentum kgm/s
v
Hcorr dHc Kalman Filter Momentum Correction kgm/s
b Appendage Relative Position Angle rad
v
rel bd Appendage Relative Angular Rate Vector rad/s
v
& rel bdd Appendage Relative Angular Acceleration rad/s
Vector
xv
Vector
v
q q Spacecraft Attitude Quaternion Vector
v
w Spacecraft Angular Rate rad/s
v
& wd Spacecraft Angular Acceleration Vector rad/s2
v
Mc Mff Feed Forward Command Torque Nm/s
v
qc qc Commanded Attitude Quaternion Vector
v
c wc Commanded Angular Rate rad/s
v
& c wdc Commanded Angular Acceleration Vector rad/s2
v
q qkf Estimated Attitude Quaternion Vector
v
wkf Estimated Angular Rate rad/s
A A Spacecraft Attitude (Direction Cosine) Matrix
A Estimated Attitude (Direction Cosine) Matrix
v
x State Vector
v x Estimated Kalman State vector
x
P P State Covariance Matrix
Q Q Plant Error Covariance Matrix
R R Measurement Error Covariance Matrix
S State Transition Matrix
K K Kalman Gain Matrix
E E Measurement Selection Matrix
T T Star Tracker to Body Alignment Matrix
RW Body to Reaction Wheel Alignment Matrix
G Body to Gyroscope Alignment Matrix
Table 1. List of Symbols
xvi
ACKNOWLEDGMENTS
The author would like to thank several people for their contributions to this effort.
Dr. F. Landis Markley of NASA Goddard Space Flight Centers Guidance, Navigation
and Control Systems Engineering branch provided me with the background for the initial
research. His vast insight and experience in all aspects of spacecraft design and control
were invaluable to my understanding of attitude determination concepts. Also critical to
the research were Girard Manke and Craig Heatwole of The Aerospace Corporation.
Their hospitality and enthusiasm in a shared interest provided the focus for this study.
The cooperative efforts and support of Marc Becraft and Don Kolve from Boeing are also
appreciated. Finally, I would like to thank Professor Brij Agrawal and Hal Titus of the
Naval Postgraduate School for their inspiration and guidance.
xvii
xviii
I. INTRODUCTION
Equation Section 1
A. ANGULAR RATE ESTIMATION
The attitude control systems of todays high tech satellites require accurate
angular rate knowledge for attitude propagation and control loop damping in order to
meet pointing and tracking requirements. Gyroscopes offer varying degrees of precision
for direct measurement of angular rates and have been the primary attitude determination
sensor used by spacecraft for many years. However, high cost and low reliability has
lead users and designers to explore other options. Using dynamic calculations,
uncertainty based estimation algorithms or a combination of these two methods, onboard
processors can estimate spacecraft angular rates without measuring them directly.
Angular rate estimation techniques can be a viable alternative for back-up or even
primary attitude determination system in many control schemes.
B. OPTIMAL ESTIMATION
An optimal estimator is a computational algorithm that processes measurements

to deduce a minimum error estimate, in accordance with stated optimality criterion, of the
state of a system by utilizing: knowledge of system measurement dynamics, assumed
statistics of system noises and measurement errors, and initial condition information.
[Ref. 1] Estimation techniques provide filtering, smoothing and prediction of state
variables in an imperfect model based on imperfect measurement update data. The most
common optimal estimator used in stochastic systems is the Kalman filter. The dynamic
model error and measurement error are assumed as zero mean Guassian white noise
processes with known covariance. Estimators are commonly used even in systems where
all state variables required by controller can be measured.
1
SYSTEM
MEASUREMENT A PRIORI
PLANT MODEL
ERROR SOURCES INFORMATION
ERRORS
SYSTEM STATE
v
STATE
v
OBSERVATION ESTIMATE
v
x y x
SYSTEM KALMAN
MEASUREMENT
FILTER
Figure 1. Basic Kalman Filter Block Diagram
Maybeck identifies three basic reasons why deterministic methods are insufficient
for describing real systems. First, no mathematical model is perfect. There are always
effects that are necessarily neglected to make the model practical and even modeled
effects are only approximations to what is physically occurring. Second, dynamic
systems are driven not only by control inputs but also by disturbances which we can
neither control nor model deterministically. Lastly sensors do not provide perfect and
complete data. [Ref. 2]
An important advantage of using an optimal estimator such as the Kalman filter is

that the attitude determination output does not affect the controller design. Therefore the
development of an optimal controller can be accomplished independently. An attitude
control algorithm optimized for ideal state inputs will remain optimized with state
estimates provided by an optimal estimator.
C. ANGULAR RATE ESTIMATION FROM VECTOR MEASUREMENTS
The concept of using attitude sensor data to produce an estimate of spacecraft

angular rate without gyroscopes is not new. Gyroless attitude and angular rate estimation
has been a prime concern for small inexpensive spacecraft that do not carry gyroscopes
but still require rate information for attitude propagation and control. Estimation
techniques also provide options for complex spacecraft that require back-up control
modes in the event of gyro failures. The problem of angular rate estimation can be
2
treated separately from attitude determination. Several reliable algorithms have been
developed that produce angular rate using on-board processors.
It is possible to extract angular rate directly from time derivatives of measured

vectors resolved in the body coordinate frame and known in an inertial frame from a
model or almanac data. This technique, however, requires at least two non-parallel
vector measurements. It also exhibits time lag and produces very noisy data since the
algorithms depend on derivatives of already noisy measurements. A better estimate of
angular rate can be achieved by treating the problem as stochastic. For this method an
estimating filter is applied that uses dynamic equations of motion to take advantage of
past data in a recursive sense. This method has the additional advantage of being able to
update rate estimates at time steps when only one vector measurement is available. The
vector derivative is treated as a noisy measurement update to a Kalman rate estimator.
A Kalman filter requires a linear state-space dynamic model but the dynamics of a
spacecraft in the rotating body coordinate frame are coupled and nonlinear. Attitude
determination schemes deal with this problem in different ways. Several methods have
been proposed for interlacing two or three linear Kalman filters together to capture the
nonlinear dynamics [Ref. 3]. It is shown that the coupled equations of motion can be
written as a linear combination of the spacecraft angular rate and two other newly defined
vectors whose components are nonlinear combinations of the angular rate vector
elements. The differential equations for each of these new vectors can be written as
linear combinations of the other and the angular rate. Adding white noise vectors to these
new equations of motion produces a set of three stochastic linear models that can be used
in separate Kalman filters. The filters are interlaced with their estimated outputs treated
as deterministic inputs to each other.
An effective method of dealing with non-linear dynamics is presented for use by

the Pseudo-linear Kalman filter (PSELIKA) and the state-dependent algebraic Riccati
equation filter (SDARE) [Ref. 4]. For these unique rate estimation techniques the
equation of motion is transformed to express the nonlinear terms in angular rate as a
product of a matrix whose elements depend on the components of the angular rate vector
and the angular rate vector itself. The pseudolinear Kalman filter or PSELIKA is
3
implemented as an ordinary discrete Kalman filter with a time dependent state transition
matrix whose angular rate dependent elements are formed from the current estimate of
angular rate. The SDARE or state-dependent algebraic Riccati equation filter is also
implemented as a discrete Kalman filter derived from the same representation of the
dynamic equation. The Kalman gain used in this filter is computed from the solution of
an algebraic Riccati equation involving the angular rate dependent matrix. This
eliminates the need to propagate and update the filter state covariance which is normally
used in the gain calculation.
The pseudolinear filter concept can be extended to use quaternion measurements.

In the attitude determination system on board the Rossi X-Ray Timing Explorer (RXTE)
satellite advanced star tracker software directly provides information in the form of
quaternions [Ref. 5]. The same propagation algorithm is used with measurement updates
in the form of quaternions. Since the device yielding the attitude measurement also
conducts a star search and identification process a time delay is introduced. Two
algorithms are presented for overcoming this delay problem.
D. FULL ORDER ESTIMATORS
It is not necessary to treat the problem of angular rate estimation separately from
attitude determination. In fact, attitude sensor data also requires filtering to smooth out
measurement noise and produce a clean attitude estimates. High precision control laws
generally require both attitude and angular rate information. These components make up
the full system state vector.
Full order estimators produce optimal state estimates in the case where imperfect
measurement data related to all of the modeled states is continuously available.
Spacecraft with a triad of rate gyroscopes and attitude sensors have the required state
measurements but the bandwidth of the attitude sensors is normally too low for the
control system. Additionally, there may be periodic gaps in sensor coverage and certain
sensors do not provide reliable attitude information in all three orthogonal dimensions.
Therefore standard full order state estimators are generally not used.
4
In most practical implementations of gyro-based attitude control systems rate
information is used in the propagation stage of an attitude estimator based on kinematic
equations. The estimator utilizes discrete vector observations, resolved in both body-
fixed coordinates and a reference frame, to estimate spacecraft attitude and gyro drift
rates.
In an ideal attitude determination system where perfect knowledge of the

spacecraft angular rates is available, the accuracy of the kinematic model for determining
attitude is limited only by the processor time step and the initial attitude state. Attitude
sensors can be used to periodically reinitialize the kinematic model but the data they
provide is discontinuous and corrupted by environmental effects and sensor design
imperfections.
In real attitude determination systems, the source of angular rate information may
be noisy as well as biased, as is the case with spacecraft gyroscopes. This results in
additional errors in the spacecraft attitude determination system. For control systems
where attitude and angular rate information is critical, optimal estimation techniques are
employed to combat the effects of model and sensor inaccuracies. As gyros degrade the
estimator can be made more robust to plant error at the expense of attitude determination
accuracy. At some point the gyro outputs may become too erratic to meet control system
requirements.
E. REDUCED ORDER ESTIMATORS
Obtaining satellite angular rate estimates without the use of rate gyroscopes or
other deterministic rate data can be accomplished with a reduced order estimator. This
class of filters produces estimates of all modeled state variables when only a subset is
directly related to the measurement data. In this case, attitude sensor information is
available but no direct rate measurement is performed. Since gyro data is not available
no update to the spacecraft rate is available between attitude measurements.
The reduced estimator Kalman filter is formulated from the state space dynamic
equations of motion linearized about the current estimate of state. Again, for spacecraft
attitude control the state vector includes both the attitude and angular rate. The
5
deterministic inputs to the model are the known externally applied moments usually
consisting only of control torques. Unmodeled disturbances and other dynamic effects
must be accounted for by a robust plant noise model. The estimates of attitude and
angular rate are propagated forward by the dynamic model with discrete corrections from
the attitude sensors.
The performance of the reduced Kalman estimator can be improved by including

higher order dynamic effects in the system plant model and including disturbance torques
as deterministic inputs. However this costs extra processing power, requires additional
sensors and cannot account for all of the unknowns. Additional dynamic complexities
are introduced in multi-body satellites that have time varying moments of inertia,
changing centers of mass and flexibility modes. Determining the magnitude and effects
of disturbances and modeling simplifications for a particular spacecraft is an important
design consideration that often requires rough calculations, simulation and engineering
judgment. In general, since disturbances are of low amplitude and low frequency
compared to control torques, their effects can be accounted for as plant error by the filter.
One approach that avoids the use of the typically uncertain dynamic model
altogether is to treat the spacecraft as a noncooperative target whose rotation must be
tracked by a stochastic estimator [Ref. 6]. This concept is borrowed from tracking
theory where it has been widely used to estimate target motion. This algorithm adds
angular acceleration to the state vector and substitutes attitude states with the integrated
rate parameters to formulate a nine state linear Kalman filter. Since this set of variables
serves as an approximate third-order attitude parameterization, the size must be
controlled by sampling interval. Instead of using dynamics, time propagation of the
estimated state variables is performed in the proposed filter by modeling the spacecraft
angular acceleration as an exponentially autocorrelated stochastic process and using a
polynomial kinematic model.
In reduced order estimators, time steps for state propagation must be kept small to
minimize the effects of dynamic simplifications and linearization. The primary limitation
of the reduced estimator becomes evident during prolonged intervals between attitude
sensor measurement updates. Gaps in attitude sensor coverage must be kept short enough
6
to ensure that estimator errors do not grow beyond required precision of the attitude
control system.
F. STAR SENSOR BASED ESTIMATORS
Bandwidth and accuracy limitations of attitude sensors have precluded the use of
reduced order estimators for spacecraft missions with high performance pointing and
tracking requirements. Aside from star trackers, most sensors do not provide accurate
and stable enough attitude data for high performance tracking requirements.
Additionally, these sensors provide discrete measurements unlike gyroscopes which
produce nearly continuous data and can be sampled at extremely high frequency limited
only by processing capabilities. However, recent advancements in sensor technology
suggest that star trackers can provide updates accurately and frequently enough to be
considered for use as primary attitude sensors for a wide range of spacecraft missions.
Key technological improvements include a wide field of view, high sensitivity, low noise
equivalent angles and high bandwidth iteration rate.
It has been proposed that advanced high-bandwidth star sensors can be used as the
primary sensors for on-orbit attitude rate determination in place of angular rate
gyroscopes [Ref. 7]. A Kalman filter algorithm, using updates from star tracker
measurements, is presented that tracks errors in the attitude and angular rate for a
kinematic based state space model. The spacecraft angular rate used for attitude
propagation is separated into a nominal component and a component due to control and
disturbance moments. The algorithm assumes that the bandwidths of the control system
and disturbance effects are at least an order of magnitude smaller than the measurement
bandwidth. This allows the rate of change in spacecraft angular rate due to unmodeled
dynamic effects to be modeled as a first-order Gauss-Markov process. Nonlinear
kinematics are used in the attitude propagation phase which is performed external to the
estimator. The errors in the this model are then tracked using an ordinary linear discrete
Kalman filter.
7
G. ANALYTICAL RATE DETERMINATION
Another method of obtaining estimates of spacecraft angular rates is through

direct calculation from the dynamic equations of motion for the system. This rate
calculation is performed in real time based on known, derived and sensed internal
parameters of the spacecraft. The software implemented dynamic model can be adjusted
to include varying levels of complexity for multi-body spacecraft that operate in different
configurations. The Aerospace Corporation has patented one methodology for this
process called the Pseudo Gyro [Ref. 8].
Information from internal sensors that detect relative orientations and angular
rates of momentum exchange devices and appendages are used to continually update the
parameters used in the dynamic calculations. These parameters determine component
contributions to the system angular momentum and inertia dyadic. Control torques and
modeled disturbances are integrated to capture external dynamic effects. Using the total
system angular momentum and inertia dyadic and the relative momentum of internal
components and appendages from the dynamic model the spacecraft angular rate can be
calculated. If the appendage masses are small or their relative motions are slow the
spacecraft inertia matrix is slow to change. The moment of inertia calculations can be
performed at a lower frequency to save processing time. Cross product effects due to the
rotating coordinate frame must also be accounted for in the dynamic model since all
measurements are taken in the body frame.
The accuracy of the calculated rate is dependent on the quality of the dynamic
model and the sensor information available. Error sources include imperfect knowledge
of component or appendage inertia matricies and mass centers as well as relative angular
position and rate data from the internal sensors. Often relative rates are not measured
directly but derived from position encoders which adds extra noise to the momentum
along the axis of rotation of the appendage. Errors are also introduced through external
disturbance effects since they cannot be perfectly modeled. All of these dynamic
influences effect the model precision in varying degrees. Additionally, due to the
finite/discrete processing capabilities of the computer performing the calculations there
will be a slight drift over time from the true state even in the case of perfect sensors and
8
input data. This drift is obviously minimized by increasing the frequency of the discrete
model calculations.
H. ATTITUDE ESTIMATION FROM CALCULATED RATE
The rate estimate produced by the dynamic calculations can be fed into a Kalman
filter in place measurement data from a gyroscope. The filter receives attitude update
information from attitude sensors to produce smoothed attitude reference and angular rate
error estimates to enhance satellite attitude determination accuracy. The plant error
introduced into the filter is a combination of all the internal sensor and modeling errors
that are used as inputs to the dynamic model. All known sensor biases must be
incorporated into the plant model since the Kalman filter assumes errors to be zero mean
Gaussian. Depending on different modes of operation, different plant noise models may
be required for acceptable filter performance. Normally, filters that receive
measurements from gyroscopes are designed to track gyro biases which remain relatively
constant in the body frame. The rate error from the dynamic estimate, however, exhibits
different characteristics. Since the dynamic model is external to the Kalman filter higher
order effects can be more easily incorporated into the dynamics but onboard processing
capability may still limit the complexity.
I. PURPOSE
The objective of this study is to develop and evaluate a practical attitude and rate
estimation scheme for a multi-body spacecraft attitude control system that incorporates
both real time angular rate calculation from the system dynamic model (dynamic gyro)
and a Kalman filter estimator with attitude sensor updates provided by star trackers. It is
hypothesized that acceptable attitude control performance can be realized by
maneuvering multi-body spacecraft without hardware gyroscope data using this
methodology.
The evaluation of this approach is performed through simulation using a model

developed in SIMULINK. A MATLAB script file is used to set up the necessary
initialization parameters and system specifications. Simulation results are presented
graphically in MATLAB plots. The performance of the developed gyroless attitude
9
determination system is compared to a conventional gyro-based system that uses the
same Kalman filter estimation algorithm and attitude updates. Results are also presented
that analyze the affects of several major error sources on the performance of this dynamic
gyro based attitude determination system.
10
II. SIMULATION OVERVIEW
Equation Section (Next)
The evaluation and analysis of attitude determination scheme proposed in this

thesis is done through modeling and simulation. The context for the development of the
attitude determination system is a multi-body attitude control system for a maneuvering
spacecraft. An attitude simulation model for the spacecraft that includes vehicle
dynamics, determination and control subsystems as well as modeled error sources is
developed using MATLAB/SIMULINK. An overview of the simulation and top level
subsystems is presented in this chapter. Subsequent chapters provide the descriptions of
the subsystems and derivation of the equations on which the model is based. In the
following chapters, actual SIMULINK diagrams are presented that show the
implementation of associated subsystems.
A. CONCEPTUAL ATTITUDE SIMULATION MODEL
A conceptual representation of the simulation developed for this study is

illustrated in Figure 2.
11
EXTERNAL
DISTURBANCE
TORQUES
ANGULAR ATTITUDE
DYNAMICS RATE QUATERNION
(EQUATIONS OF MOTION) SPACECRAFT
KINEMATICS MOTION
EXTERNAL
MOMENT OF TORQUE
INERTIA
MOMENTUM
STAR
EXCHANGE GYROS
TRACKERS
GYRO
RELATIVE
or
ANGLE DYNAMIC GYRO BIAS
RELATIVE CORRECTION
MOMENTUM
APPENDAGE ERROR
MOTION INTERNAL MOMENT OF STATE
SENSORS INERTIA KALMAN
DRIVE CONTROL FILTER
MOMENTUM
MOTOR DYNAMIC
CONTROL GYRO
TROQUE
REACTION MOMENTUM
WHEELS MAGNETIC CORRECTION
TORQUERS
ESTIMATED
CONTROLLER RATE
ATTITUDE
CONTROL PROPAGATOR
LAWS
ATTITUDE
COMMAND CORRECTION
APPENDAGE ESTIMATED
COMMAND QUATERNION
MOTION RATE ATTITUDE
FEED COMMAND DETERMINATION COMPUTER
FORWARD QUATERNION
TORQUE
COMMAND
PROFILE
Figure 2. Conceptual Attitude Simulation Block Diagram
12
The model is divided into subsystems represented in color shaded blocks. This
breakdown reduces the complexity of the overall model into manageable segments to aid
in design and analysis. Arrows indicate the dynamic coupling and flow of data between
subsystem blocks. The top level SIMULINK diagram that implements the concept is
shown in Figure 3.
Figure 3. Top Level SIMULINK Attitude Simulation Model
B. SPACECRAFT MOTION
The actual spacecraft attitude motion is simulated in the rotational dynamics and
kinematics subsystems with inputs and outputs represented by solid black lines. Relative
motion of the secondary body or appendage is treated in a separate subsystem with
dynamic effects coupled into the overall spacecraft motion through momentum exchange
directly related to the drive motor rate. Dynamic effects of reaction wheel control are
also realized through momentum exchange. Magnetic control effects are input to the
13
spacecraft dynamics as external torques along with modeled and approximated
disturbance torques. The relative appendage motion also causes changes in the
spacecrafts moments of inertia used in the dynamic model. The kinematics subsystem
propagates the actual spacecraft attitude quaternions referenced to an inertially fixed
coordinate frame.
C. CONTROL SYSTEM
Primary spacecraft attitude control is conducted through momentum exchange

with reaction wheels. The reaction wheel commands are based on control laws and up-
linked feed forward torque command profile. The control laws use errors between the
measured and commanded (desired) spacecraft angular rates and attitude quaternions.
The relative momentum generated in the wheels is subtracted from the spacecraft
momentum in the dynamics subsystem. Magnetic torques are generated based on
reaction wheel momentum build up when the system is set for momentum dumping.
These torques are directly applied in the spacecraft dynamics block.
D. ATTITUDE DETERMINATION SYSTEM
The attitude determination subsystem is of primary concern in this thesis. It

incorporates the dynamic modeling concept and an error state Kalman filter in order to
correct for attitude propagation errors. The option of using modeled mechanical gyros to
determine spacecraft angular rate is also maintained in order to allow comparison
simulations to be conducted. The concept of using a dynamic model for analytical rate
determination will be referred to in the sequel as the dynamic gyro.
The attitude determination computer simulation uses a trigger to control the

bandwidth of discrete calculations. The basic data flow within the attitude determination
computer is shown in the large rectangular subsystem block in Figure 2. The equivalent
SIMULINK diagram that implements this subsystem is shown in Figure 4.
14
Figure 4. SIMULINK Attitude Determination Subsystem
Random noise and drift rates are added to the actual angular velocity vector for
the simulated rate input of the mechanical gyro option. Bias corrections from the Kalman
estimator are applied to the angular rate before it is used in the attitude propagator. When
the dynamic gyro option is simulated the angular rate is analytically determined from
measured and known spacecraft parameters fed into a discrete dynamic model. Artificial
errors and noise are applied to measured and derived parameters though internal sensor
models. Momentum corrections to the dynamic gyro are derived from Kalman filter
updates. The calculated rate is then supplied to the attitude propagator.
The error state Kalman filter algorithm depends on updates generated by star
tracker measurements. A star tracker model is used to produce artificial noise corrupted
horizontal and vertical measurements related to the star tracker orientation. The Kalman
filter tracks rate errors for bias correction of gyro measurements or momentum
corrections for the dynamic gyro. It also produces attitude corrections that are applied to
the attitude propagator output.
15
The attitude Propagator uses discrete kinematics to convert angular rate to
quaternion attitude. These estimated parameter are then used by the attitude control
algorithm to complete the feedback loop.
16
III. DYNAMICS AND CONTROL
As three-axis stabilized spacecraft become more technologically advanced, their

operations require more slewing maneuvers and their dynamic complexities increase.
This often leads to multiple rigid or flexible components that have independent pointing
and tracking requirements. Complex satellites often consist of a primary payload that
demands strict pointing control while directional telemetry and command antennas or
secondary payloads are controlled independently for other functions such as tracking a
ground station throughout its maneuvers.
Here we consider the dynamics of a spacecraft that consists of a primary body

with momentum exchange control devices and a coupled rigid secondary body or
appendage. The secondary body rotates with one degree of freedom relative to the
spacecraft about an axis through the centers of mass of both bodies. Under these
conditions, the position of the spacecraft center of mass remains stationary during
appendage relative motion.
A. BIFOCAL RELAY MIRROR SATELLITE DYNAMICS
The example spacecraft for which the equations of motion are derived is the
Bifocal Relay Mirror satellite shown in Figure 5. This spacecraft is designed to receive
laser energy from a ground station through a receive telescope and to redirect it through
an optically coupled transmit telescope to a different point on the ground. The primary
body is the transmit telescope and the rigidly attached spacecraft bus while the smaller
receive telescope is treated as the secondary body or appendage. Weight and balance
design ensures that the centers of mass of the two bodies lie close to the coupled axis of
rotation so that the system center of mass is nearly constant during relative motion.
Pointing control of the spacecraft is accomplished with reaction wheels while a drive
motor is used to control the relative angle between the transmit and receive telescopes.
17
PRIMARY BODY
SECONDARY BODY Transmit Telescope
Receive Telescope and Bus
y X
x Relative
Y x Rotation
Axis

y z
z
Z
Laser
Beam
Up
Laser
Beam
Down
EARTH
Figure 5. Bifocal Relay Mirror Satellite
18
There are three coordinate systems defined for the development of the dynamic
equations. These coordinate frames are depicted in the system diagram shown in figure
x. The xyz coordinate system is fixed to the primary body with the origin at its center of
mass. The x-axis is oriented parallel to the axis of rotation between the primary and
secondary bodies and the z axis is parallel to the optical axis of the telescope. The y axis
is defined such that the xyz coordinate system is a right-handed mutually orthogonal
frame. The xyz coordinate system is similarly oriented to the secondary body with its
origin at the center of mass. The x and x axes remain parallel during appendage motion
while the angular displacements of the y and z axes from the y and z axes respectively
are defined by the relative rotation angle . The equations of motion are derived for the
spacecraft body coordinate system, XYZ, which is parallel to xyz frame. Its origin is at
the total spacecraft center of mass. Unit vectors along the X, Y, and Z spacecraft body
r r r
axes are given by i , j and k respectively.
1. Rigid Body Equations of Motion
The angular equations of motion are derived from the application of Newtons
second law to rotational dynamics. In the general case the equation of motion defined in
an inertial frame for a rigid body about an arbitrary point P is given by
v v&
MP = H (
v& v&
)
P + c r cdm (3.1)
m
v v
where M P is the net external torque applied to the body about P, HP is the angular
v v
momentum of the body about P, rc is the position of the center of mass relative to P, c is
the vector from the center of mass to the position of dm in the body, and dm is an
incremental unit of mass within the body
v v
If point P is coincident with the body center of mass, then c = 0 and the equation
simplifies to
v v&
M=H (3.2)
Equation (3.2) applies in inertial reference frames. To extend it to body coordinates
where it can be employed, we need to understand the concept of vector derivatives in a
19
rotating coordinate frame. If the body frame rotates relative to an inertial reference frame
v v
with angular velocity then the derivative of any vector A in inertial coordinates can
be related to the derivative in the rotating body coordinates by
v& v& v v
A =A + A B (3.3)
I B
Applying this relation to the angular momentum, the equation of motion for a rigid body
in rotating body coordinates becomes
v v& v v
M=H + H (3.4)
2. Multi-Body Equations of Motion
For the two body Bifocal Relay Mirror system with reaction wheels shown in
figure 1, the total system angular momentum in body XYZ coordinates can be written
v v v v
H = HS- + Hrel + Hw (3.5)
v v
where H is the total system angular momentum of the spacecraft, HS- is the angular
momentum of the system due to the rotation of the body coordinate frame, neglecting the
v
contribution due to relative motion of the receive telescope and reaction wheels, Hrel is
v
the angular momentum due to the relative motion of the secondary body, and H w is the
angular momentum due to the relative motion of the reaction wheels
Since the Bifocal Relay Mirror spacecraft is approximated as a system of rigid

v
bodies, we can substitute the total system angular momentum, H , into Equation (3.4) to
get the multi-body spacecraft equation of motion
v v& v& v& v v v v
M=H S- + Hrel + H w + ( HS- + Hrel + Hw ) (3.6)
v
where M is the net external torque applied to the spacecraft about its center of mass
including all control and disturbance torques.
20
3. Moments of Inertia
To determine these angular momentums we need to consider the moments of

inertia of the spacecraft and it components. We define IT to be the inertia matrix for the
primary body (transmit telescope and bus) about its center of mass in the xyz coordinate
frame.
ITxx -ITxy -ITxz

IT = -I Txy I Tyy -I Tyz (3.7)
-ITxz -ITzy ITzz

The vector from the center of mass of the primary body to the spacecraft center of mass is
given by
rT = x Ti + yT j + zT k (3.8)
The moment of inertia matrix for the receive telescope, IR, about its center of mass given
in the xyz coordinate frame is given by
I Rx'x' -IRx'y' -IRx'z'

R
IR = -IRx'y' Iy'y'
R
-I y'z' (3.9)
-IRx'z' -IRz'y' IRz'z'

The vector from the center of mass of the receive telescope to the spacecraft center of
mass is given by
rR = x Ri + y R j + z R k (3.10)
To align the receive telescope inertia matrix with the spacecraft body coordinate
frame a transformation matrix is applied. This time varying matrix is the direction cosine
matrix that defines a single rotation. The axis of rotation is parallel to the body X axis
and is of magnitude . The transformation matrix is given by
1 0 0
TX = 0 cos( ) sin()
(3.11)
0 sin( ) cos()
21
The SIMULINK subsystem used to generate this x-axis rotation matrix given an input
angle is shown in Figure 6.
Figure 6. SIMULINK Subsystem Diagram: X-Axis Rotation Matrix
To obtain the total system inertia matrix about the spacecraft center of mass in the
body coordinate frame XYZ the parallel axis theorem is applied to the inertia matricies of
each body. Additionally the rotation matrix is applied to the receive telescope inertia
matrix to align it with the body coordinates. If mT and mR are the masses of the transmit
telescope and receive telescope respectively the total system inertia matrix is given by
y 2T + z 2T x T y T x Tz T

I = IT + m T x T yT z T2 + x 2T y T z T + TXT I RTX L
x T zT y T z T x 2T + yT2

(3.12)
y 2R + z R2 x R y R xR zR

+ mR x R yR z R2 + x 2R yRzR
x R zR yR z R x 2R + y 2R

The calculation of the spacecraft system inertia matrix is implemented in the SIMULINK
subsystem shown in Figure 7.
22
Figure 7. SIMULINK Subsystem Diagram: Spacecraft Moment of Inertia Matrix
a. Rate of Change of Spacecraft Inertia Matrix
All components of the system moment of inertia matrix are constant except
TXT I R TX which varies with the relative rotation angle . Therefore the rate of change of
the system inertia matrix is given by
&I = T& T I T + T T I T
& (3.13)
X R X X R X
0 0 0 0 0 0

where T& X = 0 sin()&
cos()& and T& X = 0 sin()&
T
-cos()& . It can be
0 -cos()& sin()& 0 cos()& sin()&
shown that &I as a function of and & is given by
0

& = I& sin()
I R
x'y' I+ Rx'z'
cos() L
Ix'z' sin( ) Ix'y' cos()
R R

IRx'y' sin( ) + IRx'z' cos()
2I Ry'y' sin()cos() + 2I Ry'z' cos 2 () 2IRy'z' sin 2 () 2IRz'z' sin()cos( ) L
IRy'y' sin 2 ( ) Ry'y'
I cos2 ()+ 4IRy'z' sin()cos( ) + I Rz'z'cos 2 () Iz'z'
R
sin 2 ()
I Rx'z' sin( ) Ix'y'
R
cos()

IRy'y' sin 2 ( ) Ry'y'
I cos2 ()+ 4IRy'z' sin()cos()+ Rz'z' I cos2 () z'z'
R
I sin2 ()
) 2IRz'z' sin()cos() 2IRy'y' sin()cos() 2IRy'z' cos2 ( )
2I Ry'z' sin 2 ( +
(3.14)
The SIMULINK subsystem diagram that calculates the rate of change of the spacecraft
inertia matrix based on appendage relative angular orientation and rate is shown in Figure
8.
23
Figure 8. SIMULINK Subsystem Diagram: Rate of Change of Spacecraft Inertia
Matrix
4. Angular Velocity
The spacecraft angular velocity is defined by the angular velocity of the primary
body with respect to the inertial reference frame. This angular velocity vector is given by
x
T = y
v
(3.15)
z
The relative angular velocity between the primary and secondary bodies due to the
rotation about the body x axis is given by
&
rel = 0
v
(3.16)
0
24
The angular velocity of the receive telescope is equal to the spacecraft angular velocity
plus the relative velocity
x + &
R = T + rel = y
v v v
(3.17)
z
5. Angular Momentum
The components of the total spacecraft angular momentums given in Equation

(3.5) can now be defined. The system angular momentum neglecting the momentum due
to the relative motion of the receive telescope and the reaction wheels is given by
v v
HS- = I (3.18)
The angular momentum of the receive telescope relative to the spacecraft is given by
v v
Hrel = IR rel (3.19)
Substituting Equations (3.18) and (3.19) into Equation (3.5) we get the total system
angular momentum
v v v v
H = I + R I rel+ wH (3.20)
Substituting this relation into Equation (3.6), the spacecraft equation of motion can be
rewritten as
v &vI + v& H+ v ( vI v
I rel + wH)
v v v
M = &I + & I
+ R rel w + R (3.21)
6. Solving for Spacecraft Angular Rate

v v
To solve for the spacecraft angular rate, , we can isolate & in Equation (3.21)
and perform the integration using a computer solver. This however places an
unnecessary burden on the processor to continuously calculate the time derivative of the
v
spacecraft moment of inertia. A simpler method is accomplished by first solving for HS- .
The spacecraft equation of motion is rewritten as
25
v v& v& v& v v
M=H S- + I R rel+ wH+ H (3.22)
v&
which can be solved for H S- .
v& v v& v& v v

H S- = M I R rel H w + H (3.23)
After the integration,
( )
v v v v v v v v v
HS- = M IR &rel &Hw + H dt= ( M+ H ) dt RI rel Hw
v v
(3.24)
t t
v
the angular rate, , is obtained from
v v
= I-1HS- (3.25)
Figure 9 shows the SIMULINK diagram for implements the spacecraft dynamics for the
Bifocal Relay Mirror attitude simulation. The integration is performed using the
Dormand-Prince ode5 solver.
Figure 9. SIMULINK Subsystem Diagram: Spacecraft Dynamics
B. COMMAND
Maneuvering spacecraft often require feed forward command to maintain precise

tracking requirements throughout their maneuvers. The Bifocal Relay Mirror satellite
must maintain tightly controlled ground tracking by both the transmit and receive
26
telescopes in order to perform its mission. For this system the envisioned feed forward
command will include control torques necessary to maintain a dynamic attitude profile,
which includes the relative motion between the primary and secondary bodies, calculated
in the absence of disturbances. The spacecraft control laws will be used to correct for
errors in the calculated command profile and external disturbances based on errors in
spacecraft attitude and angular rate from the attitude determination system.
The estimated spacecraft rotational profile is predetermined prior to the execution

of a maneuver and tracking operation. The profile includes the spacecraft body attitude,
angular rate and angular acceleration as well as the relative angle, rate and acceleration
between the transmit and receive telescopes. The net external torque required to maintain
this profile in the absence of disturbances is fed forward to the control devices. This
torque can be calculated from the spacecraft equation of motion, Equation(3.21),
neglecting the reaction wheel control devices
v v v
M c = &I c + & c I
v v
( v
+ R & Irelc+ c c I
+ R
v
)
Irelc (3.26)
In this equation the subscript c is used to represent feed forward command. The
equivalent SIMULINK subsystem used in the simulation to generate the feed forward
command torque is shown in Figure 10.
Figure 10. SIMULINK Subsystem Diagram: Feed Forward Torque Command
27
C. CONTROL
Primary attitude control of the Bifocal Relay Mirror satellite is accomplished with
reaction wheels. Four reaction wheels are arranged in a pyramid constellation to achieve
redundancy in the event of a single wheel failure. Under normal operations three wheels
are operating while the forth is off line. Control torques commands are calculated in
body coordinates and then distributed to the three operating wheels.
Torque commands are generated from the feed forward command profile plus the
attitude control laws with compensation for the gyroscopic torques generated from the
spinning reaction wheels. The simple partial plus derivative type controller is chosen for
this satellite. The control laws are based on attitude quaternion and angular rate error as
calculated from the outputs of the attitude determination system where the quaternion
error is calculated by Equation (4.16) [Ref. 9]. The three body axis wheel control laws
are given by
& w1 = K q1q E1+1 KE1

H
& w2 = K q2q E2+2 KE2
H (3.27)
& = K q +K
H
w3 q3 E3 3 E3
More exotic control schemes exist but the optimal state estimates provided by the
attitude determination system are simply employed by these control laws.
The control law gains are chosen to minimize steady state errors attitude errors
while ensuring that oscillations induced by the attitude control system do not interfere
with on-orbit structural modes and payload components. To determine the optimal
control gains, attitude control simulations are conducted with a representative
maneuvering profiles and external disturbances. The state inputs used by the controller
during gain adjustment simulations are perfect attitude and rate knowledge. As explained
in Chapter I, the optimal controller can be determined independently from the attitude
determination system since its outputs are based on optimal state estimation.
The reaction wheel control subsystem used in the SIMULINK attitude simulation
model is shown in Figure 11. Saturation and rate limiting is applied to simulate nonlinear
effects in real reaction wheel control systems. The net control torque is determined and
28
then distributed to the individual reaction wheels based on their orientation within the
constellation.
Figure 11. SIMULINK Subsystem Diagram: Reaction Wheel Control
D. DISTURBANCE TORQUES
Representative disturbance torques are simulated in order to observe spacecraft

attitude performance in a realistic space environment.
1. Gravity Gradient Torque
At low altitudes the torque induced by gravity gradient on spacecraft without

matched body moments of inertia can be significant. A model for the gravity gradient
torque on a spacecraft in orbit is given by
v 3 r r
MG = 3 R 0 I R0 (3.28)
R0
r
where R0 the distance from the spacecraft center of mass to the Earths center, R 0 is the
unit vector in that direction given in body coordinates, and I is the total spacecraft
29
moment of inertia matrix. If the direction cosine matrix from the orbit to body frame,
B O
r
C , is known, R 0 is given by
0 c13
r B O c
R 0 = C 0 = 23 (3.29)
1 c33
In the attitude simulation model, the orbital reference frame direction cosine matrix
components are propagated using the spacecraft angular rate relative to the rotating orbit
frame from an initial orientation. The relative orbital rate is determined from the inertial
rate minus the rate of rotation of the orbit frame. The SIMULINK subsystem for
propagating the orbital reference is shown in Figure 12.
Figure 12. SIMULINK Subsystem Diagram: Orbital Reference Propagation

Using the orbital reference coordinates the gravity gradient torque can then be
written
c13 c13
v 3
M G = 3 c23 I c23 (3.30)
R0
c33 c33
Figure 13 shows the SIMULINK subsystem that models the gravity gradient torque based
on the orbital reference frame direction cosine matrix.
Figure 13. SIMULINK Subsystem Diagram: Gravity Gradient Torque Model

30
2. Other External Disturbance Torques
Other disturbance torques include those due to unbalanced solar pressure,

magnetic interactions, and aerodynamic drag effects. These disturbances are not as easily
modeled as gravity gradient but their characteristics and magnitudes are important
considerations in the design of spacecraft attitude control systems. In the simulation
model, secular and slowly varying periodic moments are introduced in each body axis to
account for these unknown disturbances. Worst case magnitudes are chosen based on
orbit profile and spacecraft characteristics to ensure robust control design.
E. MAGNETIC MOMENTS
The magnetic moment imparted on the spacecraft by the earths magnetic field is
v
dependent upon the magnetic field strength, B , and the spacecrafts magnetic dipole
v
vector, m . The magnetic moment is given by
v v v
Mm = m B (3.31)
1. Magnetic Field Model
Highly accurate models of the earths magnetic field have been developed but a
simple dipole model is adequate for the purposes of this simulation. This approximation
assumes a simple dipole magnetic field tilted 11 degrees from the earths spin axis. The
r
earths magnetic field is a function of the earths unit dipole vector, M , the distance from
the center of the earth to the center of mass of the spacecraft, R, and the unit vector in
r
that direction, R . The magnetic field is given by
v K r r r r
( )
B = 3 3 M gR R M
R
(3.32)
where K is the earths magnetic field constant equal to 7.943 x 1015 Nm2/a2. Using
classical orbital elements it can be shown that the components of the magnetic field
vector in orbital coordinates are given by
31
cos( ) ( cos()sin( i) sin()cos( )i cos( u) ) sin( )sin ()sin( u)
v K
BO = 3 cos()cos( i) sin()sin( )i cos( u)
R
2sin() ( cos()sin( i) sin()cos( )i cos( u) ) +2cos( ) s in()sin( u)
(3.33)
where is the true anomaly of the spacecraft, is the magnetic dipole tilt angle, i is the
orbit inclination, and u is the angle of the magnetic dipole normal to the line of ascending
nodes. The SIMULINK subsystem that propagates the earth magnetic field vector with
the orbital reference is shown in Figure 14.
Figure 14. SIMULINK Subsystem Diagram: Earth Magnetic Field Vector in Orbital
Coordinates
The magnetic field vector can be transformed to the spacecraft body coordinates
with the direction cosine matrix.
v v
B B = B C OB O (3.34)
B
where CO is the transformation matrix from the orbit coordinate frame to the body
frame. Figure 15 shows the subsystem that realizes this transformation using direction
cosine matrix components.
32
Figure 15. SIMULINK Subsystem Diagram: Magnetic Field Conversion to Body
Coordinates
2. Magnetic Control Torque
Using magnetic torquers a spacecraft magnetic dipole can be generated to react

with the earths magnetic field to produce a control torque. In this simulation magnetic
control can be used to help desaturate the reaction wheels. The control law for magnetic
dumping of reaction wheel momentum is given by
v v v
M cmd = k ( B B H w ) (3.35)
where k is the magnetic torquer gain. Figure 16 shows the SIMULINK subsystem that
simulates the magnetic dumping control torque. Saturation is added to simulate nonlinear
effects in the torque rods.
Figure 16. SIMULINK Subsystem Diagram: Magnetic Torquers for Momentum

Dumping
Magnetic torquers can also be used for attitude control. The control laws for the
v
torquers to produce a desired control moment, M CD , on the spacecraft is given by
v v
v B B M CD
M cm = v 2 (3.36)
BB
33
34
IV. ATTITUDE REPRESENTATION AND KINEMATICS
There are many ways of describing the orientation of one coordinate system
relative to another. The most common descriptors used in spacecraft attitude
determination include Euler angles, direction cosine matricies and quaternions, also
known as Euler parameters [Ref. 10]. Euler angles provide a convenient way to represent
attitude and are usually the easiest to visualize. However, singularities arise when the
relative orientation from the reference coordinate system becomes large. Therefore,
highly maneuverable spacecraft require other means of attitude representation. Direction
cosine matricies and quaternions overcome this problem. Direction cosine matricies
provide the most convenient way of transforming vectors between coordinate systems but
require significantly higher attitude processor throughputs than quaternions. As an
example, propagating an attitude matrix with angular rate data requires the integration of
nine elements while the quaternion has only four. For these reasons, the quaternion
representation is chosen for this attitude simulation model.
A. QUATERNION DEFINITION AND CHARACTERISTICS

v
The four vector quaternion based representation q R 4 can be written as
q1

v q
q = 2 (4.1)
q3

q 4
or equivalently,
v
q = q1i + q 2 j + q 3k + q 4 (4.2)
where the quaternion has a three vector imaginary part
q1
v
q I = q 2 and a scalar real part q R = q 4 .
q 3
35
1. Meaning of Quaterions
Quaternions represent the angular orientation of a body relative to a reference

coordinate frame by a single axis rotation of magnitude about the eigenvector axis
v
given by corresponding to the +1 eigenvalue of the direction cosine or attitude matrix.
The direction cosine matrix can be written in terms of these four parameters
1 0 0 0 - 3 2
vv
A = 0 1 0 cos() + 3
0 - 1 sin() + T (1 cos()) (4.3) (4.4)
0 0 1 - 2 1 0
Similarly these parameters can be represented by the four quaternions
q1 1
v
q I = q 2 = sin
2 2

q3 3
(4.5)

q R = q 4 = cos
2
which have the added property that

v
q = q12 + q22 + q32 + q42 = 1 (4.6)
2. Attitude Matrix
The equivalent attitude or direction cosine matrix can be generated from

quaternions using
q12 q22 q32 + q 24 2(q1q2 + q3 q4 ) 2(q1 q3 q 2 q4 )

A = 2(q1q 2 q 3 q4 ) q + q q + q
2
1
2
2
2
3
2
4 2(q 2q 3 + q1 q4 ) (4.7)
2(q1q 3 + q 2q4 ) 2(q 2q 3 q1 q4 ) q12 q 22 + q 32 + q 42

The SIMULINK subsystem used in this simulation model to convert a quaternion to the
equivalent attitude matrix is shown in Figure 17.
36
Figure 17. SIMULINK Subsystem Diagram: Attitude Matrix from Quaternions
3. Quaternion Multiplication
The multiplication of two quaternions to form a third defines an angular

v
orientation resulting from two eigen-axis rotations. If q represents the transformation
v
from coordinate frame A to B and q represents the transformation from frame B to C
then the transformation from frame A to C is given by
v vv
q = qq (4.8)
This multiplication can be implemented several ways. The two quaternions can
be multiplied directly using the quaternion format given in Equation (4.2)
(q1i + q2 j + q3 k + q4 ) = (q1i + q2 j + q3 k + q4 )(q1i + q2 j + q3k + q4 ) (4.9)
using the equalities
i 2 = j 2 = k 2 = 1
ij = ji = k
(4.10)
jk = kj = i
ki = ik = j
Quaternion multiplication can also be performed by treating the imaginary and real parts
separately
37
v v v v v
q I = qR qI + qR qI + qI q I
v v (4.11)
q R = qR qR qI gqI
In matrix form quaternion multiplication is given by
q1 q4 -q3 q2 q1 q1
q q q4 -q1 q2 q2
2 = 3 (4.12)
q 3 -q 2 q 1 q 4 q3 q3

q 4 -q 1 -q 2 -q 3 q4 q4
This multiplication method is used in the attitude simulation. The SIMULINK subsystem
is shown in Figure 18.
Figure 18. SIMULINK Subsystem Diagram: Quaternion Multiplication
4. Coordinate Rotations
Sequential axis rotations such as those that define the three Euler angle
representation can be realized by the successive quaternion products. Additionally, if
there are three small simultaneous rotations 1 , 2 , and 3 about the coordinate axes x, y,
and z respectively the resulting quaternion is determined as follows:

v 1 1
Let = ( + + )
2 2 2 1/2
and = 2 , then

1 2 3
3
q1
v v
q I = q 2 = sin and q R = q 4 = cos (4.13)
2 2
q 3
38
5. Quaternion Inverse and Identity
The inverse of a quaternion is its complex conjugate and is analogous to the

transpose of a direction cosine matrix. A quaternion is conjugated by reversing the sign
on the vector part.
q1 -q 1
q -q
v v
If q = 2 then q* = 2 .
q3 -q 3

q 4 q 4
The identity quaternion is found by multiplying any quaternion by its conjugate as

shown below.
0

vv * v * v 0
qq = q q = (4.14)
0

1
6. Quaternion Error
A difference quaternion can be defined between to orientations that are referenced

v
to the same coordinate frame. If q represents the transformation from coordinate frame
v
A to B and q represents the transformation from frame A to C then the transformation
from frame B to C is given by
q4 q 3 -q 2 -q1 q1

v v * v -q3 q4 q 1 -q 2 q2
q=q q = (4.15)
q2 -q1 q4 -q3 q 3

q1 q2 q3 q4 q4
This can be used to calculate the error quaternion where the target orientation is
v v
given by q and the actual spacecraft orientation is q .
v v v
q e = q*q (4.16)
39
The SIMULINK subsystem diagram for determining the error quaternion is given in
Figure 19.
Figure 19. SIMULINK Subsystem Diagram: Error Quaternion
The angular difference in radians between the original quaternions is contained in

the real part of the error quaternion.
e = 2cos 1 (q 4e ) (4.17)
7. Vector Transformations
Transforming vectors between coordinate systems requires two quaternion

v
multiplications where the vector is treated as quaternion with a real part of zero. If q
represents the orientation of coordinate frame B with respect to reference frame A then a
v
vector given in coordinates of the reference frame, v A , can be transformed to the
coordinates of frame B by
v vv v
v B = q * vA q (4.18)
B. QUATERNION KINEMATICS
Actual spacecraft motion is simulated using the continuous quaternion differential

equations. In the attitude determination computer simulation the attitude propagator uses
the discrete kinematic equations.
40
1. Continuous Kinematics
The differential equation for the quaternions of a rotating coordinate system can
be found by differentiating this equation with respect to a fixed reference frame A. If the
rotational rates for coordinate frame B are given by
x v v v
v dv A dv B dv B v v
= y , then = = 0 and = vB
dt I dt dt I
z B
The inertial derivative of equation (4.18) becomes

v v
v v dq* v v v * v dq
vB = vA q + q vA (4.19)
dt dt
v
Using the properties of quaternions, substituting the original equation for v B and
v
noting that the vector v A is arbitrary it can be shown that
v
dq 1 v v
= q (4.20)
dt 2
or in matrix differential form
q& 1 q1 0 -
3 2 1 q1
v q& q - 1 2 q 2
dq 2 1 v 2 1 3 0
= = S( ) = (4.21)
dt q& 3 2 q 3 2 2 - 1 0 3 q 3

q& 4 q 4 - 1 - 2 - 3 0 q 4
v v
v -S( ) v
where S ( ) = v T and S( ) is the skew symmetric matrix associated with the
- 0
v
vector . The kinematics subsystem diagram used in the SIMULINK model is shown in
Figure 20. The integration is performed using the Dormand-Prince ode5 solver.
41
Figure 20. SIMULINK Subsystem Diagram: Continuous Kinematics
2. Discrete Kinematics
v
The solution to the continuous differential equation when is held constant is
v 2 S()t v t v sin( ) q v
1 v
q=e q 0 = cos( ) I 4x 4 + S (4.22)
2 ( )
0

t
where = . This solution leads to the discrete implementation of the quaternion
2
kinematic equations.
The attitude determination system uses angular rate information to kinematically

propagate the spacecraft attitude quaternion in discrete time steps. The angular rate
v
vector R 3 is integrated with time step t to produce the incremental angle vector
v v
R 3 . For small time steps approximates an eigen-axis rotation in the current
body coordinate frame so it can be related to the change in the attitude quaternion by
v 1i + 2 j + 3k
q sin +cos (4.23)
2 2
The SIMULINK subsystem used in the discrete attitude propagator to obtain the
incremental quaternion step is shown in Figure 21.
42
Figure 21. SIMULINK Subsystem Diagram: Discrete Attitude Propagator
The new attitude quaternion is simply determined by

v v v
q new =q old ( q) (4.24)
where the equation implies quaternion multiplication.
43
44
V. KALMAN FILTER
The Kalman filter provides a non-deterministic means of estimating a system state

vector using a state space mathematical plant model and sensor measurement data related
to some subset of the state variables. It is a stochastic optimal estimator designed to
minimize the weighted mean square error in the state estimate. The Kalman gain matrix
determines the weighting based on the relative confidence between the past state estimate
propagated to the current time and the current partial measurement of state variables.
The error covariance matrix, the second statistical moment of the state vector, tracks the
confidence in the state estimate while a measurement error covariance matrix relates the
confidence in the measurement.
A. RECURSIVE DISCRETE KALMAN FILTER
The Kalman filter can be implemented recursively since all of the information
from past state measurements is encapsulated in the previous state estimate and error
covariance matrix with are both tracked. During state propagation the error covariance is
updated to reflect additional error added by imperfections in the plant model. A recursive
discrete Kalman filter is used in the proposed attitude determination scheme. Figue 22
illustrates the recursive nature of the Kalman filter.
45
ENTER RECEIVE
INITIAL MEASUREMENT
ESTIMATE UPDATE
AND ERROR
COVARIANCE
CORRECTION
PROPAGATE
COMPUTE
ESTIMATE AND ERROR
COVARIANCE TO NEXT KALMAN GAIN
TIME STEP
PREDICTION COMPUTE NEW

ESTIMATE AND
ERROR
COVARIANCE
Figure 22. Discrete Kalman Filter Loop
The standard equations for a recursive discrete Kalman filter are summarized in
Table 2 below. For full background development and equation derivations see Appendix
A [Ref. 1].
v v v
System/plant Model xk = x k-1 k-1w + k-1
v v
w k ~ N(0,Qk )
v v v
Measurement Model zk = H kx k + vk
v v
v k ~ N(0,R k )
v v
Initial Conditions x 0 = E[x 0 ]
v v v v
P0 = E[(x 0 x 0 )(x 0 x 0 ) T ]
v w
Assumptions E[w kv j ] = 0
(uncorrelated errors) for all j, k
Prediction v - v+
State Estimate Extrapolation xk = x k-1 k-1
k = P k-1 k-1 k-1
Q+
- + T
Error Covariance Extrapolation P k-1
Correction Kalman Gain Matrix K k = Pk- HTk (HkPk- HTk + R k )-1

v v v v
State Estimate Update x +k = x -k + K k (z k H k x k- )
Error Covariance Update Pk+ = (I K k H k )Pk-
Table 2. Discrete Kalman Filter Equations
46
B. ERROR STATE EXTENDED KALMAN FILTER IMPLEMENTATION
A linear error state extended discrete Kalman filter is implemented in the

simulated attitude determination model to estimate spacecraft attitude and angular rate.
The nonlinear attitude propagation is performed discretely outside of the Kalman filter
according to Equations (4.23) and (4.24) at the high bandwidth frequency of the attitude
processor. The angular rate estimate used in the kinematic model is provided by the
spacecraft gyros or pseudo gyro rate calculations. The Kalman filter is designed to
provide bias corrections to the gyro outputs for attitude propagation [Ref.11]. If the
pseudo gyro rate is used, the bias error is treated as a rate error from which a correction to
the system angular momentum is determined. Measurement updates are provided by star
trackers. The difference in the measured and predicted star vector is related to the
attitude error and used to provide corrections to the state estimate.
1. State Variables
The Kalman filter used in this model estimates six state variables: a three vector
of attitude errors, % R 3 , and a three vector of gyro bias errors, b% R 3 . The total state
vector is given by
v %
x= R 6 (5.1)
b%
The attitude error, % , represents the deviation in the spacecraft attitude relative to
the inertial reference frame given by a vector of three simultaneous rotations. The bias
error, b% , represents the change in bias of the angular rate data in the spacecraft body
coordinate frame.
2. Attitude Propagation Error Correction Methods
There are two ways to implement the Kalman filter corrections to the attitude
propagator. In the first method, the attitude propagator is fed with raw gyro rate
information and the Kalman filter maintains the total gyro bias tracking with the state
bias error, b% . In this method, the attitude quaternion must be corrected at each time step
47
by the filter. The attitude error vector, % , is converted to an incremental quaternion
v
rotation, q and applied to the propagated attitude quaternion as shown in Chapter IV
v
with Equations (4.23) and (4.24) or, using the small angle approximation q is given by
v % % %
q= 1 i + 2 j + 3 k + 1 (5.2)
2 2 2
Then by quaternion multiplication the propagated attitude quaternion is updated
v vv
q new =(q)qold (5.3)
The alternative method of applying the filter correction to the propagator is to use
v
the bias error state, b% , to correct the unknown gyro bias, b each time a measurement
update is obtained
v v
bnew = bold + b% (5.4)
For this method, the total unknown gyro bias is tracked separately from the filter and
added to the gyro angular rate before it is fed to the attitude propagator.
v
= ( g + b ) t
v v
(5.5)
The attitude error state correction is then only necessary at measurement updates. Using
v
this method allows the Kalman filter error vector x to be reset after each bias
measurement update. Since the filter approximates nonlinear errors with a linear model,
keeping the errors as close to zero as possible improves the estimate. This Kalman filter
correction method is used in the simulated attitude determination model.
3. Plant Model
A nominal plant model is chosen for this Kalman filter which assumes a constant
rate bias between filter updates and uncoupled (linear) attitude errors between
propagation steps. It is derived from the continuous equation of state
v
dx d% dt 03X3 A T % v
= = =f (t)x (5.6)
dt db% dt 03X3 03X3 b%
where A is the estimate of the spacecraft attitude or direction cosine matrix. This linear
48
state equation includes coupled rotational effects through the time variant estimated
attitude matrix which changes dynamically.
4. State Transition Matrix
For the discrete Kalman filter implementation of this model the state transition
matrix, k , is calculated for each time step using the matrix exponential of f(t) from the
continuous state equation
I T
At
k = e f (t)t = 3 X 3 k
(5.7)
03 X 3 I3 X 3
where ATk is the transpose of the estimated attitude matrix at the current discrete time, tk
and t is the interval between steps. The SIMULINK subsystem used to calculate the
state transition matrix is shown in Figure 23. The estimated attitude matrix is calculated
from the propagated attitude quaternion using Equation (4.7).
Figure 23. SIMULINK Subsystem Diagram: State Transition Matrix
5. Kalman Filter Prediction Equations

v
The state transition matrix allows the state estimate, x , and associated error
covariance matrix, P, to be propagated forward in the discrete Kalman filter prediction
step
v- v
x k+1 = k k x (5.8)
-
Pk+1 = k kPTk +Q (5.9)
The subscript k+1 indicates next discrete time step, tk+1, and the superscript indicates
predicted future value based only on information up to the current time step tk. The
covariance matrix relates the confidence in the associated state estimate. A larger P
49
indicates less confidence which means that the attitude determination system is assumed
to have larger errors. Measurement updates taken when P is large will greater impact on
the estimate than those taken when the confidence is high. The Kalman filter prediction
step is implemented in the SIMULINK model as shown in Figure 24.
Figure 24. SIMULINK Subsystem Diagram: Kalman Filter Prediction Step
6. Plant Noise Covariance
The plant noise covariance, Q, is a positive definite matrix that characterizes the
plant error accumulated through the time step, t , assuming that it can be modeled as
Gaussian white noise meaning normally distributed with zero mean. Q is normally taken
to be a diagonal matrix meaning that there is no known correlation between the errors of
the six state variables. In terms of the variances of the state variables the plant noise
covariance matrix is given by
1
2
0 0 0 0 0
2
0 20 0 0 0
0 0 23 0 0 0
Q= (5.10)
0 0 0 2b1 0 0
0 0 0 0 2b2 0

0 0 0 0 0 2b3
If the time step is varied or the plant error is known to be significantly changed
due to operating mode, then Q should be modified to reflect that change. Decreasing Q
effectively asserts that the plant generates less error and the predicted state estimate will
50
be weighted higher relative to the measurement updates. On the other hand, if it is
increased less confidence is placed on the estimate predicted by the model and more
significant corrections will be applied during measurement updates. In this simulation, Q
is taken to be constant.
7. Kalman Filter Initialization
In order to initiate the Kalman filter algorithm, an initial state estimate, x i and its
associated covariance, Pi, must be chosen. Least squares batch processing can be
performed on two or more star tracker measurements prior to initiating the filter or an
educated guess can be used. Since the state represents errors from truth, the selected
initial state estimate is
0 100 0 0 0 0 0
0 0 100 0 0 0 0

v 0 0 0 100 0 0 0
x i = with associated covariance Pi =n . The
0 0 0 0 1 0 0
0 0 0 0 0 1 0

0 0 0 0 0 0 1
constant, n, is chosen to reflect confidence in the estimate at initialization. The attitude
errors generated are normally about an order of magnitude up from the rate bias errors so
the variances are weighted higher.
8. Sensor Measurement Update
The state estimate and the covariance matrix propagate forward without
correction through each time step until a measurement update is produced. The state
estimate is used to kinematically propagate the attitude quaternion while the covariance
matrix builds up due to the added error through each step.
a. Measurement Vector
If an attitude measurement update is produced at the next time step, k+1,

by one or more star trackers the Kalman filter correction step is applied. The star trackers
produce horizontal and vertical outputs (H,V) corresponding to the position of the star on
51
the detector array. These outputs are internally compensated for with software and
calibrated to form a measurement vector in tracker coordinates
H v
v v s
sm = V with the normalization sm = vm (5.11)
sm
1
which is the star tracker measurement vector. For simulation modeling purposes, noise is
added to the H and V components to create an artificial star tracker measurement vector.
The measurement vector is generated in the SIMULINK model as shown in Figure 25.
Figure 25. SIMULINK Subsystem Diagram: Star Tracker Measurement
At every integration cycle, t , the attitude estimate is available from the

Kalman filter prediction step. If a star tracker observation has occurred during that time
step, the measurement vector must be propagated ahead to correspond to the current
discrete time of the attitude computer. Otherwise, the state vector and the state transition
and covariance matricies must be interpolated to accomplish the filter update. Tracker
processing latencies and transport delays must also be compensated for in the time
difference. For simplicity in this model, it is assumed that the all observations occur at an
exact discrete time step of the integration cycle. To produce star measurements for the
model, an artificial star tracker reference is chosen with representative noise applied to
the H and V outputs.
b. Predicted Measurement
The predicted vector in star tracker coordinates is needed to determine the

measurement residual for the update. This is generated using the known position of the
star in inertial space. The detected star goes through the identification process and gets
52
compensated for aberration to yield a unit vector, sI, in inertial coordinates. Using the
estimate for the inertial to body attitude matrix generated from the filter state prediction,
A-k+1 , and the calibrated body to star tracker transformation matrix, T, a predicted vector
is generated in tracker coordinates.
v -k+1 svI
sp =T A (5.12)
The inertial star vectors for the model are generated by applying the
transpose of the body to star tracker and true attitude matricies to the same star tracker
reference used to generate the measurement.
H
v T T
sI =A k+1 T V (5.13)
1
Figure 26 illustrates the SIMULINK subsystem used to produce a simulated inertial star
reference vector.
Figure 26. SIMULINK Subsystem Diagram: Inertial Star Vector
c. Measurement Residual
The measurement residual, zk+1, is formed by subtracting the predicted

from measured star vector and considering only the first two components.
v v v
z k+1 =E(sm sp ) (5.14)
1 0 0
where E= . The measurement residual is determined in the SIMULINK model
0 1 0
as shown in Figure 27.
53
Figure 27. SIMULINK Subsystem Diagram: Measurement Residual
9. Observation Matrix
Since the measurement residual is simply the difference between unit vectors, it is
easily related to the state attitude errors. The observation or feedback sensitivity matrix
defined by the relationship to the state variables
- S ( vs ) 0
H k+1 = E T A (5.15)
k+1 I 2x3
v
where S(sI ) indicates the skew symmetric matrix associated with the inertial star
reference vector. Figure 28 shows the subsystem diagram used in the SIMULINK model
to produce the observation matrix.
Figure 28. SIMULINK Subsystem Diagram: Observation Matrix
10. Kalman Gain
The Kalman gain can then be calculated for the correction step at k+1 from the
standard discrete filter equation
-
K k+1 =Pk+1H Tk+1 (R k+1 +H k+1Pk+1
- T
Hk+1 )-1 (5.16)
where R R 2 X 2 is the measurement noise covariance matrix associated with the assumed
Gaussian white noise in the H and V outputs of the star trackers. Since these errors are
54
usually similar and uncorrelated, the R matrix is normally taken to be a multiple of the
identity.
2H 0
R= 2
0 V
Noisier or less accurate star trackers will generally have a larger R matrix which
decreases the Kalman gain and thus provides less of a correction to the state and
covariance matricies. In reality, the measurement error may not be constant for each
correction step even when the same star tracker is used. Higher intensity stars, although
easier to detect and identify, produce slightly larger noise. Also stars detected toward the
field of view limits of the tracker usually have larger errors due to distortions than those
detected near the star tracker optical axis.
The Kalman gain SIMULINK subsystem is shown in Figure 29.
Figure 29. SIMULINK Subsystem Diagram: Kalman Gain
11. Kalman Filter Correction Equations
With the Kalman gain, the state and covariance matricies can be corrected with
the measurement update
v v
x= Kk+1z k+1 (5.17)
v v - v
x k+1 =xx+
k+1 (5.18)
and
Pk+1 =(I 6X6 K k+1H k+1 )Pk+1

-
(5.19)
The covariance matrix update is accomplished in the simulation model as shown in
Figure 30.
55
Figure 30. SIMULINK Subsystem Diagram: Covariance Updeate
In general, the estimated state drifts away from the true state as it is propagated
forward with the non-ideal plant model and is corrected back toward truth as
measurement updates occur. Corrections to the attitude components associated with the
body axes closely aligned with the star tracker optical axis will be corrected much less
than those that are aligned perpendicular. The spacecraft attitude error continues to grow
during periods where there are no cataloged stars in the sensor field of view or when the
tracker information is unavailable. Additionally, when only one tracker provides updates
for extended periods the angular orientation about that star trackers axis in body
coordinates goes unchecked.
C. CONTROLLER DESIGN IMPLICATIONS
Since the Kalman filter is an optimal least squares estimator, the development of
an optimal controller can be accomplished independently. Therefore, the attitude
determination output should not affect the controller design. In this model, quaternion
and rate error control laws are used as well as feed forward torque to generate reaction
wheel commands. The optimal gains for the controller as determined by simulation with
ideal deterministic attitude knowledge remain optimal with the attitude determination
based on Kalman filter state estimates.
56
VI. DYNAMIC RATE CALCULATION
The continuous dynamic equations of motion for the Bifocal Relay Mirror
spacecraft are derived in Chapter III. These equations produce the spacecraft angular rate
from external control and disturbance moments applied to the body. A similar discrete
model can be applied in the spacecraft attitude processor software to produce a real time
calculated estimate of the angular rate, referred to as the dynamic gyro. This process is
borrowed from The Aerospace Corporations Pseudo Gyro concept [Ref. 8]. At high
bandwidth processor execution the discretization of the dynamics introduces little error.
The angular rate generated by this method can be used as a substitution for conventional
gyroscope outputs. Attitude determination based on the dynamic gyro can be
implemented as a back up failure mode or a primary operating mode to increase the
expected lifetime of the satellite gyroscopes.
A. DISCRETE EQUATIONS OF MOTION
The discretized equations of motion are derived from

v v
H = M t (6.1)
v
where M is the sum of external moments applied to the spacecraft including controls,
modeled disturbances and gyroscopic stiffness. This allows the total system angular
momentum to be tracked with
v v v
Hk+1 = H k + H (6.2)
Then subtracting the relative momentum of the reaction wheels and secondary body
produces
v v v v
HS- = H Hw H rel (6.3)
The calculated spacecraft angular rate is then given by
v v
= I-1HS- (6.4)
57
B. MOMENTUM CORRECTION FROM KALMAN FILTER UPDATES
The accuracy of the dynamic angular rate calculation ultimately depends on the
tracking of the system angular momentum. If uncorrected, the numerous error sources in
the model will cause the angular rate error to grow over time. The error state Kalman
filter designed for gyro-based attitude determination systems can be used to provide the
necessary model corrections. The gyro bias error states, b% , are interpreted as spacecraft
body rate errors. Using the calculated spacecraft inertia matrix, I, a correction to the
system moment of inertia can be generated by
v
Hcorr = I b% (6.5)
The Kalman filter momentum correction is applied as if the error in the dynamic gyro is
attributable to the total spacecraft body. The relative momentum terms from the
secondary body and the reaction wheels are treated as if they are without error.
The SIMULINK dynamic gyro subsystem used in the attitude simulation model is
illustrated in Figure 31.
Figure 31. SIMULINK Subsystem Diagram: Dynamic Gyro
C. INPUTS AND ERROR SOURCES
After initial calibration, kinematic plant error in gyro-based attitude determination

systems is almost entirely attributable to a single set of imperfect gyroscope rate sensors.
As long as gyro data does not become erratic a Kalman estimator based on a slowly
58
changing rate bias plant model produces an effective attitude determination system even
with relatively noisy rate inputs.
The error in rate calculations from dynamic modeling, on the other hand, is due to
numerous factors and is much harder to characterize. There are multiple internal sensors
involved as well as dynamic modeling simplifications. Since the dynamic calculation is
produced from total system momentum tracking any error in knowledge of external
torques directly correlates to rate error. Errors in system or component moments of
inertia have the same effect. Like gyro outputs, internal position and rate sensor data are
corrupted by measurement and alignment errors. These data from all moving appendages
and momentum exchange devices are critical to the accuracy of the rate calculation.
Satellites not designed to use dynamic rate calculations for attitude determination are
usually not equipped with appendage relative rate measurement sensors. These data must
either be derived adding more error to the calculation or substituted with commanded
rates. It is important that all known biases be removed from sensor data and calculated
input errors since the Kalman filter estimator is based on the assumption of uncorrelated
zero-mean Gaussian noise. Even if all input parameters were known exactly the discrete
modeling of the spacecraft dynamics introduces some error.
1. External Control and Disturbance Torques
In the dynamic gyro, known externally applied moments are integrated in the
system angular momentum calculations. These include control moments other than those
imparted by momentum exchange devices as well as modeled disturbance torques. Since
control torques are normally of significant magnitudes, it is essential that they be
modeled correctly. Moments from magnetic torque rods depend on the imperfectly
controlled magnetic dipole and the earths magnetic field strength. The magnetic dipole
is provided by torquer current measurement. The earths magnetic field must either be
modeled or measured with magnetometer but is not precisely known. Reaction jet
moments are almost impossible to model accurately. This often result in degraded
pointing and tracking during firings operations.
The effects of external disturbance moments depends on spacecraft configuration

and orbital profile. Those disturbances that have significant effects on vehicle dynamics
59
should be modeled whenever possible to increase the accuracy of the attitude control
system. These types of errors are non-zero mean in the short term and therefore can only
be corrected for with sensor measurement updates. An extended option for The
Aerospace Corporations Pseudo Gyro includes a torque bias estimator to reduce the
effects of unknown external disturbances.
For the Bifocal Relay Mirror satellite, the gravity gradient torque is the most
significant disturbance and can be modeled as an input to increase the accuracy of the
dynamic rate calculation. The model is well understood and is generated from the inertia
matrix which is already required by the dynamic gyro and vehicle orientation with
respect to the gravity vector which can be determined from the estimated attitude and
ephemeris data. The gravity gradient model used by the dynamic gyro is equivalent to
the subsystem shown in Figure 13.
2. Reaction Wheel Relative Momentum
Instead of integrating the torques produced by momentum exchange devices their

relative momentum effects are used directly in the dynamic rate calculation. The relative
momentum of each reaction wheel is given by its orientation within the spacecraft, the
component inertia of its spinning disk and the wheel spin rate. The imperfect sensor
measurements from the reaction wheel tachometers introduce errors in system
momentum calculation. Relative orientation angles of reaction wheels are fixed and
errors can be corrected through calibration. Orientations of control moment gyros,
however, are variable. Since these devices usually carry much more momentum small
gimbal resolver errors can have a significant impact on total system momentum
calculations. In this simulation, time varying artificial alignment errors are applied to the
reaction wheel momentum measurements to observe these effects. The error corrupted
wheel momentum measurement subsystem implemented in the SIMULINK model is
shown in Figure 32.
60
Figure 32. SIMULINK Subsystem Diagram: Error Corrupted Reaction Wheel
Momentum Measurement
3. Moment of Inertia Calculations
The calculation of the total spacecraft inertia matrix is accomplished by Equation

(3.12) using the SIMULINK subsystem illustrated in Figure 7. Fixed component
moments of inertia and masses of the primary and secondary bodies are assumed to be
known as well as the relative positions of their centers of mass from the system mass
center. Imperfect knowledge of these parameters introduces errors in the angular rate
calculation. Errors generated from rotating spacecraft components are not constant in the
spacecraft body frame.
The inertia matrix also depends on internal sensor input from position encoders or
potentiometers for relative angular orientation of appendages. The model of the
potentiometer that measures the relative angle of the receive telescope includes
quantization effects and additive noise. If appendage relative motion is slow or
component moments of inertia are small, it may not be necessary to update the system
inertia dyadic at the bandwidth of the attitude processor. A trigger is added to the
SIMULINK subsystem that calculates the inertia matrix (Figure 7), so that the affects of
update rate can be evaluated.
4. Appendage Relative Momentum
Knowledge of appendage relative momentums has direct bearing on the system

momentum and therefore the dynamic rate calculation. The relative momentum of the
Bifocal Relay Mirror satellites secondary body is calculated from Equation (3.19). It
61
depends on the knowledge of the fixed secondary body inertia matrix and the relative
angular rate. The attitude simulation model assumes that there is no directed
measurement of the relative angular rate. The rate must therefore be derived from
potentiometer measurements of the relative orientation about the axis of rotation. This
sensor may have a minimum discernable incremental angle and noise corruption. Also,
since the rate is derived from position measurements it exhibits increased noise and time
lag. The simulation diagram that produces the error corrupted relative angle and rate for
the spacecraft secondary body is shown in Figure 33.
Figure 33. SIMULINK Subsystem Diagram: Appendage Relative Angle and Rate
Measurement
Appendages of significant inertia or relative rates like the receive telescope of the
Bifocal Relay Mirror satellite should be controlled with smooth gimbal drive motors in
order to minimize nonlinear relative pointing errors. Ideally, all moveable appendages
would have an associated rate sensor for each axis of rotation.
The other option for approximating appendage relative angular rates is to use the
commanded drive input. In some systems, appendage controllers can provide a smooth
relative rate through the drive motor, which significantly enhances angular momentum
tracking. Drive actuators with considerably erratic friction effects will cause errors in the
attitude control system during slew operations. Short term transients may introduce
significant settling times for error corrections.
The subsystem of the attitude determination simulation that controls all of the
internal sensor measurements and input parameter calculations for the dynamic gyro is
shown in Figure 34.
62
Figure 34. SIMULINK Subsystem Diagram: Dynamic Gyro Inputs
D. CALCULATED ANGULAR RATE ERROR CORRECTION
The error state extended Kalman filter is designed to track rate bias errors that are
relatively constant in the spacecraft body frame. This is a good approximation for
properly functioning gyro based systems. For the filter to remain effective when used
with dynamic gyro, the bias in the output rate as seen in body coordinates due to the
various error sources must be small and exhibit a bandwidth below the measurement
update rate. Since it is the case that the common error sources are not zero-mean, their
effects on the rate output must be relatively constant so that the momentum correction
supplied by the Kalman filter applies over the update interval. Transient error spiking
can require multiple star tracker measurement updates to correct. If the internal sensors
provide inconsistent measurements, dynamic gyro angular rate will be degraded.
63
64
VII. RESULTS
Simulation results demonstrate that multi-body spacecraft attitude control without

the use of rate gyroscopes can be performed using the attitude and angular rate estimation
scheme proposed in this thesis. The performance of dynamic gyro based attitude
determination system is compared with a similar gyro based system using the Bifocal
Relay Mirror attitude simulation. Simulation input parameters are varied to analyze the
effects of major error sources on the dynamic gyro model. Also evaluated are the effects
of star tracker accuracy and measurement update rates on the attitude determination
system.
A. BASELINE SIMULATION
A full set of MATLAB plots are presented in this section for baseline analysis of
the dynamic gyro based attitude determination and control scheme developed in this
thesis. This set of results validates the potential effectiveness of the proposed attitude
and angular rate estimating scheme used for multi-body spacecraft control. It also
provides a common reference for analysis and comparison with subsequent simulation
results. For other simulations, only selected plots that are required for analysis of results
will be shown.
1. Simulation Input Parameters
Table 3 shows the inputs parameters that are held constant for each simulation run
used to obtain results. These input parameters are set in the MATLAB script file that
calls the SIMULINK attitude simulation. The MATLAB code file is included as
Appendix B.
SIMULATION PARAMETERS (SIMULINK)
Simulation Time Period 500 sec
Attitude Determination Bandwidth 20 Hz
SIMULINK Solver Method ode5 (Dormand-Prince)
Solver Fixed Step Size 0.05 sec
ORBIT PARAMETERS (Circular Orbit)
Altitude 715 km
65
COMMANDED MANEUVER PROFILE
Inertial Attitude Quaternions See Figure 35
Body Axes Angular Rates See Figure 36
Secondary Body Relative Angle See Figure 37
DISTURBANCE MOMENTS
Gravity Gradient Modeled
Secular (Magnitude) 1e-4 Nm
Periodic (Magnitude, Period) 4e-4 Nm, [400,500,600] sec
Disturbance Effect on System Angular See Figure 38.
Momentum
MOMENTS OF INERTIA MATRICIES
Primary Body [2997.28025,-
3.9331,118.2824
-3.9331,3164.18285,1.1230
118.2824,1.1230,881.82105]
kgm2
Secondary Body [1721.07340,-0.0116,-
7.8530
-0.0116,1559.85414,-
12.5463
-7.8530,-
12.5463,182.89962] kgm2
CENTER OF MASS OFFSET FROM SPACECRAFT C.M.
Primary Body [0.558354158,
3.91788e-4,
0.15226902] m
Secondary Body [-1.302113918,
-9.13673e-4,
0.355100092] m
WHEEL CONTROL LAWS
Quaternion Error Gains (Kq) [3000,7000,4500]
Angular Rate Error Gains (Kw) [1000,2000,1000]
Control Law Delay for Initial 30 sec
Determination Errors
GYROSCOPE CHARACTERISTICS
Static Rate Biases 1e-4*[-1,1.5,1] rad/sec
Rate Noise Variance 1e-8
Acceleration Noise Variance (Rate 1e-12
Random Walk)
3 Gyro Alignment Aligned to Body Axes
ATTITUDE DETERMINATION INITIALIZATION ERRORS
Quaternion Errors (q1,q2,q3) [0.008,0.012,-0.008]
Angular Rate Errors [-0.001,0.001,0.002] rad/sec
KALMAN FILTER INITIALIZATION
State Estimate [0,0,0,0,0,0]
66
State Error Covariance Matrix 5e2*[100,0,0,0,0,0
0,100,0,0,0,0
0,0,100,0,0,0
0,0,0,1,0,0,0
0,0,0,0,1,0,0
0,0,0,0,0,1,0
0,0,0,0,0,0,1]
STAR TRACKER ALIGNMENT TO BODY AXES
Tracker 1 (x-rotation, y-rotation) 135 deg, 30 deg
Tracker 2 (x-rotation, y-rotation) 135 deg, -30 deg
Tracker 3 (x-rotation, y-rotation) 180 deg, 0 deg
REACTION WHEEL PARAMETERS (4 WHEEL PYRAMID
CONSTELLATION)
Number of Wheels Operating 3
Constellation Angle to xy-Plane 45 deg
Constellation Torque Saturation Limit 1 Nm
Constellation Torque Rate limits 10 Nm/sec
Table 3. Simulation Input Parameters
2. Command Attitude Profile
The 500 second maneuvering profile chosen for the Bifocal Relay Mirror attitude
simulation analysis is illustrated in Figures 35, 36 and 37. This profile resembles the
maneuver required to maintain transmit and receive telescope pointing control during an
overhead operational pass to conduct laser relay operations. The majority of the
maneuver is performed in the spacecraft pitch axis, q2, as both telescopes orient to point
at fixed ground sites. Less significant motion is required in the spacecraft roll and yaw
axes in order to ensure that the relative axis of rotation of the receive telescope is
correctly oriented during the tracking maneuver. Based on ground site separation
distance and orbital altitude, the largest relative angle require between the telescopes is
about 30 degrees during a near overhead pass between the uplink and downlink ground
sites.
67
Commanded Quaternions
1
q4
0.8 q1
q2
q3
0.6 q4
q2
0.4
0.2
q1
0
q3
-0.2
-0.4
0 100 200 300 400 500
Time (sec)
Figure 35. Baseline Commanded Attitude Profile
-3 Commanded Angular Rates

x 10
12
10 w1
w2
w2
w3
8
6
rad/s
2
w1
w3
-2
-4
0 100 200 300 400 500
Time (sec)
Figure 36. Baseline Commanded Angular Rate Profile
68
Relative Angle to Recieve Telescope (alpha)
35
30
25
20
Deg
15
10
0
0 100 200 300 400 500
Time (sec)
Figure 37. Baseline Commanded Relative Angle to Receive Telescope Profile
The magnitude of the total spacecraft angular momentum during the maneuvering
profile is plotted in Figure 38. The spacecraft attitude is completely controlled by
momentum exchange with the reaction wheels. Therefore, the changes to the angular
momentum profile are attributable the external disturbance moments. The gravity
gradient disturbance is modeled and the other disturbances are fixed so the angular
momentum profile remains essentially common in all simulation runs.
System Angular Momentum
8.1
8.05
7.95
7.9
Nms
7.85
7.8
7.75
7.7
7.65
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 38. Total Spacecraft Angular Momentum Profile
69
3. Attitude Determination System Performance Results
In this baseline simulation, updates from randomly selected star trackers are
provided to the dynamic gyro and the attitude propagator via the Kalman filter at two
second intervals. The star tracker H and V measurements are corrupted with noise
variance of 1x10-4. In order to observe the performance of the attitude determination
system without updates, a 200 second star gap is simulated starting 100 seconds into the
run. As a worst-case analysis, this star gap occurs during the peak maneuvering time of
the satellite including the rotation of the secondary body. Attempt is made to tailor the
plant and measurement error covariance matricies used in the Kalman filter. The attitude
determination system is initiated with the errors given in Table 3.
The accuracy of the angular rate calculation is entirely dependent upon the ability
of the dynamic gyro to track the total spacecraft angular momentum. The error in the
magnitude of the total system momentum compared to the simulated actual momentum is
shown in Figure 39. The steady state momentum error is held within 0.07 Nms with
consistent star tracker data but builds to 0.35 Nms after 200 seconds without stars. No
external disturbance torques are modeled as dynamic gyro inputs in this simulation run.
Angular Momentum Error
0.4
0.3
0.2
0.1
Nms
-0.1
-0.2
-0.3
-0.4
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 39. Baseline Dynamic Gyro Angular Momentum Error
After star tracker measurements are processed, the Kalman filter provides a
momentum correction to the dynamic gyro determined from the spacecraft inertia and the
70
rate bias error state estimate. The attitude quaternion is also updated and the states are
reset to zero. Figure 40 shows the magnitude of the error states determined by the
Kalman filter. No updates are provided during the 200 second star gap.
-4
x 10 State Attitude Errors
xa1
5 xa2
rad xa3
-5
0 50 100 150 200 250 300 350 400 450 500

Time (sec)
-5
x 10 State Rate "Bias" Errors
xb1
5 xb2
xb3
rad/s
-5
0 50 100 150 200 250 300 350 400 450 500

Time (sec)
Figure 40. Baseline Kalman Filter Attitude and Rate Bias Errors
Quaternion errors in the attitude determination system are plotted in Figure 41.
Steady state attitude errors are maintained within 3x10-4 during periods of continuous star
coverage but are increased an order of magnitude by the end of the 200 second star gap.
After the star gap, the attitude error build up is quickly removed through measurement
corrections.
71
-3
x 10 Quaternion Errors
3
q1
q2
2 q3
q4
-1
-2
-3
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 41. Baseline Estimated Attitude Quaternion Error
Figure 42 illustrates the nature of the angular rate estimation error produced by
the dynamic gyro based determination system. The error generated in the roll axis, w1, is
the most significant because it is aligned with the axis of rotation of the secondary body.
The relative rate is derived from imperfect potentiometer measurements. The rate errors
increase during receive telescope motion due to the potentiometer quatization effect
which produces the broken pattern of noise. The yaw axis is most coupled dynamically
to roll axis and exhibits similar errors at less magnitude. The build up of the bias in the
rate error is hard to perceive among the noise but the effects are evident in the quaternion
error plot.
72
-3
x 10 Rate Errors
1.5
w1
w2
1 w3
0.5
rad/s 0
-0.5
-1
-1.5
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 42. Baseline Estimated Angular Rate Error
4. Attitude Control System Performance Results
Attitude control is accomplished by three reaction wheels selected from a pyramid

constellation of four. Nonlinear response of these control actuators is simulated by
saturation and rate limiting the output torques as shown in Figure 11. Wheel torque
commands are generated by combining feed forward and error based control laws. The
control laws implementation is delayed at initiation to allow the attitude determination
system to converge. High gain control laws dominate the wheel torque response. Figures
43 and 44 show the torque and momentum response of the three operating reaction
wheels.
73
Wheel Torques
2
Wheel 1
Wheel 2
1.5
Wheel 3
0.5
Nm 0
-0.5
-1
-1.5
-2
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 43. Reaction Wheel Control Torques
Wheel Angular Momentums

100
Wheel 1
80 Wheel 2
Wheel 3
60
40
20
Nms
-20
-40
-60
-80
-100
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 44. Reaction Wheel Angular Momentum
Attitude control performance is dominated by attitude determination errors. The

control system is designed to limit steady state errors in the attitude quaternions to less
than 2x10-7 with ideal attitude and angular rate knowledge. With the control laws
applied, the attitude control response resembles the errors of the attitude determination
system in steady state. Figure 45 shows the attitude quaternion error response for the
baseline simulation run. Damped corrections back to steady state after the star gap are
much slower than those realized by the attitude determination system.
74
-3
x 10 Quaternion Error
4
q1
q2
3
q3
q4
2
-1
-2
-3
-4
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 45. Baseline Control Attitude Quaternion Error
The error in the controlled spacecraft angular rate is shown in Figure 46. Since
most of the attitude determination noise is removed by the controller, the uncorrected rate
error bias due to the star gap is evident. Although the errors remain small during the star
gap, the bias drives the attitude quaternion error build up.
-3
x 10 Rate Error
1.5
w1
w2
1 w3
0.5
rad/s
-0.5
-1
-1.5
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 46. Baseline Control Angular Rate Error
B. DYNAMIC GYRO VS GYRO PERFORMANCE
Results presented in this section compare a gyro based attitude determination

system with the dynamic gyro based system. Three orthogonally mounted gyroscopes
75
are modeled to provide simulated spacecraft body angular rate measurements to the
attitude propagator when the gyro option is selected. The gyro model includes static bias,
rate noise, and rate random walk supplied by integration of random noise. Gyroscope
characteristics used for the simulation results are listed in Table 3. A bias error tracking
system is also modeled in parallel with the gyro generated rates to accept bias corrections
provided by the Kalman filter. Direct comparison of the dynamic gyro and gyro based
attitude determination systems is accomplished using identical input parameters to obtain
simulation results. The SIMULINK subsystem for gyro rate measurement and bias error
tracking is illustrated in Figure 47.
Figure 47. SIMULINK Subsystem Diagram: Gyro Measurement and Bias Error
Tracker.
1. Continuous Star Tracker Coverage
With precise attitude updates from star trackers there is only little noticeable
difference between the dynamic gyro and gyro based attitude determination systems.
Figure 48 provides a side-by-side comparison of simulation results obtained with a
continuous star update interval of 2 seconds.
76
-4 -4
DG: Quaternion Errors Gyro: Quaternion Errors
x 10 x 10
4 4
q1 q1
3 q2 3 q2
q3 q3
2 q4 2 q4
1 1
0 0
-1 -1
-2 -2
-3 -3
-4 -4
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Time (sec) Time (sec)
-3 -3
DG: Rate Errors Gyro: Rate Errors
x 10 x 10
1.5 1.5
w1 w1
w2 w2
1 1
w3 w3
0.5 0.5
rad/s
rad/s
0 0
-0.5 -0.5
-1 -1
-1.5 -1.5
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 48. Results: Dynamic Gyro vs. Gyro with Continuous Star Coverage
Although the attitude quaternion errors show little difference between the two
systems the rate errors from the gyros are noticeably more accurate and better resemble
white noise. Since the characteristics of the rate errors are different between the two
systems the plant error covariance matricies used in the Kalman filter are chosen
differently. Table 4 shows the plant error covariance matricies used by the dynamic gyro
and gyro based systems. These values were determined through tuning over several
simulation runs.
77
Plant Error Covariance Matrix (Q)
Dynamic Gyro Gyro
50 0 0 0 0 0 3 0 0 0 0 0
0 50 0 0 0 0 0
0 3 0 0 0

0 0 50 0 0 0 0 0 3 0 0 0

0 0 0 0.5 0 0 0 0 0 0.03 0 0
0 0 0 0 0.5 0 0 0 0 0 0.03 0

0 0 0 0 0 0.5 0 0 0 0 0 0.03
Table 4. Plant Error Covariance Matricies
2. Gapped Star Tracker Coverage
Figure 49 shows the comparison of the dynamic gyro and gyro based attitude
determination systems when a 200 second gap in star tracker coverage is encountered
during maneuvering operations. Very little difference is evident in the rate error plots
because the bias build up during the star gap is so small compared to the noise in the
error. The effect of the bias, however, shows in the attitude quaternion error plots. At
the end of the star gap the error in the dynamic gyro based system is about five times that
of the gyro based system.
78
-3 -3
DG: Quaternion Errors Gyro: Quaternion Errors
x 10 x 10
3 3
q1 q1
q2 q2
2 2
q3 q3
q4 q4
1 1
0 0
-1 -1
-2 -2
-3 -3
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
-3 -3
DG: Rate Errors Gyro: Rate Errors
x 10 x 10
1.5 1.5
w1 w1
w2 w2
1 1
w3 w3
0.5 0.5
rad/s
rad/s
0 0
-0.5 -0.5
-1 -1
-1.5 -1.5
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 49. Results: Dynamic Gyro vs. Gyro with Gapped Star Coverage
C. DYNAMIC GYRO PLANT ERROR ANALYSIS
In this section the effects of major error sources on the performance of the
dynamic gyro are analyzed. Figure 50 shows the attitude determination plots for the
baseline simulation run. Dynamic gyro input parameters are varied and results are
compared with these plots.
79
-3 -3
Quaternion Errors Rate Errors
x 10 x 10
3 1.5
q1 w1
q2 w2
2 1
q3 w3
q4
1 0.5
rad/s
0 0
-1 -0.5
-2 -1
-3 -1.5
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Angular Momentum Error

0.4
0.3
0.2
0.1
Nms
-0.1
-0.2
-0.3
-0.4
0 50 100 150 200 250 300 350 400 450 500
Time (sec)
Figure 50. Baseline Attitude Determination Performance Results
1. Disturbance Torque Modeling
The accuracy of the dynamic gyro depends on knowledge of external disturbance

moments and internal spacecraft momentum. For the baseline simulation run no external
disturbances were modeled in the dynamic gyro. This leads to the relatively large drift in
attitude quaternions when uncorrected during the star gaps. If the spacecrafts attitude
with respect to the orbital reference frame can be determined the gravity gradient
disturbance torque can be modeled. For the Bifocal Relay Mirror satellite this
disturbance has the greatest effect. Figure 51 illustrates the increase in performance of
the dynamic gyro based attitude determination system when the gravity gradient moment
is modeled. At the end of the star gap the attitude errors are comparable to the gyro
based system. The total system angular momentum error is much smaller during the star
gap and improvement can even be seen during times of continuous star coverage.
80
-3
Quaternion Errors Angular Momentum Error
x 10
3 0.4
q1
q2 0.3
2
q3
q4 0.2
1
0.1
Nms
0 0
-0.1
-1
-0.2
-2
-0.3
-3 -0.4
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 51. Gravity Gradient Disturbance Modeled in the Dynamic Gyro
2. Rotation Axis Alignment Error
Alignment errors of momentum exchange control devices and slewing

appendages have direct effects on the dynamic gyro momentum error. Alignment errors
of fixed reaction wheels are easily compensated for by calibration but control moment
gyros have directionally variant momentum vectors with respect to the spacecraft body.
The results plots shown in Figure 52 are produced when a periodic alignment error on the
net momentum of the reaction wheels with a magnitude of approximately 0.5 degrees is
added to the simulation. A significant increase in attitude error is developed during the
star gap. The dynamic gyro does not track the system angular momentum as well even
during continuous star coverage.
-3
x 10
3 0.4
q1
q2 0.3
2
q3
q4 0.2
1
0.1
Nms
0 0
-0.1
-1
-0.2
-2
-0.3
-3 -0.4
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 52. Reaction Wheel Alignment Error Effects on the Dynamic Gyro
81
3. Potentiometer Quantization
The patterns shown in angular rate error plots for the dynamic gyro do not appear
as white noise during appendage motion because of quantization effects in the model of
the potentiometer that measures the relative orientation. The simulated quantization level
is 0.01 degrees in the baseline simulation. Since the effect alternates direction of error
the attitude quaternion errors are not affected significantly. Figure 53 shows simulation
performance results when the quantization level is decreased an order of magnitude. The
angular rate error looks more like white noise and the bias build up during the star gap
can be seen.
-3 -3
x 10 Quaternion Errors x 10 Rate Errors
3 1.5
q1 w1
q2 w2
2 1
q3 w3
q4
1 0.5
rad/s
0 0
-1 -0.5
-2 -1
-3 -1.5
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 53. Reduced Potentiometer Quantization Effect on Dynamic Gyro

Performance
4. Moment of Inertia Update Frequency
The inertia matrix for spacecraft with relatively small or slowly moving
appendages does not change quickly. In these cases processing power may be saved by
performing the system inertia calculation at a lower bandwidth than the dynamic gyro. In
the Bifocal Relay Mirror satellite the secondary body is large so even small slew rates
cause significant change in the system inertia matrix. Figure 54 shows the effects of
decreasing the inertia calculation frequency from the 20 Hz rate of the dynamic gyro to
once every 10 seconds. The quaternion error profile is significantly altered but the
magnitude of the error is only slightly increased. The momentum error in the dynamic
gyro takes longer to correct after the star gap since the Kalman filter attempts to correct
for all errors as if they were due to spacecraft momentum.
82
-3
x 10
3 0.4
q1
q2 0.3
2
q3
q4 0.2
1
0.1
Nms
0 0
-0.1
-1
-0.2
-2
-0.3
-3 -0.4
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 54. 10 Hz Moment of Inertia Calculation
D. STAR TRACKER MEASUREMENT ERROR ANALYSIS
The dynamic gyro based attitude determination system is highly reliant on

Kalman filter updates based on star tracker measurements. The quality of the Kalman
filter corrections depend on the accuracy of the trackers, the discrete measurement
interval and number of trackers providing measurements.
1. Star Tracker Accuracy
The results shown in Figure 55 are from simulations using random selection at 2
second intervals. The plot on the left is generated using a 1 sigma variance of 1x10-4 in
the horizontal and vertical measurements of the trackers. To generate the plot on the
right the variance is reduced by four times. The direct effect on attitude determination
performance is apparent.
83
-4 -4
Quaternion Errors Quaternion Errors
x 10 x 10
4 4
q1 q1
3 q2 3 q2
q3 q3
2 q4 2 q4
1 1
0 0
-1 -1
-2 -2
-3 -3
-4 -4
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 55. Star Tracker 1 Sigma Variance: 1x10 [left] vs 2.5x10-5 [right]
-4
2. Update Interval
Figure 56 shows a comparison of attitude determination performance using 2 and

6 second star measurement intervals. Shorter update intervals result in quicker
convergence and less steady state error.
-4 -4
x 10 x 10
4 4
q1 q1
3 q2 3 q2
q3 q3
2 q4 2 q4
1 1
0 0
-1 -1
-2 -2
-3 -3
-4 -4
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 56. Star Tracker Update Interval: 2 Seconds [left] vs 6 Seconds [right]
3. Star Tracker Selection
Because the star tracker measurements provide useful information along only two
axes at least two trackers must be active to maintain pointing control. Figure 57 shows
the comparison of an attitude determination simulation with updates spread evenly
between three trackers and one that uses two trackers with 95% of the updates coming
from the same sensor. There is no apparent degradation when the updates are spread
evenly between the two operational trackers.
84
-4 -4
x 10 x 10
4 4
q1 q1
3 q2 3 q2
q3 q3
2 q4 2 q4
1 1
0 0
-1 -1
-2 -2
-3 -3
-4 -4
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 57. Star Tracker Selection Comparison: Even Distribution [left] vs One
Tracker Favored [right]
85
86
VIII. CONCLUSIONS
Dynamic modeling provides an imperfect but operative means of estimating

multi-body spacecraft angular rates when mechanical gyro data are not available. An
attitude and angular rate estimation scheme is developed in this thesis that integrates the
dynamic gyro concept with an error state extended Kalman filter estimator that provides
precise attitude updates from star tracker measurements. Results indicate that the
determination system provides effective estimates for performing attitude control.
A. SUMMARY
The attitude determination system is incorporated into a multi-body spacecraft

attitude simulation for evaluation and analysis. Simulation is an extremely valuable
analysis tool for understanding the effects of error sources on system performance. It is
also a suitable mechanism for comparison with gyro based determination systems since a
mechanical gyro model can easily be inserted.
The effects of the primary error sources on the dynamic gyro performance are
investigated through simulation. It is shown that the corrections provided by a star
tracker based Kalman filter make the system robust to measurement and parameter
knowledge error sources. Significant improvement in attitude determination performance
is realized when disturbance torques are modeled. The other primary error sources
include the alignment of momentum exchange control devices and relative angle and rate
knowledge in large or quickly slewing appendages. Error effects are amplified during
star gaps when no corrections to the dynamic model are available.
The software implemented dynamic gyro essentially emulates the functions of a

set hardware gyroscopes. In a spacecraft where the mechanical gyros have failed or
become too erratic to be corrected by the Kalman filter the dynamic gyro may be a viable
replacement. Operated in tandem with mechanical gyros, either system can provide
redundant inertial angular velocity for improved attitude control.
87
This attitude determination concept is ideally suited to a spacecraft designed
specifically for its implementation with precise internal sensors and mechanisms to
monitor spacecraft parameters and integrated external torque estimation modeling. In
these systems the dynamic gyro can increase the lifetime and reliability of the spacecraft
while reducing the power requirements.
B. RECOMMENDATIONS
A strong recommendation for future work is simulation model tailoring to real

spacecraft hardware. The multi-body dynamics can be expanded to model an actual
satellites mass and inertia characteristics. Additionally, sensor parameters can be
modeled to match existing sensor error specs. With these modifications the simulation
model can be used conduct analysis and predict performance when the dynamic gyro
software is implemented on board the spacecraft. Future work may include using the
simulation to compare the predicted dynamic gyro based attitude determination
performance to actual gyro based performance recorded with telemetry data.
There are several improvements that can be made to the attitude determination
simulation model. It could easily be modified to provide the capability to mix gyroscope
and dynamic gyro data for redundant information. A parity matrix developed from the
pseudo-inverse concept can be generated to account for system observability in the over
determined case. This would allow performance analysis with selected gyro failures and
further aid in the evaluation of its utility as a back-up attitude determination scheme.
Another simulation improvement would be the incorporation of higher order dynamic
effects into the model including center of mass offsets and flexibility modes for the
appendage couplings.
Finally, improvements to the attitude determination scheme can be developed

including a torque error estimator as suggested by The Aerospace Corporation [Ref. 8].
The better the dynamics are modeled in the processing software the more effective the
dynamic gyro becomes as a replacement for the hardware gyroscope.
88
APPENDIX A: KALMAN FILTER BACKGROUND
Equation Section 1
The Kalman filter combines all available measurement data, plus prior
knowledge about the system and measuring devices, to produce an estimate of the desired
variables in such a manner that the error is minimized statistically [Ref. 2]. An
understanding of the Kalman filter requires some background in the theory of probability
of random variables and processes [Ref. 1].
A. PROBABILITY
The probability of an event, e, represents a possible outcome of a random

experiment and is written Pr(e). A random variable, X, can be thought of as a function of
the outcomes of some random experiment. The manner of specifying the probability with
which different values are taken by a random variable is the probability distribution
function, F(x)
F(x)=Pr(X x) (A.1)
This is a function over the range of possible values that shows the probability with which
the random variable takes on a value at or less than the value of the range. Its derivative
is the probability density function, f(x),
dF(x)
f(x)= (A.2)
dx
This function identifies the likelihood of a random variable assuming a particular value in
its range of possible values. Over the range of all possible values a characteristic of any
probability distribution or density function is

F()= f(u)du = 1 (A.3)
-
A joint probability density function can be defined for multiple variables. For two
random variables the joint density is given by
2 F2 (x,y)
f2 (x,y) = (A.4)
x y
89
B. EXPECTATION
The expectation of a random variable is defined as the sum of all values the
random variable may take, each weighted by the probability with which the value is
taken.

E[X]= xf(x)dx (A.5)
-
which is also called the mean value or first moment of X. This is a precisely defined
number toward which the average of a number of observations of the random variable X
tends. Since a function of a random variable is itself a random variable, the expectation
of a function of X, E[g(X)], can be expressed as

E[g(X)]= g(x)f(x)dx (A.6)
-
An important statistical parameter descriptive of X is its mean squared value defined by

E[X2 ] = x f(x)dx
2
(A.7)

which is also called the second moment of X. The root-mean-squared (rms) value of X is
E[X2 ] . The variance of a random variable is defined as the mean squared deviation of
the variable from its mean denoted by 2 .

2 = (x-E[X]) f(x)dx = E[X ]-E[X]
2 2 2
(A.8)
-
For zero mean random variables the variance is simply E[X2 ] . The square root of the
variance, , is called the standard deviation of the random variable.
A very important concept is that of statistical correlation between random

variables. A partial indication of the degree to which one variable is related to another is
given by the cross covariance, which is the expectation of the product of the deviations of
the two variables from their means,
90

E[(X-E[X])(Y-E[Y])]= (x-E[X])(y-E[Y])f 2 (x,y)dydx=E[XY]-E[X]E[Y]
- -
(A.9)
For a vector of random variables a symmetric covariance matrix can be defined

where the diagonal elements are the individual variances of the vector components and
the off diagonals are given by the cross covariances between the corresponding vector
components. The cross covariance, normalized by the standard deviations of X and Y, is
called the correlation coefficient.
E[XY]-E[X]E[Y]
= (A.10)
X Y
This is a measure of the degree of linear dependence between variables X and Y. If they
are completely independent is zero.
C. LEAST-SQUARES ESTIMATION CONCEPT
The optimality condition used in the Kalman filter is the minimization of

weighted least-squares error. The least-squares minimization problem involves a set of
v v
measurements, y , which are linearly related to the vector of state variables, x , by the
expression
v v v
y=Hx+v (A.11)
v
where v is an unknown vector of actual measurement errors with zero mean. The error
v
we seek to minimize is the measurement residual, z , given by
v v v
z=y-Hx (A.12)
v
where x is the estimate of the actual state vector. Since the sum of the squares of a
vector are generated by the inner product, the cost function, J, is
(
v v
) ( y-Hx
v v
)
T
J= y-Hx (A.13)
Minimization is obtained by differentiating and setting the result to zero
91
J v
v =0 (A.14)
x
and ensuring that the Hessian is positive semidefinite
2J
v 0 (A.15)
x 2
The resulting least-squares solution for the state estimate is
v
( ) v
1
x= H T H H Ty (A.16)
A weighted least-squares solution can be used when not all components of the
measurement residual are treated equally
v
( ) v
1
x= H T WH H TWy (A.17)
where W is a positive semidefinite matrix relating scale factors between components.

Expanding this solution to multiple estimates over time assumes all measurements are
used together in a batch processing scheme. Every time a new measurement becomes
available it is appended to the measurement vector and the estimated state vector includes
estimates corresponding to each of the accumulated measurements. The Kalman filter is
based on recursive processing where each measurement is used sequentially to generate
an optimal estimate of the current state without recomputing estimates of all previous
states. For this technique all previous information is embodied in the prior estimate and
state covariance matrix.
D. STATE ESTIMATE AND COVARIANCE
The Kalman filter estimation algorithm maintains the first two statistical moments
v
of the estimated state. The estimated state, x , is a vector of random variables whose
v
mean (first moment) is the actual state vector, x . The error in the state estimate is
defined as
v v v
x% @ x x (A.18)
and is thus assumed to be zero mean. The covariance of the state estimate error,
designated P, is given by
92
vv
% % T]
P=E[xx (A.19)
This covariance matrix associated with the state estimate error (second moment) provides
v
a statistical measure of the uncertainty in x and the correlation between the errors of its
components. For a state vector of two variables the state estimate error is given by
v% x% 1
x= x% (A.20)
2
and has the covariance matrix
x% x% 1x% 2 E[x% 12 ] E[x% 1 %x2 ]

2
P=E 1 = (A.21)
x% 1x% 2 x% 22 E[x% 1 x% 2 ] E[x% 22 ]
The diagonal elements of the state covariance matrix are the mean square errors in
the knowledge of the state variables while the off diagonal elements are indicators of
cross correlation between elements of the estimated error.
E. STOCHASTIC LINEAR DYNAMIC MODEL
The Kalman filter requires representation of system dynamics in a linearized

state-space form, a linear measurement model, and assumed characteristics of process
and measurement noise. For a continuous linear system the general state-space model
and measurement equations are given by
v& v v
x(t)=F(t)x+G(t)w(t)
v v (A.22)
y(t)=H(t)x(t)+v(t)
v v
where w(t) and v(t) are random vectors representing the unmodeled disturbance inputs
and measurement errors. These vectors are treated as unbiased (zero mean) white noise.
Note that if one of these vectors is known to have a nonzero bias the mean can be
augmented onto the state vector creating a new state space model with unbiased random
v
error. The covariance matricies associated with the process disturbance, w(t) , and the
v
measurement noise, v(t) , are
93
v v
Q = E[v(t)v(t)T ]
v v (A.23)
R = E[w(t)w(t) T ]
The equivalent discrete representation is
v v v
xk+1x= +k wk k k
v v v (A.24)
yk =Hk xk +vk
Here the subscript k represents values at time tk and subscript k+1 represents values at the
next discrete time step tk+1. Note that for time invariant systems the discrete state
transition matrix, k , is related to the continuous formulation by
k =eF(t k+1 -t k )
(A.25)
which depends only on the systems dynamics matrix, F, and the discrete time interval. If
the discrete interval is kept short enough the relation holds as an approximation. The
state transition matrix allows the calculation of state vector at the next discrete time step
in the absence of forcing functions. It obeys the differential equation
d k
= F(t) k (A.26)
dt
F. PROPAGATION OF ERRORS
The recursive Kalman filter requires the propagation estimate and error
covariance based on system dynamics. In the discrete implementation the error in the
current estimate given by
v v v
x% k = x k - x k (A.27)
has the covariance matrix representing the uncertainty in the estimate
v v
Pk = E[x% k %xTk ] (A.28)
v
The best estimate of the future state, x k+1 , is given by
v v
x k+1 = k k x (A.29)
The error in the new estimate
v v v
x%k +1 = k x% k k wk (A.30)
has the expected value
94
v v v v
E[x% k+1 ] = k E[x% k ] k E[wk ] = 0 (A.31)
v v
since x% k and w k are assumed to be unbiased. Thus the predicted estimate remains
unbiased. Note that if a deterministic input is added to the dynamic system model then
the identical quantity is added to both the actual and estimated state leaving the
estimation error unchanged. It can be shown that the associated state error covariance
matrix is given by
v v
Pk+1 = E[x% k+1x% Tk+1 ] = k Pk Tk + k Qk kT (A.32)
where Qk is the covariance of the random system disturbance. It is evident from this
equation that the size of disturbance directly impacts the uncertainty in the state estimate.
The system dynamic stability reflected in the state transition matrix also effects the
uncertainty. The covariance of neutrally stable or unstable systems will grow unbounded
over time in the absence of measurement corrections. The propagation of the state
estimate and its error covariance matrix is the prediction step of the Kalman filter
algorithm.
G. MEASUREMENT UPDATES
The correction step of the Kalman filter algorithm incorporates measurement

updates to refine the state estimate and error covariance. This step is executed only when
a measurement becomes available. When a measurement is taken we use the superscripts
- and + to denote values at a particular time before and after incorporation of the
measurement correction. Given a prior estimate of the system state at time tk we seek to
v
update our estimate based on the measurement yk in a linear, recursive form
v v v
x +k = K 'k x k- + K k y k (A.33)
where K'k and K k are matricies, as yet unspecified, that determine the relative weighting
of the prior estimate and current measurement. The error in this estimate can be shown to
be
v
( v
) v v
x% +k = K 'k + K k H k I x k + K 'k x% k- + K k v k (A.34)
v v
Since x% -k and v k are unbiased
95
v v
E[x% -k ]=0
v v (A.35)
E[v k ] = 0
this estimate can only remain unbiased for all given states if
K 'k = I K k H k (A.36)
With this requirement the estimator takes the form
v v
(v v
x +k = x -k + K k y k H k x -k ) (A.37)
and has corresponding error

v v v
x% +k = ( I K kH k ) x% -k + K k v k (A.38)
The state error covariance must also be updated.

v v
Pk+ = E[x% +k x% +T
k ]
v v vT T v v% -T v T
= E{(I Kk H k)x% -k [x% -T
k (I K k Hk ) + vk Kk ] + K k v k [x k (I K k Hk ) + v k Kk ]}
T T
(A.39)
Using the definitions
v v
Pk- = E[x% -k x% -T
k ] (A.40)
and
v v
R k = E[v k vTk ] (A.41)
and noting that the measurement errors are uncorrelated
v v v v
E[x% -k vTk ] = E[vk x% -T
k ]= 0 (A.42)
The covariance can be simplified to
Pk+ = (I K k H k )Pk- (I K k H k )T + K k R kK Tk (A.43)
The criterion for choosing the optimal Kk is to minimize a weighted scalar sum of
the diagonal elements of the error covariance matrix, Pk+ . Thus the cost function is
v v% +
J k = E[x% +T
k Sx k ] (A.44)
where S is any positive semidefinite matrix. Choosing S = I,
96
J k = trace[Pk+ ] (A.45)
which is equivalent to minimizing the length of the estimation error vector. To find the
value of Kk which yields a minimum, it is necessary to take the derivative of Jk with
respect to Kk and set it equal to zero. The result is
2(I Kk H k )Pk- HTk + 2Kk R k = 0 (A.46)

Solving for Kk,
K k = Pk- HTk [H k Pk- HTk ]1 (A.47)

which is referred to as the Kalman gain matrix. Since the Hessian of Jk is positive
semidefinite so Kk does indeed produce a minimum. Substituting into the equation for
the updated error covariance gives
Pk+ = Pk- Pk- H Tk [H kPk- HTk + R k ]1H k Pk- = ( I K kH k ) Pk- (A.48)
97
98
APPENDIX B: MATLAB CODE
A. INPUT PARAMETERS AND SIMULATION CALL
%Bifocal Relay Mirror Spacecraft Attitude Control System

%Input parameters and control options for attitude simulation model
clear
tic
stoptime=500; %Simulation stoptime

dt=.05; %fixed step size and computer time step
m1=2267.6059; %Mass of X-mit telescope (kg)

m2=972.3628; %Mass of RCV telescope (kg)
%Orbital parameters (circular)

Re=6378.1363; %Earth radius (km)
mu=3.986004415e5; %gravitational constant (km^3/s^2)
h=715; %orbit altitude (km)
r=Re+h; %orbit radius (km)
we=7.2921158533e-5; %Earth rotation rate (r/s)
w0=sqrt(mu/r^3); %orbital rate based on circular orbit/altitude (r/s)
nu0=0; %initial true anomaly (r)
u0=0; %initial angle of the magnetic dipole normal to the
%line of ascending nodes (r)
e1=11*pi/180; %magnetic dipole tilt angle
inc=40*pi/180; %inclination
Km=7.943e15; %magnetic constant (Nm^2/a)
Bm=Km/(r*1000)^3; %magnetic field strength (N/am)
%Disturbance torques
Mds=1e-4*[1;1;1]; %secular disturbance torques (Nm)
Mdp=4e-4*[1,1,1]; %periodic disturbances torque magnitudes (Nm)
pdper=[400,500,600]; %periodic disturbance torque periods (s)
pdfreq=2*pi./pdper; %periodic disturbance torque frequencies (r/s)
pdphase=[pi/4,pi/2,pi]; %periodic disturbance torque phase (r)
%Magnetic dumping control

magoo=0; %magnetic dumping on/off control (0=off, 1=on)
magsat=180; %maximum torque rod output
magk=5e5; %magnetic torquer gain
%Appendage (RCV scope) motion relative to body x-axis (sinusoidal)

99
bper=300; %period of oscillation (s)
bstart=.25*bper; %start motion (s)
bmotion=1*bper; %length of motion (s)
bstop=bstart+bmotion; %stop motion (s)
bfreq=2*pi/bper; %frequency of motion (rad/s)
bamp=15*pi/180; %amplitude/half maximum angular offset (rad)
bi=0; %initial angular offset (rad)
bdi=0; %initial relative angular rate (rad/s)
bacc=bamp*bfreq^2; %maximum realtive angular acceleration (rad/s^2)
bnv=1e-5^2; %relative angular noise variance from potentiometer (rad)
bns=0; %initial seed of potentiometer noise
bq=.01*pi/180; %quantization of potentiometer readings (rad)
%Command profile inertial to body accelerations (sinusoidal)

wdcfreq=2*pi*[1/800;1/500;1/600]; %frequency of acceleration profile variation
wdc=[.05*wdcfreq(1)^2;.5*wdcfreq(2)^2;-.1*wdcfreq(3)^2];
%amplitude of acceleration profile
wdcphase=[pi/8;-pi/12;0]; %phase of acceleration profile
cpoo=1; %turn commanded profile on/off (0=off, 1=on)
I1=[2997.28025,-3.9331,118.2824;
-3.9331,3164.18285,1.1230;
118.2824,1.1230,881.82105];
%moment of inertia matrix of primary body (XMIT scope) about its cm
I2=[1721.07340,-.0116,-7.8530;
-.0116,1559.85414,-12.5463;
-7.8530,-12.5463,182.89962];
%moment of inertia matrix of appendage (RCV scope) about its cm
dx1=.558354158; %Distance of X-mit telescope cm from S/C cm in x dir (m)

dy1=3.91788e-4; %Distance of X-mit telescope cm from S/C cm in y dir (m)
dz1=-.15226902; %Distance of X-mit telescope cm from S/C cm in z dir (m)
dx2=-1.302113918; %Distance of RCV telescope cm from S/C cm in x dir (m)
dy2=-9.13673e-4; %Distance of RCV telescope cm from S/C cm in y dir (m)
dz2=.355100092; %Distance of RCV telescope cm from S/C cm in z dir (m)
IC1=I1+m1*[dy1^2+dz1^2,dx1*dy1,dx1*dz1;
dx1*dy1,dx1^2+dz1^2,dy1*dz1;
dx1*dz1,dy1*dz1,dx1^2+dy1^2];
%moment of inertia matrix of primary body (XMIT scope) about S/C cm
roti=[1,0,0;0,cos(bi),sin(bi);0,-sin(bi),cos(bi)];
%rotation matrix corresponding to initial relative angular position
%of appendage (RCV scope)
%x-axis rotation of magnitude bi
Ic2i=roti'*I2*roti; %initial rotated moment of inertia matrix of appendage
%(RCV scope)
100
I2m=m2*[dy2^2+dz2^2,dx2*dy2,dx2*dz2;
dx2*dy2,dx2^2+dz2^2,dy2*dz2;
dx2*dz2,dy2*dz2,dx2^2+dy2^2];
%moment of inertia of appendage (RCV scope) about S/C C.M.
%due to mass offset
IC2i=roti'*I2*roti+I2m; %initial moment of inertia matrix of appendage
%(RCV scope) about S/C cm
Ii=IC1+IC2i; %initial S/C moment of inertia about S/C cm
%Wheel control law gains

%u=-Kqqe-Kwwe=-Kq(qc*q)-Kw(w-wc)
cldel=30; %delay control laws for attitude determination
Kq1=3000;
Kq2=7000;
Kq3=4500;
Kw1=1000;
Kw2=2000;
Kw3=1000;
%Initial body with respect to inertial quaternions

q1i=0;q2i=0;q3i=0;
qi=[q1i,q2i,q3i,sqrt(1-q1i^2-q2i^2-q3i^2)];
%Initial body with respect to orbit quaternions

q1oi=0;q2oi=0;q3oi=0;
qoi=[q1oi,q2oi,q3oi,sqrt(1-q1oi^2-q2oi^2-q3oi^2)];
%Initial body with respect to inertial direction cosine matrix

dcmi(1,1)=qi(1)^2-qi(2)^2-qi(3)^2+qi(4)^2;
dcmi(1,2)=2*(qi(1)*qi(2)+qi(3)*qi(4));
dcmi(1,3)=2*(qi(1)*qi(3)-qi(2)*qi(4));
dcmi(2,1)=2*(qi(1)*qi(2)-qi(3)*qi(4));
dcmi(2,2)=-qi(1)^2+qi(2)^2-qi(3)^2+qi(4)^2;
dcmi(2,3)=2*(qi(2)*qi(3)+qi(1)*qi(4));
dcmi(3,1)=2*(qi(1)*qi(3)+qi(2)*qi(4));
dcmi(3,2)=2*(qi(2)*qi(3)-qi(1)*qi(4));
dcmi(3,3)=-qi(1)^2-qi(2)^2+qi(3)^2+qi(4)^2;
Ci=[dcmi(1,1);dcmi(1,2);dcmi(1,3);
dcmi(2,1);dcmi(2,2);dcmi(2,3);
dcmi(3,1);dcmi(3,2);dcmi(3,3)];
%Direction cosine component vector
Ai=[dcmi(1,1),dcmi(1,2),dcmi(1,3);
dcmi(2,1),dcmi(2,2),dcmi(2,3);
dcmi(3,1),dcmi(3,2),dcmi(3,3)];
%Direction cosine/attitude matrix
101
%Initialize S/C momentum with wi at orbital rate
wi=[0;-w0;0]; %Initial angular rate of S/C primary body due to orbital rate
Hi=Ii*wi; %initial angular momentum of S/C due to orbital rate
%Gyro characteristics
grsb=1e-4*[-1,1.5,1]; %arbitrary gyro static bias
grv=1e-9*[1,1,1]; %variance of gyro rate noise
grn=[0,1,2]; %initial seed of gyro rate noise
gasb=[0,0,0]; %gyro acceleration static bias (zeros)
gav=1e-12*[1,1,1]; %variance of gyro acceleration noise
gan=[3,4,5]; %initial seed of gyro acceleration noise
grrwi=[0,0,0]; %initial rate random walk (zeros)
wbi=[0;0;0]; %initial gyro bias correction
gx=0*pi/180; %body to gyro(IRU) x-axis rotation
gy=0*pi/180; %body to gyro(IRU) y-axis rotation
G=[cos(gy),sin(gy)*sin(gx),-sin(gy)*cos(gx);
0,cos(gx),sin(gx);
sin(gy),-cos(gy)*sin(gx),cos(gy)*cos(gx)];
%Gyro alignment matrix - body to IRU
%Dynamic Gyro input options

dg=1; %choose Dynamic Gyro or regular gyros (1=DG, 0=gyros)
dggg=0; %calculate gravity gradient torque for DG (1=yes(requires or=1), 0=no)
dgsd=0; %secular disturbance torque known for DG (1=yes, 0=no)
dgpd=0; %periodic disturbance torque known for DG (1=yes, 0=no)
moiu=1; %steps per periodic moment of inertia calculation (1=every step)
or=0; %perform orbital reference calculations
%Attitude determination initialization with error

qe1i=q1i+.008;qe2i=q2i+.012;qe3i=q3i-.008;
qei=[qe1i,qe2i,qe3i,sqrt(1-qe1i^2-qe2i^2-qe3i^2)];
dcmei(1,1)=qei(1)^2-qei(2)^2-qei(3)^2+qei(4)^2;
dcmei(1,2)=2*(qei(1)*qei(2)+qei(3)*qei(4));
dcmei(1,3)=2*(qei(1)*qei(3)-qei(2)*qei(4));
dcmei(2,2)=-qei(1)^2+qei(2)^2-qei(3)^2+qei(4)^2;
dcmei(3,3)=-qei(1)^2-qei(2)^2+qei(3)^2+qei(4)^2;
Cei=[dcmei(1,1);dcmei(1,2);
dcmei(1,3);dcmei(2,1);dcmei(2,2);dcmei(2,3);
dcmei(3,1);dcmei(3,2);dcmei(3,3)];
Aei=[dcmei(1,1),dcmei(1,2),dcmei(1,3);
dcmei(2,1),dcmei(2,2),dcmei(2,3);
dcmei(3,1),dcmei(3,2),dcmei(3,3)];
102
qe1oi=q1oi+.02;qe2oi=q2oi+.03;qe3oi=q3oi-.02;
qeoi=[qe1oi,qe2oi,qe3oi,sqrt(1-qe1oi^2-qe2oi^2-qe3oi^2)];
wei=wi+1e-3*[-1;1;2];
Hei=Ii*wei;
%Initialize Kalman filter

xi=0*[ones(3,1)*eps*1e10;ones(3,1)*eps*1e8]; %initial state vector
Pi=[100*eye(3),zeros(3);zeros(3),eye(3)]*500; %initial covariance matrix
if dg==1
Qn=.5;
else %plant noise covariance constant DG vs gyro
Qn=.03;
end
Q=[100*eye(3),zeros(3);zeros(3),eye(3)]*Qn; %plant noise covariance
%small Q = good plant, biq Q = bad plant
R=10000*eye(2); %measurement noise covariance (10000)
E=[1,0,0;0,1,0]; %choose horizontal and vertical components only
%Star tracker parameters

corr=2*1/dt; %computer steps between stars
gapstart=100; %start gap in star tracker data (s)
stargap=200; %length of star gap (s)
gapstop=gapstart+stargap; %end gap in star tracker data (s)
t1x=135*pi/180; %body to star tracker one x-axis rotation
t1y=30*pi/180; %body to star tracker one y-axis rotation
T1=[cos(t1y),sin(t1y)*sin(t1x),-sin(t1y)*cos(t1x);
0,cos(t1x),sin(t1x);
sin(t1y),-cos(t1y)*sin(t1x),cos(t1y)*cos(t1x)];
%body to star tracker one rotation matrix
t2x=135*pi/180; %body to star tracker two x-axis rotation
t2y=-30*pi/180; %body to star tracker two y-axis rotation
%body to star tracker two rotation matrix
t3x=180*pi/180; %body to star tracker two x-axis rotation
t3y=0*pi/180; %body to star tracker two y-axis rotation
%body to star tracker two rotation matrix
mixT3=1; %choose whether or not to mix in star tracker 3 data (1=yes, 0=no)
T3th=.3; %choose threshold for random mix of star tracker 3 data
%(0=50%, higher=less, lower=more)
T12b=0; %choose to bias random selection between tracker 1 and 2
%(0=50/50, lower=favor1, higher=favor2)
103
stnb=[0,0]; %unknown star tracker bias error
stnv=1e-4^2*[1,1]; %variance of star tracker noise
stns=[0,1]; %initial seed of star tracker noise
%3 of 4 reaction wheels in pyramid constellation

rwsat=1*[1,1,1]; %constellation wheel torque saturation limits (Nm)
rwrl=10*[1,1,1]; %constellation wheel torque rate limiter (Nm/s)
hwi=-I2*[bdi;0;0]; %initialize wheels to cancel appendage momentum
rwa=45*pi/180; %angle of reaction wheels to x-y plane
RW=inv([-cos(rwa)*cos(pi/4),-cos(rwa)*cos(pi/4),cos(rwa)*cos(pi/4);
-cos(rwa)*cos(pi/4),cos(rwa)*cos(pi/4),cos(rwa)*cos(pi/4);
-sin(rwa),-sin(rwa),-sin(rwa)]); %body to wheel transform matrix
hwsi=RW*hwi; %distribute initial wheel momentum
rwnm=[0;0;0]; %mean of wheel noise
rwnv=1e-2^2*[1;1;1]; %variance of wheel noise
rwns=[7;8;9]; %initial seed of wheel noise
rwev=1e-8*[1,1,1]; %variance of wheel alignment walk error noise
rwes=[10,11,12]; %initial seed of gyro acceleration noise
rwam=1*5e-3*[1,2,-1]; %reaction wheel constellation alignment error magnitudes
rwaf=2*pi*[100,50,20]; %reaction wheel constellation alignment error frequencies
rwap=[pi/4,pi/3,pi/2]; %reaction wheel constellation alignment error phase
sim('acs_sim509') %Call to SIMULINK simulation
%Call ploting files for simulation results analysis

profileplots %command profile
ADplots509 %attitude determination performance
ACplots509 %attitude control performance
toc
B. SIMULATION RESULTS PLOTS
1. Commanded Profile Plots
%Laser Realy Spacecraft command profile plots
figure(1)
clf
plot(tout,qc)
title('Commanded Quaternions')
legend('q1','q2','q3','q4')
xlabel('Time (sec)')
grid on
104
figure(2)
clf
plot(tout,wc)
title('Commanded Angular Rates')
legend('w1','w2','w3')
ylabel('rad/s')
grid on
figure(3)
clf
plot(tout,b*180/pi)
ylabel('Deg')
title('Relative Angle to Recieve Telescope (alpha)')
grid on
figure(4)
clf
Hm=sqrt(H(:,1).^2+H(:,2).^2+H(:,3).^2);
plot(tout,Hm)
title('System Angular Momentum')
ylabel('Nms')
grid on
2. Attitude Determination Performance Results
%Bifocal Realy Mirror satellite attitude determination plots

%analyze performance of dynamic gyro and Kalman filter
%DG momentum error or gyro bias

if dg==1
figure(5)
clf
Hm=sqrt(H(:,1).^2+H(:,2).^2+H(:,3).^2);
Hmdg=sqrt(Hdg(:,1).^2+Hdg(:,2).^2+Hdg(:,3).^2);
plot(tout,Hm-Hmdg)
title('Angular Momentum Error')
ylabel('Nms')
axis([0,stoptime,-.4,.4])
grid on
else
figure(5)
clf
plot(tout,wb)
title('Gyro Biases')
105
legend('wb1','wb2','wb3')
ylabel('rad/s')
grid on
end
%KF states
figure(6)
clf
subplot(211),plot(tout,x(:,1:3))
title('State Attitude Errors')
legend('xa1','xa2','xa3')
ylabel('rad')
axis([0,stoptime,-.0008,.0008])
grid on
subplot(212),plot(tout,x(:,4:6))
title('State Rate "Bias" Errors')
legend('xb1','xb2','xb3')
ylabel('rad/s')
axis([0,stoptime,-.00008,.00008])
grid on
%AD errors
figure(7)
clf
plot(tout,(q-qkf))
title('Quaternion Errors')
legend('q1','q2','q3','q4')
axis([0,stoptime,-3e-3,3e-3])
grid on
figure(8)
clf
plot(tout,w-wkf)
title('Rate Errors')
ylabel('rad/s')
axis([0,stoptime,-1.5e-3,1.5e-3])
grid on
3. Attitude Control Performance Results
%Bifocal Realy Mirror satellite attitude control plots

106
%analyze commanded versus controled attitude
figure(9)
clf
plot(tout,Mws)
title('Wheel Torques')
ylabel('Nm')
legend('Wheel 1','Wheel 2','Wheel 3')
grid on
figure(10)
clf
plot(tout,hws)
title('Wheel Angular Momentums')
ylabel('Nms')
legend('Wheel 1','Wheel 2','Wheel 3')
grid on
figure(11)
clf
plot(tout,(q-qc))
title('Quaternion Error')
legend('q1','q2','q3','q4')
axis([0,stoptime,-4e-3,4e-3])
grid on
figure(12)
clf
plot(tout,w-wc)
title('Rate Error')
ylabel('rad/s')
axis([0,stoptime,-1.5e-3,1.5e-3])
grid on
107
108
LIST OF REFERENCES
1. Gelb, Arthur (ed.), Applied Optimal Estimation, The M.I.T. Press, Cambridge,
MA, 1974
2. Maybeck, Peter S., Stochastic Models, Estimation, and Control, Academic Press,
New York, NY, 1979.
3. Azor, Ruth, Bar-Itzhack, Itzhack Y., Harman, Richard R., Satellite Angular Rate
Estimation from Vector Measurements, Journal of Guidance, Control, and
Dynamics, Vol. 21, No. 3, May-June, 1998.
4. Harman, Richard R., Bar-Itzhack, Itzhack Y., Pseudolinear and State-Dependent

Riccati Equation Filters for Angular Rate Estimation, Journal of Guidance and
Control, Vol. 22, No. 5, 1999, pp. 723-725.
5. Azor, Ruth, Bar-Itzhack, Itzhack Y., Deutschmann, Julie K., Harman, Richard R.,
Angular-Rate Estimation Using Star Tracker Measurements, NASA Technical
Paper 1999.
6. Oshman, Yaakov, Makeley, F. Landis, Sequential Attitude and Attitude-Rate

Estimation Using Integrated_Rate Parameters, Journal of Guidance, Control,
and Dynamics, Vol. 22, No. 3, May-June 1999.
7. Gai, E., Daly, K., Harrison, J., Lemos, L. K., Star-Sensor-Based

Attitude/Attitude Rate Estimator, Journal of Guidance and Control, Vol. 8, No.
5, Sept.-Oct. 1985, pp. 560-565.
8. Herman, Louis K., Manke, Girard M., Heatwole, Craig M., Hamada, Brian T.,
McCain, Ian M., Pseudo Gyro, U.S. Patent, Docket No.: D365.
9. Wie, Bong, Space Vehicle Dynamics and Control, American Institute of

Aeronautics and Astronautics, inc., Temple, AZ, 1998.
10. Kolve, D. I., Describing an Attitude, American Astronautical Society Paper,

AAS 93-037, 6-10 February 1993.
11. Gray, C. W., Generic Star Tracker/IRU Attitude Determination Filter,

Proceedings of the American Astronautical Societys Guidance and Control
Conference, February 1999.
12. Stengel, Robert F., Stochastic Optimal Control: Theory and Application, John
Wiley-& Sons, inc., New York, NY, 1986.
109
13. Brown, Robert G., Introduction to Random Signals and Applied Kalman Filtering
3rd ed., John Wiley & Sons, Inc., New York, NY, 1997.
14. Markeley, F. L., Nelson, J. D., Zero-Gyro Safemode Controller for the Hubble
Space Telescope, Journal of Guidance, Control, and Dynamics, Vol. 17, No. 4,
pp. 815-822, July-Aug. 1994.
15. Crassidis, John L., Markeley, F. Landis, Predictive Filtering for Attitude
Estimation Without Rate Sensors, Journal of Guidance, Control, and Dynamics,
Vol. 20, No. 3, May-June 1997.
16. Voth, Chris, Herman, Louis K., Heatwole, Craig M., Manke, Girard M., Attitude
Determination for the STEX Spacecraft Using Virtual Gyros, Technical Paper,
October 15, 1999.
17. Kim, H. Y., Juang, J. N., Junkins, J. L., Recursive attitude prediction, American
Astronautical Society Paper, AAS 00-272, 20-21 March 2000.
18. Gray, C. W., Herman, L. K., Kolve, D. I., Westerlund, G. W., On-Orbit Attitude
Reference Alignment and Calibration, American Astronautical Society Paper,
AAS 00-272, 20-21 March 2000.
110
INITIAL DISTRIBUTION LIST
1. Defense Technical Information Center....................................................................2

8725 John J. Kingman Road, Suite 0944
Ft. Belvoir, VA 22060-6218
2. Dudley Knox Library ...............................................................................................2

Naval Postgraduate School
411 Dyer Road
Monterey, CA 93943-5101
3. Department Chairman, Code AA.............................................................................1

699 Dyer Road
4. Professor Brij N. Agrawal, Code AA/Ag ................................................................2

699 Dyer Road
5. Professor Hal A. Titus, Code EC/Ts........................................................................1

Department of Electrical and Computer Engineering
6. Professor Barry Leonard, Code AA/Ln ...................................................................1

699 Dyer Road
7. SRDC Research Library, Code AA .........................................................................1

699 Dyer Road
111
8. Donald I. Kolve, Associate Technical Fellow .........................................................1
Information, Space & Defense Systems
The Boeing Company
P.O. Box 3999 MC 8F-05
Seattle, WA 98124-2499
9. F. Landis Markley, Aerospace Engineer..................................................................1

Guidance, Navigation and Control Systems Engineering Branch
CODE 571 NASA Goddard Space Flight Center
Greenbelt, MD 20771
10. Girard M. Manke, Senior Engineering Specialist....................................................1

The Aerospace Corporation
2350 East El Segundo Blvd
El Segundo, CA 90245-4691
11. LT William J. Palermo, USN...................................................................................2

402 Kemp Lane
Chesapeake, VA 23325
112

01jun Palermo

Uploaded by

Copyright:

Available Formats

01jun Palermo

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

01jun Palermo

Uploaded by

Copyright:

Available Formats

Calhoun: The NPS Institutional Archive

Theses and Dissertations Thesis and Dissertation Collection

Angular rate estimation for multi-body

Downloaded from NPS Archive: Calhoun

Thesis Advisor: Brij N. Agrawal

Approved for public release; distribution is unlimited

16. PRICE CODE

ANGULAR RATE ESTIMATION FOR MULTI-BODY SPACECRAFT

Submitted in partial fulfillment of the

MASTER OF SCIENCE IN ASTRONAUTICAL ENGINEERING

NAVAL POSTGRADUATE SCHOOL

Approved by: ___________________________________________

Spacecraft with high performance attitude control systems requirements have

Figure 1. Basic Kalman Filter Block Diagram..................................................................2

Table 1. List of Symbols ...............................................................................................xvi

were invaluable to my understanding of attitude determination concepts. Also critical to

Naval Postgraduate School for their inspiration and guidance.

A. ANGULAR RATE ESTIMATION

An optimal estimator is a computational algorithm that processes measurements

Figure 1. Basic Kalman Filter Block Diagram

An important advantage of using an optimal estimator such as the Kalman filter is

C. ANGULAR RATE ESTIMATION FROM VECTOR MEASUREMENTS

The concept of using attitude sensor data to produce an estimate of spacecraft

It is possible to extract angular rate directly from time derivatives of measured

An effective method of dealing with non-linear dynamics is presented for use by

The pseudolinear filter concept can be extended to use quaternion measurements.

D. FULL ORDER ESTIMATORS

In an ideal attitude determination system where perfect knowledge of the

E. REDUCED ORDER ESTIMATORS

The performance of the reduced Kalman estimator can be improved by including

F. STAR SENSOR BASED ESTIMATORS

Another method of obtaining estimates of spacecraft angular rates is through

H. ATTITUDE ESTIMATION FROM CALCULATED RATE

The evaluation of this approach is performed through simulation using a model

The evaluation and analysis of attitude determination scheme proposed in this

A. CONCEPTUAL ATTITUDE SIMULATION MODEL

A conceptual representation of the simulation developed for this study is

Figure 2. Conceptual Attitude Simulation Block Diagram

Figure 3. Top Level SIMULINK Attitude Simulation Model

Primary spacecraft attitude control is conducted through momentum exchange

D. ATTITUDE DETERMINATION SYSTEM

The attitude determination subsystem is of primary concern in this thesis. It

The attitude determination computer simulation uses a trigger to control the

As three-axis stabilized spacecraft become more technologically advanced, their

Here we consider the dynamics of a spacecraft that consists of a primary body

A. BIFOCAL RELAY MIRROR SATELLITE DYNAMICS

Figure 5. Bifocal Relay Mirror Satellite

1. Rigid Body Equations of Motion

2. Multi-Body Equations of Motion

angular momentum due to the relative motion of the reaction wheels

Since the Bifocal Relay Mirror spacecraft is approximated as a system of rigid

To determine these angular momentums we need to consider the moments of

ITxx -ITxy -ITxz

I Rx'x' -IRx'y' -IRx'z'

Figure 6. SIMULINK Subsystem Diagram: X-Axis Rotation Matrix

a. Rate of Change of Spacecraft Inertia Matrix

shown that &I as a function of and & is given by

The components of the total spacecraft angular momentums given in Equation