Simulation of Spacecraft Attitude and Orbit Dynamics

Pasi Riihimäki, Jean-Peter Ylén

Control Engineering Laboratory
Helsinki University of Technology
PL-5500, 02015 TKK
E-mail: [email protected], [email protected]

KEYWORDS MODEL STRUCTURE AND MATHEMATICS

Simulation Model, Satellite, FDIR, Quaternion
Model Structure
ABSTRACT The spacecraft simulation model is organized like any
actual control loops (See Figures 1.). The interfaces of
In this paper, the simulation model of satellite attitude the components are defined and modeled in such a way
and orbit dynamics is discussed. The satellite attitude that the simulation model would be as modular as
model has been represented in term of a quaternion and possible. Modifications to the simulation model are easy
a ordinary differential equation is used to describe the to do and one part of the model can be easily replaced
satellite orbital motion. The different actuators and with another.
sensors have been modeled with suitable faults and
failures. The simulation model enables us to consider The model is initialized and controlled from the
the satellite motion under different environmental coordination level. This means that the model
perturbations (for example aerodynamic drag, external parameters and possible faults and failures in the FDIR
celestial body etc.) and failure in actuators and sensors. simulation case are defined.
The simulation model is utilized in the development of
attitude and orbit control algorithms or fault detection,
isolation and recovery (FDIR) technologies. Simulation
results are also given.

INTRODUCTION
During the last decades modeling, simulation, and wider
computational science and engineering have become
more and more important tools in the research and
development projects. The design phase has to be
reduced in time and cost when the use of new ideas and
tools becomes possible. This is also the trend in space
application in which the real tests are not possible or at
least they are expensive. New demands on the aerospace
and control engineering have become up and they have
Figures 1. Spacecraft simulation model structure.
to be able to answer to requirements.
Coordinate Systems
Spacecraft simulators or simulators in general, are
software tools that can be used by researchers, Three different coordinate systems are defined in the
engineers, students or everybody to analyze and assess simulator:
system operations, behaviors, and to answer to the 1. Inertial Coordinate System (ICS),
questions regarding phenomenon or product. The 2. Orbit Coordinate System (OCS), and
simulations are essential tools in the mission and 3. Body Coordinate System (BCS).
spacecraft control design. For example, the scientific
missions are unique and the instrumentation of a The inertial coordinate system is usually defined such
spacecraft is designed only for this specific mission. that the center of mass of the Earth (cm) acts as origin
There are not any ready-to-use platforms that can be and the direction of the axes are fixed to the solar
used. Hence, it is not possible to verify the operation of system. This kind of coordinate system is not exactly
control algorithms and strategies in real process but the inertial but it is enough for all engineering purposes
simulation environments can be used. (Sidi 1997). The Z-axis of the ICS is the rotation axis of
the Earth in a positive direction and the X-Y plane is the
There are plenty of companies that offer their simulation equatorial plane of the Earth, which is perpendicular to
services to the research institutes and space companies. the Earth’s rotation axis. The vernal equinox vector ϒ is
selected to be the X-axis of the ICS. Finally, the Y-axis

Proceedings 19th European Conference on Modelling and Simulation

Yuri Merkuryev, Richard Zobel, Eugène Kerckhoffs © ECMS, 2005
ISBN 1-84233-112-4 (Set) / ISBN 1-84233-113-2 (CD)
has been chosen in such a way that the ICS is right- ª q4' q3' −q2' q1' º ª q1 º
handed orthogonal coordinate system. « ' » « »
−q q4' q1' q2' » « q2 »
q '' = qq ' = « ' 3 ⋅ (3)
Orbit coordinate system is also a right-handed « q2 − q1' q4' q3' » « q3 »
« ' » « »
orthogonal coordinate system with origin in the center ¬« −q1 − q2' − q3' q4' ¼» ¬« q4 ¼»
of the satellite mass. The Z-axis is pointing towards the
center of the Earth; X-axis to the direction of satellite
perpendicular to the Z-axis, and Y-axis completes the SIMULATION MODEL
coordinate system such that it is right-handed and
The simulation model is realized in the MATLAB/
orthogonal. The third coordinate system, which has
SIMULINK-environment.
been fixed to the moving and rotating spacecraft,
defines satellite orientation.
Orbit Model
Rotation The motion of a celestial body is based on the quite
elementary principles of celestial mechanics. In the 17th
The attitude transformation in space can be executed by
century, J. Kepler provided three basic empirical laws
using various different aspects. In the simulation model,
that describe the motion of planet in unperturbed
the quaternion technique is used. The main feature of
planetary orbit. The orbital dynamics of a satellite is
quaternions is that they provide a convenient product
extensively explained in many books, for example (Sidi
rule for successive rotations and they have simple form
1997) and (Wertz 1978).
of kinematics (Wertz 1978, Wis´niewski 1996).
If we consider a system of two particles P1 and P2 of
The basic definition of the quaternion is a consequence
masses m1 and m2 and apply Newton’s second law and
of the property of the direction cosine matrix A that it
the law of gravity to the two-body system, we can get
has at least one eigenvalue of unity. This means that
the fundamental equation (4) of the motion of the two-
there is an eigenvector e (Euler axis) that is unchanged
in every rotation. The quaternion is defined as a vector body system where the symbol µ = G(m1+m2) and G is
the universal constant of gravitation. This equation
(1) where qi∈R, i, j and k satisfy the Hamilton’s rule
describes the motion of the particle P1 relative to the
(2) and where the length of the quaternion is unity. (Sidi
second mass P2.
1997)

q = q4 + q1i + q2 j + q3k µ
(1) r+
 r=0 (4)
r3
i 2 = j2 = k 2 = ijk = −1 In general, if a particle P moves in a force field F, the
ij = − ji = k momentum of the force F about origin O is
(2)
jk = −kj = i
M=r×F
ki = −ik = j
where r is the position vector of the particle P. The
When the Euler axis e of the rotation is known the angular momentum about origin is
connection between quaternion and the rotation Euler
axis is h=m(r×v)= r×p

­ q1 = e1 sin (α 2 ) where p is the linear momentum of the particle. Thus,

° the time rate of the angular momentum h is equal to the
° q2 = e2 sin (α 2 )
® moment of the force F.
° q3 = e3 sin (α 2 )
°q = cos (α 2 ) dh d
¯ 4 = ( r × mv ) = 0 + r × F = M (5)
dt dt
where ei is a component of Euler axis and α is the
magnitude of the rotation. The equation (5) states the fundamental fact that the
momentum acting on a particle is equal to the time rate
The final combined rotation of two successive rotations of the change of its angular momentum.
can be performed as a matrix-vector multiplication (3)
where q and q´ are the individual rotations. In space science it is common to describe the satellite
orbit by five numbers, known as orbital elements or
classical orbital elements (COE). A sixth element is
added to determine the location of the satellite in its
orbit (Wertz 1978 and Sidi 1997). The classical orbital
elements have been described in the Table 1. Because the derivative of classical orbital elements when the
these elements are poorly defined if e and/or i is equal perturbing force is conservative or non-conservative.
to zero, so-called equinoctial orbital elements (EOE) Knowing the initial condition of COE the Gauss
have been defined in terms of the classical orbital equations can be integrated to calculate the evolution of
elements. The equinoctial orbital elements have been the elements. Gauss equation is represented in equation
defined in Table 2. (7), where θ is the angle between satellite location
vector and the vector pointing towards perigee (See
Table 1. The classical orbital elements. (Sidi 1997) Figures 3.), p = a(1-e2), n = sqrt(µ/a3), r = p/(1+e
Symbol cos(θ)), and fr, fθ and fz are the components of the
a the semi major axis perturbing force along the radius vector direction r, the
e the eccentricity transverse orbit direction and the direction of the normal
i the inclination to the orbit plane, respectively.
Ω the right ascension of the ascending node
ω the argument of perigee To avoid the singularity due to the poorly defined
M the mean anomaly parameters, the Gauss equations can be rewritten in the
terms of the equinoctial elements as in (8), where
b=a⋅sqrt(1-P12-P12), h=nab, p/r=1+ P1sin(L)+P2cos(L),
L=ω+Ω+θ, and K=ω+Ω+E. (We and Roithmayr, C.M.
2001)

Figures 2. The definitions of the elements.

Table 2. The definitions of the equinoctial orbital

elements (EOE). (We and Roithmayr, C.M. 2001)
EOE
a a Figures 3. The spacecraft orbit and auxiliary circle.
P1 e sin ( Ω + ω )
Attitude Model
P2 e cos ( Ω + ω )
§i· Dynamics
Q1 tan ¨ ¸ sin ( Ω ) From equation (5) we get that the torque acting on the
©2¹
satellite body is equal to the derivative of the angular
§i·
Q2 tan ¨ ¸ cos ( Ω ) momentum of the spacecraft in the inertial coordinate
©2¹ system. Hence, in the rotating body coordinate system
l Ω+ω + M
¦ T = h I = h + Ȧ × h .
In Keplerian orbit the derivative of the first five orbital
elements are equal to zero. If the satellite orbital If momentum exchange devices are used in the control,
elements are known the satellite location r and the the angular momentum vector h = hB+hw where hB is the
velocity vector v can be calculated, and vise versa. angular momentum of satellite rigid body and hw is the
Algorithms to do this can be found in any textbook angular momentum of the momentum exchange
concerning orbital dynamics, for example (Sidi 1997), devices. Hence, the time rate of angular velocity of the
(Wertz 1978). satellite body is like in equation (6).

In the general case, in which any kind of perturbing

−1 ª º
force can exist, the equation of orbital motion is  = ( I s ) ⋅ « −I s Ȧ + ¦ Ti − h w − Ȧ × ( I s Ȧ ) − Ȧ × h w » (6)
Ȧ
¬ i ¼
µ
r+
 r = fp Kinematics
r3 The spacecraft attitude has been modeled as a
quaternion representation q=(q1, q2, q3, q4). Hence, the
with initial condition and where fp is the perturbing
equation (9), where ωi is the satellite angular velocity
force per unit mass. Due to the perturbing acceleration
about satellite body axis i, gives the derivative of
the orbital elements are not constants. Hence, so-called
quaternion vector.
Gauss form of Lagrange’s planetary equations describes
 2a 2
a = ª f r e sin (θ ) + fθ (1 + e cos (θ ) ) º¼
µp ¬ Tt = rt × Ft (10)
pª § e + cos (θ ) · º
e = « f r sin (θ ) + fθ ¨¨ cos (θ ) + ¸» The idea of reaction wheels (or momentum exchange
µ «¬ © 1 + e cos (θ ) ¸¹ »¼
devices) is to transfer the angular momentum of the
rf z cos (ω + θ ) whole system between different parts of the spacecraft
i =
µp without changing its overall internal angular
momentum. The achieved torque level is of the order of
 = rf z sin (ω + θ )
Ω (7) 0.05 – 2 Nm. (Sidi 1997)
µ p sin ( i )
f z r sin (ω + θ ) cos ( i ) The reaction-wheel is modeled as equation (11).
ω = −
µ p sin ( i )
°­Ȧ  w = f ( Tdem ) − f µ
1 pª § r· º (11)
− « f r cos (θ ) − fθ ¨1 + ¸ sin (θ ) » ®
°̄ h w = I w Ȧ w
e µ¬ © p¹ ¼
2rf 1 − e ª 2
§ r· º
M = n − 2r + « f r cos (θ ) − fθ ¨1 + ¸ sin (θ ) » In magnetotorquer, the control torque Tmag is generated
na nae ¬ © p ¹ ¼ by an interaction of the Earth’s geomagnetic field B(t)
with the magnetic dipole momentum m(t) (See
2a 2 ª p º
a =
h ¬ « ( P2 sin ( L ) − P1 cos ( L ) ) f r + fθ »
r ¼
equations (12) and (13)) where ncoil is the number of
coil, icoil(t) is the magnetotorquer current, Acoil loop area,
rªp ª § p· º and n̂ is the unit normal vector to the plane of the loop.
P1 = « cos ( L ) f r + « P1 + ¨1 + ¸ sin ( L ) » fθ
h¬r ¬ © r¹ ¼
Tmag ( t ) = m(t ) × B(t ) (12)
− P2 ( Q1 cos ( L ) − Q2 sin ( L ) ) f z º¼
m(t ) = ncoil icoil (t ) Acoil ⋅ nˆ (13)
rªp ª § p· º
P2 = « sin ( L ) f r + « P2 + ¨1 + ¸ cos ( L ) » fθ
h¬r ¬ © r¹ ¼ Sensors
+ P1 ( Q1 cos ( L ) − Q2 sin ( L ) ) f z º¼
(8) In the simulation model the modeled sensors are:
Q1 =
r
(1 + Q12 + Q22 ) sin ( L ) f z • coarse Earth and Sun Sensor (CESS),
2h • star tracker,
Q1 =
r
(1 + Q12 + Q22 ) cos ( L ) f z • magnetometer,
2h • gyro, and
r ª a p§ 2b · • GPS.
l = n − « ¨ P1 sin ( L ) + P2 cos ( L ) + ¸ f r
h ¬a + b r © a ¹
a § p· The CESS is modeled as a component that gives the
+ ¨1 + ¸ ( P1 cos ( L ) − P2 sin ( L ) ) fθ direction of Sun and Earth in the body coordinate
a+b© r¹
system.
+ ( Q1 cos ( L ) − Q2 sin ( L ) ) f z º¼
The star tracker is modeled as a component that gives
ª q1 º ª 0 ω3 −ω2 ω1 º ª q1 º the satellite attitude contaminated with an uncertainty
«q » « ω2 »» «« q2 »» that depends on the satellite angular speed.
d « 2 » 1 « −ω3 0 ω1
= (9) Magnetometer is modeled as a component that gives the
dt « q3 » 2 « ω2 −ω1 0 ω3 » « q3 »
« » « »« » magnitude and direction of the Earth’s magnetic field.
«¬ q4 »¼ «¬ −ω1 −ω2 −ω3 0 »¼ «¬ q4 »¼ WMM magnetic model is used in the simulator.

Actuators A gyroscope is modeled as an instrument that measures

The actuators are used to produce the control torques the angular speed of the spacecraft. The actual angular
and forces for the satellite attitude control. The modeled speed is contaminated with relatively small Gaussian
actuators are: random uncertainty. A GPS is modeled as an instrument
• thruster, that gives the satellite location in the inertial Earth
centered coordinate system.
• reaction-wheels, and
• magnetotorquer.
Faults
The thruster has been modeled as a thrust force vector One of the main aims of the spacecraft simulation
Ft affecting the satellite in position rt. Hence, the torque model is that it can be used in the FDIR-simulation.
about the center of the mass of the spacecraft is the Hence, the faults and failures have to be taken into
cross product between the position and the force vectors account already in the design and modeling phase. The
(Equation (10)). faults can occur in any part of the model and any kind of
faults are possible. Usually, the faults are either additive For spherical Earth, the gravitational force dfi acting on
or parametric but also a total blackout of a component is a s/c mass element dmi located at a position Ri is
possible.
− µ R i dmi
Perhaps, the most prevalent fault is ice building on the df i = .
Ri3
surface of any optical instrument increasing the
inaccuracy of this element.
Hence, the torque about the satellite geometric center
Environmental Torques due to a force dfi, at position ri, is

The main sources of the environmental torques are dTi = ri × dfi = ( ȡ + r 'i ) × dfi
represented in the Table 3.

Table 3. The main environmental torques (Wertz 1978). where ρ is the vector from the geometric center to the
Source Dependence Dominant cm and ri’ from cm to the mass element dmi. Assuming
that the cm and the geometric center of the s/c lie in the
Aerodynamic e-αr below ~ 500 km
Magnetic 1/r3 ~ 500 - 35 000 km
same point the gravity-gradient is
Gravity Gradient 1/r3 ~ 500 - 35 000 km
Solar Radiation Independent Interplanetary space 3µ ˆ ˆ º
Tgg = ªR s × I s R
Rs3 ¬
above synchronous s¼

altitude
Micrometeorites independent Normally negligible
where Rˆ s is a unit vector along R s and Is is the
The aerodynamic drag is one of the main environmental spacecraft inertial matrix. (Wertz 1978)
torques for the spacecraft in low orbit. The aerodynamic
drag model has been explained extensively, for A SIMULATION CASE
example, in the book (Wertz 1978). In this section, some simulation results obtained by the
above-described simulation model are presented. The
The force dfa on the surface elements dA is given by simulation case is simple and fancied.
equation (14) where N̂ is a outward normal of the
surface element dA, V̂ a unit vector of the translational The orbit of the simulation case is circular with an
velocity, ρ is the air density and CD is the drag altitude of 450 km and inclination 87°. The moments of
coefficient of the surface. In real terms, the drag inertia of the satellite are Ixx=36, Iyy=17, Izz=26, and Ixy=
coefficient CD is a function of the surface structure and Ixy= Ixz = Iyz =0 kgm2. The aim is that the attitude control
the local angle of attachment and its value is usually system shall ensure a three-axis stabilization of the
between 1 and 2. For all practical applications, the value satellite. The satellite attitude is measured by GPS
CD=2 can be used. (Wertz 1978) sensor and only three reaction wheels are used in the
control. The reaction wheels are mounted orthogonally
such that the rotation axes are along X, Y and Z-axis of
1
df aero = − CD ρV0 2 N
2
(
ˆ ⋅V
ˆ V
ˆ dA ) (14) the satellite body. Three PID-controllers are used to
calculate the control torques. The simulation results
In the simulation model, the satellite structure has been have been represented in Figures 4-8.
approximated by a collection of simple geometrical 5
Roll angle, [deg]

figures. Hence, the total aerodynamic torque is the sum

over the torques acting on individual parts of the 0

spacecraft. -5
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
5
Yaw angle, [deg]

Taero = ¦ ri × fi
i 0

=
1
2
CD ρV0 2 ¦ Ai N
i
(
ˆ ⋅V
i
ˆ V
ˆ ×r)
i
-5
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
5
Pitch angle, [deg]

Any nonsymmetrical body in orbit is subject to a 0

gravitational torque because of the variation in the -5

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Earth’s gravitational force over the object. Usually in Time, [s]
the literature, the gravity-gradient is only derived for the Figures 4. Attitude angles.
unrealistic spherical Earth model. Due to the non-
sphericity and the non-homogenous mass distribution
of the Earth the real gravitational field is granular.
 -3
x 10
0.8 q1
1 q2
0.6
ωx

0 q3
-1 0.4 q4

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
-3
0.2
x 10
1
0

qi
0
ωy

-1
-2 -0.2
-3
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 -0.4
-3
x 10
-0.6
5
0
ωz

-5 -0.8
-10
-15 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Time, [s]
Time, [s]
Figures 8. The attitude quaternions.
Figures 5. The satellite angular rates.
6
x 10 REFERENCES
6.97
6.96
Naasz, B.J. 2002. “Classical Element Feedback Control for
a

0.1
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Spacecraft Orbit Maneuvers” Thesis, Master of Science.
Virginian Polytechnic Institute and State University,
e

0.0998

2.5664
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Aerospace Engineering. 90 p.
2.5664
Sidi, M.J. 1997. “Spacecraft Dynamics & Control – A
i

2.5664
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Practical Engineering Approach” Cambridge University
4
3.8
3.6 Press. 409 p. ISBN 0-521-78780-7
Ω

3.4
3.2
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 We, B. Roithmayr, C.M. 2001. “Integrated Orbit, Attitude,
3.142
3.14 and Structural Control Systems Design for Space Solar
ω

3.138
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Power Satellites” NASA/TM-2001-210854. [10.3.2005]
2
0
http://techreports.larc.nasa.gov/ltrs/PDF/2001/tm/NASA-
M

-2 2001-tm210854.pdf
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Time, [s] Wertz, J.R. 1978. “Spacecraft Attitude Determination and
Figures 6. The classical orbital elements in simulation. Control” D. Reidel Publishing Company, Dordecht,
Holland.858 p. ISBN 90-277-0959-9
CONCLUSION Wis´niewski, R. 1996. “Satellite Attitude Control Using Only
Electromagnetic Actuation” Ph. D. Thesis. Denmark,
The satellite attitude and orbit simulation model with Aalborg University, Department of Control Engineering.
the most common actuators and sensors have been 175 p.
introduced in this paper. The simulation model can be
utilized both in the control algorithm designs and in the AUTHOR BIBLIOGRAPHIES
development of FDIR methods.
PASI RIIHIMÄKI was born in Ähtäri,
The simulation models have been implemented in the Finland and went to Helsinki
MATLAB/SIMULINK environment. University of Technology, where he
received a master’s degree in
30 automation and system science.
20 Nowadays, he is a postgraduate student
10 in the Control Engineering Laboratory in Helsinki
0 University of Technology.
-10

-20
JEAN-PETER YLÉN was born in
ωw

-30
Helsinki, Finland and went to Helsinki
-40
University of Technology, where he
received a master’s degree in chemical
-50
ω1
engineering and a doctor’s degree in
-60 ω2
automation and system science.
-70 ω3

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Time, [s]

Figures 7. The angular rates of reaction wheels.