Downloaded from orbit.dtu.

dk on: Dec 24, 2022

Satellite Dynamics and Control in a Quaternion Formulation (2nd edition)

Blanke, Mogens; Larsen, Martin Birkelund

Publication date:
2010

Document Version
Publisher's PDF, also known as Version of record

Link back to DTU Orbit

Citation (APA):
Blanke, M., & Larsen, M. B. (2010). Satellite Dynamics and Control in a Quaternion Formulation (2nd edition).
Technical University of Denmark, Department of Electrical Engineering.

General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright
owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.

 Users may download and print one copy of any publication from the public portal for the purpose of private study or research.
 You may not further distribute the material or use it for any profit-making activity or commercial gain
 You may freely distribute the URL identifying the publication in the public portal

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately
and investigate your claim.
Satellite Dynamics and Control in a Quaternion
Formulation
-
Lecture note for course 31365
Spacecraft Dynamics and Control
at DTU

Mogens Blanke and Martin Birkelund Larsen

Automation and Control Group
DTU Electrical Engineering Department
Technical University of Denmark

Version 2f - September, 2010

 Abstract

This lecture note treats modelling and attitude control design using a quater-
nion description of attitude for a rigid body in space. Dynamics and kine-
matics of a satellite is formulated as a non-linear model from Euler’s moment
equations and a description of kinematics using the attitude quaternion to
represent rotation. A general linearised model is derived such that the user
can specify an arbitrary point of operation in angular velocity and wheel
angular momentum, specifying the a inertia matrix for a rigid satellite. A
set of Simulinkr models that simulate the satellite’s nonlinear behaviour
are described and a Matlabr 1 function is described that has been written
to calculate the linear model in an arbitrary point of operation.

Copyright ⃝ c 2005-2010 Mogens Blanke, Martin B. Larsen and Technical

University of Denmark.
This document may be used for academic purposes with reference to this
document.

1
Matlab and Simulink are trademarks of MathWorks Inc. USA.
Contents

1 Introduction 3

2 Models in State-space Form 5

2.0.1 Solution of an LTI state-space equation . . . . . . . . 6
2.0.2 Laplace transform of an LTI system . . . . . . . . . . 6
2.0.3 Stability of an LTI model . . . . . . . . . . . . . . . . 7
2.0.4 Dealing with LTI models in Matlab . . . . . . . . . . 7
2.1 Nonlinear models and linearization . . . . . . . . . . . . . . . 8
2.1.1 Scalar case . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1.2 Linearisation in scalar case - example . . . . . . . . . 8
2.1.3 Vector case . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.4 Linearisation vector case - example . . . . . . . . . . . 10

3 Satellite Equations of Motion 11

3.1 Dynamics of satellite with reaction wheels . . . . . . . . . . . 11
3.2 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 Reaction Wheels as Actuators . . . . . . . . . . . . . . . . . . 14
3.3.1 Wheel Aligned with Satellite Coordinate System . . . 14
3.3.2 Wheel Aligned in Other Direction . . . . . . . . . . . 14
3.4 Actuator Dynamics - the Wheel Control Loop . . . . . . . . . 15
3.4.1 Combined Model . . . . . . . . . . . . . . . . . . . . . 15
3.5 Nonlinear model . . . . . . . . . . . . . . . . . . . . . . . . . 15

4 Nonlinear Simulation Model 19

5 Linear Model for Analysis 23

5.1 Linear Satellite Model with Actuator Dynamics . . . . . . . . 26

6 Matlab Function to get Linearised Model 27

6.1 Linear model without wheel dynamics . . . . . . . . . . . . . 27
6.2 Linear Models with wheel dynamics . . . . . . . . . . . . . . 28
6.3 Linear model for sampled system . . . . . . . . . . . . . . . . 28

1
CONTENTS 2

7 Attitude Control 29
7.1 Reference in both angular rate and attitude . . . . . . . . . . 30
7.2 Reference in Attitude . . . . . . . . . . . . . . . . . . . . . . 32
7.3 Stability Analysis and Simple Controller Design . . . . . . . . 33

A Programme Listings and Contents CD 39

A.1 linSatModNogg . . . . . . . . . . . . . . . . . . . . . . . . . . 40
A.2 romerEx . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
A.3 SATsimstart . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
A.4 RWAsim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
A.5 SATsim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Chapter 1

Introduction

The purpose of this lecture note is to derive and provide the complete non-
linear model of a satellite with reaction wheels, using quaternion represen-
tation of attitude, and to derive the linearised equations of motion of such
satellite for use in linear analysis and control design. The note extends and
supplements the treatment of attitude models that are provided in section
4.8 and later in the book by Marcel J. Sidi [5], which is based on Euler
angle description of attitude. Sidi analyse a number of feedback schemes
in chapter 6, based on an Euler angle model and assuming small deviations
from a reference. In chapter 7, he introduces attitude control based on
quaternion error feedback. However, the complete quaternion based model
was not derived by [5]. A linearised model was needed in the quaternion
formulation as well.
When we wish to make an analysis of the quaternion feedback scheme
similar to that done in chapter 6, a small signal model is needed for the
satellite described with attitude represented as a quaternion.
This note is aimed at students at a first course in spacecraft dynamics
and control, at the first year of the MSc level, where participants have
different backgrounds. Therefore, this note starts with some prerequisites.
In chapter 2, we define models of linear time invariant (LTI) systems as a
set of coupled first-order differential equations, referred to as the state-space
formulation in the automatic control literature, and analysis of stability
properties are discussed. It is then shown how a continuous-time nonlinear
system is linearised, using the Jacobian to obtain the state-space model in
vector-matrix notation.
Chapter 3 considers the dynamics of a spacecraft with reaction wheels
and derives a model of spacecraft attitude dynamics and kinematics in a
quaternion representation. The dynamics associated with torque control of
reaction wheels1 are included in the model. Gravity gradient effects are not
included in this version.
1
The terms reaction wheel and momentum wheel are used as synonyms in this note.

3
CHAPTER 1. INTRODUCTION 4

Chapter 4 derives a linearised model for analysis in an arbitrary point of

operation, i.e. angular momentum from reaction wheels and angular velocity
vector.
Chapter 5 introduces ways to design and implement feedback control and
discusses the issue of how to command reference pointing in various control
configurations. Finally, stability analysis is demonstrated on an example.
The appendix comprises Matlabr code used in examples in the text.
A comprehensive introduction to rotation sequences and quaternions is
given by Kuipers [4]. Some advanced topics in attitude control are presented
by Wie [7]. A comprehensive coverage of satellite dynamics is also available
in the now classical text by Huges [2].
Chapter 2

Models in State-space Form

Systems described in continuous-time by differential equations are often re-

written to a set of first-order coupled differential equations because these
are more easily solved. A differential equation in the variable x(t) and input
u(t)
...
x (t) + a3 ẍ (t) + a2 ẋ (t) + a1 x (t) = b1 u (t) , (2.1)
can be re-written if we introduce
x1 = x (t) , x2 = ẋ (t) , x3 = ẍ (t) ,
to    
x x2
d  1   
x2 = x3
dt
x3 −a1 x1 − a2 x2 − a3 x3 + b1 u
In a matrix-vector notation, when we write in component form,
      
x1 0 1 0 x1 0
d 
x2  =  0 0 1   x2  +  0 u
dt
x3 −a1 −a2 −a3 x3 b1
dx
In condensed matrix-vector notation, and using ẋ = dt ,
ẋ(t)= A(t)x + B(t)u.
If an output is included, we get the general form, referred to as a state-space
model
ẋ(t)= A(t)x(t) + B(t)u(t),
(2.2)
y(t)= C(t)x(t) + D(t)u(t).
A particular and important case is when the matrices A, B, C, D are
time-invariant. We then have a linear time-invariant (LTI) system:
ẋ(t)= Ax(t) + Bu(t),
(2.3)
y(t)= Cx(t) + Du(t).
The LTI system is used throughout classical control for analysis of various
properties, most important stability, and for design of linear control systems.

5
CHAPTER 2. MODELS IN STATE-SPACE FORM 6

2.0.1 Solution of an LTI state-space equation

The Laplace transform is a transform from the time domain to a domain
in the complex frequency s. The Laplace transform £ of an LTI system has
the following properties,

£(x(t)) = x(s), (2.4)

£(ẋ(t)) = sx(s) − x(t0 ).

The final value theorem Eq. (2.5) shows how one finds the steady state value
of a bounded function x(t), from its Laplace transform x(s) as:

lim (x(t)) = lim (sx(s)). (2.5)

t→∞ s→0

As an example, let x(t) be a unit step at time zero,

{
1 t≥0 1
x(t) = ⇒ £(x(t)) = . (2.6)
0 t<0 s

For a vector, x(t) and a time-invariant matrix A,

£(x(t)) = x(s), (2.7)

£(Ax(t)) = Ax(s).

Rigorous treatments of the Laplace transform is available in most introduc-

tory texts in automatic control and signal processing.

2.0.2 Laplace transform of an LTI system

Applying the Laplace transform to the LTI state-space model, Eq. (2.3),
gives:
sx(s)= Ax(s) + Bu(s) + x(t0 )
(2.8)
y(s)= Cx(s) + Du(s).
Let I ∈ Rn×n , be a diagonal unity matrix of dimension n×n , this is rewritten
as follows:

sIx(s)=Ax(s) + Bu(s) + x(t0 ) (2.9)

⇔ (sI − A) x(s) = Bu(s) + x(t0 )
⇔ x(s) = (sI − A)−1 (Bu(s) + x(t0 ))
( )
y(s)= C (sI − A)−1 B + D u(s) + C (sI − A)−1 Bx(t0 ).

The transfer function Hyu (s) from input u(s) to output y(s) is
( )
y(s)= C (sI − A)−1 B + D u(s) ≡ Hyu (s)u(s), (2.10)
CHAPTER 2. MODELS IN STATE-SPACE FORM 7

where Hyu (s) has k rows (same as number of elements in y) and m columns
(same as number of inputs), Hyu (s) ∈ Ck×m . Each element in Hyu (s) is a
single-input single-output (SISO) transfer function from a particular input
to a particular output.
Since Hyu (s) = C (sI − A)−1 B + D, all transfer functions in Hyu (s)
have the same denominator, namely the polynomial det (sI − A) .The nu-
merators will in general be different from each other.

2.0.3 Stability of an LTI model

The stability of a transfer function is determined by the roots of its denom-
inator, also referred to as the poles of the system. These poles are the roots
of the equation
det (sI − A) = 0, (2.11)
which are also the eigenvalues of the matrix A.
All eigenvalues of A must have a strictly negative real part for the system
to be asymptotically stable. Eigenvalues with zero real and zero imaginary
parts indicate pure integrations in the system. Pairs of eigenvalues symmet-
ric on the positive and negative imaginary axis indicate that an undamped
oscillating mode may exist in the system. Eigenvalues with positive real
part show that the linear system is unstable.
Eigenvalues of a system in state-space form are generally found using
numerical methods. Analytic solutions may be achieved in special cases.

2.0.4 Dealing with LTI models in Matlab

Matlab has commands dedicated to show features of LTI systems.
One need first to define a LTI system in state space form for Matlab:
sys = ss(A,B,C,D)
[Wn,z,p] = damp(sys) returns vectors Wn, z, p containing the
natural frequencies, damping factors and poles of the LTI model sys.
p = pole(sys) computes the poles p of the LTI model sys

pzmap(sys) computes the poles and (transmission) zeros of the LTI

model sys and plots them in the complex plane. The poles are plotted
as x’s and the zeros are plotted as o’s. Transmission zeros are the complex
frequencies s = z at which the rank of Hyu (z) is less than the normal rank
of Hyu (s). These are outside the scope of this note and will not be treated
further in this context [6, (section 4.5)].
The Matlab LTI toolbox documentation should be consulted for further
details.
CHAPTER 2. MODELS IN STATE-SPACE FORM 8

2.1 Nonlinear models and linearization

2.1.1 Scalar case
Let a system be described by the nonlinear scalar function

ẋ = g(ẋ, u),
(2.12)
y = h(x, y),
where x ∈ R1 , u ∈ R1 , y ∈ R1 . We wish to obtain a model that describes op-
˙ x̄, ū. Introduce deviations from the trajectory
eration around a trajectory x̄,
by

x̃ = x − x̄,
ũ = u − ū, (2.13)
ỹ = y − ȳ.
Then a Taylor expansion of Eq. (2.12) gives
∂g(x̄,ū) ∂g(x̄,ū)
x̃˙ + x̄˙ = g(x̄, ū) + ∂x x̃ + ∂u ũ + ...
(2.14)
∂h(x̄,ū) ∂h(x̄,ū)
ỹ + ȳ = h(x̄, ū) + ∂x x̃ + ∂u ũ + ...

Since x̄˙ = g(x̄, ū) and ȳ = h(x̄, ū) by definition of the ”bar” variables, and
truncation to first order of the Taylor series gives
∂g(x̄,ū) ∂g(x̄,ū)
x̃˙ = ∂x x̃ + ∂u ũ,
(2.15)
∂h(x̄,ū) ∂h(x̄,ū)
ỹ = ∂x x̃ + ∂u ũ.

This expression is linear in the tilted variables, i.e., we have obtained a linear
model. The parameters can be functions of the trajectory around which we
linearise.

2.1.2 Linearisation in scalar case - example

Let a nonlinear system be

ẋ = a0 + a1 x + a2 x2 + b1 u + b2 xu,
(2.16)
y = c1 x + c3 x3 .

The model linearised about the (x̄, ū) trajectory is

x̃˙ = (a1 + 2a(2 x̄ + b2 ū) x̃) + (b2 x̄) ũ,

(2.17)
ỹ = c1 + 3c3 x̄2 x̃.

A special case is that (x̄, ū) is time-invariant. In this case we obtain a linear
time-invariant system (LTI).
CHAPTER 2. MODELS IN STATE-SPACE FORM 9

2.1.3 Vector case

Let a system be described by the nonlinear vector equation

ẋ(t)= g(x(t), u(t))

(2.18)
y(t) = h(x(t), y(t))
where x ∈ Rn is the state vector, u ∈ Rm the control input, y ∈ Rk the
output. We wish to obtain a model that describes operation around a tra-
˙ x̄, ū. Introduce therefore the (small) deviation from the
jectory defined by x̄,
trajectory by
x̃ = x − x̄
ũ = u − ū (2.19)
ỹ = y − ȳ
Taylor expansion of Eq. 2.18 gives Eq. 2.20, where the argument (t) has
been omitted for brevity,

x̃˙ + x̄˙ = g(x̄, ū, d̄)+ ∂g(x̄,ū,

∂x
d̄)
x̃ + ∂g(x̄,ū,d̄)
∂u ũ + ...
(2.20)
ỹ + ȳ = h(x̄, ū, d̄)+ ∂h(x̄,ū,
∂x
d̄)
x̃ + ∂h(x̄,ū,d̄)
∂u ũ + ...
Using that x̄˙ = g(x̄, ū), ȳ = h(x̄, ū) and truncating the Taylor series
after the first order,

x̃˙ = ∂g(x̄,ū,
∂x
d̄)
x̃ + ∂g(x̄,ū,d̄)
∂u ũ
(2.21)
ỹ = ∂h(x̄,ū,
∂x
d̄)
x̃ + ∂h(x̄,ū,d̄)
∂u ũ
This defines a linear system
x̃˙ = Ax̃ + Bu ũ
(2.22)
ỹ = Cx̃ + Du ũ
The matrices A, Bu , are defined by the derivatives of the vector valued
function g with respect to vectors x, and u (and C, Du equivalently by
differentiation of h):

A = ∂g(x̄,ū,
∂x
d̄)
C = ∂h(x̄,ū,
∂x
d̄)
(2.23)
∂g(x̄,ū,d̄) ∂h(x̄,ū,d̄)
Bu = ∂u Du = ∂u
The derivative of a vector with respect to a vector is referred to as the
Jacobian,
 ∂g1 ∂g1 ∂g1 
∂u1 ∂u2 ··· ∂um
 
 ∂g2 ∂g2
··· ∂g2 
∂g(x̄, ū, d̄)  ∂u1 ∂u2 ∂um 
=
 .. .. ..

 (2.24)
∂u  .. 
 . . . . 
∂gn ∂gn ∂gn
∂u1 ∂u2 ··· ∂um
CHAPTER 2. MODELS IN STATE-SPACE FORM 10

is a matrix with n rows and m columns.

2.1.4 Linearisation vector case - example

Let a system be given by
[ ]
a0 + a1 x1 x2 + a2 x21 + b1 u1 + b2 x2 u1
ẋ = g(x, u) =
a3 + a4 x21 x22 + b1 u1 + b2 x21 u1
[ ]
y = h(x, u) = c1 x1 + c3 x1 x22 ,

hence
ẋ1 = a0 + a1 x1 x2 + a2 x21 + b1 u1 + b2 x2 u1
ẋ2 = a4 x21 x22 + b1 u1 + b3 x21 u1
y = c1 x1 + c3 x1 x22 ,
and application of differentiation as defined in Eq. 2.24 gives:
[ ]
[ ] x̃1
x̃˙ 1 = a1 x̄2 + 2a2 x̄1 a1 x̄1 + [b1 + b2 x̄2 ] ũ1
x̃2[ ]
[ ] x̃1 [ ]
x̃˙ 2 = 2a4 x̄1 x̄2 + 2b3 x̄1 ū1 2a4 x̄1 x̄2
2 2 + b1 + b3 x̄21 ũ1 (2.25)
x̃
[2 ]
[ 2
] x̃1
ỹ = c1 + c3 x̄2 2c3 x̄1 x̄2 .
x̃2

The resulting state-space equation is

[ ] [ ][ ] [ ]
x̃˙ 1 a1 x̄2 + 2a2 x̄1 a1 x̄1 x̃1 b1 + b2 x̄2
= + ũ1
x̃˙ 2 2a4 x̄1 x̄22 + 2b3 x̄1 ū1 2a4 x̄21 x̄2 x̃2 b1 + b3 x̄21
[ ]
[ ] x̃1
ỹ = c1 + c3 x̄22 2c3 x̄1 x̄2
x̃2
(2.26)
Chapter 3

Satellite Equations of Motion

3.1 Dynamics of satellite with reaction wheels

Recall from Euler’s moment equation [5, (Eq. 4.8.1) (p.107)] , [3, (p.50)]

ḣtot = Ne − ω × htot (3.1)

∑N
where htot = j (Ij ωj ) is the summation over all (N ) parts within the satel-
lite and the kinematics is described by the attitude quaternion [5, (p.104)]
, [1, (p.511)]
 
0 ω3 −ω2 ω1
1 1  −ω3 0 ω1 ω2 
q̇ = Ωq ≡  
q (3.2)
2 2 ω2 −ω1 0 ω3 
−ω1 −ω2 −ω3 0

A rotation is described by the rotation matrix [5, Eq. A.4.15] where the four
elements of the quaternion are the parameters used in the rotation matrix.
The standard form of satellite dynamics, with satellite inertia matrix Is ,
angular velocity of the satellite ω and hw is the total angular momentum
of wheels, all seen in the satellite’s coordinate system, Ne is total external
torque. With this notation, the dynamic equation reads
d
(Is ω) + ḣw = Ne − ω × Is ω − ω × hw (3.3)
dt
or
ω̇ = −I−1 −1 −1 −1
s (ω × Is ω) −Is ω × hw − Is ḣw + Is Ne (3.4)
Writing the cross product as a matrix operation using S() makes the
notation simpler in the sequel
 
0 −ω3 ω2
S(ω) =  ω3 0 −ω1  (3.5)
−ω2 ω1 0

11
CHAPTER 3. SATELLITE EQUATIONS OF MOTION 12

gives
ω̇ = −I−1 −1 −1 −1
s S(ω)Is ω − Is S(ω)hw − Is ḣ + Is Ne (3.6)
The combined nonlinear dynamic model of the satellite is completed by
noting that control torque in the body coordinate system is Nc and gives
the rate of change of the total angular momentum from wheels,

ḣ = −Nc , (3.7)

and hence the dynamics of the satellite actuated by reaction wheels

ω̇ = −I−1 −1 −1 −1
s S(ω)Is ω − Is S(ω)h + Is Nc + Is Ne . (3.8)

3.2 Kinematics
The kinematics of the satellite, is described using the attitude quaternion q
= [q1 , q2 , q3 , q4 ]T to represent a rotation. The orientation of the satellite is
obtained rotating from the inertial frame to the satellite frame by the direc-
tional cosine (rotation) matrix ASI . The rotation matrix is parameterized
by the quaternion, ASI (qSI ) or just ASI (q) for brevity. A vector measured
in the inertial frame is aI, the vector measured in the satellites body frame
has coordinates aS , then aS = ASI (q)aI . The kinematics of the satellite is
then     
q1 0 ω3 −ω2 ω1 q1
d   
 q2  = 1  −ω3 0 ω1 ω2   
  q2  (3.9)
dt  q3  2  ω2 −ω1 0 ω3   q3 
q4 −ω1 −ω2 −ω3 0 q4
In a compact notation,
1
q̇ = Ω(ω)q (3.10)
2
This compact notation is subsequently referred to in block diagrams. An-
other form of the detailed kinematics is also useful in the further analysis.
Equation 3.9 is easily rewritten to separate terms with q4 from other ele-
ments of the quaternion. Define a vector part of the first three components
of the quaternion and denote this by g, which is referred to as the Gibbs
vector  
q1
g =  q2  , (3.11)
q3
then 3.10 reads
1 1
ġ = − ω × g+ q4 ω (3.12)
2 2
1 T
q̇4 = − ω g
2
CHAPTER 3. SATELLITE EQUATIONS OF MOTION 13

Ne
Satellite dynamics Satellite kinematics

³ ³
Nc
+ Is -1

-1

hw
+ S( ) ½

Is
Is

Figure 3.1: Satellite model with nonlinear dynamics and kinematics, Equa-
tions 3.4 and 3.10, hw is the angular momentum of wheels, Ne an external
torque and Nc the control torque from reaction wheels, all given in the SBC
system. Output is q and ω

or written out in component form,

   
q1 0 −ω3 ω2  
    q1
d  q2  1 ω3 0 −ω1  
  =−   q2  (3.13)
dt q 3 2 −ω 2 ω1 0 
q3
q4 ω1 ω2 ω3
 
q4 0 0  
ω1
1 0 q4 0  
+  ω2  . (3.14)
2  0 0 q4 
ω3
0 0 0
[ ]
We can then write the kinematics equation using q = gT , q4 where g
is the first three components of the quaternion,
[ ] [ ] [ ]
d g 1 −S(ω) 1 I3×3
= g+ q4 ω (3.15)
dt q4 2 −ω T 2 0

The unit length of the quaternion is preserved in this operation. The model
of the satellite comprise the dynamics from Equation 3.4 and the kinematics
from Equation 3.10. This model is illustrated in Figure 3.1
CHAPTER 3. SATELLITE EQUATIONS OF MOTION 14

3.3 Reaction Wheels as Actuators

A reaction wheel is an electromechanical device composed of an electric
motor that drives a disc or wheel designed to have a large inertia for the
mass it has. The motor can be integrated in the wheel structure. With
permanent magnets incorporated in the circumference of the rotor and coils
wound in the stator part, a brush-less DC motor would result.

3.3.1 Wheel Aligned with Satellite Coordinate System

assume there is a single reaction wheel in the satellite. First, assume the
wheel’s axes is aligned with the satellite body coordinate system. Denote
the wheel inertia by J, the angular velocity of the wheel by ωw and the
angular momentum of the wheel by hw . The torque from the stator to the
rotor of the wheel as Nw , then
d
ḣ = (Jωw ) = Nw (3.16)
dt
For the satellite and the wheel as an entirety, we observe from Equation 3.4
that the control torque on the satellite is
d
Nc = −ḣ = − (Jωw ) (3.17)
dt

When a reaction wheel is used as actuator, ḣ of the wheel provides the (nega-
tive) control torque on the satellite and the angular momentum of the wheel
hw (which was measured in satellite coordinates) gives a cross-coupling term
in the dynamic equations of motion of the satellite, see Equation 3.4.

3.3.2 Wheel Aligned in Other Direction

When the wheel is not aligned with an axis of the satellite body coordinate
system (SBC), the shaft axis of the wheel, is rotated with respect to the
SBC. The angular momentum vector of the wheel, seen from the satellite
(SBC) is then needed, using a rotation from wheel coordinates to SBC. The
components of this angular momentum in the SBC system is a projection
of the wheel’s angular momentum along its axis of rotation. The projection
is a column vector, the directional cosines of which are
     
h1 ew · e1s a1w
 h2  =  ew · e2s  [|h|] =  a2w  [h] (3.18)
w w
h3 sbc ew · e3s a3w

If the wheel assembly consists of several wheels, say four, we have the
directional cosines for each wheel in the columns of the direction matrix for
CHAPTER 3. SATELLITE EQUATIONS OF MOTION 15

the wheels,
 
    hw1
h1 a11 a12 a13 a14  
 h2  =  a21 a12 a13 a14   hw2  (3.19)
 hw3 
h3 sbc a31 a12 a13 a14
hw4
Abbreviated to vector form, this is
h = Aw hw (3.20)
The matrix Aw has three rows and a number of columns equal to the
number of reaction wheels in the satellite.

3.4 Actuator Dynamics - the Wheel Control Loop

The above model assumed an ideal wheel servo, i.e. the wheel can accelerate
as fast as demanded. In practice, however, a torque control loop can not
be ideal but has a limited, although high bandwidth. The block diagram in
Fig. 3.2 shows the schematics of a wheel control loop.
Referring to the idealized case (A) in Figure 3.2
Nc (t) = −Nw (t) (3.21)
In case (B) of Figure 3.2, the speed control loop around ḣ gives rise to a
low-pass dynamic relation between Nw and Nc .
τw ḧ = −ḣ + Nw
Nc = −ḣ (3.22)
∫ t
h(t) = ḣ(τ )dτ + h(0)
0
In the Laplace domain, torque on the satellite reads, from Eq. 3.22
1
Nc (s) = − Nw (s) (3.23)
1 + sτw

3.4.1 Combined Model

The combined model of satellite and wheels is shown in Figure 3.3

3.5 Nonlinear model

The combined dynamic and kinematics equations Eqs. 3.4 and 3.15give the
general nonlinear model for the satellite’s angular motion,
   −1 
ω Is (−S(ω)Is ω − S(ω)h + Nc +Ndist )
d   
 g = − 21 S(ω)g+ 12 q4 I3×3 ω 

    (3.24)
dt q4 −2ω g
1 T

h −Nc
CHAPTER 3. SATELLITE EQUATIONS OF MOTION 16

-1
Nc
A: Idealized wheel model
dhw/dt

Ndw hw

-1
B: Model with wheel-servo Nc
dynamics included
dhw/dt

Ndw hw

³ ³
+ 1
Ww
-

Figure 3.2: Schematic wheel control loop. The upper part (A) shows the
idealized model where ḣ is identical to the demand torque. The bottom
part (B) shows a more realistic configuration where a speed control servo
loop is closed around the wheel. This gives the generation of ḣ a low-pass
behaviour. The bandwidth of the ḣ loop is equal to ωw .
CHAPTER 3. SATELLITE EQUATIONS OF MOTION 17

Reaction wheel assembly

-1 Aw
dhw/dt

hw
Nd

³ ³
1
AwT + Aw
Ww

-1

Ne Satellite dynamics Satellite kinematics

³ ³
Nc
+ Is -1

-1

+ S( ) ½

Is
Is

Figure 3.3: Satellite model with reaction wheel assembly as actuator, satel-
lite dynamics and kinematics.
CHAPTER 3. SATELLITE EQUATIONS OF MOTION 18

Limitations of the model

It is noted that the effects of flexibility of the satellite’s structure or ap-
pendages is not included in this model.

Remark 1 Change of reference coordinate system

As a remark, it is noted that the dynamic model has a structure that is in-
variant with the choice of coordinate system. As long as all quantities are
referred to the same coordinate system, the structure of the equation is un-
changed. The values of elements in the inertia tensor do change, needles to
say. This is seen as follows.
The effects caused by a rotation of both measured values and reference coordi-
nates for the inertia tensor to a satellite geometry defined coordinate system
is achieved by rotating the various quantities in the dynamic equation,
d
(Is ωp ) = −ω p × Is ω p −ωp × hp − ḣp + Ndist,p (3.25)
dt
then, rotating to the body system using ω p = Apb ω b , hp = Apb hb , Np =
Apb Nb , Is = Apb Isb ATpb gives

d ( )
Apb Isb ATpb Apb ω b = −Apb ω b × Apb Isb ATpb Apb ω b (3.26)
dt
−Apb ω b × Apb hpb − Apb ḣpb + Apb Ne,b + Apb Ndist,b

since Aω × Ah = A(ω × h), this simplifies to

d
(Isb ωb ) = −ω b × Isb ω b −ω b ×hb − ḣb + Ne,b + Ndist,b (3.27)
dt
which is again the dynamics in its original structure but now expressed in
body axes instead of principal axes.
Chapter 4

Nonlinear Simulation Model

The complete nonlinear model of satellite dynamics and kinematics, and

wheel actuation was implemented in Simulinkr in the program sat3 sim.mdl
(by the author).

The torque on the satellite is Nc = −Aw Ndw where Ndw is the wheel
torque demand and Aw the distribution matrix for wheel torque that is
defined in the Matlab program RWA sim.m. The distribution matrix has
its number of columns equal to the number of wheels. The number of rows
is always 3, since the Nc torque is the torque the three axes of the satellite
coordinate system.
The spacecraft parameters used are, as an example:

• Principal moments of inertia Ixx = 5, Iyy = 6, Izz = 7 kgm2 .

• The maximum angular momentum of the reaction wheel 0.12 Nms.

• The maximum wheel speed 280 rad/sec.

• moments of inertia of the reaction wheel 0.00043 kgm2 .

• The nominal wheel speed 100 rad/sec.

• The closed loop bandwidth the reaction wheel 10.0 rad/s.

• Maximum torque produced by the reaction wheel 0.0075 N.

• The reaction wheel assembly is defined by the number of wheels and

their alignment.

19
CHAPTER 4. NONLINEAR SIMULATION MODEL 20

Nw_dem
q_sat
To Workspace5
To Workspace

1
Nw_dem Nw_dem
Qeci2scb 1
Qics2sbc
2 w_w_ini
w_w_ini scbWbody 2
q_sat_init wbody_sbc

w_wheels 3
(0.0 0.0 0.0 1.0) w_sat_init
w_wheels
q_sat_init Rigid Body

(0.0 0.01 0.01) w_wheel w_sat

To Workspace6 To Workspace1
w_sat_ini

12:34

Digital Clock Display

Figure 4.1: Simulink mainwindow for sat3 sim.mdl. Input are control torque
demand, initial values for attitude quaternion of satellite, angular rate of
satellite, initial val be present in the workspace before you start a simulation.

3
q_sat_init

4
w_sat_init In q4 notation

Initial quaternion
initial rate

Qeci2scb 1
[0,0,0] scbNdist Qeci2scb
Ndist = 0 scbWbody scbWbody
scbHw

Kinematics
scbNw
2
Satellite body dynamics scbWbody

sbcHw
1 Nw_dem
Nw_dem
sbcNw

2 w_w_ini
w_wheels 3
w_w_ini
w_wheels
Momentum wheels

Figure 4.2: Level 2 of sat3sim with three components: satellite dynamics,

satellite kinematics and wheel model.
CHAPTER 4. NONLINEAR SIMULATION MODEL 21

scbWbody

scbNcc

3 scbHw
scbHw
Spacecraft Body

4 −1 K*u
scbNw 1
−1 Gain Iscb^−1 s 1
xo scbWbody
2
scbNdist

Initils
1
initial rate

Figure 4.3: Satellite dynamics is a direct implementation of Eulers moment

equation.

signal1
f(u) signal1
2 signal2
scbWbody Fcn1

f(u) signal2

Fcn2 1
xo s Q Norm 1
f(u) signal3 Qeci2scb
Normalization
Fcn3

f(u) signal4

Fcn4
1
Initial quaternion Initials

Figure 4.4: Kinematics implemented using quaternion parameterization.

CHAPTER 4. NONLINEAR SIMULATION MODEL 22

1 K*u 1
Nw_dem xos
Jest^−1
Integrator1

3
w_wheels

K*u K*u 1
xos K*u 1
Matrix Saturation J^−1 sbcHw
Gain Integrator2 J*A_w

K*u 2
sbcNw
A_w
2
w_w_ini

Figure 4.5: Simulation of reaction wheel assembly. The number of wheels

are determined by the dimension of the number of columns in the A war
matrix, defined in the Matlab function RWA sim.

• An example for n w = 4 (four wheels) is the thetraeder directions of

the wheels and the distribution matrix:
√ √ √ √ 
1 1
− 1
− 1
√ 3 √3 3 3
 2 
Aw =  3 − 3 2
0 0 . (4.1)
 √ √ 
0 0 − 23 2
3

• An example with three orthogonal wheels aligned with the axes of the
satellite gives the distribution matrix
 
1 0 0
Aw = 0 1 0 . (4.2)
0 0 1
Chapter 5

Linear Model for Analysis

This non-linear equation of motion is now linearised in an arbitrary point

of operation (ω̄, ḡ, q̄4 , h̄) in order to arrive at a set of linear state space
equations. This is needed to enable a stringent stability and performance
analysis for the linear control systems. The deviation from steady state
(point of linearization) is denoted by a tilde above the variables, ω = ω̄ + ω̃,
but the quaternion representation of attitude poses a specific problem. With
dg denoting the orientation at time t + dt relative to the attitude at time t,
then, since
   ( )   1 
[ ] dg1 e1 sin ( 12 ωdt) 2 ω1 dt
g̃  dg2   e2 sin 1 ωdt   1 ω2 dt 
=   (2 )   2
 dg3  =  e3 sin 1 ωdt  ≃  1 ω3 dt 
 (5.1)
q̃4 (12 ) 2
dq4 cos 2 ωdt 1
d
then dt q4 = 0 and S(ω)dg = 0 hence

d 1 1 1
g̃ = − S(ω)g̃+ q̃4 I3×3 ω = I3×3 ω (5.2)
dt 2 2 2
d d
h = h̄ + h̃; h = h̃ (5.3)
dt dt
The desired form of the linear equation of motion has a state vector

x = (ω̃1 , ω̃2 , ω̃3 , g̃1 , g̃2 , g̃3 , h1 , h2 , h3 )T

and has control input u = Nc . The external and disturbance torques have
zero mean. The state space equation is then

ẋ(t)= A(t)x(t)+Bu (t)Nc (t) + Bd (t)Ndist (t) (5.4)

where
∂fi ∂fi
Aij = ; Bij = (5.5)
∂xj ∂uj

23
CHAPTER 5. LINEAR MODEL FOR ANALYSIS 24

and  −1 
−I−1 −1
s S(ω)Is ω − Is S(ω)h + Is (Nc + Ndist )
 − 12 S(ω)g+ 12 q4 I3×3 ω 
f =


 (5.6)
−2ω g
1 T

−Nc
Using symbolic manipulation to calculate the Jacobians Eq. 5.5, we obtain

 
I−1
s Aω,ω 0 I−1
s Aω,h
A =  12 I3×3 0 0 ; (5.7)
0 0 0
 −1   −1 
Is Is
Bu =  0  ; Bd =  0  ;
−I3×3 0

with
Aω,ω = [Aω,1 , Aω,2 , Aω,3 ] (5.8)
and the three columns of Aω,ω are
 
ω2 I31 −ω 3 I21
Aω,1 =  −2I 31 ω1 −I 32 ω2 −ω 3 I33 +ω 3 I11 +h3  (5.9)
I21 ω1 +ω 2 I22 +I 23 ω3 −ω 2 I11 −h2
 
I31 ω1 +2I 32 ω2 +ω 3 I33 −ω 3 I22 −h3
Aω,2 = ω3 I12 −ω 1 I32  (5.10)
−ω 1 I11 −2I 12 ω2 −I 13 ω3 +ω 1 I22 +h1
 
−I 21 ω1 −ω 2 I22 −2I 23 ω3 +ω 2 I33 +h2
Aω,3 =  ω1 I11 +I 12 ω2 +2I 13 ω3 −ω 1 I33 −h1  (5.11)
ω1 I23 −ω 2 I13
Further,  
0 ω3 −ω2
Aω,h =  −ω3 0 ω1  (5.12)
ω2 −ω1 0
Equation 5.4 with matrices defined in Equations 5.7, 5.8, 5.9, 5.10, 5.11
and 5.12 constitute the linear model of the satellite in an arbitrary nominal
condition (point of operation). The nominal condition is expressed through
the parameters ω = ω̄, the average angular rate of the satellite, and h = h̄,
the average stored angular momentum expressed in satellite body coordi-
nates.
This linear model accept an arbitrary moment of inertia tensor, which
enables subsequent use for both controller design and analysis of sensitivity
(robustness) properties. Uncertainties include magnitude and rotation of
the inertia tensor and alignment of wheels. The basic dynamic properties
CHAPTER 5. LINEAR MODEL FOR ANALYSIS 25

change with the resulting angular momentum of the wheels. The changes
in linear dynamics could be analyzed should the satellite be demanded to
rotate along one of its axes, e.g. during maneuvers. The linear model is
then available should quantitative analysis be desired. Nonlinear analysis is
limited to mainly stating the more fundamental question of stability.

Example 1 Inertial Pointing Satellite

Let the following parameters apply for a satellite in inertial pointing con-
dition. Assume the principal axes coincide with the satellite body coordinate
system.

ω = (0, 0, 0)[rad/s]; hi ∈ [0, hmax ] [N ms]; (5.13)

 
Ixx 0 0
Is =  0 Iyy 0  [kgm2 ] (5.14)
0 0 Izz
The nominal model for small signals (linear model) case is hence, when
the satellite has ω = 0,

 
I−1
s Aω,ω 0 0
 1
0 0 ;
A= 2 I3×3 (5.15)
0 0 0
 −1   
Is I−1
s
Bu =  0 ; Bd =  0  ; (5.16)
−I3×3 0

and Aω,ω of Equations 5.8, 5.9, 5.10 and 5.11 simplifies to

 
0 −h3 h2
Aω,ω =  h3 0 −h1  (5.17)
−h2 h1 0

Example 2 Earth pointing satellite An Earth-pointing satellite in a circular

orbit has orbit period To . The orbit rate is ωo = 2π
To . The nominal angular
rate of the satellite is then ω=(0,-ωo ,0) [rad/s]. The negative sign arrives
from the definition of the orbit coordinate system with its x-axis along the
satellite’s velocity vector and the z-axis pointing towards Nadir (Center of
the Earth). Investigation of properties of the satellite should be done by
linearizing the system at ω=(0,-ωo ,0).
CHAPTER 5. LINEAR MODEL FOR ANALYSIS 26

5.1 Linear Satellite Model with Actuator Dynam-

ics
When the dynamic model of satellite Eqs. 5.4 and 5.4 gets its input torque
from wheels with actuator dynamics described in Eq. 3.22, the wheel state ḣ
need be added to the state vector. The linear dynamic model of the satellite
is hence extended with three states h˙x , h˙y and h˙z for ḣ in each of the three
axes of the satellite body coordinate system, denoted sbc. The state vector
( )T
for the aggregated model is x = ω, g, hsbc , ḣsbc . Input to the state
equation is Nd and Ndist (t)

ẋ(t) = Ax(t) + Bu Nd (t) + Bd Ndist (t) (5.18)

 
I−1
s Aω,ω 0 I−1s Aω,h −I−1
s
 1 I3×3 0 0 0 
A= 
2 , (5.19)
0 0 0 I3×3 
0 0 0 −kw I3×3
   −1 
0 Is
 0   0 
Bu = 

 , Bu = 
  0 

0
−kw I3×3 0

Note that torque demand and the wheel angular momentum vectors are in
satellite body coordinates in this equation.
The wheel dynamics should be included when a high performance satel-
lite control loop is designed. For more basic designs, it is certainly sufficient
to choose gains in the satellite attitude control such that the poles of the
closed loop satellite control are 10-15 times lower than the wheel bandwidth
(w bw in the parameter list below). The poles for the satellite with closed
loop attitude control are found using either of the following Matlab com-
mands
damp(satmodcl) or
[p,z]=pzmap(satmodcl)
on the LTI structure satmodcl for the closed loop system. The latter com-
mand provides both system poles and system zeros.
Chapter 6

Matlab Function to get

Linearised Model

The different versions of the models described above are calculated by the
Matlab function lin sat wheels mod. The call is

[satmod c, satmod dist] = lin sat wheels mod(Is, ws, hs, w bw, Ts)

Input and output parameters are described in the table:

variable unit type explanation

input
Is kgm2 matrix (3,3) Satellite’s inertia matrix
ws rad/s vector (3,1) Satellite’s mean angular velocity
hs N ms vector (n,1) Mean ang. wheel ang. momentum in sbc
w bw rad/s scalar Wheel bandwidth - same for all wheels
Ts s scalar Sampling time for a discrete time model
output use Ts = 0 to get continuous time model
satmod c Matlab LTI model using Ndw as input
satmod dist Matlab LTI model using Ndist as input
The last two parameters w bw and Ts can be omitted from the call.

6.1 Linear model without wheel dynamics

To obtain the model of Eq. 5.4 for the case with idealized wheels, use the call

[satmod c, satmod dist] = lin sat wheels mod(Is, ws, hs)

The result is a set of two Matlab LTI (Linear Time Invariant) structures,
satmod c, satmod dist, that comprises the system description in state space

27
CHAPTER 6. MATLAB FUNCTION TO GET LINEARISED MODEL28

form. The first comprise the model using wheel torque as input. The second
describes the system using disturbance torque as input. The matrices of the
state space model can be retained using the command

[A, Bu, C, D] = ssdata( satmod c )

[A, Bd, C, D] = ssdata( satmod dist )

6.2 Linear Models with wheel dynamics

To obtain the model of Eq. 5.19 with wheel dynamics, use the call
[satmod c, satmod dist] = lin sat wheels mod(Is, ws, hs, w bw)

6.3 Linear model for sampled system

To obtain a sampled (discrete time) model for the case with wheel dynamics
included, use the call

[satmod c, satmod dist] = lin sat wheels mod(Is, ws, hs, w bw, Ts)

where Ts is the sampling time in seconds.

Chapter 7

Attitude Control

Attitude control means to use the measurements and a reference to calculate

a torque demand that will make the measured state equal to the reference.
Following the argument in [5, (pp154-156)], that the satellite body coor-
dinate system and a target coordinate system will be aligned when zero
rotation is needed to align the satellite body system with the axes of the
target system. This can be done using a control law

u = −Kp qe − Kpd ω (7.1)

where Kp and Kpd are controller gain matrices, ω the angular velocity vector
of the satellite and qe is the quaternion error,
    
q1e q4r q3r −q2r −q1r q1s
 q2e   −q3r q4r q1r q2r   
 =   q2s  (7.2)
 q3e   q2r −q1r q4r −q3r   q3s 
q4e q1r q2r q3r q4r q4s
i.e. a quaternion that describes the rotation needed to align the satellite
body coordinate axes with those of the target coordinate system. The
quaternion qs describes the rotation to go from inertial axes to statelite
body coordinate system axes, qr the rotation of inertial axes to the refer-
ence (target) coordinate system axes.
The feedback in Eq. 7.1 has qt as reference for pointing, but no reference
in ω. When we want to make a reorientation manoeuvre with a spacecraft, it
could be convenient to introduce a reference also in ω, to make a maneuvre
that can follow an angular rate, e.g. when a satellite shall track another
object or a point on Earth. A control law with reference in both quaternion
and in angular velocity is:

u = −Kp qe − Kpd ωe (7.3)

where
ωe = ωs − ωr (7.4)

29
CHAPTER 7. ATTITUDE CONTROL 30

Continuous time Ndist

closed loop control Bu
state-space model

³
(gref, wref)T
Ndem
K Bu C
+
reference
+
(gs, ws, hw )T
T
(gs, ws)

A (gs, ws)T

measurements

Figure 7.1: Closed loop with output feedback and reference to the loop from
both ω and g in the linear state space representation of the satellite.

We will analyse both cases in the sequel, Eq. 7.3 with reference values
in both angular velocity and attitude, Eq. 7.1 with reference quaternion but
zero angular velocity as reference.

7.1 Reference in both angular rate and attitude

Attitude control means to use the measurements and a reference from both
angular rate and direction (quaternion) to obtain a closed loop control as
illustrated in Figure 7.1

The feedback law is taken to be

   
q1e q4e ω1e
Nd = −Kp  q2e q4e  − Kpd  ω2e  (7.5)
q3e q4e ω3e
where the quaternion error is defined as

qe = q−1
s qr (7.6)

and the angular velocity error is

ωe = ωs − ωr (7.7)

For small angles, and taking q4 positive,

 
g1r − g1s
 g2r − g2s 
qe ≃  
 g3r − g3s  . (7.8)
1
CHAPTER 7. ATTITUDE CONTROL 31

and the control law can be written as

   
g1r − g1s ω1r − ω1s
Nd = Kp  g2r − g2s  + Kpd  ω2r − ω2s 
g3r − g3s ω3r − ω3s
 
ω1r − ω1s
 ω2r − ω2s 
 
[ ]  ω3r − ω3s 
= Kpd Kp    (7.9)
g − g 
 1r 1s 
 g2r − g2s 
g3r − g3s

The linear model has the state vector x = [ω1 , ω2 , ω3 , g1 , g2 , g3, h1 , h2 , h3 ]T

and the measurement is y = [ω1 , ω2 , ω3 , g1 , g2 , g]T . With reference in both
ω and g, yref = [[ω1r , ω2r , ω3r , g1r , g1r , g3r ]T y this is the same as

[ ]
Nd = Kpd Kp (yr − y)
[ ] [ ]
= − Kpd Kp Cx+ Kpd Kp yr (7.10)
= −Ky Cx + Kr yr
[ ]
where ] lines define the gain matrices Ky = Kpd Kp and
[ the last two
Kr = Kpd Kp .
The transfer function from reference to output with the closed feedback
is then
y(s) = C(sI9,9 − (A − Bu Ky C))−1 Bu Kr yr (s) (7.11)
The form of this expression is

y(s) = C(sI9,9 − Acl )−1 Bcl yr (s) (7.12)

with

Acl = (A − Bu Ky C) (7.13)
[ ]
Bcl = Bu Kr = Bu Kpd Kp (7.14)
Ccl = C and Dcl = 06,6 (7.15)

The closed loop transfer functions 7.12 are calculated using the commonplace
Matlab functions (see below), using the LTI system definition

sys cl = ss(Acl , Bcl , Ccl , Dcl ). (7.16)

CHAPTER 7. ATTITUDE CONTROL 32

Continuous time Ndist

closed loop control Bu
state-space model

Ndem +

³
[ gref ]T [ws, gs]T
+ +
Kr Bu C

+ + +
Reference
for g only [ws, gs, hw ]T

- Ky
measurements

[ws, gs]T

- Ky

Figure 7.2: Closed loop with output feedback and reference to the loop from
g only. Linear state space model of the satellite.

7.2 Reference in Attitude

When the reference is not desired in angular velocity, the block diagram 7.2
is applicable.
With reference from g only, yr = [g1r , g2r , g3r ]T . This gives

[ ]
Nd = − Kpd Kp y + Kp gr
[ ]
=− Kpd Kp Cx + Kp gr (7.17)
= −Ky Cx + Kr gr

where the last two lines define the gain matrices Ky and Kr The transfer
function from reference to output is therefore

y(s) = C(sI9,9 − (A − Bu Ky C))−1 Bu Kr gr (s) (7.18)

Using the abbreviations

Acl = (A − Bu Ky C) (7.19)
Bcl = Bu Kr = Bu Kp (7.20)
Ccl = C and Dcl = 06,3 (7.21)

we can again analyse the closed loop system using the form

y(s) = C(sI9,9 − Acl )−1 Bcl yr (s) (7.22)

CHAPTER 7. ATTITUDE CONTROL 33

7.3 Stability Analysis and Simple Controller De-

sign
The essential properties of a closed loop control system include closed loop
eigenvalues to analyse stability, bode plot from reference input to output,
step response to unit step at the reference input, and response to torque
disturbances.
Design of closed loop control systems can be achieved in several elegant
ways, including optimal and robust control methods. Even with a limited
background in control theory, one can design well performing controllers
using the stepwise procedure below.

Algorithm 1 Stepwise procedure for simple control system design

Commands are Matlabr control systems toolbox and the supplied func-
tions for the course.

1. Determine the linear model without wheel dynamics using

[satmod c, satmod dist] = lin sat wheels mod(Is, ws, hw)

2. Eigenvalues in the open loop (no control) system using

damp(satmod c)
pzmap(satmod c)

3. Get the A, Bu , C, D matrices using

[A, Bu , C, D] = ssdata(satmod c)

4. Get the A, Bdist , C, D matrices using

[A, Bd , C, D] = ssdata(satmod dist)

5. Define the closed loop by

Dcl = zeros(size(C,1),size(Kr ,2))
[satmod c closed] = ss((A − Bu Ky C, Bu Kr , Ccl , Dcl ))

6. Investigate eigenvalues of the closed loop by

damp(satmod c closed)
pzmap(satmod c closed) and
[P,Z]=pzmap(satmod c closed)

7. Once a suitable stable controller has been found, determine frequency

response properties by
bode(satmod c closed)

8. and step responses by

step(satmod c closed)
CHAPTER 7. ATTITUDE CONTROL 34

9. Investigate disturbance rejection by

[satmod dist closed] = ss((A − Bu Ky C, Bd , C, D))
bode(ssatmod dist closed)
step(satmod dist closed)

10. iterate in the above procedure until specifications are met.

Example 3 Rømer-like Satellite

Let the Rømer-like satellite have

ω = (0, 0, 0)T [rad/s]; hi = [0, 0.5, 0.2]T [N ms]; (7.23)

 
14.3 0 0
Is =  0 13.6 0  [kgm2 ] (7.24)
0 0 4.6
Define the desired bandwidth of the closed-loop system by ωbw cls = 0.2
[rad/s]. A simple design approach is to select gains to be
2
kp (j) = Is (j, j)ωbw cls
kpd (j) = 2Is (j, j)ωbw cls (7.25)

Kp = diag(kp (1), kp (2), kp (3))

Kpd = diag(kpd (1), kpd (2), kpd (3)) (7.26)

We investigate the properties of this simplistic design idea by using the

steps of algorithm 1. The gain factors used are the above given Kp and
Kpd . Combinations of two sets of nominal values are used, {ωs = [0; 0; 0],
ωs = [0; −0.003; 0] }and {hs = [0; 0; 0],hs = [0; 0; 0.5] }The results are shown
in Figures 7.3, 7.5 and 7.7.
CHAPTER 7. ATTITUDE CONTROL 35

pole zero plot open loop pole zero plot closed loop
0.1 0.1

0.08 0.08

0.06 0.06

0.04 0.04

0.02 0.02
imag

imag
0 0

−0.02 −0.02

−0.04 −0.04

−0.06 −0.06

−0.08 −0.08

−0.1 −0.1
−0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05
real real

(a) open loop (b) closed loop

Figure 7.3: R loop (b). Angular velocity is ωs = [0; 0; 0] rad/s and angular
momentum is hs = [0; 0; 0] Nms. The closed loop eigenvalues are all in -0.2
as intended. The small imaginary part is due to numerics. Three eigenvalues
at the origin are cancelled by three zeros. The origin of this is the physical
integration from ḣ to h.

pole zero plot open loop pole zero plot closed loop
0.01 0.01

0.008 0.008

0.006 0.006

0.004 0.004

0.002 0.002
imag

imag

0 0

−0.002 −0.002

−0.004 −0.004

−0.006 −0.006

−0.008 −0.008

−0.01 −0.01
−0.01 −0.008 −0.006 −0.004 −0.002 0 0.002 0.004 0.006 0.008 0.01 −0.01 −0.008 −0.006 −0.004 −0.002 0 0.002 0.004 0.006 0.008 0.01
real real

(a) open loop (b) closed loop

Figure 7.4: Same as Figure 7.3 but zoomed

CHAPTER 7. ATTITUDE CONTROL 36

pole zero plot open loop pole zero plot closed loop
0.1 0.1

0.08 0.08

0.06 0.06

0.04 0.04

0.02 0.02
imag

imag
0 0

−0.02 −0.02

−0.04 −0.04

−0.06 −0.06

−0.08 −0.08

−0.1 −0.1
−0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05
real real

(a) open loop (b) closed loop

Figure 7.5: Roemer-like satellite pole-zero plot for open (a) and closed loop
(b). Angular velocity is ωs = [0; −0.003; 0] rad/s and angular momentum
hs = [0; 0; 0.5] Nms. The angular velocity around the y-axis, and the mo-
ment of inertia with Ix < Iy < Iz , cause the open loop system to be unstable.
The closed loop is stable with closed loop poles only slightly affected by the
perturbation of the open loop system.

pole zero plot open loop pole zero plot closed loop
0.01 0.01

0.008 0.008

0.006 0.006

0.004 0.004

0.002 0.002
imag

imag

0 0

−0.002 −0.002

−0.004 −0.004

−0.006 −0.006

−0.008 −0.008

−0.01 −0.01
−0.01 −0.008 −0.006 −0.004 −0.002 0 0.002 0.004 0.006 0.008 0.01 −0.01 −0.008 −0.006 −0.004 −0.002 0 0.002 0.004 0.006 0.008 0.01
real real

(a) open loop (b) closed loop

Figure 7.6: Same as Figure 7.5 but zoomed

CHAPTER 7. ATTITUDE CONTROL 37

pole zero plot open loop pole zero plot closed loop
0.1 0.1

0.08 0.08

0.06 0.06

0.04 0.04

0.02 0.02
imag

imag
0 0

−0.02 −0.02

−0.04 −0.04

−0.06 −0.06

−0.08 −0.08

−0.1 −0.1
−0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05
real real

(a) open loop (b) closed loop

Figure 7.7: Roemer-like satellite pole-zero plot for open (a) and closed loop
(b). Angular velocity is ωs = [0; −0.003; 0] rad/s and angular momentum
hs = [0; 0; 0.5] Nms. The bias in angular momentum has moved the open
loop pole-zero pattern significantly. Despite the large perturbation, the
closed loop poles show that the closed loop system is still nicely damped.

pole zero plot open loop pole zero plot closed loop
0.01 0.01

0.008 0.008

0.006 0.006

0.004 0.004

0.002 0.002
imag

imag

0 0

−0.002 −0.002

−0.004 −0.004

−0.006 −0.006

−0.008 −0.008

−0.01 −0.01
−0.01 −0.008 −0.006 −0.004 −0.002 0 0.002 0.004 0.006 0.008 0.01 −0.01 −0.008 −0.006 −0.004 −0.002 0 0.002 0.004 0.006 0.008 0.01
real real

(a) open loop (b) closed loop

Figure 7.8: Same as Figure 7.7 but zoomed

Bibliography

[1] James A. Wertz (ed). Spacecraft Attitude Determination and Control.

Kluwer Academic Publishers, 7th. reprint edition, 1984.

[2] Peter C. Huges. Spacecraft Attitude Dynamics. John Wiley and Sons,
1986.

[3] Marshall H. Kaplan. Modern Spacecraft Dynamics and Control. John

Wiley and Sons, 1976.

[4] Jack B. Kuipers. Quaternions and Rotation Sequences. Princeton Uni-

versity Pres, 2002.

[5] Marcel J. Sidi. Spacecraft Dynamics and Control - A Practical Engineer-

ing Approach. Cambridge University Press, 2001.

[6] Sigurd Skogestad and Ian Postlethwaite. Multivatiable Feedback Control:

Analysis and Design. Wiley, 2005.

[7] Bong Wie. Space Vehicle Dynamics and Control (2nd edition). American
Institute of Aeronautics and Astronautics, 2008.

38
Appendix A

Programme Listings and

Contents CD

The CD enclosed with the printed version of this report comprise the func-
tions listed below. The files are available for participants in DTU course
31365: Spacecraft Dynamics and Control via file-sharing at the DTU Cam-
pusnet.

• linSatModel.m. Calculate A,B,C,D matrices of linear model.

• roemerEx.m. Script file that runs the Roemer-like satellite design

example.

• SATsimstart.m. Script to set parameters used in the SIMULINK

simulation.

• RWAsim.m. Script to set reaction wheel assembly parameters used

in the SIMULINK simulation.

• SATsim.m. Script file to set satellite main body parameters in the

simulation.

• sat3sim.mdl. SIMULINK model with nonlinear simulation of satel-

lite with reaction wheel assembly as actuator. Gravity gradient is not
included in this version.

39
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 40

A.1 linSatModNogg
function [satmod_c, satmod_dist] = lin_sat_wheels_mod(Is, ws, hw,
tau_w, Ts)
%
% procedure to get system matrices of satellite with
momentum wheels
% and no gravity gradient
%
% Terms:
% SBC Satellite
Body Coordinate system
%
% (C) copyright Mogens Blanke 2001-2005
%
version 2.01 November 2005
%
% \qquad input:
% Is satellite inertia
matrix in SBC real(3,3)
% ws satellite’s steaty state angular
velocity (SBC frame) real(3,1)
% hw wheel momentum in SBC frame
real(3,1)
%
% tau_w wheel loop timeconstant scalar(1,1)
[s]
% Include w_bw if you want to calculate the model with wheel
%
dynamics included
% Ts optional sampling time to get a
discrete model
% \qquad\qquad\qquad\qquad the discrete model is calculated
if the Ts parameter is specified
%
% The last two parameters can be ommitted
from the call.
%
% \qquad output:
% \qquad satmod_c statespace model
(Matlab LTI format):
% \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
dx(t)/dt =\qquad A*x(t) + B_u*N_c(t)
% \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 41

y(t) \qquad\qquad=\qquad C*x(t) + D_u*N_c(t) \qquad\qquad

%
% \qquad satmod_c
%
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad dx(t)/dt =\qquad A*x(t)
+ B_d*N_d(t)
% \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad y(t) \qquad\qquad=\qquad
C*x(t) + D_d*N_d(t) \qquad\qquad
%
% where N_c is control torque demand
to reaction wheel(s)
% and \qquad N_d is external disturbance torque on
satellite

Aw(:,1) = [0, Is(3,1), Is(2,1); ...

-2*Is(3,1), -Is(3,2),
-Is(3,3)+Is(1,1);...
2*Is(2,1), Is(2,2)-Is(1,1), Is(2,3)]*ws;

Aw(:,2)
= [Is(3,1), 2*Is(3,2), (Is(3,3)-Is(2,2));...
-Is(3,2), 0, Is(1,2);...

Is(2,2)-Is(1,1), 2*Is(1,2),-Is(1,3)]*ws;

Aw(:,3)= [-Is(2,1),Is(3,3)-Is(2,2),-2*Is(2,3);...

Is(1,1)-Is(3,3), Is(1,2), 2*Is(1,3);...

Is(2,3), -Is(1,3), 0]*ws
;

invrtIs = inv(Is);

Ahw = [0, -hw(3), hw(2);...

hw(3), 0, -hw(1);
...
-hw(2), hw(1), 0];

Awh = [0, ws(3), -ws(2);...

-ws(3), 0,
ws(1); ...
ws(2), -ws(1), 0];
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 42

Aww = Aw + Ahw;

A = [invrtIs*Aww,
zeros(3,3), invrtIs*Awh;...
0.5*eye(3,3), zeros(3,3), zeros(3,3);...

zeros(3,3), zeros(3,3), zeros(3,3)];

Bu = [invrtIs; zeros(3,3); -eye(3,3)];

Bd = [invrtIs; zeros(3,3);
zeros(3,3)];

C = [eye(6,6), zeros(6,3)]; D
= zeros(6,3)

satmod_c = ss(A,Bu,C,D); satmod_dist = ss(A,Bd,C,D);

%
case with wheel dynamics included (12 states in state vector)
if (nargin
== 4 | nargin == 5)
A = [invrtIs*Aww, zeros(3,3), invrtIs*Awh, -invrtIs;...

0.5*eye(3,3), zeros(3,3), zeros(3,3), zeros(3,3);...

zeros(3,3),
zeros(3,3), zeros(3,3), eye(3,3);...
zeros(3,3), zeros(3,3), zeros(3,3),
-eye(3,3)*1/tau_w];

Bu = [zeros(3,3); zeros(3,3); zeros(3,3); -eye(3,3)*1/tau_w];

Bd =
[invrtIs; zeros(3,3); zeros(3,3); zeros(3,3)];

C = [eye(6,6), zeros(6,6)];
D = zeros(6,3);

satmod_c = ss(A,Bu,C,D); satmod_dist = ss(A,Bd,C,D);

end

if
(nargin == 5)
satmod_c = SS(A, Bu, C, D, Ts);
end
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 43

A.2 romerEx
function romer_ex
% roemer-like example
% MB 9. nov 2009
%

Is = [14.6 0 0; 0 13.6 0; 0 0 4.6];

hs0 = [0; 0; 0];
hs1 = [0; 0; 0.5];
ws0 = [0; 0; 0];
ws1 = [0; -0.0030; 0];

%
% parameters for plot
re_min = -0.25; re_max = 0.05; im_min=-0.1; im_max = 0.1;
re_axis = 0.9*[re_min, re_max];
im_axis = 0.9*[im_min, im_max];
re_min_zoom = -0.01; re_max_zoom = 0.01;
im_min_zoom=-0.01; im_max_zoom = 0.01;
re_axis_zoom = 0.9*[re_min_zoom, re_max_zoom];
im_axis_zoom = 0.9*[im_min_zoom, im_max_zoom];

limplotvec = [re_min re_max im_min im_max];

limplotveczoom = [re_min_zoom re_max_zoom im_min_zoom im_max_zoom];

% run open and closed loop for ang. velocity wsy and zero angular
% momentum
[satmod_Nw2y_open, satmod_Nd2y_open] = lin_sat_wheels_mod(Is, ws0, hs0);
[P,Z] = pzmap(satmod_Nw2y_open)
figure(1)
pole_zero_plot(P,Z,re_axis, im_axis, limplotvec, ...
’pole zero plot open loop’, ’romex_ws0_hs0_open.eps’)
pole_zero_plot(P,Z,re_axis_zoom, im_axis_zoom,limplotveczoom, ...
’pole zero plot open loop’, ’romex_ws0_hs0_open_zoom.eps’)

[A,B,C,D] = ssdata(satmod_Nw2y_open);
om_bw = 0.2;
Kpd = 2*Is*om_bw;
Kp = Is*om_bw^2;
Ky = [Kpd, 2*Kp];
Kr = -2*[Kp];
D_cl = zeros(size(C,1),size(Kr,2));
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 44

satmod_ref2y_closed = ss(A-B*Ky*C, B*Kr, C, D_cl );

[P,Z] = pzmap(satmod_ref2y_closed)
figure(2)
pole_zero_plot(P,Z,re_axis, im_axis, limplotvec, ...
’pole zero plot closed loop’, ’romex_ws0_hs0_closed.eps’)
pole_zero_plot(P,Z,re_axis_zoom, im_axis_zoom,limplotveczoom, ...
’pole zero plot closed loop’, ’romex_ws0_hs0_closed_zoom.eps’)

% run open and closed loop for ang. velocity wsy and zero angular
% momentum
[satmod_Nw2y_open, satmod_Nd2y_open] = lin_sat_wheels_mod(Is, ws1, hs0);
[P,Z] = pzmap(satmod_Nw2y_open)
figure(3)
pole_zero_plot(P,Z,re_axis, im_axis, limplotvec, ...
’pole zero plot open loop’, ’romex_ws1_hs0_open.eps’)
pole_zero_plot(P,Z,re_axis_zoom, im_axis_zoom,limplotveczoom, ...
’pole zero plot open loop’, ’romex_ws1_hs0_open_zoom.eps’)

[A,B,C,D] = ssdata(satmod_Nw2y_open);
om_bw = 0.2;
Kpd = 2*Is*om_bw;
Kp = Is*om_bw^2;
Ky = [Kpd, 2*Kp];
Kr = -2*[Kp];
D_cl = zeros(size(C,1),size(Kr,2));
satmod_ref2y_closed = ss(A-B*Ky*C, B*Kr, C, D_cl );
[P,Z] = pzmap(satmod_ref2y_closed)
figure(4)
pole_zero_plot(P,Z,re_axis, im_axis, limplotvec, ...
’pole zero plot closed loop’, ’romex_ws1_hs0_closed.eps’)
pole_zero_plot(P,Z,re_axis_zoom, im_axis_zoom,limplotveczoom, ...
’pole zero plot closed loop’, ’romex_ws1_hs0_closed_zoom.eps’)

% run the example for open and closed loop with wsy and hs1 (angular
% momentum bias)
[satmod_Nw2y_open, satmod_Nd2y_open] = lin_sat_wheels_mod(Is, ws1, hs1);
[P,Z] = pzmap(satmod_Nw2y_open)
figure(5)
pole_zero_plot(P,Z,re_axis, im_axis, limplotvec, ...
’pole zero plot open loop’, ’romex_ws1_hs1_open.eps’)
pole_zero_plot(P,Z,re_axis_zoom, im_axis_zoom,limplotveczoom, ...
’pole zero plot open loop’, ’romex_ws1_hs1_open_zoom.eps’)

om_bw = 0.2;
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 45

Kpd = 2*Is*om_bw;
Kp = Is*om_bw^2;
Ky = [Kpd, 2*Kp];
Kr = -2*[Kp];
D_cl = zeros(size(C,1),size(Kr,2));
satmod_ref2y_closed = ss(A-B*Ky*C, B*Kr, C, D_cl );
[P,Z] = pzmap(satmod_ref2y_closed)
figure(6)
pole_zero_plot(P,Z,re_axis, im_axis, limplotvec, ...
’pole zero plot closed loop’, ’romex_ws1_hs1_closed.eps’)
pole_zero_plot(P,Z,re_axis_zoom, im_axis_zoom,limplotveczoom, ...
’pole zero plot closed loop’, ’romex_ws1_hs1_closed_zoom.eps’)

end

function pole_zero_plot(P, Z, re_axis, im_axis, limplotvec, plottitle, saveplotas)

zeroataxis = [0, 0];

plot( re_axis,zeroataxis,’k-’,zeroataxis,im_axis,’k-’,’LineWidth’,2), hold on
plot(real(P),imag(P),’r+’,real(Z),imag(Z),’bo’,’LineWidth’,2,’MarkerSize’,16), grid;
axis(limplotvec); hold off
title(plottitle,’FontSize’,14),
xlabel(’real’,’FontSize’,12), ylabel(’imag’,’FontSize’,12)
eval([’print -depsc2 -tiff ’,saveplotas])
end
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 46

A.3 SATsimstart
% script file to test Simulink satellite simulation
%
(C) copyright Mogens Blanke 2001-2003
% version 03a December 2003
%
global
fpc_gain;
global fpc_flag;
global fpc_param;
global fpc_state;
global time;

%
initialize all simulation parameters
SAT_sim; RWA_sim;

% FPC_sim;
%
AD_sim;
% start simulation % in front of the line you don’t want !

sat3_sim

A.4 RWAsim
% RWAsim
% This script initializes your wheel model.
%
It enables different number of wheels and orientations by
% selecting
the parameters Awar, w-w_ini and N_dw.
% Number of columns in Awar MUST
equal number of wheels
% Wheel initial velocity is specified in w_w_ini
%
N_dw specifies the torque demand of each wheel
% n_w specifies the number
of wheels
%
% Note: # columns in Awar == dim(w_w_ini) == dim(N_dw) ==
n_w
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 47

%
% (C) Copyright 2001-2003 Mogens Blanke, DTU
% version 03a December
2003
%
% Parameters for each wheel

wheel_inertia = 0.0010;% kg m^2

tau_w
= 0.25;
wheel_loop_bandwidth = 1/tau_w; % rad/s
wheel_infspeed = 500;%
rad/s
wheel_inftorque = 0.050; %max torque on wheel
wheel_kqn = wheel_loop_bandwidth
* wheel_inertia;

% Example 1: 4 wheels in tetraeder orientation

n_w =
4; % number of wheels
a = sqrt(1/3); b =sqrt(2/3);
Awar = [a b 0; a -b
0; -a 0 -b; -a 0 b]’;
w_w_ini = [100, 100, 100, 100];

% Example 2: 3 wheels
in orthogonal orientation
n_w = 3; a = 1;
Awar = [a 0 0; 0 a 0; 0 0 a ]’;
w_w_ini
= [100, 100, 100];

% Example 2: 1 wheel in direction 45 deg of satellite’s

x-y plane
n_w = 1;
a = sqrt(1/2);
Awar = [a a 0]’;
w_w_ini = [100];

%
Remember to put % in front of lines for the cases you don’t use
APPENDIX A. PROGRAMME LISTINGS AND CONTENTS CD 48

A.5 SATsim
% SATsim
% script to initialize satellite main body
parameters

Iscb = diag([3,14,15]);
Tsamp = 1.0;