Textbook PDF

Chapter 1
Rigid Body dynamics
In order to describe the attitude of a rigid body and to determine its evolution as a
function of its initial angular velocity and applied torques, Euler’s angles and Euler’s
equations of motion need to be introduced. The transformation matrix between different
reference frames will be recalled and the concept of inertia tensor will also be briefly
discussed.
1.1 Frames of reference and transformation matrices

Assuming that a satellite is a rigid body is a reasonable initial model for attitude dy-
namics and control. However, in practice, this assumption can only be used as a first
approximation. For satellites with large deployable solar arrays the structure can be
quite flexible. The elastic modes in the structure can be excited through attitude con-
trol thrusters firings. This leads to vibrations which reduce the pointing accuracy of the
payload. In addition, fuel consumption and fuel slosh in propellant tanks can cause the
inertia properties of the satellite to be time varying, leading to a more complex control
problem. But if we assume that our spacecraft is a rigid body, we can attach to it a body
frame, FB , described by a set of unit vectors (ê1 , ê2 , ê3 ). The position of FB with respect
to an inertial reference frame FI , identified by the unit vectors (Ê 1 , Ê 2 , Ê 3 ), completely
describes the attitude of our spacecraft.
Assuming that ~v is a vector quantity, it is possible to write it as
~v = xê1 + yê2 + zê3

or, equivalently, ~v = X Ê 1 + Y Ê 2 + Z Ê 3
The column vectors v B = (x, y, z)T and v I = (X, Y, Z)T provide the component repre-
sentations of the same vector quantity ~v in the reference frames FB and FI , respectively.
If we now consider the components êiI = (e1,i , e2,i , e3,i )T of the i–th unit vector êi in
FI , that is
êi = e1,i Ê 1 + e2,i Ê 2 + e3,i Ê 3
we can write
~v = xê1 + yê2 + zê3

= x(e1,1 Ê 1 + e2,1 Ê 2 + e3,1 Ê 3 ) +
+ y(e1,2 Ê 1 + e2,2 Ê 2 + e3,2 Ê 3 ) +
+ z(e1,3 Ê 1 + e2,3 Ê 2 + e3,3 Ê 3 )
1
2 CHAPTER 1. RIGID BODY DYNAMICS
~v = (e1,1 x + e1,2 y + e1,3 z)Ê 1 +

+ (e2,1 x + e2,2 y + e2,3 z)Ê 2 +
+ (e3,1 x + e3,2 y + e3,3 z)Ê 3
This means that the components of ~v in FI can be expressed as a function of those in

FB as follows:
X=e1,1 x + e1,2 y + e1,3 z
Y =e2,1 x + e2,2 y + e2,3 z
Z=e3,1 x + e3,2 y + e3,3 z
or, in compact matrix form,
v I = LIB v B
where the transformation matrix LIB is given by
 
e1,1 e1,2 e1,3 · ¸
  .. ..
LIB = e2,1 e2,2 e2,3 = ê1I .ê2I .ê3I
e3,1 e3,2 e3,3
LIB is made up by the components of the unit vectors êi as expressed in FI . Any matrix
made up by mutually orthogonal row or column unit vectors is an orthogonal matrix and
is characterized by several properties, among which we only recall that:
• the inverse of an orthogonal matrix L is given by its transpose: L−1 = LT ;
• the determinant of an orthogonal matrix is det(L) = ±1 (and that of a rotation

matrix is 1);
• if L1 and L2 are orthogonal matrices, their product L1 L2 is an orthogonal matrix.
Thanks to the first property, it is possible to write the inverse coordinate transforma-
tion as
v B = LBI v I = (LIB )−1 v I = (LIB )T v I
which means that it is also  
êT1I
 ... 
 
LBI =
 êT2I 

 ... 
êT3I
As an exercise, demonstrate that the dual relations
 T 
Ê 1B
· ¸  ... 
. .  
 T 
LBI = Ê 1B ..Ê 2B ..Ê 3B ; LIB =  Ê 2B 
 
 ... 
T
Ê 3B
also hold.
1.2. EULER’S ANGLES 3
1.2 Euler’s angles

It is possible to use the coordinate transformation matrix LBI to describe the attitude
of the spacecraft through the unit vectors êi of the body frame attached to it, coming
out with a total of 9 parameters. As a matter of fact, these 9 parameters are not free to
vary at will, inasmuch as they must satisfy 6 constraints, expressed by the orthonormality
condition, that is ½
0 if i 6= j
êi · êj = δi,j = (1.1)
1 if i = j
Roughly speaking, only 9 − 6 = 3 parameters should be sufficient to describe the attitude
of FB w.r.t. FI .
One of the set of three parameters most widely used to describe the attitude of a rigid
body (or equivalently the attitude of the body frame attached to it) w.r.t. a fixed frame
are the Euler’s angles, a sequence of three rotations that take the fixed frame and make
it coincide with the body frame. The original sequence of rotations proposed by Euler to
superimpose FI onto FB is the sequences 3-1-3:
1. the first rotation is about the third axis of the initial frame, that is Ê3 , in our case,
and takes the first axis Ê 1 to the direction ê01 perpendicular to the plane determined
by the unit vectors Ê 3 and ê3 ; Ê 2 is rotated onto ê02 ; the rotation angle is called
precession angle Ψ;
2. the second rotation is about the first axis transformed after the first rotation, ê01 ,
and takes the axis ê03 into the position of ê3 ; ê02 is moved onto ê002 ; the rotation angle
is called nutation angle Θ;
3. the third and final rotation is about ê3 and brings ê001 = ê01 and ê002 to their final
positions, ê1 and ê2 , respectively; the rotation angle is called spin angle Φ.
The three angles, representing the amplitude of the three, successive rotations Ψ, Θ, Φ,
respectively about the third, the first, and again the third axis, can be used to represent
the attitude of the frame FB : The nutation angle represents the inclination of the third
body axis ê3 w.r.t. the local vertical Ê 3 ; The precession angle represents the angle between
the first inertial axis Ê 1 and the line of the nodes ξ, i.e. the intersection between the
planes perpendicular to ê3 and Ê 3 ; The spin angle is the rotation about the third body
axis.
The transformation matrix LBI can be expressed as a function of these three angles,
in terms of three elementary rotation matrices, as will be derived in the sequel.
Other sequences
It must be remembered that the sequence of rotations here described is not the only
possible choice for rotating FI onto FB . Many other sequences are available and equally
useful. In atmospheric flight mechanics the most widely used sequence of rotations is the
3–2–1, also known as the Bryant’s angles.
In this case the first rotation is about the third axis, Ê 3 , and its amplitude is called
yaw angle ψ. The second rotation about the second axis ê02 is the pitch angle, θ, and
takes the first axis onto its final position. The third rotation about ê1 is the roll angle,
φ. This set of angles is used also in space flight dynamics, to describe the attitude of a
spacecraft with respect to the Local Horizontal – Local Vertical (LHLV) reference frame.
Ê 3 ê03 ≡ Ê 3
ê02
Ψ
Ψ
Ê 1 Ê 2 Ê 1 Ê 2
ê01
ê2
Ê 3 ê002 Ê 3 ê002
ê003 ê3 ≡ ê003
Θ Θ Φ
Θ ê02 ê02
Φ ê1
Ψ Ψ
Ê 1 Ê 2 Ê 1 Ê 2
ê01 ≡ ê001 ê01
Figure 1.1: Euler’s angles
In many textbooks also this latter set of rotations is often referred to as Euler’s angles,
and this fact may lead to some confusion.
As a final observation, the order of the rotation sequence is important: Rotations do
not commute! This means that the rotation sequence 1–2–3, performed with the same
angles about the same axis, will take the initial frame to another one. The rotation
sequence 1–2–3 is known as Cardan angles.
Singularity
There is another problem with the representation of rotations in a three dimensional space,
that is the singularity of all the descriptions in terms of three parameters. This means
that there will always be positions of the two frames that can be described in different
ways, once a particular sequence of rotations is chosen. As an example, if the original
Euler’s angle sequence is employed, the case in which Θ = 0 is singular, inasmuch as the
precession and spin rotations will be about the very same axis, Ê 3 ≡ ê3 . This means
that all the triplets (Ψ, 0, Φ) for which Ψ + Φ is constant represent the same change of
reference frame.
Similarly, when the Bryant’s angles are used, the case θ = ±π/2 is singular, as in this
case all the triplets (ψ, ±π/2, φ) for which ψ − φ is constant will provide the same final
attitude for FB .
The problem of coordinate transformation singularity has some unpleasant mathemat-
ical consequence that will be underlined in the following paragraphs.
y 1
y 2
y1
v
x 2
x1 x 1
Figure 1.2: Planar rotation of amplitude α.
1.2.1 Building the coordinate transformation matrix from ele-

mentary rotations
Planar rotations
Consider the sketch of Fig. 1.2, where two planar reference frames F1 and F2 , with the
same origin O are represented. The angle α, assumed positive for counter–clockwise
rotations, allows one to identify univocally the position of the axes X2 –Y2 of F2 w.r.t. the
frame F1 defined by the axes X1 and Y1 .
Given the components x1 and y1 of a vector ~v expressed in F1 , the components x2 e
y2 can be expressed as a function of the angle α. The following relations can be easily
inferred from Fig. 1.2:
x2 = x1 cos(α) + y1 sin(α)
y2 = −x1 sin(α) + y1 cos(α)
or, in matrix form,

½ ¾ · ¸½ ¾
x2 cos(α) sin(α) x1
=
y2 − sin(α) cos(α) y1
This relation expresses the coordinate transformation that takes the component of a
vector quantity expressed in F1 into those of a reference frame F2 rotated w.r.t. F1 of an
angle α.
In compact notation we can write
v 2 = L21 v 1 = R(α)v 1
where the subscript near the vector indicates the frame in which the components of the
vector quantity are considered, the matrix L21 is the coordinate transformation matrix
from F1 to F2 that in the two dimensional case coincides with the elementary rotation
matrix R(α).
The inverse transformation from F2 to F1 is given by
v 1 = L12 v 2 = (R(α))−1 v 2
Recalling the properties of orthogonal matrices, it is
L12 = (R(α))−1 = (R(α))T = R(−α)
Elementary rotations for the sequence 3–1–3

Each one of the Euler’s rotations can be considered an elementary rotation about a given
axis, that remains unchanged during the transformation. It is still possible to apply
the relations derived for the planar case, adding a further equation that states that the
coordinate relative to the rotation axis does not vary.
The coordinate transformation during the first rotation is given by
x0 = X cos(Ψ) + Y sin(Ψ)
y 0 = −X sin(Ψ) + Y cos(Ψ)
z0 = Z
that, in matrix form, can be written as:

 0    
 x  cos(Ψ) sin(Ψ) 0  X 
y0 =  − sin(Ψ) cos(Ψ) 0  Y
 0   
z 0 0 1 Z
In an analogous way it is possible to demonstrate that, during the second rotation

about ê01 , the coordinate transformation is given by
x00 = x0
y 00 = y 0 cos(Θ) + z 0 sin(Θ)
z 00 = −y 0 sin(Θ) + z 0 cos(Θ)
that in matrix form becomes:

 00    0 
 x  1 0 0  x 
y 00 
= 0 cos(Θ) sin(Θ)  y0
 00   0 
z 0 − sin(Θ) cos(Θ) z
Finally, the third rotation about ê003 is represented by the transformation
x = x00 cos(Φ) + y 00 sin(Φ)

y = −x sin(Φ) + y 00 cos(Φ)
z = z 00
or, in matrix form:

    
 x  cos(Φ) sin(Φ) 0  x00 
y =  − sin(Φ) cos(Φ) 0  y 00
   00 
z 0 0 1 z
The three elementary rotation matrices of the Euler’s sequence 3–1–3 can thus be
defined as
   
cos(Ψ) sin(Ψ) 0 1 0 0

R3 (Ψ) = − sin(Ψ) cos(Ψ) 0  
; R1 (Θ) = 0 cos(Θ) sin(Θ)  ;
0 0 1 0 − sin(Θ) cos(Θ)
 
cos(Φ) sin(Φ) 0

R3 (Φ) = − sin(Φ) cos(Φ) 0 
0 0 1
where the subscript near the rotation matrix symbol R indicates the axis around which
the rotation is performed, while the argument indicates the amplitude of the rotation.
Summing up...
When passing from the inertial frame FI to the body frame FB using Euler’s sequence, the
coordinate transformation of vector quantities can be obtained combining in the correct
order the elementary rotation matrices, as follows:
v 0 = R3 (Ψ)v I
v 00 = R1 (Θ)v 0
v B = R3 (Φ)v 00
that is
v B = R3 (Φ)R1 (Θ)R3 (Ψ)v I
This means that
LBI = R3 (Φ)R1 (Θ)R3 (Ψ)
Performing the row–column products, the following expression for LBI is obtained:
 
cos Φ cos Ψ sin Φ cos Θ cos Ψ sin Φ sin Θ
 − sin Φ cos Θ sin Ψ + cos Φ sin Ψ 
 
 
LBI = 
 − cos Φ cos Θ sin Ψ cos Φ cos Θ cos Ψ cos Φ sin Θ 

 − sin Φ cos Ψ − sin Φ sin Ψ 
 
sin Θ sin Ψ − sin Θ cos Ψ cos Θ
As the product of orthogonal matrices is an orthogonal matrix, the inverse of which is

equal to its transpose, the inverse coordinate transformation matrix LIB is simply given
by
 
cos Φ cos Ψ − cos Φ cos Θ sin Ψ sin Θ sin Ψ
 − sin Φ cos Θ sin Ψ − sin Φ cos Ψ 
 
 
LIB = LBI −1 = LBI T = 
 sin Φ cos Θ cos Ψ cos Φ cos Θ cos Ψ − sin Θ cos Ψ 

 + cos Φ sin Ψ − sin Φ sin Ψ 
 
cos Φ sin Θ cos Φ sin Θ cos Θ
Elementary rotations for the sequence 3–2–1

It is left as an exercise to the reader the composition of elementary rotation matrices for
the sequence 3–2–1. Adopting the same notation used above, it is
   
cos(ψ) sin(ψ) 0 cos(θ) 0 − sin(θ)
R3 (ψ) =  − sin(ψ) cos(ψ) 0  ; R2 (θ) =  0 1 0  ;
0 0 1 sin(θ) 0 cos(θ)
 
1 0 0
R1 (φ) =  0 cos(φ) sin(φ) 
0 − sin(φ) cos(φ)
and the final result is given by

 
cos θ cos ψ cos θ sin ψ − sin θ
 
 sin φ sin θ cos ψ sin φ sin θ sin ψ sin φ cos θ 
 
LBI = 
 − cos φ sin ψ + cos φ cos ψ 

 
 cos φ sin θ cos ψ cos φ sin θ sin ψ cos φ cos θ 
+ sin φ sin ψ − sin φ cos ψ
A first consequence of the Euler’s angle singularity

When Θ = 0, the coordinate transformation matrix does not depend on Ψ and Φ sepa-
rately, but only on their sum. In such a case, it is
 
cos Φ cos Ψ sin Φ cos Ψ 0  
 − sin Φ sin Ψ + cos Φ sin Ψ  cos(Φ + Ψ) sin(Φ + Ψ) 0
 
   
LBI =   
 − cos Φ sin Ψ cos Φ cos Ψ 0  =  − sin(Φ + Ψ) cos(Φ + Ψ) 0 

 − sin Φ cos Ψ − sin Φ sin Ψ 
  0 0 1
0 0 1
As an exercise, demonstrate that LBI depends on ψ − φ only, when Bryant’s rotation

sequence is employed and θ = ±π/2.
How to build elementary rotation matrices

There is a simple way to build mnemonically the elementary rotation matrices. The
matrices are 3 by 3. If a rotation of an angle α about the i-th axis is being considered,
place 1 in position i, i, and fill the remaining elements of the i–th row and i–th column
with zeroes. All the other elements of the principal diagonal are cos α and the last two
outside the diagonal are sin α. The sin element above the row with the 1 must have a
minus sign. As an example, let us consider a rotation θ about the second axis (like in the
second rotation of the Bryant’s sequence). We start filling the matrix with 0s along the
second row and column, with a one in position 2,2:
 
· 0 ·
 0 1 0 
· 0 ·
ê2
Ê 3
ê002
Ψ̇
ê3 Φ̇
ê02
ê1
Ê 1 Ê 2
Θ̇
ê01
Figure 1.3: Angular velocity as a funciton of Euler’s angle rates.
Then we fill the diagonal with cos θ:

 
cos θ 0 ·
 0 1 0 
· 0 cos θ
and we put sin θ in the remaining places, with a minus sign in the row above the 1:
 
cos θ 0 − sin θ
 0 1 0 
sin θ 0 cos θ
This is R2 (θ). When a rotation about the first axis is considered, the 1 is on the first row
and apparently there is no row above it. But it is sufficient to cycle and start again from
the bottom: In this case the minus sign is on the sin in the third row.
1.2.2 Angular velocity and the evolution of Euler’s angles

~ is given by
The angular velocity ω
~ = ω1 ê1 + ω2 ê2 + ω3 ê3
ω
but it is also (see Fig. 1.4)
~ = Ψ̇Ê 3 + Θ̇ê01 + Φ̇ê3
ω
The components of the unit vector Ê 3 in FB are given by the third column of the
matrix LBI , that is E 3B = (sin Φ sin Θ, cos Φ sin Θ, cos Θ)T , while the components of ê01
are (cos Φ, − sin Φ, 0)T . Thus
~ =
ω Ψ̇(sin Φ sin Θê1 + cos Φ sin Θê2 + cos Θê3 ) +
+ Θ̇(cos Φê1 − sin Φê2 ) +
+ Φ̇ê3
= (Ψ̇ sin Φ sin Θ + Θ̇ cos Φ)ê1 +
+ (Ψ̇ cos Φ sin Θ − Θ̇ sin Φ)ê2
+ (Ψ̇ cos Θ + Φ̇)ê3
or, in matrix form

    
 ω1  sin Φ sin Θ cos Φ 0  Ψ̇ 
ω2 =  cos Φ sin Θ − sin Φ 0  Θ̇
   
ω3 cos Θ 0 1 Φ̇
Inverting the 3 × 3 matrix, one obtains the law of evolution of Euler’s angles as a
function of angular velocity components in body axis, that is
    
 Ψ̇  sin Φ/ sin Θ cos Φ/ sin Θ 0  ω1 
Θ̇ = cos Φ − sin Φ 0  ω2
   
Φ̇ − sin Φ/ tan Θ − cos Φ/ tan Θ 1 ω3
or, in explicit form,
Ψ̇ = (ω1 sin Φ + ω2 cos Φ)/ sin Θ

Θ̇ = ω1 cos Φ − ω2 sin Φ
Φ̇ = (−ω1 sin Φ − ω2 cos Φ)/ tan Θ + ω3
These equations can be integrated to obtain the evolution of the Euler angles, if the
angular velocity is known. But they also show an unpleasant feature of Euler’s angle
singularity, that is the spin and precession rates go to infinity when Θ approaches 0. This
fact has some serious consequences on the problem of attitude representation, inasmuch
as it is not possible to accept a set of attitude parameters the evolution of which cannot
always be described in an accurate way.
If the Bryant’s angles are used, the reader can demonstrate that
    
 ω1  0 − sin φ cos θ cos φ  ψ̇ 
ω2 =  0 cos φ cos θ sin φ  θ̇
   
ω3 1 0 − sin θ φ̇
and the inverse relation is
φ̇ = ω1 + (ω2 sin φ + ω2 cos φ) tan θ

θ̇ = ω2 cos φ − ω3 sin φ
ψ̇ = (ω2 sin φ + ω3 cos φ)/ cos θ
Again, in the neighborhood of the singular condition θ = ±π/2 the rate of change of the
roll and yaw angles goes to infinity.
1.2.3 The quaternions

Euler’s eigenaxis rotation theorem states that it is possible to rotate a fixed frame FI
onto any arbitrary frame FB with a simple rotation around an axis â that is fixed in both
frames, called the Euler’s rotation axis or eigenaxis, the direction cosines of which are
the same in the two considered frame.
A very simple algebraic demonstration of Euler’s theorem can be obtained from the
following considerations:
• the eigenvalues of any (real) orthogonal matrix L have unit modulus; indicating
with H the Hermitian conjugate, which, for a real matrix is coincident with the
transpose, one has
La = λa ⇒ aH LT La = λ̄λaH a ⇒ (1 − λ̄λ)aH a = 0
that for any nontrivial eigenvector a implies that
λ̄λ = 1 ⇒ |λ| = 1
• (at least) one eigenvalue is λ = 1; any n × n real matrix has at least one real
eigenvalue if n is an odd number, which means that a 3 × 3 orthogonal matrix must
have at least one eigenvalue which is λ1 = ±1. The other couple of eigenvalues will
be, in the most general case, complex conjugate numbers of unit modulus, which
can be cast in the form λ2,3 = exp(±iφ). The determinant is equal to the product
of the eigenvalues, which is one, for an orthogonal matrix, so that
λ1 λ2 λ3 = 1 ⇒ λ1 = 1
The eigenvector relative to the first eigenvalue satisfies the relation
La = 1 · a
This means that there is a direction a which is not changed under the action of transfor-
mation matrix L. If L represents a coordinate change, the vector a will be represented
by the same components in both the considered reference frames,
â = a1 ê1 + a2 ê2 + a3 ê3

= a1 Ê 1 + a2 Ê 2 + a3 Ê 3
For this reason, the transformation that takes the initial frame onto the final one can be
considered as a single rotation α about the Euler axis â.
In order to express the coordinate transformation matrix LBI as a function of α and
â it is sufficient to consider the following sequence of rotations:1
1. take the unit vector Ê 1 onto â, so that the new frame F 0 is given by the unit vectors
â, ê02 , ê03 ; call this rotation R̄;
2. rotate both frames FI and F 0 about the eigenaxis of the rotation angle α; because
of the definition of Euler rotation, FI goes onto FB , while F 0 will rotate into a
new frame F 00 given by the unit vectors â, ê002 , ê003 ; this rotation is represented by the
elementary rotation matrix R1 (α);
¯ that takes F 00 onto F has the
3. at this point it should be noted that the rotation R̄ B
¯ = R̄T .
same magnitude of R̄, but it is performed in the opposite direction so that R̄
Summing up it is
T
LBI = R̄ R1 (α)R̄
1
This derivation is taken from B. Wie, Space Vehicle Dynamics and Control, AIAA Education Series,
Reston (VA), USA, 1998, Chap. 5, pp. 312–315 and 318–320.
Ê 3 Ê 3 ê03
(a) ê3 (b)

â â
ê2
ê02
Ê 1 α Ê 2 Ê 1 Ê 2
ê1
ê003 Ê 3 ê03 ê003
(c) ê3 ê002 (d) ê3 ê002

â â
ê2 ê2
ê02
Ê 1 α Ê 2
ê1 ê1
where  
a1 a2 a3
R̄ =  R21 R22 R23 
R31 R32 R33
Carrying out the calculations, one get, for the components Lij of the coordinate transfor-
mation matrix LBI the following expressions:
L11 = a21 + (R21

2 2
+ R31 ) cos α
L12 = a1 a2 + (R21 R22 + R31 R32 ) cos α + (R21 R32 − R31 R22 ) sin α
L22 = a22 + (R22
2 2
+ R32 ) cos α
2
L33 = a23 + (R23 2
+ R33 ) cos α
Taking into account the orthogonality conditions for R̂, one gets
a21 + R21
2 2
+ R31 =1 ⇒ 2
R21 2
+ R31 = 1 − a21
2 2 2
a2 + R22 + R32 = 1 ⇒ 2 2
R22 + R32 = 1 − a22
a23 + R23
2 2
+ R33 =1 ⇒ 2
R23 2
+ R33 = 1 − a23
a1 a2 + R21 R22 + R31 R32 = 0 ⇒ R21 R22 + R31 R32 = −a1 a2
a2 a3 + R22 R23 + R32 R33 = 0 ⇒ R22 R23 + R32 R33 = −a2 a3
a3 a1 + R21 R23 + R31 R33 = 0 ⇒ R21 R23 + R31 R33 = −a1 a2
while remembering that the first row of an orthogonal matrix is given by the cross product
of the second and the third ones, it is
a1 = R22 R33 − R23 R32
a2 = R23 R31 − R21 R33
a3 = R21 R32 − R22 R31
Substituting these results into the expressions of the coefficients Lij the following
expression for LBI is obtained:
 
cos α + a21 (1 − cos α) a1 a2 (1 − cos α) + a3 sin α a1 a3 (1 − cos α) − a2 sin α
LBI =  a2 a1 (1 − cos α) − a3 sin α cos α + a22 (1 − cos α) a2 a3 (1 − cos α) + a1 sin α 
a3 a1 (1 − cos α) + a2 sin α a3 a2 (1 − cos α) − a1 sin α cos α + a23 (1 − cos α)
(1.2)
or, in compact matrix form
LBI = cos α1 + (1 − cos α)aaT − sin αÃ
where 1 is the 3 × 3 identity matrix and Ã is the cross product equivalent matrix form
 
0 −a3 a2
Ã =  a3 0 −a1 
−a2 a1 0
such that a × b = Ãb.
We now define the Euler parameters or quaternions as
q0 = cos(α/2)
q1 = a1 sin(α/2)
q2 = a2 sin(α/2)
q3 = a3 sin(α/2)
By letting q = â sin(α/2) and rembering that cos α = cos2 (α/2)−sin2 (α/2) = q02 −q·q and
sin α = 2 cos(α/2) sin(α/2) = 2q0 sin(α/2), it is easy to demonstrate that the coordinate
transformation matrix is given by
LBI = (q02 − q · q)1 + 2qq T − 2q0 Q̃
where the ˜ indicates again the cross product matrix equivalent
 
0 −q3 q2
Q̃ =  q3 0 −q1 
−q2 q1 0
1.2.4 Evolution of the quaternions

The evolution of the quaternions is described by the set of linear differential equations,
represented in matrix form as2
    

 q̇ 0 
 0 −ω1 −ω2 −ω3   q0 

  1  
q̇1 ω 0 ω3 −ω2 
 q1
= 

1

 q̇2  2 ω2 −ω3 0 ω1  
 q2 

   
q̇3 ω3 ω2 −ω1 0 q3
The equivalent matrix form is given by

1
q̇0 = − ω · q
2
1
q̇ = (q0 ω − ω × q)
2
1.2.5 Quaternions vs Euler angles

The quaternions allow for representing the attitude of a rigid body with several advantages
over Euler’s angles, above all the absence of inherent geometric singularity. Moreover, the
linear equation to be integrated in time in order to determine their evolution as a function
of angular velocity components is less computationally expensive than that derived for
the Euler’s angles. The price to pay is that 4 parameters are used, instead of only three,
that are not independent, inasmuch as they must satisfy the constraint
q02 + q · q = 1
Moreover, their geometric interpretation during an evolution is less immediate than that
of the Euler’s angles, the geometric meaning of which is intuitive. For this reason the
attitude of a satellite is often integrated in strapdown attitude determination systems in
terms of quaternions but then represented in terms of Euler angles.
1.3 Time derivative of vector quantities

If we consider a vector quantity in an inertially fixed reference frame FI ,
~v = X Ê 1 + Y Ê 2 + Z Ê 3
its time derivative is given simply by

d~v
= Ẋ Ê 1 + Ẏ Ê 2 + Ż Ê 3
dt
that is · ¸ n oT
dv
= Ẋ, Ẏ , Ż = v̇ I
dt I
When the same vector quantity ~v is expressed in terms of components in a moving ref-
erence frame FB , rotating with angular velocity ω~ BI = ω
~ with respect to FI , the time
derivatives of
~v = xê1 + yê2 + zê3
2
See above, pp. 326–328
1.4. EULER’S EQUATIONS OF MOTION OF A RIGID BODY 15
is given by3
d~v
~ × (xê1 ) +
= ẋê1 + ω
dt
~ × (yê2 ) +
+ ẏê2 + ω
~ × (zê3 )
+ żê3 + ω
This means that, in terms of vector components in FB , it is
· ¸
dv
= v̇ B + ω B × v B
dt B
where
v̇ B = {ẋ, ẏ, ż}T
ω B = {ω1 , ω2 , ω3 }T
1.4 Euler’s equations of motion of a rigid body

1.4.1 The inertia tensor
~ of a mass element δm, moving with velocity ~v is
The angular momentum δ h
~ = ~r × (δm~v )
δh
where ~r is the position vector of the mass, with respect of the pole used for the evaluation
of moments of vector quantities.
For an extended rigid body (Fig. 1.4), the total angular momentum is given by
Z
~ = (~r × ~v )δm
h
B
If the body is rotating around its center of mass, the velocity of every mass element is
~v = ω
~ × ~r
so that Z
~ =
h ω × ~r )] δm
[~r × (~
B
Expressing the vector quantities in body components as ω B = (ω1 , ω2 , ω3 )T and r B =
(x, y, z)T , the vector product ω
~ × ~r is given by
~ × ~r = (ω2 z − ω3 y)ê1 + (ω3 x − ω1 z)ê2 + (ω1 y − ω2 x)ê3
ω
Carrying on the calculations, the product ~r × (~ ω × ~r ) is
£ ¤
~r × (~
ω × ~r ) = (y 2 + z 2 )ω1 − (xy)ω2 − (xz)ω3 ê1 +
£ ¤
+ −(xy)ω1 + (x2 + z 2 )ω2 − (yz)ω3 ê2 +
£ ¤
+ −(xz)ω1 − (yz)ω2 + (x2 + y 2 )ω2 ê3
3
Remeber the Poisson formula for the time derivative of a unit vector,
dêi
~ × êi
=ω
dt
Ê 3
ê3
~v
dm ê2
~r
~r 0 CM
~
h
Ê 1 Ê 2
ê1 ~
ω
Figure 1.4: A rotating rigid body.
ω × ~r )] δm is strictly function of the mass

The integration over the body B of [~r × (~
distribution only, as angular velocity components are independent of body shape and
location. This means that, letting
~ = h1 ê1 + h2 ê2 + h3 ê3

h
it is
h1 = Ix ω1 − Ixy ω2 − Ixz ω3
h2 = −Ixy ω1 + Iy ω2 − Iyz ω3
h3 = −Ixz ω1 − Iyz ω2 + Iz ω3
where the moments of inertia Ix , Iy , Iz , and the products of inertia Ixy , Ixz , Iyz , are
Z Z Z
¡ 2 2
¢ ¡ 2 2
¢ ¡ 2 ¢
Ix = y + z δm ; Iy = x + z δm ; Iz = x + y 2 δm
B B
Z Z ZB
Ixy = (xy) δm ; Ixz = (xz) δm ; Iyz = (yz) δm
B B B
In matrix form the angular momentum components in body axes are given by
hB = Iω B
where the symmetric matrix

 
Ix −Ixy −Ixz
I =  −Ixy Iy −Iyz 
−Ixz −Iyz Iz
is the inertia matrix that represent the inertia tensor in body axes.
The same relations can be derived directly in a more compact vector form remembering
that, for the double vector product, the following relation holds:
~ × (~y × ~z ) = (~x · ~z ) ~y − (~x · ~y ) ~z

x
so that, in the present case, it is
~r × (~
ω × ~r ) = (~r · ~r ) ω
~ − (~r · ω
~ ) ~r
Taking into account the definition of the dyadic tensor
(~x~y ) ~z = (~y · ~z )~x
and the fact that the angular velocity vector is constant and can be taken out of the
integral, it is
Z
~ =
h [~r × (~ ω × ~r )] δm
B
µZ ¶
= [(~r · ~r ) − (~r~r )] δm ω~
B
~
= Ĩ ω
where Ĩ is the inertia tensor. Expressing ~r in terms of body components and integrating
over B the previous expression for the inertia matrix I is obtained.
Principal axes of inertia

The matrix I is real and symmetric, so its eigenvalues are real4 and its eigenvectors are
mutually orthogonal.5 This means that there exists a body reference frame FP such that
4
From the definition of eigenvalue and eigenvector of a complex matrix A, it is easy to derive the
following equation,
xH Ax
Ax = λx ⇒ λ = H
x x
If A is Hermitian (i.e. the linear operator represented by A is self–adjoint), it is
y H Ax = (Ay)H x
so that, remembering the properties of the hermitian inner product, such that xH y = (y H x), it is
xH Ax (Ax)H x xH Ax
λ̄ = = = =λ
xH x xH x xH x
which means that the eigenvalue λ is equal to its conjugate, i.e. it must be real.
5
Given two distinct eigenvalues λi 6= λj and their respective eigenvectors xi and xj , the following
relations hold:
Axi = λi xi
Axj = λj xj
Multiplying the first equation by xH H

j and the second one by xi , taking the complex conjugate of the
second and subtracting it from the first, one gets
xH H T T
j Axi − (xi Axj ) = λi xj xi − (λj xi xj )
Remembering that the eigenvalues are real, it is
xH H H
j Axi − (Axj ) xi = (λi − λj )xj xi
the inertia matrix is diagonal,

 
Jx 0 0
I =  0 Jy 0 
0 0 Jz
where the principal moment of inertia Jx , Jy , and Jz are the eigenvalues of I. The
eigenvectors are called principal axes.
Symmetries
If the mass distribution of the body B is characterized by symmetries, this property
reflects onto the inertia matrix I. As an example, if B has a plane of symmetry, one of the
principal axis will be perpendicular to the plane and the other two will lie on that plane.
If this case, the products of inertia relative to the axis perpendicular to the symmetry
plane will be zero. This case is typical of fixed wing aircraft, where the longitudinal
plane is approximately a symmetry plane. The y body axis, directed perpendicular to the
symmetry plane, is characterized by zero products of inertia, so that the inertia matrix
of an aircraft is typically equal to
 
Ix 0 −Ixz
I= 0 Iy 0 
−Ixz 0 Iz
If the body is axi-symmetric (or simply has a regular polygonal mass distribution
w.r.t. an axis σ), the symmetry axis σ is a principal axis of inertia, while any couple of
perpendicular axes on the plane Σ normal to σ will complete the set of principal axes.
In this case the principal moments of inertia relative to the axes perpendicular to the
symmetry axis will be equal. Assuming that σ = ê3 is a symmetry axis, the inertia
tensor becomes  
Jt 0 0
I =  0 Jt 0 
0 0 Js
where the subscripts t and s indicate the transverse and spin (or axial) moments of
inertia, respectively.
1.4.2 Kinetic energy

The rotational kinetic energy of a rigid body is given by
Z
1
T = (~v · ~v ) δm
2 B
Remembering that ~v = ω ~ × ~r , the argument of the integral becomes (~ω × ~r ) · (~

ω × ~r ).
Also, taking into account the permutation property of the scalar triple product
~ · (~y × ~z ) = ~y · (~z × x
x ~ ) = ~z · (~x × ~y )
Taking into account the definition of Hermitian operator the first term of the last equation is zero and so
the Hermitian product xH j xi is zero if λi 6= λj . Since both xj and xi are real, their Hermitian product
coincides with the scalar product, so that distinct eigenvectors are real and perpendicular.
one obtains the equivalence
ω × ~r ) · (~
(~ ω × ~r ) = ω
~ · [~r × (~
ω × ~r )]
Substituting this expression into the integral, and taking the (constant) angular velocity
out of the integration symbol, one gets
Z
1
T = ω ~ · [~r × (~ ω × ~r )] δm
2 B
that, remembering the definition of the inertia tensor, brings

1 ³ ´ 1
T = ω ~ = ω TB (Iω B )
~ · Ĩ ω
2 2
or, equivalently,
1 ~ = 1 ω T hB
T = ω~ ·h
2 2 B
1.4.3 Euler’s equation of motion

The second fundamental law of rigid body dynamics states that the time derivative of
the angular momentum is equal to the external torque applied to the body B. In vector
form, it is
dh~
=M ~
dt
Expressing the vector quantities in body axis components one gets
ḣB + ω B × hB = M B
If the inertia matrix I is constant, it is
I ω̇ B + ω B × (Iω B ) = M B
When a set of principal axes is chosen as the body axes, the inertia tensor is diagonal
and the Euler’s equation of motion for a rigid body are obtained:
Jx ω̇1 + (Jz − Jy )ω2 ω3 = M1

Jy ω̇2 + (Jx − Jz )ω3 ω1 = M2
Jz ω̇3 + (Jy − Jx )ω1 ω2 = M3
These equations can be integrated as a function of the applied torque to obtain the
time history of the angular velocity components. These, in turn, can be used to determine
the variation with time of the Euler’s angles (or of the quanternions), thus describing the
evolution of the rigid body attitude.
1.4.4 Conservation of angular momentum

Writing Euler’s equations in a set of principal axes such that (without loss of generality)
Jx > Jy > Jz , torque–free motion is described by the following set of ODEs,
Jx ω̇1 + (Jz − Jy )ω2 ω3 = 0

Jy ω̇2 + (Jx − Jz )ω3 ω1 = 0
Jz ω̇3 + (Jy − Jx )ω1 ω2 = 0
It is easy to demonstrate that the magnitude h of the angular momentum vector
~ = h1 ê1 + h2 ê2 + h3 ê3

h
= Jx ω1 ê1 + Jy ω2 ê2 + Jz ω3 ê3
is constant when a torque–free motion is considered. A first intuitive derivation of this

property is that if the applied torque vanishes the angular momentum vector is constant in
FI , and its magnitude is independent of the considered reference system. It is also possible
to demonstrate analytically that h = ||h|| is constant, by taking the time derivatives of
h2 = (Jx ω1 )2 + (Jy ω2 )2 + (Jz ω3 )2
in the hypothesis of torque–free motion (M1 = M2 = M3 = 0),
dh2 ¡ ¢
= 2 Jx2 ω1 ω̇1 + Jy2 ω2 ω̇2 + Jz2 ω3 ω̇3
dt
= 2 [Jx ω1 (Jy − Jz )ω2 ω3 + Jy ω2 (Jz − Jx )ω1 ω3 + Jz ω3 (Jx − Jy )ω1 ω2 ]
= 2 (Jx Jy − Jx Jz + Jy Jz − Jy Jx + Jz Jx − Jz Jy ) ω1 ω2 ω3 = 0
Geometrically, the angular velocity vector must lie on an ellipsoid (the angular momentum
ellipsoid ) in FB , the equation of which takes the form
ω12 ω22 ω32

+ + =1
(h/Jx )2 (h/Jy )2 (h/Jz )2
1.4.5 Conservation of kinetic energy

In a similar way, it is also possible to demonstrate that the kinetic energy of a rigid body
is constant if no external torque is applied. Again, taking the time derivative of
1 1¡ ¢
T = ωB · h = Jx ω12 + Jy ω22 + Jz ω32
2 2
it is
dT
= Jx ω1 ω̇1 + Jy ω2 ω̇2 + Jz ω3 ω̇3
dt
= ω1 (Jy − Jz )ω2 ω3 + ω2 (Jz − Jx )ω1 ω3 + ω3 (Jx − Jy )ω1 ω2
= 2 (Jy − Jz + Jz − Jx + Jx − Jy ) ω1 ω2 ω3 = 0
This means that the angular velocity vector must satisfy also the equation
ω12 ω22 ω32

+ + =1
(2T /Jx ) (2T /Jy ) (2T /Jz )
that is, it must point the surface of the kinetic energy (or Poinsot) ellipsoid in FB . The
combination of these two last results demonstrate that the angular velocity vector describe
the curve given by the intersection of the angular momentum ellipsoid and the kinetic
energy ellipsoid, which is called the polhode.
1.5. GENERALIZED EULER EQUATIONS 21
1.5 Generalized Euler equations

In their original formulation, Euler equations are written in a body–fixed reference frame
FB with the origin in the center of mass O of the body B. On one side, the expression
employed for the angular momentum h of B requires that B is rigid, and centering B in
O greatly simplify the expression. At the same time, linear and angular momentum con-
servation laws do apply to any mechanical system (under sufficiently mild assumptions),
so that it is possible to obtain a generalized formulation for the equation of motion in a
frame which is non centered in CM.
The classical equations, referred to the center of mass O, are
m~aO = F~
~O
dh
= M~O
dt
where, m is the mass of B and F~ is the external force, producing an acceleration ~aO of
the center of mass. Considering an arbitrary point A with arbitrary motion, where ~r A is
the position vector of O with respect to A, momentum and torque can be referred to A,
~ O + m~r A × d~r A
~A = h
h
dt
~ ~ ~
M A = M O + ~r A × F
while the acceleration of O can be rewritten as a function of the absolute acceleration of

A
d2~r A
~aO = ~aA +
dt2
By substituting the above expressions in the angular momentum equation, one gets
µ ¶ · µ 2
¶¸
d ~ d~r A ~ A − ~r A × m ~aA + d ~r A
hA − m~r A × =M
dt dt dt2
which can be rewritten as

~A
dh ~ A × ~aA = M
+S ~ A
dt
~ A = m~r A is the static moment of the body with respect to A.
where S
Chapter 2
Passive Stabilisation of Rigid

Spacecraft
Spin stabilisation in a simple, low cost and effective means of attitude stabilisation. Prior
to, or just after deployment, the satellite in spun up about its axis of symmetry. For
this reason, spin stabilised satellites are usually short cylinders. The angular momen-
tum accumulated about the spin axis the provides “gyroscopic stability” against external
disturbance torques.
Although simple and reliable, spin stabilised satellites are inefficient for power gener-
ation. Since the satellite is continually spinning, the entire surface of the satellite must
be covered with solar cells. In addition, payload efficiency is particularly low when only
one direction is fixed in space and maneuvering of the spin axis complex.
When the requirement on pointing accuracy is weak (of the order of some tenths of
a degree) gravity–gradient torque may be used for stabilizing an Earth pointing satellite,
while the use of a dual–spin system allows one to despin a part of the satellite, while
providing gyroscopic stability to the whole spacecraft.
2.1 Torque–free motion of axi–symmetric satellites

The principal moments of inertia of an axi–symmetric satellite will be given by
Jt = Jx = Jy
Js = Jz
where the subscripts t and s stand for transverse and spin, respectively, and we assume
that the symmetry axis coincides with the third (ê3 ) axis of the body frame FB .
During the spin–up manoeuver, the satellite will accumulate angular momentum about
the spin axis, but because of various perturbations or imperfections, such as thruster
misalignment, the final condition at the end of the spin–up will hardly be a pure spin
about the spin axis ê3 . The imperfections will cause some (hopefully residual) nutation.
For torque–free motions
M1 = M2 = M3 = 0
of axial symmetric spacecraft, Euler’s equations take the following form:
Jt ω̇1 + (Js − Jt )ω2 ω3 = 0
Jt ω̇2 + (Jt − Js )ω1 ω3 = 0
Js ω̇3 = 0
23
24 CHAPTER 2. PASSIVE STABILISATION OF RIGID SPACECRAFT
The first two equations are coupled, while the third one is independent of the other
two. This means that the latter one can be integrated on its own. The resulting (trivial)
solution is given by
ω3 = Ω
where Ω is the (constant) spin rate about the spin axis.
Letting
Js − Jt
λ= Ω
Jt
the first two equations can be rewritten as
ω̇1 + λω2 = 0
ω̇2 − λω1 = 0
Multiplying the first equation by ω1 and the second by ω2 and summnig up, one gets
ω̇1 ω1 + ω̇2 ω2 = 0
that is
ω12 + ω22 = ω12
2
= constant
~ that lies in the ê1 − ê2
This means that the component of the angular velocity vector ω
plane, namely
ω~ 12 = ω1 ê1 + ω2 ê2
has a constant magnitude. As also ω3 is constant, we get that
ω || = ω12 + ω22 + ω32 = ω12

||~ 2
+ Ω2 = constant
The first two equations of motion,
ω̇1 + λω2 = 0
ω̇2 − λω1 = 0
can be easily integrated. In fact, deriving the first one with respect to the time t, one gets
ω̈1 + λω̇2 = 0
that, taking into account the second equation, becomes
ω̈1 + λ2 ω1 = 0
which is formally identical to the well known equation of the linear harmonic oscillator.
The general solution
ω1 (t) = A cos(λt) + B sin(λt)
for initial conditions
ω1 (t = 0) = ω1,0 ; ω̇1 (t = 0) = ω̇1,0
becomes
ω̇1 (0)
ω1 (t) = ω1 (0) cos(λt) + sin(λt)
λ
= ω12 sin[λ(t − t0 )]
2.1. TORQUE–FREE MOTION OF AXI–SYMMETRIC SATELLITES 25
0.1
ω1 [rad s−1]
0.05
0
−0.05
−0.1
0.1
ω2 [rad s−1]
0.05
0
−0.05
−0.1
0.15
ω3 [rad s−1]
0.1
0.05
0
0 0.5 1
t [103 s]
Figure 2.1: Time–history of angular velocity components for a free spin condition.
Deriving the solution for ω1 w.r.t. t and substituting into the first equation, it is
ω̇1 ω̇1 (0)

ω2 = − = ω1 (0) sin(λt) − cos(λt)
λ λ
= −ω12 cos[λ(t − t0 )]
The evolution of ω1 and ω2 shows that ω~12 whirls around ê3 with angular velocity λ.
Writing the angular velocity as
~ =ω
ω ~ 12 + Ωê3
~ describes a cone around the axis of symmetry eˆ3 of the spinning
during the evolution, ω
body, which is called the body cone.
The angular momentum vector can be written in the form
~ = J1 ω1 ê1 + J2 ω2 ê2 + J3 ω3 ê3
h
= Jt (ω1 ê1 + ω2 ê2 ) + Js ω3 ê3
~ 12 + Js Ωê3
= Jt ω
~ and ω
It can be observed that h ~ are both a linear combination of the vectors ω~ 12 and
~ ~ and ê3 lie in the same
ê3 . Thus, during the motion of the spinning body, the vectors h, ω
plane Π, that rotates around ê3 , if we look at the motion from FB .
In the most general case h~ and ω~ are not aligned. They are aligned only if ω12 = 0,
that is, if we have a pure spinning motion about the symmetry axis. This is the desired
spin condition, where a sensor placed on the satellite on its symmetry axis points a fixed
direction in space. If ω12 6= 0, the motion of the spinning body is more complex (and the
pointing less accurate). A geometric description of this motion will now be derived.
The assumption of torque–free motion, M ~ = 0, implies that
~
dh ~ = constant
~ = 0 =⇒ h
=M
dt
θ PP
P
P
q
¡
¡
ª¡ )³
³
γ PP
q
AA B
~ A
h ~B
ω
ê3
~ 12
h
~ 12
ω
ê1 ê2
Figure 2.2: Torque–free spinning of an axi–symmetric satellite as seen in FB .
in the inertial frame, and Π rotates around h, ~ if we look at the motion from FI .
It is possible to define two angles, θ and γ, that remains constant in time,

h12 Jt ω12 
tan θ = = 

h3 Js ω 3 Jt
=⇒ tan θ = tan γ
ω12 
 Js
tan γ = 
ω3
~ θ is the (constant) nutation angle and

If we choose an inertial frame such that Ê 3 k h,
it defines the orientation of the symmetry axis ê3 in the inertial space. The angle γ is the
semi–aperture of the body–cone. This means that the angle between h ~ and ω ~ , equal to
|θ − γ|, is also constant, and ω~ describes a cone around h, ~ fixed in the inertial frame, the
space cone.
The body cone is attached to the body axes, but it is not fixed in space. On the
converse, the space cone, attached to the vector h, ~ that is constant in the inertial frame,
is fixed in FI . The two cones remains tangent along ω ~ , that is the axis of instantaneous
rotation of the body, and the motion of the satellite can be represented by the body cone
rolling along the surface of the space cone.
As a final observation, it should be noted that, when Js > Jt (oblate body, that is
a disk–shaped body), it is γ > θ, and the space cone is inside the body cone. On the
converse, when Js < Jt (prolate body, that is a rod–shaped body), it is γ < θ, and the
space cone is outside the body cone.
For what concern the attitude resulting from such a motion, substituting the expres-
2.1. TORQUE–FREE MOTION OF AXI–SYMMETRIC SATELLITES 27
~
Ê 3 k h
~
ω
ê3
Ê 1 Ê 2
Figure 2.3: Torque–free spinning of an axi–symmetric satellite as seen in FI .
sions for the angular velocity components determined above in the following equation,
ω1 = Ψ̇ sin Φ sin Θ + Θ̇ cos Φ

ω2 = Ψ̇ cos Φ sin Θ − Θ̇ sin Φ
ω3 = Ψ̇ cos Θ + Φ̇
and remembering that Θ ≡ θ = constant ⇒ Θ̇ = 0, one gets
Ψ̇ sin Φ sin Θ = ω12 sin[λ(t − t0 )]

Ψ̇ cos Φ sin Θ = −ω12 cos[λ(t − t0 )]
Ψ̇ cos Θ + Φ̇ = Ω
By squaring and summing the first two equations, it is evident that
Ψ̇2 sin2 Θ = (ω12 )2
so that Ψ̇ is constant; Ψ̇ is called precession rate or coning speed, and it is the angular
velocity of the line of the nodes on the horizontal plane.
Dividing the first equation by the second, the spin angle Φ is determined,
tan Φ = − tan[λ(t − t0 )] ⇒ Φ = −λ(t − t0 )
Deriving w.r.t. time, the inertial spin rate is obtained
Φ̇ = −λ
that is also constant. At this point, it is possible to evaluate the precession rate from the
third equation,
Ω − Φ̇
Ψ̇ cos Θ + Φ̇ = Ω ⇒ Ψ̇ =
cos Θ
But from the definition of λ,

Js − Jt Jt Jt
λ= Ω ⇒ Ω= λ= Φ̇
Jt Js − Jt Jt − Js
one gets
Js Φ̇
Ψ̇ =
Jt − Js cos Θ
If Jt > Js , that is we have a prolate body, Ψ̇ and Φ̇ have the same sign and we have
direct precession, that is the precession rate is in the same direction of the spin rate. On
the converse, if Js > Jt and an oblate body is dealt with, Ψ̇ and Φ̇ have different signs
and we have retrograde precession, the precession rate being in the opposite direction with
respect to the spin rate.
An observation
It is important to note that the derivation presented in this paragraph are valid for any
rigid body which has two equal principal moment of inertia. This is a category much wider
than that of axi–symmetric bodies, including any prism with a basis made of a regular
polygon, but also any other body of irregular shape such that there exists a set of principal
axes of inertia such that Jx = Jy 6= Jz .
2.2 Torque–free motion of tri–inertial satellites

An analytical solution of the motion of tri–inertial rigid body can be derived in terms of
Jacobi elliptic functions. Luckily there is also a geometric description of the same motion,
due to Poinsot, which is much simpler, nonetheless extremely useful for the description
of the dynamics of an arbitrary rigid body.
Remembering that the kinetic energy
1 ~ =d
T = ω~ ·h
2
and the angular momentum vector h are constant, it is possible to consider the (constant)
dot product
~
h 2T
~ · =
ω (2.1)
h h
as the length d of a (constant) segment ON along the direction of h. ~ The invariable plane
σ, which is the plane perpendicular to the direction of h, ~ placed at a distance d from the
body center of mass O, is thus fixed in FI , and it represents the locus of all the possible
~ that satisfy Eq. (2.1). Remembering that the Poinsot ellipsoid is the locus of all the
ω
possible ω ~ that satisfy the kinetic energy equation, the intersection between the ellipsoid
and the invariable plane must contain the angular velocity vector.
It is easy to demonstrate that such an intersection is a single point P , i.e. the Poinsot
ellipsoid and the invariable plane are tangent. Since the time derivative of the kinetic
energy is zero,
dT 1 d~
ω ~
= ·h
dt 2 dt
the increment d~ ω and h~ are perpendicular, thus d~ ~ + d~
ω lie on σ. But since the vector ω ω
must be also on the Poinsot ellipsoid, d~ ω must be tangent to its surface. These two
2.2. TORQUE–FREE MOTION OF TRI–INERTIAL SATELLITES 29
~
h
σ N
~
ω
Figure 2.4: The invariable plane.
~
h
σ N
ω +́´
d~
~
ω
Figure 2.5: The Poinsot ellipsoid rolls on the invariable plane.

conditions can be satisfied only if the Poinsot ellipsoid is always tangent to σ. Moreover,
the Poinsot ellipsoid is fixed in the body frame FB , so that d~ ω is the same in both FI
and FB . This means that the Poinsot ellipsoid rolls without slipping on σ.
2.2.1 Drawing the polhode curves

The locus of all the possible ω~ s on the Poinsot ellipsoid is given by the polhode curve,
which is the intersection between the Poinsot ellipsoid and the angular momentum ellip-
soid. Thus, during the rolling motion, the tangent point moves on the Poinsot ellipsoid
along the polhode. The corresponding curve on σ is the herpolhode. When the body is
axisymmetric, both the polhode and the herpolhode are circles and the situation can be
depicted in terms of space and body cones. In general the herpolhode is not a closed
curve, but the polhode must be a closed curve on the Poinsot ellipsoid, inasmuch as after
a revolution around the spin axis ω ~ must attain again the same value, in order to satisfy
conservation of both kinetic energy and angular momentum.
In order to draw the shape of the polhodes on the Poinsot ellipsoid, it is sufficient to
recall the equations of the Poinsot ellipsoid and the angular momentum ellipsoid, that are
ω12 ω22 ω32
+ + = 1
(2T /Jx ) (2T /Jy ) (2T /Jz )
ω12 ω22 ω32
+ + = 1
(h/Jx )2 (h/Jy )2 (h/Jz )2
Subtracting the first equation from the second and multiplying the result by h2 , one gets
µ ¶ µ ¶ µ ¶
h2 2 h2 2 h2
Jx Jx − ω1 + Jy Jy − ω2 + Jz Jz − ω32 = 0
2T 2T 2T
the polhode equation. In order to have real solutions for the above equation, the three
coefficients cannot have the same sign. For this reason the parameter J ∗ = h2 /(2T ) must
lie between the maximum and the minimum moment of inertia. Assuming, without loss
of generality, that Jx > Jy > Jz , it is
Jz ≤ J ∗ ≤ Jx
The easiest way to determine the shape of the polhodes is to consider their projection
onto the planes of the three–dimensional space ω1 –ω2 –ω3 . Eliminating ω3 between the
equations of the two ellipsoids brings the equation
Jx (Jx − Jz )ω12 + Jy (Jy − Jz )ω22 = 2T (J ∗ − Jz )
which represents an ellipse, since all the coefficients are positive. Similarly, eliminating
ω1 one gets
Jy (Jy − Jx )ω22 + Jz (Jz − Jx )ω32 = 2T (J ∗ − Jx )
which is again the equation of an ellipse, inasmuch as all the coefficients are negative. On
the converse, eliminating ω2 , which is the angular velocity component with respect to the
intermediate axis, brings
Jx (Jx − Jy )ω12 + Jz (Jz − Jy )ω32 = 2T (J ∗ − Jy )
which represent a hyperbola, the coefficients of the left–hand side being of different sign.
It should be noted that, depending on the sign of J ∗ − Jy , which can be either positive
2.2. TORQUE–FREE MOTION OF TRI–INERTIAL SATELLITES 31
ω3
−1
−2
2
1 2
0 1
ω2 0
−1
−1 ω1
−2 −2
Figure 2.6: Polhode curves on the Poinsot Ellipsoid.
or negative, the axis of the hyperbola will be vertical or horizontal, in the ω1 –ω3 plane.
When J ∗ = Jy , the polhode equation degenerates into the form
Jx (Jx − J ∗ ) ω12 + Jz (Jz − J ∗ ) ω32 = 0
which represents the separatrices, the boundaries of motion about the axis of maximum
and minimum inertia.
2.2.2 Stability of torque–free motion about principal axes

Spinning about any of the principal axis is an equilibrium condition for a rigid body of
arbitrary inertias. The shape of the polhodes already provide an information about the
stability of these equilibria, the axes of maximum and minimum inertia being centers
surrounded by finite size orbits, while the intermediate axis is a saddle point, that is a
small perturbation will take the angular velocity vector “far” from the initial point in the
neighborhood of the saddle.
These facts can be demonstrated analytically. Let us consider a spinning motion about
the z axis of the principal frame, such that ω ~ 0 = Ωê3 . Assuming that ω~ =ω ~ 0 + ∆~ ω,
T
where the body frame components of ∆~ ω , given by ∆ω B = {ω1 , ω2 , ω3 } , are small
perturbations with respect to Ω, Euler’s equations can be rewritten as follows
Jx ω̇1 + (Jz − Jy )Ωω2 = 0
Jy ω̇2 + (Jx − Jz )Ωω1 = 0
Jz ω̇3 = 0
where second and higher order terms were neglected. The third equation is decoupled,
thus stating that a perturbation on the spinning axis does not affect the other two. The
first two equations can be rewritten in matrix form
½ ¾ · ¸½ ¾
ω̇1 0 (Jy − Jz )Ω/Jx ω1
=
ω̇2 (Jz − Jx )Ω/Jy 0 ω2
The characteristic equation is thus given by

(Jy − Jz )(Jx − Jz )
λ2 + Ω2 =0
Jx Jy
The roots are pure imaginary,
s
λ = ±iΩ
Jx Jy
if either Jy < Jz and Jx < Jz , or Jy > Jz and Jx > Jz , that is if the spin axis is either
the axis of maximum or minimum moment of inertia. In such a case, the solution of
the equation is confined in the neighborhood of the spin condition which is Ljapunov
(although not asymptotically) stable. On the converse, if the spin axis is the intermediate
one, the product (Jy − Jz )(Jx − Jz ) is negative, inasmuch as one of the two factor is
positive and the other is negative. In this case, the roots are both real,
s
λ = ±Ω −
Jx Jy
and one of the two eigenvalues is positive. This means that spinning around the interme-
diate axis is an unstable equilibrium for the spinning rigid body.
2.3 Nutation Damping

Nutation damping is a simple yet effective way to restore a state of pure spin, if a nu-
tation angle different from zero should be induced by some external cause. We will now
investigate how it is possible to exploit the effects of energy dissipation in order to make
a pure spin condition asymptotically stable.
2.3.1 Effects of energy dissipation

The rigid body is an abstraction. Usually flexible appendages and/or fuel sloshing induce
some energy dissipation the effects of which can be easily determined. Remembering that
the kinetic energy of a rotating rigid body is given by
1 1
T = Jt (ω12 + ω22 ) + Js Ω2
|2 {z } |2 {z }
transverse spin
and its angular momentum is
~ = Jt ω1 ê1 + Jt ω2 ê2 + Js ω3 ê3 ⇒ h2 = J 2 (ω 2 + ω 2 ) + J 2 Ω2
h t 1 2 s
the quantity (2Js T − h2 ) is

(2Js T − h2 ) = Js Jt (ω12 + ω22 ) + Js2 Ω2 − Jt2 (ω12 + ω22 ) − Js2 Ω2
= Jt (ω12 + ω22 )(Js − Jt )
But assuming that the (possible) variation of the nutation angle is slow enough to make
its rate negligible with respect to Ψ̇, it is also
ω12 + ω22 ≈ Ψ̇2 sin2 Θ
2.3. NUTATION DAMPING 33
0.1
ω1 [rad s−1]
0.05
0
−0.05
−0.1
0.1
ω2 [rad s−1]
0.05
0
−0.05
−0.1
0.15
ω3 [rad s−1]
0.1
0.05
0
0 0.5 1
t [103 s]
Figure 2.7: Time–history of angular velocity components with nutation damping.
As a consequence
2Js T = h2 + Jt Ψ̇2 sin2 Θ(Js − Jt )
If we differentiate with respect to time the term 2Js T , taking into account that h and
Ψ̇ are constant, we get
2Js Ṫ = 2Jt Ψ̇2 (Js − Jt ) sin Θ cos ΘΘ̇

Jt 2
Ṫ = Ψ̇ (Js − Jt ) sin Θ cos ΘΘ̇
Js
In a real system there are several causes that induce energy dissipation, such as fuel
sloshing. In such a case, it is Ṫ < 0. So, when energy is dissipated, the quantity
(Js − Jt ) sin(2Θ)Θ̇ is also negative, and a nutation rate is induced by the dissipation.
For a prolate body, when Js < Jt , the resulting nutation rate, for small nutation angles
Θ, is positive. This means that the nutation angle increases with time and the satellite
tends towards a flat spin condition, with Θ = π/2. On the converse, if the satellite is an
oblate body and Js > Jt , Θ̇ is negative, for small nutation angles, and tends to decrease
with time, because of dissipation. This means that the spin around the symmetry axis is
stable, as Θ → 0.
Most satellites use thruster to remove nutation and maintain the desired pure spin
condition. But for oblate bodies a nutation damper provides a very simple and useful
means for passively stabilizing the pure spin motion. An example of damper is the ball–
in–tube: a plastic tube is filled with viscous fluid and a ball bearing is free to move in the
fluid, thus dissipating energy. The only drawback of this technique is that it may take
several minutes (if not hours) to remove nutation.
As ω1 and ω2 are damped out, ω3 increases because of conservation of angular mo-
mentum. Letting Ωf be the final spin rate, it is
h2 = Jt2 (ω12 + ω22 ) + Js2 Ω2 = Js2 Ω2f

ê3
ê2
A
(md /m)ξ
´
u+́
O
~b
ê1
md
ξ
~
n
Figure 2.8: A satellite equipped with a “ball–in–tube” damper.
Therefore,
Jt2 2
Ω2f = Ω2 + 2
(ω1 + ω22 )
Js
and the final spin rate results higher than the initial one, as expected.
2.3.2 Equations of motion of a rigid satellite with damper

In order to derive the equations of motion of a satellite equipped with a spring–mass–
dashpot damper, it is necessary to consider the generalized form of the Euler equation,
~
dh
+S~ × ~aA = M ~
dt
written in a reference frame FA attached to the satellite platform P, with origin in the
point A, which is the center of mass of the satellite when the damper mass md is in the
undeformed position P (ξ = 0) for the spring. All the “moments” are taken with respect
to A.
Assuming that I 0 is the inertia tensor when ξ = 0, the inertia tensor in the general
case is given by h i
I = I 0 + md (ρ2 − b2 )1 − (~ρ~ρ − ~b~b)
where the position vector ~ρ of the mass is
~ρ = ~b + ξ n̂
n̂ being the orientation of the damper and ~b the position vector of P . The overall angular
momentum is
~ = Iω
h ˙ ρ × n̂)
~ + md ξ(~
˙ ~b × n̂)
~ + md ξ(
= Iω
~
~ is the angular velocity of P. It is then possible to write the time derivative of h
where ω
in FA as
Ã !
~
dh
= ḣA + ω A × hA
dt
A
h i
¨ A × nA ) + ω A × Iω A + md ξ(b
= İω A + I ω̇ A + md ξ(b ˙ A × nA )
The time derivative of the inertia matrix I is

· 2 ¸
dρ T T
İ = md 1 − (ρA ρ̇A + ρ̇A ρA )
dt
˙ T nA + b2 , it is
where, being ρ2 = ξ 2 + 2ξbA
dρ2 ˙ + bT n A )
= 2ξ(ξ A
dt
while
˙ A ⇒ ρA ρ̇T + ρ̇A ρT = ξ(b
ρ̇A = ξn ˙ A nT + nA bT + 2ξnA nT )
A A A A A
Summing up, it is
£ ¤
İ = md ξ˙ 2(ξ + bTA nA )1 − bA nTA − nA bTA − 2ξnA nTA
At this point, all the elements required for writing the time derivative of the angular
momentum in FA are explicitly given. In order to express the term S ~ × ~aA , it is now
necessary to derive the actual position of the center of mass O, given by
md
~r A = ξ n̂
m
where m is the overall satellite mass. The static moment can thus be written as
~ = m~r A = md ξ n̂ ⇒ S A = mr A = md ξnA
S
while assuming that the external force F ~ is zero, so that the absolute acceleration of the
center of mass O is zero, the resulting (absolute) acceleration of A is
d2~r A
~aA = −
dt2
By applying two times the rule for the time derivative of a vector quantity in FA , it is
µ 2 ¶
d ~r A
= r̈ A + 2ω A × ṙ A + ω̇ A × r A + ω A × (ω A × r A )
dt2 A
where
md ˙ md ¨
ṙ A = ξ n̂ ; r̈ A = ξ n̂
m m
As a consequence, it is
³m ´ h i
~
(S × ~aA )A = −md ξnA ×
d ¨ ˙
ξnA + 2ξω A × nA + ξ ω̇ A × nA + ξω A × (ω A × nA )
µ 2¶h m
md i
= − ˙ A × (ω A × nA ) + ξ 2 nA × (ω̇ A × nA ) + ξ 2 (ω T nA )(nA × ω A )
2ξ ξn A
m
In order to complete the model, it is necessary to derive the equation of motion for
the damper mass md . The position of md with respect to the (inertially fixed) center of
gravity O is given by
~r m = −~r A + ~b + ξ n̂
so that the Newton’s second law
d2~r m ~m
md =F
dt2
~ m is the overall force acting on md , can be written in FA as
where F
µ 2 ¶
d r mA ~m
md =F
dt2 A
A
Working out the expression of r mA , it is

µ ¶
dr mA ˙ A + ω A × (−r A + bA + ξnA )
= −ṙ A + ξn
dt A
= (1 − md /m)ξn ˙ A + ω A × [(1 − md /m)ξnA + bA ]
µ 2 ¶
d r mA ¨ A + ω̇ A × [(1 − md /m)ξnA + bA ] +
= (1 − md /m)ξn
dt2 A
+ 2(1 − md /m)ξω ˙ A × nA + ω A × {ω A × [(1 − md /m)ξnA + bA ]}
Since md has only one degree of freedom along the direction n̂, it is sufficient to consider
only the component of motion along n̂ itself, which is obtained by taking the scalar
product of Newton’s law times nA . Since all the terms like nTA (ω A × nA ) are zero,
this simplify the expression of the absolute acceleration component of md along n̂. The
equation of motion of md thus becomes1
³
−cξ − kξ = md (1 − md /m)ξ¨ + nTA (ω̇ A × bA ) +
˙
© ª £ ¤¢
+ (nTA ω A ) ω TA [(1 − md /m)ξnA + bA ] − ω 2 (1 − md /m)ξ + nTA bA
where the only forces acting on md along n̂ are the viscous damping term and the elastic
force of the spring.
Collecting all the terms, it is now possible to write the equation of motion in FA of a
rigid satellite equipped with a damper placed in an arbitrary position ~b, with arbitrary
orientation n̂:
h i
¨ A × nA ) + ω A × Iω A + md ξ(b
İω A + I ω̇ A + md ξ(b ˙ A × nA ) +
µ 2¶h i
md ˙ 2 2 T
− 2ξ ξnA × (ω A × nA ) + ξ nA × (ω̇ A × nA ) + ξ (ω A nA )(nA × ω A ) = M A
m
³ © ª
md (1 − md /m)ξ¨ + nTA (ω̇ A × bA ) + (nTA ω A ) ω TA [(1 − md /m)ξnA + bA ] +
£ ¤¢
−ω 2 (1 − md /m)ξ + nTA bA + cξ˙ + kξ = 0
1
To complete the transformaiton, it is necessary to apply the double vector product rule
~ × (~y × ~z ) = (~
x x · ~z )~y − (~
x · ~y )~z
to the last term.

2.3.3 A practical case: the axial damper

The equation of motion for the torque–free motion (M ~ = 0) of a rigid satellite equipped
with a damper will now be specialized to the case where the damper is mounted parallel
to the z–axis, that is n̂ = k̂. Also, a set of principal axes of inertia for the configuration
with ξ = z = 0 is chosen as the reference frame FA and it is assumed that ~b = bî. This
means that
nA = (0, 0, 1)T ; bA = (b, 0, 0)T ; ρA = (b, 0, z)T
Under these assumptions the inertia matrix is
 
Jx + md z 2 0 −md bz
I= 0 Jy + m d z 2 0 
−md bz 0 Jz
and its time derivative is

 
2md z ż 0 −md bż
İ =  0 2md z ż 0 
−md bż 0 0
Also,
bA × nA = biA × kA = −bj A = (0, −b, 0)T
The additional term of the generalized Euler equation is given by
 
µ 2 ¶  2z żω1 + z 2 ω̇1 − z 2 ω2 ω3 
~ × ~aA )A = − md
(S 2z żω2 + z 2 ω̇2 + z 2 ω1 ω3
m  
0
Collecting all the terms, one gets the following system of ordinary differential equations,
Jx ω̇1 + (Jz − Jy )ω2 ω3 + md (1 − md /m)z 2 ω̇1 − md bz ω̇3 + 2md (1 − md /m)z żω1 +

−md bzω1 ω2 − md (1 − md /m)z 2 ω2 ω3 = 0
Jy ω̇2 + (Jx − Jz )ω1 ω3 + md (1 − md /m)z 2 ω̇2 − md bz̈ + 2md (1 − md /m)z żω2 +

+md bz(ω12 − ω32 ) + md (1 − md /m)z 2 ω1 ω3 = 0
Jz ω̇3 + (Jy − Jx )ω1 ω2 − md bz ω̇1 − 2md bżω1 + md bzω2 ω3 = 0
As for the motion of the damper mass, it is
md (1 − md /m)z̈ − md bω̇2 + md bω1 ω3 − md (1 − md /m)z(ω12 + ω22 ) + cż + kz = 0
Stability of a spinning satellite with damper

It is easy to demonstrate that a condition of spin about the z–axis, such that ω B =
(0, 0, ωS )T and z = ż = 0 is an equilibrium for a rigid satellite equipped with a damper such
that n̂ = ê3 . It is then possible to linearize the equation of motion in the neighborhood of
such a spin condition in order to determine its stability. Letting ω B = (ω1 , ω2 , ωS + ω3 )T
and assuming that ω1 , ω2 , ω3 , z and ż are first order perturbation, a set of linear ordinary
differential equation that describe the evolution of the perturbation in the neighborhood
of the nominal spin condition is obtained when higher order terms are neglected. The
linearized form of the model of a spinning satellite with damper is
Jx ω̇1 + (Jz − Jy )ωS ω2 = 0

Jy ω̇2 + (Jx − Jz )ωS ω1 − md bz̈ − md bωS2 z = 0
Jz ω̇3 = 0
md (1 − md /m)z̈ − md bω̇2 + md bωS ω1 + cż + kz = 0
In order to simplify the above expression the following parameters can be defined:
Jz − Jy Jz − Jx k c z md b2
λ1 = ωS ; λ2 = ωS ; µ = md /m ; p2 = ; β= ; ζ= ; δ=
Jx Jy md md b Jy
so that the linear system can be rewritten in the form
ω̇1 + λ1 ω2 = 0
ω̇2 − λ2 ω1 − δ ζ̈ − δωS2 ζ = 0
(1 − µ)ζ̈ − ω̇2 + ωS ω1 + β ζ̇ + p2 ζ = 0
where the third equation (relative to perturbations of the spin condition) is not included,
as it is decoupled from the others.
Taking the Laplace transform L of the time–varying states, such that L{ωi (t)} = ω̄i (s)
and L{ζ(t)} = ζ̄(s), one gets the following transformed system:
sω̄1 + λ1 ω̄2 = 0
sω̄2 − λ2 ω̄1 − s2 δ ζ̄ − δωS2 ζ̄ = 0
s2 (1 − µ)ζ̄ − sω̄2 + ωS ω̄1 + sβ ζ̄ + p2 ζ̄ = 0
which becomes, in matrix form

    
s λ1 0  ω̄1   0 
 −λ2 s −(s2 + ωS2 )δ  ω̄2 = 0
2 2    
ωS −s s (1 − µ) + sβ + p ζ̄ 0
The characteristic equation

¯ ¯
¯ s λ1 0 ¯
¯ ¯
¯ −λ2 s −(s2
+ ωS2 )δ ¯=0
¯ ¯
¯ ωS −s s2 (1 − µ) + sβ + p2 ¯
can be thus written in polynomial form as
(1 − µ − δ)s4 + βs3 + [p2 − δλ1 ωS − δωS2 + λ1 λ2 (1 − µ)]s2 +

+ λ1 λ2 βs + (λ1 λ2 p2 − λ1 δωS3 ) = 0
Application of the Routh’s criteria

Once the coefficients of the polynomial characteristic equation are known, nowadays classic
numerical root finding techniques can be easily applied, using the particular set of satellite
parameters of interest. The drawback of this approach is that it does not offer any physical
insight into the influence of each parameter on the global stability properties of the system.
Routh’s criteria for asymptotic stability can be easily applied, and when the system
order is sufficiently low, it is possible to determine analytically the stability condition
and, as a consequence, stability boundaries in the parameter space. Given the polynomial
characteristic equation
an sn + an−1 sn−1 + . . . + a2 s2 + a1 s + a0 = 0
Routh’s criteria states that a necessary condition for stability is that all the coefficients
have the same sign. Assuming that the coefficient are normalized with respect to an , so
that an = 1, the condition a0 > 0 is the generalized static stability requirement.
In order to evaluate necessary and sufficient condition for asymptotic stability, it is
necessary to build the Routhian matrix, which has the form
an an−2 an−4 ...
an−1 an−3 an−5 ...
b3,1 b3,2 b3,3 ...
b4,1 b4,2 b4,3 ...
where the first two rows are build using the coefficients of the polynomial, while those of
the following ones are evaluated as
bi−1,j bi−2,j+1 − bi−2,j bi−1,j+1
bi,j =
bi−1,j
The number of non zero elements of each row decreases of one unit as the row index i is
increased by two, until the last two rows are reached such that only one element on each
row is non-zero. The following row will be mare entirely of zeroes.
The Routh’s condition for asymptotic stability guarantees that the real part of all the
roots of the characteristic equation are negative if and only if all the elements in the first
column of the Routhian matrix are not null and have the same sign. If the polynomial
equation is normalized with respect to an this condition requires that all the terms bi,1
are positive, with b1,1 = 1 and b2,1 = an−1 .
Application of Routh’s criteria to the present case, i.e. the spinning satellite with
axial damper, brings to the following Routhian matrix:
(1 − µ − δ) [p2 − δλ1 ωS − δωS2 + λ1 λ2 (1 − µ)] (λ1 λ2 p2 − λ1 δωS3 ) 0
β λ1 λ2 β 0 0
(p2 − δωS2 − λ1 δωS + δλ1 λ2 ) (λ1 λ2 p2 − λ1 δωS3 ) 0 0

δβλ1 (λ2 − ωS )(λ1 λ2 − ωS2 )
0 0 0
p2 − δωS2 − λ1 δωS + δλ1 λ2
λ1 (λ2 p2 − δωS3 ) 0 0 0
Since β > 0 for all dissipative dampers, the conditions for asymptotic stability of the
spin condition about the z-axis are
1−µ−δ > 0
p2 − δ(ωS2 + λ1 ωS − λ1 λ2 ) > 0
δβλ1 (λ2 − ωS )(λ1 λ2 − ωS2 )
> 0
p2 − δωS2 − λ1 δωS + δλ1 λ2
λ1 (λ2 p2 − δωS3 ) > 0
Since usually md ¿ m, so that µ ¿ 1 and δ ¿ 1, we can consider the limiting case µ → 0

and δ → 0, so that most of the inequalities become trivial. Taking into account that p2
and β are strictly positive for a real damper, the only conditions left are
λ1 (λ2 − ωS )(λ1 λ2 − ωS2 ) > 0

λ1 λ2 > 0
Taking into account the definitions of λ1 and λ2 , the first inequality can be rewritten
in the form µ ¶2
4 Jz − Jx − Jy
ωS (Jz − Jx ) > 0
Jx Jy
which is equivalent to the condition
Jz > Jx
As a consequence, because of the second condition, it must also be
Jz > J y
so that only spinning about the axis of maximum inertia is a stable condition for a satellite
equipped with a damper. It should be noted that the presence of the damper makes such
a spin condition asymptotically stable.
2.4 Attitude Maneuvers of a spinning satellite

It is possible to use a nutation angle and the precession rate to rotate the spin axis of
a spinning spacecraft. If we consider a thruster generating a force F at position r in
the body frame centered in the center of mass CM, the resulting moment acting on the
satellite is
M =r×F
If we have a firing time od ∆t, starting at t0 , we have an overall change of angular
momentum equal to Z t0 +∆t
∆h = (r × F )dt
t0
but if the firing time is sufficiently small with respect to the period of rotation of the
spinning spacecraft, we can approximate ∆h as
∆h = (r × F )∆t
This thrust pulse changes the direction of the spacecraft’s angular momentum vector,
thus creating a nutation angle θ. At this point the spacecraft is no longer spinning around
its symmetry axis, and a precession rate builds up, equal to
Js Φ̇
Ψ̇ =
Jt − Js cos Θ
It must be remembered the the direction of the precession motion depends on whether
an oblate or a prolate body is considered, but in the latter case the instability of the
spin condition around the axis of minimum inertia prevents the use of this technique for
changing the spin axis direction.
2.4. ATTITUDE MANEUVERS OF A SPINNING SATELLITE 41
h1 h1 ∆h
∆h1 2
∆h1 h2
h0 h0 θ»» θ
θ»» : ¾X z
X
: ´
+́
´
+́
F F
CM r CM r
Figure 2.9: Changing the spin axis direction with nutation angle and precession rate.
Since the precession rate is constant, the precession angle will vary by
Js Φ̇
∆Ψ = ∆t
Jt − Js cos Θ
After a time interval equal to
(Jt − Js ) cos Θ
T =π
Js Φ̇
the precession angle is varied by π rad, and the direction of the symmetry axis of the
spinning satellite has changed of 2θ, so that T is the duration of the manoeuvre. At this
point a second thruster firing is required to stop the precession and achieving again a pure
spin condition about the new re-oriented axis.
From simple trigonometric consideration, it is
||∆h1 || = ||∆h2 || = h0 tan θ = Js Ω tan θ
where h0 is the magnitude of the initial angular momentum.
But is is also
||∆h|| = F `∆t
where ` is the distance between the thruster axis and the symmetry axis of the satellite
(in the figure above this distance is equal to the radius of the spinning body). The total
required firing time for the manoeuver is thus equal to
Js Ω tan θ
τ = 2∆t = 2
F`
In order to estimate the amount of propellant necessary to perform the reorientation, it
must be remembered that the specific impulse of an engine is given by
F
Isp =
∆mg/τ
that is the thrust delivered divided by the weight of propellant burned in a unit time.
The specific impulse is one of the most important characteristics of rocket engines of any
size, and from its value it is possible to determine
Fτ
∆m =
gIsp
The total propellant mass required to rotate the spin axis of an angle equal to 2θ is
Js Ω tan θ
∆m = 2
gIsp `
Propellant mass ∆ m [kg]
0.05 5
Manoeuvre time T [s]

∆ m ∝ tanθ
Nman = 4
0.04 4
0.03 ∆m∝θ 3
0.02 2 Nman = 2
0.01 1 Nman = 1
0 0
0 20 40 60 0 20 40 60
2 θ [deg] 2 θ [deg]
3 400
single pulse ∆ t [ms]
Total firing time [s]
2.5
300
2
1.5 200
1
100
0.5
0 0
0 20 40 60 0 20 40 60
2 θ [deg] 2 θ [deg]
Figure 2.10: Spinning satellite manoeuvre: Fuel consumption (a), Total manoeuvre time
(b), Total firing time (c), single pulse duration (d).
For small reorientation angles it is tan θ ≈ θ, and the propellant consumption is

roughly proportional to the reorientation angle itself. But as θ increases, the fuel con-
sumption gets larger, so that it is more convenient to perform the overall reorientation
in several smaller steps. Moreover, for a given Isp , the pulse duration may become a not
negligible fraction of the spin period, thus reducing the change in angular momentum.
The drawback of a large reorientation performed in N steps is that the manoeuvre time
is increased and the overall manoeuvre is usually less accurate.
Figure 2.10 shows the results for reorientation 2θ angles between 0 and 70 deg, for a
satellite with Js = 20 kg m2 , Jt = 15 kg m2 , spinning at Ω = 2 rad s−1 , that is controlled
by thrusters with a 1 m moment arm, which produce 20 N of thrust with a specific impulse
Isp = 120 s.
The minimum pulse duration considered is 10 ms, that puts a limit on the resolution
of the reorientation, the minimum manoeuvre angle being 0.573 deg. On the other side,
a thrust firing longer than one tenth of the spin period was too long for being considered
a pulse, and this limits the maximum reorientation angle that can be achieved by a single
nutation manoeuvre, that is 17.9 deg. Using 2 manoeuvres, the propellant savings at low
2.5. GRAVITY–GRADIENT STABILIZATION 43
angles are negligible (less than 1%), but it is possible to achieve spin axis reorientation
up to 35 deg. With 4 manoeuvres the maximum reorientation angle becomes as high as
71 deg.
2.5 Gravity–gradient stabilization

Spin stabilization allows one to provide the satellite with gyroscopic stability, which can
be made asymptotic, by a proper use of dissipative devices. At the same time, some
significant disadvantages affects the operational use of spin stabilized satellites, that are
• the presence of only one inertially fixed direction along which to point the payload;
• the inefficient production of electric power, inasmuch as the entire surface of the
satellite must be covered with photovoltaic cells, only half of which faces the sun
during the spinning motion (with less than one fourth at an angle greater than 45
deg with respect to the direction of the sun rays);
• the complexity, limited accuracy and cost in terms of fuel of attitude maneuvers for
reorientation of the spin axis.
A very interesting possibility for stabilizing an Earth facing satellite in low Earth
orbit is to exploit the gravitational torque that acts on any rigid body of finite size.
The differential gravitational acceleration across the satellite provides a small, but non
negligible torque. Simply put, the side of the satellite closest to the Earth experiences a
slightly greater gravitational acceleration that the opposite side. The resulting torque will
tend to align the satellite with the local vertical direction. For some mission applications
the gravity gradient torque provides a disturbance which must be countered. However, for
an Earth facing satellite, the gravity gradient provides a simple, low cost means of attitude
stabilization, albeit rather inaccurate. The great advantage is that no fuel nor energy
dissipation is required to maintain the desired attitude, if a proper stability condition is
met.
2.5.1 Origin of the gravity–gradient torque

The gravity acceleration acting on a mass element dm is given by
~ + ~r
R
~g = −GM
~ + ~r ||3
||R
where G is the universal gravity constant, M is the mass of the primary body (the Earth,
for an Earth orbiting satellite), while R~ and ~r are the position vectors of the center of
mass of the satellite, CM, w.r.t. the center O of the primary and the position vector of
the mass element dm w.r.t. CM, respectively. Letting µ = GM , the total moment w.r.t.
CM is given by
Z
~
G = (~r × ~g ) dm
B
 ³ ´
Z ~r × R ~ + ~r
= −µ   dm
B ~ + ~r ||3
||R
ô2
ê3
ô1
dm
Ê 3 ê2
~g ~r
ô3
CM
~
R
Ê 1 Ê 2
ê1
Figure 2.11: Gravity on a finite–size rigid body.
Since ~r × ~r is zero, it is
Z Ã ~
!
~ = −µ ~r × R
G dm
B ~ + ~r ||3
||R
~ ·R
As a first observation, it is possible to note that G ~ = 0, i.e. the gravity torque around
the local vertical is always zero. In the above expression, it is possible to write
h³ ´ ³ í−3/2 h i−3/2
~ −3 ~
||R + ~r || = R + ~r · R + ~r ~ 2 ~
= R + 2~r · R + r 2
where r and R are the magnitude of ~r and R, ~ respectively. The latter expression can be
rewritten as µ ¶−3/2
~ −3 −3 2 ~ r2
||R + ~r || = R 1 + 2 ~r · R + 2
R R
Being r ¿ R, the third term inside the parentheses can be neglected, and the second one
is a first order perturbation compared to 1. Moreover, we can expand a term like (1 + ε)n
in the neighbourhood of ε = 0 as
(1 + ε)n = 1 + nε + O(ε2 )
that in the present case gives
µ ¶−3/2 Ã !
2 r 2
~
r · ~
R
~ +
R−3 1 + 2 ~r · R ≈ R−3 1 − 3 2
R R2 R
Using this latter expression in the definition of the gravity torque, one gets
Z "Ã !
~ ³ ´
#
~ = − µ ~
r · R ~
G 1−3 2 ~r × R dm
R3 B R
µZ ¶ ·Z ³ ´ ¸
µ ~ 3µ ~ ~
= − 3 ~r dm × R + 5 ~r · R ~r dm × R
R B R B
2.5. GRAVITY–GRADIENT STABILIZATION 45
The first term is zero from the definition of center of mass. Remembering the definition
of dyadic,2 the second term can be rewritten as
·µZ ¶ ¸
~ = 3µ ~ ×R ~
G ~r~r dm R
R5 B
Remembering that the inertia tensor can be written as

Z
¡ 2 ¢
I= r 1 − ~r~r dm
B
and noting that ·µZ ¶ ¸

2 ~ ×R
r 1dm R ~ =0
B
it is possible to write the gravity gradient torque as
~ = 3µ ô3 × (I ô3 )
G
R3
~
where ô3 = −R/R is the unit vector along theplocal vertical passing through the satellite
center of mass. Since the orbit rate is ω0 = µ/R3 , it is possible to rewrite the above
expression as
G~ = 3ω 2 ô3 × (I ô3 )
0
2.5.2 Equilibria of a rigid body under GG torque

The equation of motion of a satellite under the action of the gravity–gradient torque can
thus be written as
dh~
− 3ω02 ô3 × (I ô3 ) = 0
dt
It can be observed that the gravity gradient torque is zero whenever
ô3 × (I ô3 ) = 0
which requires that one of the principal inertia axes is oriented along the local vertical ô3 .
Such a condition can be maintained on a circular orbit if the satellite is purely spinning
about an axis perpendicular to the orbit plane at an angular velocity which is equal to
~ = −ω0 ô3 . In such a case a relative equilibrium with respect to
the orbit rate, so that ω
the orbit frame is obtained, with the satellite principal axes aligned with the orbit frame.
Expressing the equation in terms of body frame components gives the following equation
of motion
I ω̇ B + ω B × (Iω B ) − 3ω02 o3B × (Io3B ) = 0
2 ~ ~y is a tensor such that
The dyadic x
x~y )~z = (~y · ~z )~

(~ x
Expressing the vector components in a frame FF , it is
x~y )F = xF y TF
(~
If the attitude of the satellite is expressed in terms of Bryan angles (roll, pitch and
yaw) with respect to the orbit frame FO = {ô1 ; ô2 ; ô3 }, the body frame components of
ô3 are given by the third column of the matrix LBO , i.e.
o3B = (− sin θ, cos θ sin φ, cos θ cos φ)T
which, assuming that all the angles remain “small”, becomes
o3B = (−θ, φ, 1)T
The (absolute) angular velocity components are
ω B = ω rB − ω0 o2B
where ω rB = (ω1 , ω2 , ω3 )T is the angular velocity of the satellite with respect to FO and
o2B is given by
o2B = (ψ, 1, −φ)T
for small perturbation. As a result, angular velocity components in body frame can be
expressed as
ω B = (φ̇ − ω0 ψ, ϑ̇ − ω0 , ψ̇ − ω0 φ)T
Substituting this perturbation expressions for ω B and o3B in the equation of motion and
neglecting second order terms, one gets
Jx φ̈ − (Jz + Jx − Jy )ω0 ψ̇ + 4(Jy − Jz )ω02 φ = 0

Jy ϑ̈ + 3(Jx − Jz )ω02 ϑ = 0
Jz ψ̈ + (Jz + Jx − Jy )ω0 φ̇ + (Jy − Jx )ω02 ψ = 0
First of all, it can be noted how the pitch motion is decoupled from the roll–yaw
motion. The equation of motion of pitch (harmonic) oscillation can be rewritten in the
form
ϑ̈ − aϑ = 0
where the quantity
Jz − Jx 2
a=3 ω0
Jy
must be negative, in order to have finite size oscillations, the amplitude of which depends
on the initial conditions ϑ(t = 0) and ϑ̇(t = 0). The stability condition for pitch motion
simply requires that
Jz < Jx
If this condition is not met, the coeffient
√ a becomes positive and the characteristic equation
2
λ = a has two real solutions, λ = ± a, so that the pitch angle will exponentially grow
after a perturbation from the equilibrium condition ϑ = ϑ̇ = 0.
As for the roll–yaw equation, taking the Laplace transform, we have, in matrix form,
the following equation
· ¸½ ¾ ½ ¾
Jx s2 + 4(Jy − Jz )ω02 −(Jz + Jx − Jy )ω0 s φ̄(s) 0
=
(Jz + Jx − Jy )ω0 s Jz s2 + (Jy − Jx )ω02 ψ̄(s) 0
so that the characteristic polynomial p(s) = det(A) can be writtein in the form
p(s) = b0 s4 + b1 s2 + b2
2.6. DUAL–SPIN SATELLITES 47
where
b 0 = Jx Jz
¡ ¢
b1 = Jy2 − 3Jz2 − Jx Jy + 2Jy Jz + 2Jx Jz ω02
¡ ¢
b2 = 4 Jy2 − Jx Jy − Jy Jz + Jx Jz ω04
Let p
−b1 ± b21 − 4b0 b2
x1,2 =
2b0
denote the roots of the equation b0 x2 + b1 x + b2 = 0. If x1,2 is a pair of
√complex conjugate
roots
√ x = α ± iβ, the poles of the system can be written as s1,2 = α − iβ, and s3,4 =
α + iβ. In both cases, two of the resulting roots of the characteristic equation have a
positive real part. So, the discriminant of the characteristic equation
∆ = b21 − 4b0 b2
must be positive for stability. Moreover, if one (or both) of the solutions xi is (are) real
and positive, this means that one (or two) pole(s) of the system is (are) real and positive
and the equilibrium is unstable.
p As b0 is the product of two positive quantities, x1,2
2
are negative only if −b1 ± b1 − 4b0 b2 < 0. If b1 is negative, this condition cannot be
satisfied, and at least one of the solutions xi will be positive. If b1 > 0, both solutions
will be negative if and only if b2 is also positive, in which case ∆ < b21 . Summarizing, the
set of stability conditions can be written in the following form:
b1 > 0
b2 > 0
∆ = b21 − 4b0 b2 < 0
Letting
Jy − Jz Jy − Jx
k1 = ; k2 =
Jx Jz
it is possible to write the stability condition as a function of the two inertia parameters
k1 and k2 only, as
1 + 3k1 + k1 k3 > 0
k1 k3 > 0
2
(1 + 3k1 + k1 k3 ) − 16k1 k3 > 0
The requirement for pitch stability, Jx > Jz can be rewritten as k1 > k2 , and the resulting
stability region is plotted in Fig. 2.12
2.6 Dual–spin satellites

The limit of the use of gyroscopically stabilized satellites is that there is only one direction
which is fixed in space and along which it is possible to point a sensor or a communi-
cation antenna. Moreover, maneuvering the satellite in order to reorient the spin axis is
both difficult and expensive in terms of fuel consumption, the required propellant being
proportional to the amount of angular momentum stored in the spinning satellite.
1
0.5
k2
0
−0.5
−1
−1 −0.5 0 0.5 1
k1
Figure 2.12: Stability region of the Eart–pointing attitude for a rigid satellite under
gravity gradient torque.
In order to overcome this negative features, an alternative satellite configuration was

proposed and used in the past, where the payload and the communication hardware is
mounted on a platform P which has a relative rotational degree of freedom with respect
to the spinning part of the satellite, the rotor R. During orbit injection the two elements
are locked and the satellite behaves like a standard rigid body. After the operational
orbit is reached, the spin–up motor accelerates the rotor, with respect to the platform,
in order to transfer on the rotor the entire vehicle angular momentum, thus despinning
the platform. The de–spin maneuver is ended when the platform is fixed in space, while
the rotor provides gyroscopic stability to the entire vehicle. The main advantage of such
a solution is that, pivoting the payload mounted on the platform allows aiming it in any
direction in space, with great flexibility, without changing the direction of the rotor spin
axis.
From the configuration point of view, the first dual–spin satellites where made of a
large rotor, containing most of the satellite equipment, with solar cells mounted on the
surface of the rotor, like in any other spinning satellite. A small platform was attached
to it, as represented in Fig. 2.13. Later, the platform became larger (about one half of
the satellite). The Intelsat IV was the first prolate dual spin satellite (2.14). Finally, the
rotor was substituted by a momentum wheel, which is a small (and light) wheel spinning
at a very high speed, placed inside the platform. In this latter case, the satellite is called
a momentum–bias satellite. Both configurations (the dual–spin and the momentum bias
satellite) can be referred to as gyrostat.
Figure 2.13: Dual spin satellite: configuration with large rotor (e.g. GOES-7 satellite).
Figure 2.14: Dual spin satellite: configuration with equivalent rotor and platform (e.g.
Intelsat IV satellite).
ê3
ê2
Ω ê1
W
â P
Figure 2.15: Sketch of a satellite equipped with a momentum wheel.
2.6.1 Mathematical model of a gyrostat

The presence of an inertially fixed component and a spinning part has an effect on the
mathematical model of the satellite dynamics. First of all, it is necessary to derive from
the usual angular momentum balance the equation of motion of the satellite. Figure 2.15
shows the generic arrangement of a momentum wheel W spinning about the spin axis â
at an angular rate Ω with respect to the satellite platform P. It should be noted that,
although there is a rotating mass, the mass distribution does not change because of the
relative degree of freedom between W and P. For this reason, it is possible to define in
this case a body frame FB , fixed with respect to the platform, the inertia tensor I of the
gyrostat P + W being constant in FB .
The overall angular momentum is given by
~ = Iω
h ~ + hr â
where, letting Is be the moment of inertia of W about â, the quantity hr = Is Ω is the
relative angular momentum. The absolute angular momentum stored in the spinning
rotor
~
ha = hr + Is â · ω
is made up of two contributions, the relative angular momentum due to the spin motion
with respect to P, and the angular momentum related to the motion of the platform.
Writing the components in body frame, it is
hB = Iω B + Is ΩaB
¡ ¢
ha = Is Ω + aTB ω B
The equation of motion of the gyrostat can be written as
ḣB = −ω B × hB + M B
ḣa = g
where MB is the external torque acting on the satellite, while ga is the torque applied by
the de–spin engine to the rotor. It should be noted how ga , which is an internal torque,
does not appear in the equation of the time evolution of the overall angular momentum.
The kinematic problem for the attitude of the satellite will be solved expressing ω B
as a function of the total and wheel angular momentum, that is
ω B = I −1 (hB − Is ΩaB )
Integrating the Euler angle rate equation (or the quaternion evolution equation), as usual,
one gets the time evolution of satellite attitude. At the same time, the rotor angular speed
can be expressed as
Ω = ha /Is − aTB ω B
It is also possible to express the gyrostat dynamics directly in terms of the angu-
lar velocity vector and rotor spin rate. In this case, it is possible to write the angular
momentum vector as
¡ ¢
hB = Iω B + ha − Is aTB ω B aB
¡ ¢
= I − Is aB aTB ω B + ha aB
The time derivative of the angular momentum vector components are thus expressed as
¡ ¢
h˙B = I − Is aB aTB ω̇ B + ga aB
so that the equations of motion of the gyrostat become

−1
h ³ ´ i
ω̇ B = Ĩ M B − ω B × Ĩω B + ha aB − ga
ω̇ = ga /Is − aTB ω̇ B
¡ ¢
where Ĩ = I − Is aB aTB is the pseudo–inertia tensor.
These equations of motion can be used for describing the gyrostat spin–up dynamics,
during platform de–spin, that is, a spin axis torque ga is applied until the platform angular
velocity drops to zero.
Assuming that FB is a set of principal axes of inertia and that the spin axis is parallel
to ê3 (axial rotor), the satellite and rotor absolute angular momentum can be written,
respectively, as  
 Jx ω 1 
hB = Jy ω2 ; ha = Is (Ω + ω3 )
 
Jz ω3 + Is Ω
The equation of motion of the axial gyrostat can be written out in full as follows:
Jx ω̇1 + (Jz − Jy )ω2 ω3 + Is Ωω2 = M1

Jy ω̇2 + (Jx − Jz )ω1 ω3 − Is Ωω1 = M2
Jz ω̇1 + (Jy − Jz )ω1 ω2 + Is Ω̇ = M3
Is (Ω̇ + ω̇3 ) = ga
2.6.2 The apparent gyrostat and the Kelvin gyrostat

For the torque–free case (M B = 0) it is possible to define two cases which can be treated
using the same set of equations, provided that the symbols are defined in a proper way.
If the rotor relative spin rate is assumed to be constant, Ω(t) = const, the Kelvin
gyrostat is obtained. The variable hr becomes a system parameter and the platform
equation of motion decouples from rotor spin dynamics. In such a case, it is
I ω̇ B = −ω B × (Iω B + hr aB ) + M B
In order to keep Ω constant, it is necesary to apply a certain spin torque to the rotor,
which, from the equation of motion of rotor spin dynamics can be evaluated as
ga = Is aTB ω̇ B
In real applications it is not possible to keep the rotor speed exactly constant and a
simple feed–back loop is employed for determining the required spin torque. A tachometer
is mounted on the wheel shaft, which measure the actual rotor spin rate Ω and the error
with respect to the nominal (desired) rotor speed Ωdes is the input for the rpm regulator.
As an example, the spin torque can be set equal to
ga = K(Ωdes − Ω)
where K is the gain of the rotor control system. In this way the rotor is accelerated
(ga > 0) whenever the angular speed is lower than the prescribed one, and decelerated
when Ω > Ωdes . Of course, if the spin wheel dynamics is modeled including this simple
control scheme, it is no longer possible to treat Ω as a constant in the platform angular
velocity equations, and the full system equations must be accounted for.
If the rotor spin–up torque is zero (ga = ḣa = 0), the rotor absolute angular momentum
becomes constant and the apparent gyrostat model is dealt with. In this case ha is a system
parameter and the equation of motion can be rewritte as
³ ´
Ĩ ω̇ B = −ω B × Ĩω B + ha aB + M B
with the usual meaning for the symbols.
2.6.3 Stability of the de–spun condition

When Ω = 0, the dual–spin satellite behaves like a standard rigid body. This situation is
called all-spun condition. After platform de–spin, the satellite angular velocity must be
zero and all the satellite angular momentum must be stored in the rotor. The stability
of the equilibrium for the de–spun condition of a Kelvin gyrostat with axial rotor can be
analyzed as usual assuming the components ω B = (ω1 , ω2 , ω3 ) as small and neglecting
second order terms in the equation of motion, which can be written for the torque–free
case as
Jx ω̇1 + Is Ωω2 = 0
Jy ω̇2 − Is Ωω1 = 0
Jz ω̇3 = 0
The third equation is decoupled and the other two can be written in matrix form as
½ ¾ · ¸½ ¾
ω̇1 0 −(Is /Jx )Ω ω1
=
ω̇2 (Is /Jy )Ω 0 ω2
It is easy to show that the eigenvalues of the state matrix are a pair of conjugate imag-
inary numbers, so that the variation of ω1 and ω2 remains confined in the neighborhood
of the origin, provided that Ω 6= 0. The rotor provides the required gyroscopic stability
to the entire satellite.
Chapter 3
Active Stabilisation and Control of

Spacecraft
For many applications where accurate payload pointing is required full three-axis control
is used. The satellite does not spin (except for slew manoeuvres) and is stabilised by
reaction wheels (gyros) and/or thruster relative to an Earth (or other direction) facing
attitude. Since the satellite is not spinning, large Sun facing solar panels can be used for
power generation. The satellite is also capable of fast reorientation manoeuvres (slews) in
any arbitrary direction. However, due to the added complexity, 3-axis stabilised satellites
are far more expensive than simple spin–stabilised ones. The application and its accuracy
requirements must justify the increase in costs and the reduced life expectation of the
vehicle.
In order to actively stabilize the satellite and allow for autonomous maneuvering, the
satellite must be equipped with a suitable set of sensors, which provide the on–board
computer with the necessary information on the current attitude. The precision of the
control action depends on the accuracy of the attitude measurement. Moreover, the
control law itself depends on the type of actuators available. As a consequence, a brief
outline on the more common attitude sensors and control actuator is given, prior to the
discussion of attitude active control.
3.1 Attitude sensors

In order to implement some feedback control law for actively stabilizing a spacecraft, it
is necessary to have an information on the current satellite attitude. This information
is provided by sensors that can be based on significantly different principles, depending
on the satellite application and the related pointing accuracy requirements. Roughly
speaking, it is possible to group the attitude determination hardware into 5 categories:
• Earth sensors;
• Sun sensors;
• star sensors (and star trackers);
• gyroscopic sensors;
• magnetometers.
55
56 CHAPTER 3. ACTIVE STABILISATION AND CONTROL OF SPACECRAFT
3.1.1 Infrared Earth sensors (IRES)
Earth sensors are capable of determining the attitude of the satellite with respect to
the Earth, by sensing the Earth’s horizon. For this reason they are also referred to as
“horizon sensors”. Two techniques can be used: (i) the dynamic crossing of the horizon,
with precise determination of the crossing points, and (ii) the static determination of the
horizon inside the filed of view of the instrument.
In principle the sensor could be based on any radiation source coming from the Earth,
but the radiation in the visual wavelength region depends too strongly on the character-
istics of the reflecting surface. Infrared (IR) radiation provides a range of wavelengths
in which the energy is more uniform. Still some problems remains on how to model the
Earth’s surface and how to take into account the thickness of the atmosphere, which also
gives a sizable contribution to the IR radiation of the Earth.
Horizon–crossing sensor (IRHCES)
Static Horizon sensor (IRSES)
3.1.2 Sun sensors
Single–axis analogic sun sensor
Two–axis analogic sun sensor
Digital sun sensor
3.1.3 Star trackers
3.1.4 Rate and rate–integrating gyroscopes
3.1.5 Magnetometers
3.2 Actuators
3.2.1 Reaction control system
Cold gas jets
Chemical propulsion
Mono–propellant propulsion system

Bi–propellant propulsion system
3.3. LINEAR MODEL OF SATELLITE ATTITUDE MOTION 57
Electric propulsion
3.2.2 Momentum exchange devices

Reaction wheels
Bias–momentum wheel
Double–gimbal bias–momentum wheel
Control Moment Gyroscopes
3.2.3 Other attitude control devices
3.3 Linear model of satellite attitude motion

For a general three axis motion, Euler’s equation states that
ḣB + ω B × hB = M B (3.1)
If there are no internal moving parts (e.g. no spinning rotors) this equation can be written
in the form
I ω̇ B + ω B × (Iω B ) = M B
Assuming that the rotation is slow, the angular velocity can be considered as a first order
perturbation of the steady state ω B = 0, so that ||ω B || ¿ 1 and ω B × (Iω B ) is a second
order (negligible) term in the equations of motion. Euler’s equation thus becomes
I ω̇ B = M B
When dealing with a stabilization problem, that is small perturbations with respect to
a prescribed nominal attitude are considered, also the angular displacement is “small”. In
such a case, it is possible to demonstrate that, letting the variation of the Euler’s angles
be equal respectively to
∆Ψ = θ3 ; ∆Θ = θ2 ; ∆Φ = θ1
and defining θ = {θ1 , θ2 , θ3 }T , it is

θ̇ = ω B
Euler’s equation in linear form become simply
I θ̈ = M B (3.2)
Finally, assuming that a set of principal axes of inertia is chosen as the body frame, the
inertia matrix is in diagonal form, I = diag(J1 , J2 , J3 ), and three independent second
order linear ordinary differential equations are obtained, one for each axis,
Ji θ̈i = Mi , i = 1, 2, 3
It should be noted that for single axis rotations, a scalar equation in the same form as
Eq. (3.2) is obtained, regardless of the angular speed and angular displacement, that is
J θ̈ = M
This equation will be the starting point for the analysis of attitude slew maneuvers.
If the satellite is equipped with a momentum wheel with a moment of inertia Is
spinning about the axis â with a relative spin rate Ω, the total satellite angular momentum
is
hB = Iω B + hr aB
where hr = Is Ω is the relative angular momentum; Ω can be written as
Ω = Ω0 + ∆Ω
being ∆Ω ¿ Ω0 , i.e. ∆Ω is a “small perturbation” of the reference wheel spin rate Ω0 .
Equation (3.1) thus becomes
I ω̇ B + Is ∆Ω̇aB + ω B × [Iω B + Is (Ω0 + ∆Ω) aB ] = M B
Neglecting higher order term one gets
I ω̇ B + Is ∆Ω̇aB + Is Ω0 ω B × aB = M B
The evolution of the perturnations of Euler’s angles is again
θ̇ = ω B
so that the equations of motion becomes can be expressed in scalar form as follows:
J1 θ̈1 + Is ∆Ω̇a1 + Is Ω0 (θ̇2 a3 − θ̇3 a2 ) = M1
J2 θ̈2 + Is ∆Ω̇a2 + Is Ω0 (θ̇3 a1 − θ̇1 a3 ) = M2
J3 θ̈3 + Is ∆Ω̇a3 + Is Ω0 (θ̇1 a2 − θ̇2 a1 ) = M3
This set of equations is coupled with the dynamics of the spin–wheel, described by the
(linear) equation
ga
Ω̇ ≡ ∆Ω̇ = − ω̇ TB aB
Is
When the rotor spin axis is parallel to the i–th principal axis of inertia (â ≡ êi ), the
i–th equation of motion becomes decoupled from the others, because aj = ak = 0, for
j, k 6= i. As an example, let us assume that â ≡ ê2 , so that a1 = a3 = 0 and a2 = 1. The
equations of motion become
J1 θ̈1 − Is Ω0 θ̇3 = M1
J2 θ̈2 + Is ∆Ω̇ = M2
J3 θ̈3 + Is Ω0 θ̇1 = M3
The pitch dynamics is decoupled from the roll and yaw motions, but is coupled with the
spin–wheel dynamics
ga
∆Ω̇ = − ω̇2
Is
Summing up, the torque–free pitch dynamics (M2 = 0) is described by the equations
(J2 − Is ) θ̈ = −ga
J2
Ω = Ω0 − θ̇
Is
On the converse, roll and yaw motions are coupled by the gyroscopic term. Moreover,
if the rotor spin speed is time–varying and ∆Ω is not negligible, the roll and yaw coupling
terms depend on the pitch dynamics through Ω.
3.4. USE OF THRUSTERS FOR ATTITUDE CONTROL 59
3.4 Use of thrusters for attitude control

3.4.1 Single axis slews (open loop)
If we consider a constant torque acting around the generis i–th axis, the angular acceler-
ation is
αi = θï = Mi /Ji
Integrating from the initial condition θ(t0 ) = θ0 , θ̇(t0 ) = θ̇0 (and dropping the subscript
i), one has
θ̇(t) = θ̇0 + α(t − t0 )

1
θ(t) = θ0 + θ̇0 (t − t0 ) + α(t − t0 )2
2
Writing
dθ̇ dθ
θ̈ =
dθ dt
it is
dθ̇
α = θ̇
dθ
Integrating the previous equation in θ, one gets
Z θ̇ Z θ
ϑ̇dϑ = αdϑ
θ̇0 θ0
that, for constant applied torque, becomes

1 2
(θ̇ − θ̇02 = α(θ − θ0 )
2
or,
θ̇2 = θ̇02 + 2α(θ − θ0 )
Starting from a rest position in the nominal attitude (θ̇ = 0 and θ = 0), one gets
√
θ̇2 = 2αθ =⇒ θ̇(t) = 2αθ
Representing the manoeuvre in the phase plane θ − θ̇, the state variables evolve as
reported in Fig. 3.1: There is an initial acceleration phase and the spacecraft achieves
the desired slew rate θ̇d for the reorientation, until, at a proper time instant, it is stopped
by a second impulse from the thrusters.
For an approximate evaluation of the total time of the manoeuvre and the propellant
necessary for performing it, it is sufficient to neglect the duration of the thruster firings
for accelerating and decelerating the angular motion of the vehicle, taking into account
only the drift phase. Assuming that the drift angular velocity is θ̇d , the time required for
a manoeuvre of amplitude equal to θf is simply given by T = θf /θ̇d .
The configuration of the thrusters used for spinning the satellite up to the drift angular
velocity can be different. As reported in Fig. 3.2, a single, unbalanced thruster can be
used, firing which a force F ~ and a torque of magnitude F ` will act on the satellite. If
a balanced configuration is employed, the simultaneous firing of the two thrusters will
produce only a torqe F `, with a null net force, so that the orbit motion of the spacecraft
is not perturbed.
.
θ
acceleration drift deceleration
.
θd
θ0 θ θf
Figure 3.1: Evolution of θ and θ̇ for a rest–to–rest manoeuvre.
Assuming a sharp thrust profile, it is

F`
θ̇d = ∆t
I
and the pulse width (thruster on–time) for the spin–up is given by
I θ̇d
∆t1 =
F`
Again, from the definition of the thruster specific impulse
F ∆t
Isp =
∆mg
(that is the ratio between the momentum gained and the weight of the propellant used)
one gets
F ∆t1 I θ̇d
∆m1 = =
gIsp gÌsp
But, from the evaluation of the manoeuvre time, it is also θ̇d = θf /T , so that the
propellant mass necessary for the spin–up becomes equal to
Iθf
∆m1 =
gÌsp T
At the end of the drift interval, a second pulse of equal duration in the opposite direction
will stop the motion of the satellite and the new attitude will be acquired. The total
propellant mass is thus equal to
Iθf
∆mtot = ∆m1 + ∆m2 = 2
gÌsp T
²¯ ²¯
X
X
»
~
-F X
X
» - F~/2
» ±° » ±°
` `
j j
`
²¯
»
»
F~/2 ¾ X
X
±°
Figure 3.2: Unbalanced and balanced thruster configurations.
It can be observed that allowing larger manoeuvre times, the same angle variation
can be achieved with a significantly smaller amount of fuel. Moreover, low energy thrust
pulses are less likely to excite vibrations in the structure and/or in flexible appendages
attached to the spacecraft (such as antennas or solar panels). For large values of T , smaller
on–time for the thruster are required, and the initial assumption that the acceleration and
deceleration phases have a negligible duration with respect to the overall manoeuvre time
is confirmed.
On the other side, for agile spacecraft, that achieve high angular velocities, it is pos-
sible to improve the accuracy of the estimate of manoeuvre time taking into account the
duration of the initial and terminal phases. The time for accelerating the satellite up to
θ̇d is, again
I θ̇d
∆t1 =
F`
but this time ∆t1 is not negligible with respect to the total manoeuvre time T .
and it is equal to the time ∆tdec necessary for bringing it back to rest at the end of the
maneuver. During the acceleration, the variation of θ is
1 1 θ̇d2
θ1 = α∆t21 =
2 2α
and the total variation during the manoeuvre is given by
θf = θ̇d td + θ1 + θ2
where the angle variation during the acceleration and deceleration phases, θ1 and θ2 , are
equal. The drift time is td = T − 2∆t1 , where T is the total manoeuvre time, so that
θf = θ̇d (T − 2∆t1 ) + α∆t21
Ã !
θ̇d θ̇2
= θ̇d T − 2 + d
α α
The total manoevure time becomes thus
θf θ̇d
T = +
θ̇d α
In the limit, the minimum–time bang–bang control law brings to zero the drift time
td . The reader can demonstrate that in this limit case it is
q
T = 2 θf /α
In both cases, the faster manoeuvres are performed at the expenses of a significant
increase in propellant consumption.
3.4.2 Closed–loop control

The ideal case: thrust modulation
If the satellite is equipped with a set of sensors that provides the necessary information
on the attitude motion, it is possible to implement a control law which provides a torque
command that drives the satellite towards the prescribed attitude.
If the torque command is simply proportional to the attitude error, e = θdes − θ, the
closed–loop equation of motion becomes
¾
θ̈ = M/J K
=⇒ θ̈ = (θdes − θ) =⇒ θ̈ + p2 θ = p2 θdes
M = K(θdes − θ) = Ke J
where p2 = K/J. The solution of this second order equation for an initial condition of
rest for θ = 0 is
θ(t) = θdes [1 − cos (pt)]
i.e. an unacceptable undamped oscillation of amplitude θdes about the desired condition.
In order to damp the oscillation and asymptotically reach the desired condition, it is
necessary to add a damping term in the control law, proportional to the angular rate θ̇.
The command torque is thus
M = K(θdes − θ) + Kd θ̇ (3.3)
and the closed–loop equation of motion becomes
θ̈ + cθ̇ + p2 θ = p2 θdes
where the coefficient of the damping term c = 2ζp = −Kd /J is positive (damped oscilla-
tions) if the gain associated to the angular rate is negative. In this case the time–history
of the angular motion is
θ(t) = θdes [1 − exp(−ζpt) cos (pd t)]

p
where pd = p 1 − ζ 2 , for
p 0 < ζ < 1. If ζ > 1, the characteristic equation has two real
solutions λ1,2 = p(−ζ ± ζ 2 − 1), which are both negative. In this case the rotation angle
evolves as follows:
· ¸
λ2 λ1
θ(t) = θdes 1 − exp(λ1 t) + exp(λ2 t)
λ2 − λ1 λ2 − λ1
Two important issues need to be considered, at this point. The first one is the choice
of the gains. The case with ζ < 0 is not considered at all, inasmuch as one of the poles of
the closed loop system would be a positive real number, that is, the closed loop system
would be unstable. When ζ < 1 (sub–critical damping), the time τ to damp out 99% of
the initial oscillation amplitude increases as ζ gets smaller according to the equation
exp(−ζpτ ) = 0.01 =⇒ τ = − log(0.01)/(ζp)
At the same time, if ζ > 1 (super-critical damping) the time constant of the slowest
mode becomes larger, being approximately equal to ζ/p, for large values of ζ. The best
performance in terms of maneuver agility are obtained for the critical damping ζ = 1,
which guarantees the fastest convergence time to the desired position. The over–damped
2
1 1
.5s+1
Input Output
Modulator
filter
Schmidt trigger
Figure 3.3: Pulse–Width/Pulse–Frequency (PWPF) modulator.
case may be of interest in those cases when fuel consumption is more a stringent concern
than maneuver time. On the contrary, the under–damped case, which exhibits some
overshoot, is of no practical interest, since it is characterized by a longer settling time
and an increased fuel consumption, because of the thrusting system firing alternatively in
both directions.
The second issue concerns the actual implementation of a control law in the form (3.3).
Such a control law would be realistic only if the control torque provided by the thruster
could be modulated. Unfortunately direct thrust modulation would require a complex
hardware and, at the same time, would be extremely inefficient in large portions of the
required thrust range. For these reason, thruster are always on–off devices. As a conse-
quence, some form of pulse modulation is required to provide a feasible implementation of
the control torque demand. Two techniques will be presented, the Pulse–Width/Pulse–
Frequency (PWPF) modulator, which is an analogic device, and the Pulse–Width Mod-
ulator (PWM), which lends itself to a discrete time implementation.
Pulse–Width/Pulse–Frequency modulation
Figure 3.4 shows the structure of a Pulse–Width/Pulse–Frequency (PWPF) modulator.
The main element of the modulator is the Schmidt trigger, which consist of a double relay
with hysteresis, separated by a dead band. In order to provide a quasi–linear steady–state
response, a modulation filter is added, the input of which is the signal em = ẽ − ym , i.e.
the difference between the feed–back signal error and the modulator output. The output
of the filter v is the activation signal for the Schmidt trigger.
The output y of the modulator remains zero until the activation signal v remains below
the threshold Uon , that is until ẽ = em < Uon /K. The dynamics of the modulator filter is
1
v̇ = (Kem − v)
τ
The response of the modulator for a steady input ẽ is
v(t) = v0 exp (−t/τ ) + Kem [1 − exp (−t/τ )]
until y = 0. When v reaches Uon the filter output switches from 0 to Um and the thruster
is activated by the trigger, while the modulator input changes to ẽ − Um .
Calling ton the time in which the thruster is switched on and assuming that ẽ is
unaffected by y, integration of the filter linear dynamics starting from the initial condition
v(ton ) = Uon brings to the following expression for the activation function
µ ¶ · µ ¶¸
t − ton t − ton
v(t) = Uon exp − + K(em − Um ) 1 − exp −
τ τ
µ ¶
t − ton
= K(em − Um ) + [Uon − K(em − Um )] exp −
τ
The thruster is switched off when the activation function becomes equal to Uoff , that is
when
µ ¶
toff − ton Uoff − K(em − Um )
exp − =
τ Uon − K(em − Um )
Uon − Uoff
= 1− (3.4)
Uon − K(em − Um )
It should be noted that, if K(em − Um ) > Uoff , i.e. em > Uoff /K + Um , the above equation
has no real solution, that is the input is so great that the modulator calls for a continuous
thruster engagement. This level constitutes the saturation level of the modulator. On
the converse, if em > Uoff /K + Um , the thruster are switched off when t = toff . Calling
∆ton = toff − ton the on–time of the thruster, and assuming that it is sufficiently small, it
is possible to approximate exp (−∆ton /τ ) = 1 − ∆ton /τ , so that
Uon − Uoff
∆ton ≈ τ
Uon − K(em − Um )
After the thruster pulse, the modulator input becomes again em = ẽ and the response
of the filter for an initial condition given by v(to f f ) = Uoff is
µ ¶ · µ ¶¸
t − toff t − toff
v(t) = Uoff exp − + Kem 1 − exp −
τ τ
The thruster will be activated again when v = Uon , that is when

µ ¶
t − toff Uon − Kem
exp − =
τ Uoff − Kem
Assuming again that the off time is also sufficiently small, one gets
Uon − Uoff
∆toff = τ
Kem − Uoff
The average output of the thruster is thus given by
∆ton
ȳ = Um
∆toff + ∆ton
for Uon /K < em < Uoff /K + Um . It is zero below the activation threshold, equal to the
thruster force beyond the saturation value. The proper sizing of the trigger and filter
characteristics allows for realizing a modulator which takes into account the physical
characteristics of the thrusters, such as the minimum pulse time and the minimum time
between pulses. Figure ?? shows the response of a satellite controlled by a set of thruster
driven by a PWPF modulator (continuous line). The response is compared with that
obtained by thrust modulation (dashed line).
DC
1.0
0.5
0.0
0 2 4 6 8 10
ẽ
Figure 3.4: Duty cycle of the PWPF modulator.

y 6
Um
0 ¾
- ¾ -¾ -
-
ts ∆ton ∆toff t
¾ ∆tdc -
Figure 3.5: Pulse–Width Modulation.
Pulse–Width Modulation
The PWPF modulator is a heritage of the analogic era, although it is still widely used. The
main advantage of PWPF modulation remains the possibility of commanding a sequence
of pulses with a quasi–linear response, in terms of averaged output. Nonetheless, many
actuators are now commanded in the framework of a digital implementation of the control
logic. It is well known that a digital controller obtained from the discretization of an
analogic one can attain at most the performance of the original analogic counterpart. If a
digital controller is designed directly in the discrete time domain, taking into account the
sampling frequency and A/D and D/A conversions, better performance can be obtained.
Pulse–Width modulation provide a means for designing a thruster control law where
the implementation of the switching logic is inherently digital. As an advantage with
respect to PWPF (analogic) modulation, it is much easier to satisfy restrictions on the
on–time ∆ton and the time–between–pulses ∆toff .
Figure 3.5 shows a sequence of pulses. If ts is the sampling period of the digital
controller, a pulse can be commanded for each duty cycle, the duration of which is ∆tdc =
N ts . The minimum impulse bit (MIB) Non,min is chosen so that Non,min ts > ∆ton,min .
At saturation, the on–time is extended to the entire duty cycle. Below saturation, the
ȳ
Um
1
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1 1.2 e
esat
Modulator data: ts = 10 ms; Ndc = 20; ∆ton,min = 12 ms ⇒ Non,min = 2;
MTBP = 17 ms ⇒ Non,max = 18.
Figure 3.6: Response of a PW Modulator.
maximum pulse must satisfy the constraint on minimum–time–between–pulses (MTBP),

that is ∆toff = (N − Non )ts >MTBP.
The average output of a thruster controlled in PWM is
∆ton Non
ȳ = Um = Um
∆toff + ∆ton N
The response of the PWM to a constant input error signal ẽ is depicted in Fig. 3.6,
assuming that a quasi–linear average output ȳ = kẽ is desired between a deadband edb
and the saturation value esat , when the on–time equals the duty cycle and ȳ = Um . It
should be noted that, by a proper software implementation of the activation signal, it is
possible to obtain arbitrary pulse–width modulation curves.
3.4.3 Fine pointing control

After slew to target, the control system switches to fine pointing control, in order to
maintain the payload aimed at the desired target with a prescribed tolerance, in presence
of external disturbances. This can be done using momentum exchange devices and/or gas
jets. A simple way to achieve such a fine pointing is to use cold gas jets, that produce
thrust in the range between 0.05 and 20 N, with short pulse times (of the order of 10−2 s)
for fine control. As far as the satellite will be bouncing back and forth between the edges
of a limit cycle, it is important to estimate the amount of fuel necessary as a function of
the prescribed pointing accuracy and thruster characteristics.
Torque–free control
Let us consider the toruqe–free case, that is no external disturbances act on the spacecraft,
where the payload must be pointed in a given direction, within some deadband, which in
0.01
drift
θ [deg s ]
−1
.
thruster firing
thruster firing
0
drift
deadband
−0.01
−3 0 θ [10−4 deg] 3
Figure 3.7: Limit cycle for gas–jet fine pointing control.
the case of astronomy payloads can become as small as 1/10 arcsec = 0.000278 deg!
The satellite is allowed to drift slowly inside the deadband, with drift velocity θ̇0 .
When it reaches the limit of the deadband, θdb , the thruster is fired so as to revert the
drift rate in the opposite direction. The satellite will drift towards the other side of the
deadband, where the thruster are fired again, thus capturing the satellite into a limit cycle
of width approximatley equal to 2θdb .
The thruster generates a moment
M = F ` = J θ̈
For a thruster pulse width ∆t, the change of the attitude rate is
F`
∆θ̇ = ∆t
J
and, in order to revert the satellite motion after the firing, it must be ∆θ̇ = 2θ̇0 . Remem-
bering that
F ∆t
Isp =
∆mg
one gets, for each thruster firing, a fuel consumption of
F ∆t1 J∆θ̇
∆m1 = =
gIsp gÌsp
Considering both sides of the deadband, the total fuel consumption for each cycle is
J∆θ̇
∆m = 2
gÌsp
Ignoring the (usually negligible) thruster pulse time, the duration of the limit cycle is
θdb θdb
τ =4 =8
θ̇0 ∆θ̇
and the mean fuel consumption (that is the fuel burned per unit time) is
∆m J∆θ̇2
ṁ = =
τ 4`gIsp θdb
which can be rewritten as

(F `∆t)2
ṁ =
4`gIsp θdb
This latter equation shows clearly that in order to achieve a great accuracy (small value
of θdb ) without paying too a great penalty in terms of fuel consumption, the characteristics
of the thrusters are of paramount importance, inasmuch as ṁ remains small onfly if high
specific impulse thruster are employed, capable of delivering the force F for a very short
pulse time ∆t.
Bias torque compensation

In most cases, the payload must be aimed at a target (be it on Earth or in the open
space) in presence of external disturbances. If the resulting torque is roughly constant,
such as that due to the gravity gradient, it is possible to use this torque for reversing the
motion on one side of the deadband. In this case, the environmental torque will bring the
spacecraft drifting towards the top of the tolerance interval. At this point the thruster are
fired, to reverse the motion, that will be slowed down by the adverse disturbance torque.
Rather than allowing the spacecraft to reach the other side of the deadband and fire a
second thruster pulse, one allows the bias torque to slow and reverse the drift inside the
deadband interval. The limit cycle will be made of two parabolic arcs, one relative to the
thruster firing phase, similar to that of the previous case, that reverse the satellite motion
on one side of the deadband, and the second arc, relative to the slow deceleration and
motion reversal due to the external torque.
If the disturbance torque generates an angular acceleration
θ̈ = αD = MD /J
it is
θ̇(t) = θ̇0 + αD ∆t
θ(t) = θ0 + θ̇0 ∆t + αD ∆t2 /2
where the subscript 0 indicates the initial condition at time t0 . From the first equation,
it is ∆t = (θ̇(t) − θ̇0 )/αD , that, substituted in the second one, after some simple algebra
gives
θ̇2 = θ̇02 + 2αD (θ(t) − θ0 )
If the initial condition is on the deadband limit, it is θ0 = θdb . In order to keep the limit
cycle inside the dead band, the turning point at which θ̇ drops to zero and reverts its sign
must be crossed on the other side of the deadband, when θ0 = −θdb . In this case, it is
θ̇2 = 0 = θ̇02 + 2αD (−2θdb )
that is
J θ̇0
θdb =
4MD
3.5. MOMENTUM EXCHANGE DEVICES FOR ATTITUDE CONTROL 69
0.01
drift
θ [deg s−1]
.
thruster firing
0
deadband
−0.01
−3 3
0 θ [10−4 deg]
Figure 3.8: Limit cycle for gas–jet fine pointing control with constant disturbance torque.
In order to close the limit cycle, the thruster pulse must be such that the total ∆θ̇ is
twice the drift angular velocity θ̇0 at deadband crossing,
∆θ̇ = 2θ̇0 = F `∆t/J
Substituting the latter expression in the equation of θdb one gets
(F `∆t)2
θdb =
16MD J
and again the requirement on the thruster is to deliver a small duration pulse for keeping
the deadband small enough to satisfy fine pointing requirements.
3.5 Momentum exchange devices for attitude control

Gas jets can deliver very high torque to manoeuvre a spacecraft and change its orientation,
but they have a limited life, due to the finite propellant mass budget available. For
this reason, it is important to rely on some different kind of actuators for stabilising
and controlling satellite attitude during normal operations, without using propellant.
Reaction wheels (RWs) can exchange angular momentum with the spacecraft bus using
only electrical power, that is obtained from the solar panels for the whole lifetime of the
satellite.
Each RW is attached to the satellite structure trough an electric motor, that can
be used to accelerate and decelerate the wheel, relatively to the satellite. If we assume
that initially both the satellite and the wheel are still in the inertial frame (zero angular
momentum initial condition), spinning up the wheel in one direction will cause the bus to
rotate in the opposite direction, because of conservation of overall angular momentum.
Let us consider a single axis rotation about the principal axis of inertia êi . The
spacecraft has a principal moment of inertia Ji about êi , while the wheel spins around
the axis â = êi , having a spin moment of inertia Jsw about â and angular velocity Ω,
relative to the spacecraft bus.
The overall angular momentum about êi is
hi = Ji ωi + Jsw (Ω − ωi )
where the first term is the angular momentum of the satellite bus and the second one is
the angular momentum stored in the wheel spinning at an absolute angular velocity Ω−ωi
with respect to FI . If initially ωi = Ω = 0, so that hi = 0, and taking into account that
the total angular momentum is conserved, inasmuch as the wheel is spun up or slowed
through the exchange of internal torques applied by the electric motor, one gets
Jsw
ωi = − Ω (3.5)
Ji + Jsw
Usually it is Jsw ¿ Ji , so that ωi ¿ Ω, that is the angular velocity achieved by the
spacecraft is only a small fraction of the relative spin angular velocity of the wheel.
In order to manoeuvre the spacecraft about 3 axes, it is necessary to equip the satellite
with at least 3 RWs, one for each body axis. The cluster of RWs usually contains a
fourth spare wheel, the axis of which is oriented in such a way so as to exchange angular
momentum components to or from any of the three axes, if one of the main wheels fails.
One of the wheels can be used to de–spin the satellite after orbit acquisition. In such
a case, this wheel shall be oriented as the spin axis of the satellite and apogee kick rocket
motor stack and it is called a momentum wheel. Its angular momentum is unlikely to
drop back to zero and the control torques about its axis will be obtained accelerating and
decelerating the wheel.
3.5.1 Sizing a reaction wheel for single axis slews

Assuming for simplicity torque–free motion, we will size a reaction wheel with respect to
requirements in single axis slew manoeuvre. Dropping the i subscript, the total angular
momentum about one of the principal axis is
h = hs + hw
where hs = Jω is the spacecraft angular momentum, while hw = Jsw (Ω − ω) is the wheel
angular momentum. Since there are no external toruqes, h =const and
ḣ = ḣs + ḣw =⇒ ḣs = −ḣw
One pulse is used for spinning up the reaction wheel, achieving a certain spacecraft drift
angular velocity, according to Eq. (3.5), while a second pulse will spin down both the
spacecraft and the wheel, achieving at the end of the manoeuvre the desired attitude θf .
After the first pulse, the spacecraft gains an angular momentum equal to
ḣs (t2 − t1 ) = J θ̇ = −ḣw (t2 − t1 )
But it is also
ḣs = J θ̈ = Jα
where α is the angular acceleration of the spacecraft, so that
1
θ(t2 ) = α(t2 − t1 )2
2
1 ḣs
= (t2 − t1 )2
2J
6T
ḣs
θ = θf
t3 t4 -
t1 t2 t
θ=0
−ḣs
Figure 3.9: Torque profile for a single axis slew manoeuvre.
For a minimum–time rotation, a bang–bang control must be used, so that t2 ≡ t3 and

the final rotation angle is
ḣs
θf = 2θ(t2 ) = (t2 − t1 )2
J
But assuming t1 = 0, it is also tf = 2(t2 − t1 ), for bang–bang control, and one gets
1 ḣs 2
θf = t
4J f
Because of conservation of angular momentum, we recall that ḣs = −ḣw , so that the
maximum torque that the electric rotor must produce is
Jθf
ḣw = 4
t2f
But for t = t2 , it is also
hs (t2 ) = −hw (t2 ) = −ḣw (t2 − t1 ) = −ḣw tf /2
when the wheel achieves the maximum angular momentum. This means that the momen-
tum capacity of the wheel must be
θ2
hmax
w ≥ 2J
tf
Example
We require a spacecraft with J = 100 kg m2 to perform a 0.2 rad slew in 10 sec. This means
that the electric motor must apply torques in both direction equal to
Jθf 4 × 100 × 0.2

ḣw = 4 = 0.8 Nm
t2f 102
and the momentum capacity of the RW must be at least
θ2 2 × 100 × 0.2
hmax
w = 2J = = 4 kg m2 s−1
tf 10
Assuming a wheel inertia Jsw = 0.03 kg m2 , the maximum wheel angular rate with respect to
the satellite bus is
J + Jsw J + Jsw hs
Ω = − ω=−
Jsw Jsw J
J + Jsw hw
= hw ≈ = 133 rad s−1 = 21.2rpm
JJsw Jsw
3.5.2 Closed–loop control using a reaction wheel for single axis

slews
With the usual meaning of the symbols, the equation of motion for single–axis rotation is
d ³ ´
J θ̇ + Is Ω = Mext
dt
so that it is
J θ̈ + Is Ω̇ = Mext (3.6)
The reaction wheel dynamics is described by the first order equation
³ ´
Is Ω̇ + θ̈ = ga ⇒ Is Ω̇ = ga − Is θ̈
Substituting the latter expression in Eq. (3.6) and assuming for simplicity that the external
torque and the momentum bias are both zero, it is possible to rearrange the dynamics
equation as
(J − Is )θ̈ = −ga
By choosing
ga = KP (θ − θdes ) + KD θ̇ (3.7)
the closed–loop dynamics is described by the second order differential equation
J ∗ θ̈ + KD θ̇ + KP (θ − θdes )
where J ∗ = J − Is is the moment of inertia of the platform without the spin–wheel. By

defining the closed–loop system bandwidth as p2 = KP /J ∗ and a damping parameter ζ
such that 2ζp = KD /J ∗ , the equation of motion can be rewritten in the form
θ̈ + 2ζpθ̇ + p2 θ = p2 θdes
which can be solved analytically and clearly presents a nice way of choosing the gains KP
and KD .
3.5.3 Bias torque and reaction wheel saturation

Let us assume for simplicity that the satellite is subject to a constant disturbance torque
MD . It is easy to see that the angular displacement θ evolves according to the differential
equation
θ̈ + 2ζpθ̇ + p2 θ = p2 θdes + MD /J ∗
and asymptotically converge to the steady state
θss = θdes + MD /KP 6= θdes

This means that the use of the control law of Eq. (3.7) may not guarantee a sufficient
pointing accuracy. The steady state error can be made smaller by increasing the propor-
tional gain KP , but it will never vanish. Moreover, increasing the gain will likely degrade
the performance of the real system, since the sensor noise will also be multiplied by a high
gain. It is possible to asymptotically track the prescribed position by adding an integral
term to the control law, i.e.
ga = KP (θ − θdes ) + KD θ̇ + KI y (3.8)
where Z
ẏ = θ − θdes ⇒ y = (θ − θdes )dt
It should be noted that a steady–state can be achieved only for ẏ = 0, i.e. θ = θdes ,
with a value of y such that
KI y = MD ⇒ ga,ss = MD
This means that if the spacecraft is subject to a constantly not null external torque, the
average value of which is not zero (aerodynamic torque, gravity gradient), the platform
can be maintained aimed at a fixed position only if an equivalent motor torque is delivered
to the reaction wheel. As a consequence, the wheel will absorb these disturbances by spin–
up and fix the orientation, but its angular velocity will constantly increases. Should the
wheel reach the maximum spinning speed, the attitude control system is saturated.
In order to avoid such a problem, before reaching saturation it is necessary to dump
the angular momentum excess. This can be done simply braking the wheel and firing the
thrusters in the opposite direction. In this way the momentum accumulated in the wheel
is brought back to zero, because of the external torque produced by the thrusters. During
wheel desaturation the pointing accuracy of the control system is usually reduced.
The use of thrusters for despinning the wheels shows that also when using RWs for
attitude control it is necessary to carefully design the propellant mass budget. Thruster
firings will be also in this case the main source of fuel consumption during the station
keeping phase, and will be the main limit to the satellite operative life.
3.5.4 Roll–yaw axes control in presence of momentum bias using

reaction wheels
In order to derive the (linearized) equation of motion of a satellite equipped with a mo-
mentum wheel parallel to the pitch axis and two reaction wheels for roll and yaw control,
the angular momentum can be expressed as
hB = Iω B + H B
where
H B = Ib (Ω0 + ∆Ω)e2B + Ir Ω1 e1B + Ir Ω3 e3B
is the angular momentum stored in the spinning wheel, Ib being the moment of inertia of
the pitch–axis bias–momentum wheel while Ir is the moment of inertia of the two reaction
wheels. Writing H B as
H B = H 0B +∆H B = (0, h0 , 0)T +(∆h1 , ∆h2 , ∆h3 )T = (0, Ib Ω0 , 0)T +(Ir Ω1 , Ib ∆Ω, Ir Ω3 )T
and assuming that ||Iω B ||, ||∆H B || ¿ h0 , the equation of motion can be written as
I ω̇ B + Ib ∆Ω̇e2B + Ir Ω̇1 e1B + Ir Ω̇3 e3B + ω B × Ib Ω0 = M B
where higher order terms have been neglected. Each wheel is characterized by a dynamics
in the form
I Ω̇ = ga − I ω̇
Assuming finally a small angular displacement, so that ω̇ B = θ̈, the pitch equation is
decoupled from the others and takes the form already discussed in the previous paragraph,
while the coupled roll–yaw dynamics is described by the following system of ordinary
differential equation:
J1 θ̈1 + Ir Ω̇1 − Ib Ω0 θ̇3 = M1

J3 θ̈1 + Ir Ω̇3 + Ib Ω0 θ̇1 = M3
Coupling these equation with the dynamics of the reaction–wheels
Ω̇1 = ga1 /Ir − θ̈1

Ω̇3 = ga3 /Ir − θ̈3
one gets the system of (linear) equations
(J1 − Ir )θ̈1 − Ib Ω0 θ̇3 = −ga1 + M1

(J3 − Ir )θ̈3 + Ib Ω0 θ̇1 = −ga3 + M3
Choosing for both reaction wheels spin–torque command laws in the form
gai = KPi θi + KDi θ̇i
and taking the Laplace transform of the system of linear ordinary differential equation,
one gets for the torque–free case (M1 = M3 = 0) the following algebraic equation in the
Laplace variable s:
· 2 ¸µ ¶ µ ¶
s + 2ζ1 p1 s + p21 −γ1 s θ̄1 (s) 0
2 2 =
γ2 s s + 2ζ2 p2 s + p2 θ̄3 (s) 0
where p2i = KPi /(Ji − Ir ) and 2ζi pi = KD /(Ji − Ir ), while the gyroscopic coupling term is
represented by the coefficient γi = Ib Ω0 /(Ji − Ir ). With a proper choice of the gains, it is
possible to have p1 = p2 and ζ1 = ζ2 , which means that perturbations about roll and yaw
axes are damped with similar time behaviour. Moreover, in such a case the number of
design parameter is reduced to two and a simple parametric study is possible. It should be
noted that, independently of the choice of the controller parameters, the proposed control
law requires an explicit estimate of the yaw angle, which may not be a trivial task. A
direct measure of the roll angle is usually available from a Earth horizon sensor, while the
yaw angle can be measured directly only with the aid of an expensive star tracker. As an
alternative, a rate–integrating gyro may be used, where care must be taken on the closed
loop dynamics of the effects of gyro drift and integration of the sensor noise over large
period of time.
Figure 3.10 shows this study for γ = 1. As it can be observed, for ζ = 0, the closed
loop dynamics presents undamped oscillascions for all the values of p, which means that,
as usual, a feedback on angular displacement only is not capable of providing asymptotic
8
Re(λ) 0
6
Im(λ) ζ=0 −5
4
−10
2
−15
0
0 1 2 3 4 5 p
−2
−4 Im(λ) 5
−6 0
−8 −5
−15 −10 −5 0
Re(λ) 0 1 2 3 4 5 p
8
Re(λ) 0
−5
6
Im(λ) ζ = 0.5
4 −10
2 −15
0 1 2 3 4 5 p
0
−2 Im(λ) 5
−4 0
−6
−5
−8
−15 −10 −5 0 0 1 2 3 4 5 p
Re(λ)
8
Re(λ) 0
6
Im(λ) ζ=1 −5
4
−10
2
−15
0
0 1 2 3 4 5 p
−2
−4 Im(λ) 5
−6 0
−8 −5
−15 −10 −5 0
Re(λ) 0 1 2 3 4 5 p
8
Re(λ) 0
−5
6
Im(λ) ζ = 1.5
4 −10
2 −15
0 1 2 3 4 5 p
0
−2 Im(λ) 5
−4 0
−6
−5
−8
−15 −10 −5 0 0 1 2 3 4 5 p
Re(λ)
Figure 3.10: Root locus as a function of p for different values of ζ.
stability to the origin. For subcritical damping (ζ ≤ 1), the eigenstructure is characterized
by two pairs of complex conjugate eigenvalues, for low values of the proportional gain
KP . Beyond a certain critical value (corresponding to pcr = 1.65, for ζ = 0.5), a couple of
complex conjugate roots collapses and two real negative eigenvalues are found for higher
values of p. One of them tends towards the origins, as p is increased, so that extremely
high proportional gains are not convenient. The critical value at which a pair of complex
conjugate roots is replaced by two real eigenvalues decreases as ζ is increased. In the
supercritical range, for ζ > 1, the root locus changes its structure and regions of p can be
found where all the eigenvalues are real.
3.5.5 Roll–yaw axes control using a double–gimbal momentum

wheel
An approach very similar to that used for the case of a pitch axis momentum bias wheel
plus two reaction wheels can be used of deriving the equation of motion of a satellite
equipped with a double–gimbal momentum wheel. The angular momentum can be ex-
pressed again as
hB = Iω B + H B
where
H B = Ib (Ω)aB
is the angular momentum stored in the spinning wheel, Ib being the moment of inertia
of the pitch–axis momentum wheel, while the spin axis can be deflected by the gimbals
from its nominal position parallel to the pitch axis, so that
aB = (− cos δ2 sin δ1 , cos δ1 cos δ2 , sin δ2 )T
Assuming that the gimbal rotation angles δ1 and δ2 are small, it is aB ≈ (−δ1 , 1, δ2 )T . In
this case H B can be written in the form
H B = H 0B + ∆H B = (0, h0 , 0)T + (∆h1 , ∆h2 , ∆h3 )T

= (0, Ib Ω0 , 0)T + (−Ib Ω0 δ1 , Ib ∆Ω, Ib Ω0 δ2 )T
Assuming again that ||Iω B ||, ||∆H B || ¿ h0 , the linearized equation of motion for small
amplitude angular displacement is obtained neglecting higher order terms as
I ω̇ B + ∆ḢB + ω B × Ib Ω0 = M B
which becomes, choosing the set of principal axes of inertia as the body frame,
J1 θ̈1 + ∆Ḣ1 − h0 θ̇3 = M1

J2 θ̈1 + ∆Ḣ2 = M2
J3 θ̈1 + ∆Ḣ3 + h0 θ̇1 = M3
As usual, the pitch dynamics is decoupled and is not repeated here, while roll and yaw
equation of motion are coupled by the gyroscopic term. Letting ur = −∆ḣ1 and uy =
−∆ḣ3 be the control variables, the coupled roll–yaw dynamics is described by the following
set of linear second order equations:
J1 θ̈1 − h0 θ̇3 = ur + M1
J3 θ̈3 + h0 θ̇1 = uy + M3
3.6. QUATERNION FEEDBACK CONTROL 77
where the gimbal rotation command is
δ̇1 = ur /h0
δ̇2 = −uy /h0
It is clear that the same control law derived for the cluster of momentum bias wheel and
reaction wheels can be implemented with a double–gimbal hardware, provided that the
control power made available by the rotation of the gimbals is sufficient for the prescribed
stabilization and control purposes.
3.6 Quaternion feedback control

Many satellites are reoriented performing successive rotation about the control axis, in
order to achieve the desired attitude. Unfortunately, this strategy, although very simple to
be implemented, is not time–optimal nor optimal in terms of fuel or energy consumption.
The overall angular path is far in excess the minimum one described by the Euler rotation
about the Euler’s axis. The quaternion feedback control provides a means for obtaining
a nearly–optimal reorientation, with a control logic only marginally more complex.
The rigid body dynamic equations written in vector form are given by
I ω̇ B + ω B × (Iω B ) = M B
while the equation that describe the evolution of the quaternions is

    

 q̇0 
 0 −ω1 −ω2 −ω3   q0 

  1   
q̇1 ω1 0 ω3 −ω2  q1
=  

 q̇2 
 2 ω2 −ω3 0 ω1   q2 

   
q̇3 ω3 ω2 −ω1 0 q3
Quaternions can be easily computed by modern attitude determination systems and for
this reason their use is nowadays very common. As a consequence, a simple feedback
control law based on the information obtained from the attitude sensors would be very
easily implementable on–board for autonomous manoeuvre management.
It is possible to demonstrate that a feedback law in the form
M B = −Kq e − Cω B
where q e = (q1e , q2e , q3e )T is the attitude error quaternion vector, is globally asymptotically
stabilizing for a wide choice of the gain matrices K and C.
In general, the error quaternion vector is given by
    

 q 0e 
 q 0c q 1c −q2c −q 3c 
 q 0 

   −q1c   q1 
q1e  q 0c −q3c q 2c 
=

 q2e  q2c q3c q0c q1c   q2 
   
q3e q3c −q2c −q1c q0c q3
where (q0c , q1c , q2c , q3c )T represents the commanded (desired) value of the quaternions and
(q0 , q1 , q2 , q3 )T their current value. Without loss of generality, it is possible to assume as
the “absolute” reference frame the desired attitude (and in case switch to another absolute
frame if a new attitude needs to be acquired), so that (q0c , q1c , q2c , q3c )T = (1, 0, 0, 0)T and
q e ≡ q. In this case the control logic can be rewritten as
M B = −Kq − Cω B
It should be noted that, if the desired attitude is assumed to be (q0c , q1c , q2c , q3c )T =
(−1, 0, 0, 0)T (which is the same attitude, but represented by the complementary quater-
nion vector1 ), the quaternion feedback control laws needs to be rewritten as
M B = +Kq − Cω B
Among the simplest possible globally stabilizing controllers, we have
Controller 1: K = k1, C = diag(c1 , c2 , c3 )

k
Controller 2: K= 1, C = diag(c1 , c2 , c3 )
q43
Controller 3: K = k sign(q4 )1, C = diag(c1 , c2 , c3 )
Controller 4: K = [αI + β1]−1 , C = diag(c1 , c2 , c3 )
where k and ci are positive scalar constants, 1 is the 3 × 3 identity matrix, sign(q4 ) is the
sign function,2 and α and β are non–negative scalars.
Controller 1 is a special case of controller 4, when α = 0. It is possible to choose β = 0
only if α 6= 0. Controllers 2 and 3 approach the origin (1, 0, 0, 0)T taking the shorter
angular path, that is, tracking an Euler rotation.
All these control laws can be implemented under the assumption that there is a set
of actuators that can deliver the required amount of external torque M B , that is a set
of thruster. The problem of the jet selection logic and torque saturation are beyond the
scope of the present discussion.
3.7 Control Moment Gyroscopes

When using momentum management devices for attitude stabilisation and control it can
be useful to write the attitude equation of motion in terms of angular momentum, instead
than working with the angular velocity. Although this is true in general, it becomes almost
necessary for modeling a satellite equipped with a cluster of Control Moment Gyroscopes
(CMG).
A CMG is made of a spinning wheel mounted on a pivoting gimbal, the axis of which is
perpendicular to the wheel spin axis. Instead than accelerating and decelerating the wheel
for obtaining the proper reaction from the spacecraft platform, the momentum exchange
between the wheel and the bus is achieved rotating the wheel spin axis about the gimbal.
In this way a gyroscopic torque is obtained in a direction perpendicular to both the spin
and the gimbal axes.
1
Remember that the quaternions (q0 , q1 , q2 , q3 )T and (−q0 , −q1 , −q2 , −q3 )T represent the same atti-
tude.
2
The sign function is 1 when the argument is positive, -1 when it is negative.
3.7. CONTROL MOMENT GYROSCOPES 79
When the dynamic state of the spacecraft is expressed in terms of the (conserved)
angular momentum vector, the equation of motion is rewritten as
ḣB + ω B × hB = M B
where hB is the overall angular momentum (expressed in body–frame components). If

the spacecraft is a rigid body, it is simply
hB = Iω B
and the equation can be trivially rewritten in the previous formulation. If on the converse
the spacecraft has n rotating parts, it is
n
X
hB = Iω B + hiB
i=1
where hiB is the relative angular momentum vector stored inside the i–th spinning wheel.3
Letting
X n
HB = hi
i=1
it is
hB = Iω B + H B
and the equation of motion becomes
I ω̇ B + ω B × (Iω B ) + ḣB + ω B × hB = M B
where M B is the external torque acting on the spacecraft. Assuming for the sake of
simplicity a torque–free environment and a zero initial angular momentum, it is possible
to write the last equation as
I ω̇ B + ω B × (Iω B ) = u
where
u = −Ḣ B − ω B × H B
is the control torque generated by a proper action on the spinning wheels inside the
spacecraft bus. The control torque demand can be determined on the basis of some
feedback control logic, such as the quaternion feebback control presented in the previous
section. Once the control torque demand u is known, the CMG steering logic can be
determined as
Ḣ B = −u − ω B × H B (3.9)
For a cluster of n CMG, the internal angular momentum vector H B is a function of
the rotation angles δ1 , δ2 , . . . , δn of each wheel about its gimbal axis, that is
H B = H(δ) (3.10)
3
The inertia tensor I is assumed to take into consideration the contribution of the spin–wheels to the
mass distribution as if they were still. For this reason only the relative contribution to angular momentum
is added.
Taking the time derivative of this latter equation, one gets

dH B ∂H dδ
=
dt ∂δ dt
that can be rewritten in compact form as
Ḣ B = A(δ)δ̇ (3.11)
where A is the jacobian matrix

∂H
A=
∂δ
The problem is now that of finding a proper inversion of the equation Ḣ B = Aδ̇, that
provides a gimbal rate command δ̇ capable of delivering the required control torque u for
the current gimbal position δ and angular velocity ω B .
The case of a cluster of 4 CMG mounted with the gimbal axis perpendicular to the
faces of a pyramid the sides of which are inclined of a skew angle β is the most popular
in the literature. This CMG configuration, where the angular momentum vectors of each
wheel span one of the faces of the pyramid, is particularly popular as it delivers a nearly
spherical momentum envelope, when β = 54.73 deg and each wheel spins at the same
angular velocity relative to the spacecraft bus. Assuming that the x and y body axes are
perpendicular to the sides of the base of the pyramid and that the z axis is parallel to the
pyramid height, it is
n
X
H B = H(δ) = hiB
i=1
       
−cβ sin δ1 − cos δ2 cβ sin δ3 cos δ4
= h  cos δ1  +  −cβ sin δ2  +  − cos δ3  +  cβ sin δ4 
sβ sin δ1 sβ sin δ2 sβ sin δ3 sβ sin δ4
where h is the angular momentum stored inside each CMG, while sβ = sin β and cβ =
cos β.
In order to perform the reorientation, it is necessary to drive the gimbal according to
Eq. (3.11), where, for the case of the pyramid mounting, it is
 
−cβ cos δ1 sin δ2 cβ cos δ3 − sin δ2
A =  − sin δ1 −cβ cos δ2 sin δ3 cβ cos δ4 
sβ cos δ1 sβ cos δ2 sβ cos δ3 sβ cos δ4
and the CMG angular momentum rate is given by Eq. (3.9). A possible inversion of Eq.
(3.11) is given by the Moore–Penrose pseudo–inverse matrix
δ̇ = A# Ḣ B (3.12)
where
A# = AT (AAT )−1
Although most of the CMG steering logic proposed in the literature are variants of this
form, it should be noted that the pseudo–inverse does not exists when the matrix AAT
becomes singular. The singular points where the determinant det(AAT ) = 0 are charac-
terized by a gimbal configuration such that the CMG cluster cannot deliver torque in one
3.7. CONTROL MOMENT GYROSCOPES 81
direction, called the singular direction. These points are distributed along hypersurfaces
in the δ–space that satisfy the singularity condition
det(AAT ) = 0
or, upon transformation according to Eq. (3.10), in the angular momentum space.
The presence of singularities in the δ space has to major consequences: (i) When
approaching a singular point the gimbal rate command goes to infinity because of the
increasing norm of the matrix (AAT )−1 , and (ii) it is not possible to implement the
steering logic with an arbitrary torque demand u, because of the presence of the singular
direction along which the CMG cluster is not capable of producing a torque component.
Several singularity avoidance strategies have been developed and discussed in the
literature. The most common approach, derived from the solution derived for problem
of similar nature in robotic manipulators, is based on the general solution to Eq. (3.11),
which can be written in the form
δ̇ = AT (AAT )−1 Ḣ B + γn
where the null vector n satisfies the relation
An = 0
that is, n spans the null space of the jacobian matrix A. In this way, the steering logic
(3.12) is complemented with some null motion, where the parameter γ represents the
amount of null motion that must be added in order to “stay away” from the singularities.
The null vector can be expressed as
n = [I − AT (AAT )−1 A]d
where d is an arbitrary non–zero n–dimensional vector.

Although there are several possible approaches for determining a proper value for γ,
depending on the current value of the gimbal angles, in general it is necessary to add null
motion only when approaching a singularity. In fact, the Moore–Penrose pseudo–inverse
represents the minimum norm solution for Eq. (3.11), and as such it can be considered an
optimal solution, when singularities are not an issue. On the converse, a proper amount
of null motion needs to be added when approaching a singular gimbal configuration, in
order to avoid the increase of ||δ||. In this way, on one side the control torque demand is
always satisfied, and at the same time the control effort is increased, with respect to the
pseudo–inverse steering logic, only when this is necessary to avoid a singularity. For this
reason the parameter γ is usually chosen as a function of the singularity measure
q
m = det(AAT )
where m quantifies the distance from a singular state. As an example, it is possible to set
γ = m6 , for m ≥ 1
γ = m−6 , for m < 1
A variety of heuristic approaches have also been proposed, in order to solve the problem
of the presence of singular states, while keeping the computational burden of the steering
logic to a minimum. A heuristic approach is represented by the so–called singularity

robust pseudo–inverse, written in the form
A∗ = AT (AAT + λ1)−1
where 1 is the identity matrix and λ is a positive scalar, that is tuned as a function of m
as
λ = 0, for m ≥ m0
2
λ = λ0 (1 − m/m0 ) , for m < m0
where λ0 and m0 are small positive constants to be properly selected. In this case, the
pseudo–inverse steering logic is “altered” when approaching the singular state in such a
way so as to avoid the singularity of the matrix (AAT )−1 simply by the addition of a
small contribution to its principal diagonal.

Textbook PDF

Uploaded by

Copyright:

Available Formats

Textbook PDF

Uploaded by

Document Information

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

Textbook PDF

Uploaded by

Copyright:

Available Formats

Chapter 1

Rigid Body dynamics

1.1 Frames of reference and transformation matrices

~v = xê1 + yê2 + zê3

~v = xê1 + yê2 + zê3

~v = (e1,1 x + e1,2 y + e1,3 z)Ê 1 +

This means that the components of ~v in FI can be expressed as a function of those in

• the inverse of an orthogonal matrix L is given by its transpose: L−1 = LT ;

• the determinant of an orthogonal matrix is det(L) = ±1 (and that of a rotation

• if L1 and L2 are orthogonal matrices, their product L1 L2 is an orthogonal matrix.

1.2 Euler’s angles

Figure 1.1: Euler’s angles

Figure 1.2: Planar rotation of amplitude α.

1.2.1 Building the coordinate transformation matrix from ele-

or, in matrix form,

Recalling the properties of orthogonal matrices, it is

L12 = (R(α))−1 = (R(α))T = R(−α)

Elementary rotations for the sequence 3–1–3

that, in matrix form, can be written as:

In an analogous way it is possible to demonstrate that, during the second rotation

that in matrix form becomes:

Finally, the third rotation about ê003 is represented by the transformation

x = x00 cos(Φ) + y 00 sin(Φ)

or, in matrix form:

LBI = R3 (Φ)R1 (Θ)R3 (Ψ)

As the product of orthogonal matrices is an orthogonal matrix, the inverse of which is

Elementary rotations for the sequence 3–2–1

and the final result is given by

A first consequence of the Euler’s angle singularity

As an exercise, demonstrate that LBI depends on ψ − φ only, when Bryant’s rotation

How to build elementary rotation matrices

Figure 1.3: Angular velocity as a funciton of Euler’s angle rates.

Then we fill the diagonal with cos θ:

1.2.2 Angular velocity and the evolution of Euler’s angles

or, in matrix form

or, in explicit form,

Ψ̇ = (ω1 sin Φ + ω2 cos Φ)/ sin Θ

and the inverse relation is

φ̇ = ω1 + (ω2 sin φ + ω2 cos φ) tan θ

1.2.3 The quaternions

that for any nontrivial eigenvector a implies that

The eigenvector relative to the first eigenvalue satisfies the relation

â = a1 ê1 + a2 ê2 + a3 ê3

(a) ê3 (b)

ê003 Ê 3 ê03 ê003

(c) ê3 ê002 (d) ê3 ê002

L11 = a21 + (R21

1.2.4 Evolution of the quaternions

The equivalent matrix form is given by

1.2.5 Quaternions vs Euler angles

1.3 Time derivative of vector quantities

its time derivative is given simply by

1.4 Euler’s equations of motion of a rigid body

Figure 1.4: A rotating rigid body.

ω × ~r )] δm is strictly function of the mass

~ = h1 ê1 + h2 ê2 + h3 ê3

where the symmetric matrix

~ × (~y × ~z ) = (~x · ~z ) ~y − (~x · ~y ) ~z

so that, in the present case, it is