AER506H1F

Spacecraft Dynamics and Control I

A course presented at the
University of Toronto
Notes prepared by
C. J. Damaren
 1

1 Kinematics of Spacecraft Motion

Spacecraft are free bodies; that is they are free to translate and rotate. Mathematically,
they have six degrees of freedom: three in translation and three in rotation. The study of
spacecraft translational motion is described within the framework of orbital dynamics. The
study of its rotational motion is the subject of attitude dynamics. We shall see later on in
the course that, for all intents and purposes, the two classes of motion are uncoupled and
can be treated separately.
Before writing down the equations that describe the dynamics of attitude motion, we must
examine ways of describing it. In other words, we shall study the kinematics of rotational
motion. For example, the position of a point mass can be described with a vector consisting
of three components. Rotational motion is described by expressing the orientation of one
reference frame with respect to another.

Fb
~ 6
' $
u -
3 c.m.





rc 
& %
FI ~  
spacecraft
~ 
6 



 -
O

In the context of spacecraft dynamics, one reference frame is fixed to the spacecraft (F b )
and we keep track of the frame’s motion with respect to a nonspinning, inertial reference
frame (F I ). Orbital dynamics is the study of how the position of the spacecraft centre of
mass (rc ) evolves in time. Attitude dynamics is concerned with the orientation of F b with
~
respect to F I . This course will deal with both subjects. We endeavour to provide a complete
description of a rigid body moving in the gravitational field of the earth.
Kinematics of Spacecraft Motion 2

1.1 Reference Frames and Rotation Matrices

We shall take a vector to be a quantity r. This vector can be expressed in a reference frame
~
as

r = r1 11 + r2 12 + r3 13
13
 r ~ ~ ~ ~
~6
~
11
F1



~ 
1 
 

- = [r1 r2 r3 ] 
 ~2 
12
13
 
~
11 ~
T
~ = r F1 (1)
~

The quantity  
r1
 r2 
r= 

r3
is a column matrix containing the components of r. The quantity
~
 
11
~ 
1
 
F1 = 
 ~2

~

13
 
~
is a column containing the basis (or unit) vectors forming the reference frame F 1 . We shall
refer to F 1 as a vectrix.
~
Note that (1) can also be written as
 
 r1

r = 11 12 13  r2 
 
~ ~ ~ ~ r3
= F T1 r
~

Dot Product
Consider two vectors r and s expressed in the same reference frame F 1 :
~ ~ ~
 
11  
~   s1

1
 
r = [r1 r2 r3 ] 
 ~2 , s = 11 12 13  s2 
  
~ ~ ~ ~ ~


13
 s3
~
Kinematics of Spacecraft Motion 3

The dot product is given by

 
11  
~    s1
1 
 
 ~ 2  · 11 12 13  s2
r · s = [r1 r2 r3 ] 
 

~ ~  ~ ~ ~ s3
13


 ~ 
11 · 11 11 · 12 11 · 13  
~ ~ ~ ~ ~ ~  s1
1 2 · 11 12 · 12 12 · 13
 
= [r1 r2 r3 ] 
~ ~
 s2 
~ ~ ~ ~  

13 · 11 13 · 12 13 · 13
 s3
~ ~ ~ ~ ~ ~
But
11 · 11 = 12 · 12 = 13 · 13 = 1
~ ~ ~ ~ ~ ~
and
11 · 12 = 12 · 13 = 13 · 11 = 0
~ ~ ~ ~ ~ ~
Therefore,
r · s = rT 1s = rT s = r1 s1 + r2 s2 + r3 s3
~ ~
The notation 1 will be used to designate the identity matrix. Its dimension can be inferred
from context.
Cross Product
The cross product of two vectors expressed in the same reference frame is given by:
 
11 × 11 11 × 12 11 × 13  
~ ~ ~ ~ ~ ~  s1
1 ×1 1 ×1 1 ×1 
 
r × s = [r1  ~ 2 ~ 1 ~ 2 ~ 2 ~ 2 ~ 3   s2 
r2 r3 ] 

~ ~  s
13 × 11 13 × 12 13 × 13

3
 ~ ~ ~ ~ ~ ~
0 13 −12  
 ~ ~ ~  s1
−1 0 11 
 
= [r1 r2 r3 ] 
 ~3  s2 

~ ~  

12 −11 0 s3

~ ~ ~  
  0 −r3 r2 s1
= 11 12 13  r3

0 −r1    s2 
 
~ ~ ~ −r2 r1 0 s3
= F T1 r× s
~
Hence, if r and s are expressed in the same reference frame, the 3 × 3 matrix
~ ~
 
0 −r3 r2
× ∆ 
r =  r3 0 −r1 

−r2 r1 0
Kinematics of Spacecraft Motion 4

can be used to construct the components of the cross product. This matrix is skew-
symmetric, that is,
(r× )T = −r×
You should convince yourself that
r× r = 0
where 0 is a column matrix of zeros and

r× s = −s× r

Rotation Matrices
Let us consider two reference frames F 1 and F 2 and express r in each frame:
~ ~ ~
~


r
13 23
r = F T1 r1 = F T2 r2
~ ~
~ ~ ~
BMB 6

B
22
B
 *~

B
 - 12
A ~
A
AA
U
11 21
~ ~

We seek to discover a relationship between the components in F 1 , r1 and those in F 2 , r2 .

Proceed as follows:

F T2 r2 = F T1 r1 (2)
~ ~
F 2 · F 2 r2 = F 2 · F T1 r1
T
~ ~ ~ ~
r2 = C21 r1 (3)

We have defined

C21 = F 2 · F T1 (4)
~ ~
21
~   
2
 
=  2  · 11 12 13
 
~  ~ ~ ~
23
 ~ 
21 · 11 21 · 12 21 · 13
~ ~ ~ ~ ~ ~ 
2 2 · 11 22 · 12 22 · 13
 
= 
~ ~
 (5)
~ ~ ~ ~ 
23 · 11 23 · 12 23 · 13
 
~ ~ ~ ~ ~ ~
The matrix C21 is called a rotation matrix. It is sometimes referred to as a “direction cosine
matrix” since the dot product of two unit vectors is just the angle between them.
Kinematics of Spacecraft Motion 5

The unit vectors in F 2 can be related to those in F 1 by inserting (3) into (2):

F T1 = F T2 C21 (6)
~ ~

Rotation matrices possess some special properties. From Eq. (3),

r1 = C−1
21 r2 = C12 r2

But, from (5), CT21 = C12 . Hence,

C12 = C−1 T
21 = C21 (7)

We say that C21 is an orthogonal matrix because its inverse is equal to its transpose.
Consider three reference frames F 1 , F 2 , and F 3 . The components of a vector r in these
~ ~ ~ ~
three frames are r1 , r2 , and r3 . Now,

r3 = C32 r2 = C32 C21 r1

But, r3 = C31 r1 , and therefore

C31 = C32 C21

Principal Rotations
The most important rotations of one frame with respect to another are those about one of
the coordinate axes. The situation where F 2 has been rotated from F 1 through a rotation
about the 3-axis is shown below. The rotation matrix in this case is

13 ,23 
cos θ3 sin θ3 0

~ 6~
C3 =  − sin θ3 cos θ3 0  (8)
 
22
*~

 0 0 1
 6θ3 - 12

-A ~
θ3 A
AA
U
11 21
~ ~

For a rotation about the 2-axis:

 
13 cos θ2 0 − sin θ2
23 ~ C2 =  0

1 0 
(9)
~ 6 
AAθ2
K sin θ2 0 cos θ2

A
A - 12 ,22
 ~ ~

θ AU
A
2
11  21
~ ~
Kinematics of Spacecraft Motion 6

The rotation matrix for a rotation about the 1-axis is given by

 
13 1 0 0
23 ~ C1 =  0 cos θ1 sin θ1  (10)
 
~ 6
AAθ1 22
K
 *
 ~ 0 − sin θ1 cos θ1
A 
A 6θ1 - 12
~

11 ,21
~ ~

1.2 Angular Velocity

Let frame F 2 rotate with respect to frame F 1 . The angular velocity of frame 2 with respect
to frame 1 is denoted by ω21 . The angular velocity of frame 1 with respect to 2 is ω12 = −ω21 .
~ ~ ~
13 23
 ω
~
21
~ ~
BM 6

B
22
B
 *~
B
 - 12
A ~
A
AU
11 21
~ ~
r
The magnitude of ω21 , | ω21 |= (ω21 · ω21 ), is the rate of rotation. The direction of ω21 , i.e.,
~ ~ ~ ~ ~
the unit vector in the direction of ω21 (| ω21 |−1 ω21 ) is the instantaneous axis of rotation.
~ ~ ~
Observers in the frames F 2 and F 1 do not see the same motion because of their own relative
motions. Let us denote the vector time derivative as seen in F 1 by (˙) and that seen in F 2
by ( ◦ ). Therefore,
◦
Ḟ 1 = 0 , F 2 = 0
~ ~ ~ ~
It can be shown that

2̇1 = ω21 × 21 , 2̇2 = ω21 × 22 , 2̇3 = ω21 × 23

~ ~ ~ ~ ~ ~ ~ ~ ~
or equivalently    
2˙1 2˙2 2˙3 = ω21 × 21 22 23
~ ~ ~ ~ ~ ~ ~
or
Ḟ T2 = ω21 × F T2 (11)
~ ~ ~
We want to determine the time derivative of an arbitrary vector expressed in both frames:

r = F T1 r1 = F T2 r2
~ ~ ~
Kinematics of Spacecraft Motion 7

Therefore, the time derivative as seen in F 1 is

ṙ = Ḟ T1 r1 + F T1 ṙ1 = F T1 ṙ1 (12)
~ ~ ~ ~
In a similar way,
◦ ◦ ◦ ◦
r = F T2 r2 + F T2 r2 = F T2 r2 = F T2 ṙ2 (13)
~ ~ ~ ~ ~
◦
(Note that for nonvectors, (˙) = ( ◦ ), i.e., r2 = ṙ2 ).
Alternatively, the time derivative as seen in F 1 , but expressed in F 2 , is
ṙ = F T2 ṙ2 + Ḟ T2 r2
~ ~ ~
T
= F 2 ṙ2 + ω21 × F T2 r2
~ ~ ~
◦
= r + ω21 × r (14)
~ ~ ~
The above is true for any vector r. The most important application occurs when r denotes
~ ~
position, F 1 is a nonrotating inertial reference frame, and F 2 is a frame that rotates with
a body, particle, etc. In this case, Eq. (14) expresses the velocity in the inertial frame in
terms of the motion in the second frame.
Now, express the angular velocity in F 2 :
ω21 = F T2 ω 21 (15)
~ ~
Therefore, using (12) and (15), (14) becomes
ṙ = F T1 ṙ1 = F T2 ṙ2 + ω21 × r
~ ~ ~ ~ ~
T T ×
= F 2 ṙ2 + F 2 ω 21 r2
~ ~
= F 2 (ṙ2 + ω ×
T
21 r2 ) (16)
~
If we want to express the ‘inertial time derivative’ (that seen in F 1 ) in F 1 , then we can use
the rotation matrix C12 :
ṙ1 = C12 (ṙ2 + ω ×
21 r2 ) (17)

Acceleration
Let us denote the velocity by
◦
v = ṙ = r + ω21 × r
~ ~ ~ ~ ~
The acceleration can be calculated by applying (14) to v:
~
◦
r̈ = v̇ = v + ω21 × v
~ ~ ~ ~ ~
◦◦ ◦ ◦ ◦
= ( r + ω21 × r + ω21 × r) + (ω21 × r + ω21 × ω21 × r)
~ ~ ~ ~ ~ ~ ~ ~ ~ ~
◦◦ ◦ ◦
= r + 2ω21 × r + ω21 × r + ω21 × ω21 × r (18)
~ ~ ~ ~ ~ ~ ~ ~
Kinematics of Spacecraft Motion 8

The matrix equivalent in terms of components can be had by using (13), (15) and by making
the following substitutions:
◦◦ ◦
r̈ = F T1 r̈1 , r = F T2 r̈2 , ω21 = F T2 ω̇ 21
~ ~ ~ ~ ~ ~
The result for the components is
h i
r̈1 = C12 r̈2 + 2ω × × × ×
21 ṙ2 + ω̇ 21 r2 + ω 21 ω 21 r2 (19)

The various terms in the expression for the acceleration have been given special names:
◦◦
r acceleration w.r.t. F 2
~ ◦
2ω21 × r coriolis accleration
◦
~ ~
ω21 × r angular acceleration
~ ~
ω21 × ω21 × r centripetal acceleration
~ ~ ~

Angular Velocity Given Rotation Matrix

Begin with Eq. (6) which relates relates two reference frames via the rotation matrix:
F T1 = F T2 C21
~ ~
Now take the time derivative of both sides as seen in F 1 :
0 = Ḟ T2 C21 + F T2 Ċ21
~ ~ ~
Substitute Eq. (11) for Ḟ 2 :
0 = ω21 × F T2 C21 + F T2 Ċ21
~ ~ ~ ~
Now use (15) to get
0 = ω T21 F 2 × F T2 C21 + F T2 Ċ21
~ ~ ~ ~
T ×
= F 2 ω 21 C21 + Ċ21
~
Therefore, we conclude that
Ċ21 = −ω ×
21 C21 (20)
Given the angular velocity as measured in the frame F 2 , the rotation matrix relating F 1
to F 2 can be determined by integrating the above expression. This is termed ‘strapdown
navigation’ because the sensors that measure ω 21 are strapped down in the rotating frame
F 2.
Eqn. (20) can be written as an explicit function of ω 21 :
ω× −1
21 = −Ċ21 C21
= −Ċ21 CT21 (21)
which gives the angular velocity when the rotation matrix is known as a function of time.
Kinematics of Spacecraft Motion 9

1.3 Attitude Parameterization

To this point we have seen one way of discussing the orientation of one reference frame
with respect to another: the rotation matrix. This requires 9 parameters (they are not
independent). There are a number of other alternatives.
Euler Angles
The orientation of one reference frame with respect to another can also be specified by a
sequence of three principal rotations. One possible sequence is as follows:

1. A rotation ψ about the original 3-axis.

2. A rotation γ about the intermediate 1-axis.

3. A rotation θ about the transformed 3-axis.

13 ,I3 I3 T3 ,23
~6 ~ T3 ~ ~ 6~
~ 6
I2 KAAγ T2 22
~ *~ *~


* A  

 6ψ - 12 A 6γ - I2  6θ - T2
-A ~ ~ -A ~
ψA θ A
AA AAU
U
11 I1 I1 ,T1 T1 21
~ ~ ~ ~ ~ ~

This tranformation is called a 3-1-3 transformation and is the sequence originally used by
Euler. In the classical mechanics literature, the angles are referred to by the following names:

θ = spin angle
γ = nutation angle
ψ = precession angle

The rotation matrix from frame 1 to frame 2 is given by:

C21 (θ, γ, ψ) = C2T CT I CI1

= C3 (θ)C1 (γ)C3 (ψ)
 
cθ cψ − s θ cγ s ψ s ψ cθ + cγ s θ cψ s γ s θ
 −cψ sθ − cθ cγ sψ −sψ sθ + cθ cγ cψ sγ cθ 
=  (22)


sψ sγ −sγ cψ cγ

We have made the abbreviations s = sin, c = cos.

Another possible sequence that can be used is as follows:

1. A rotation θ3 about the original 3-axis (‘yaw’ rotation).

Kinematics of Spacecraft Motion 10

2. A rotation θ2 about the intermediate 2-axis (‘pitch’ rotation).

3. A rotation θ1 about the transformed 1-axis (‘roll’ rotation).

This sequence which is very common in aerospace applications is called the 3-2-1 attitude
sequence or the ‘yaw-pitch-roll’ convention.
In this case, the rotation matrix from frame 1 to frame 2 is given by:

C21 (θ1 , θ2 , θ3 ) = C1 (θ1 )C2 (θ2 )C3 (θ3 )

 
c2 c3 c2 s 3 −s2
=  s 1 s 2 c3 − c1 s 3 s 1 s 2 s 3 + c1 c3 s 1 c2  (23)
 

c1 s 2 c3 + s 1 s 3 c1 s 2 s 3 − s 1 c3 c1 c2
where si = sin θi , ci = cos θi .
The above transformations have singularities. If γ = 0 for the 3-1-3, then the angles θ and
ψ become associated with the same degree of freedom and cannot be uniquely determined.
For the 3-2-1, a singularity occurs when θ2 = π/2. In this case:
 
0 0 −1
π
 sin(θ1 − θ3 )
C21 (θ1 , , θ3 ) =  cos(θ1 − θ3 ) 0 
2

cos(θ1 − θ3 ) − sin(θ1 − θ3 ) 0
Therefore, θ1 and θ3 are associated with the same rotation.
Infinitesimal Rotations
Consider the 3-2-1 transformation, (23), when the angles satisfy |θ1 |, |θ2 |, |θ3 |  1, i.e.,
. .
small angles. In this case, we make the approximations ci = 1, si = θi and neglect products
.
of small angles, θi θj = 0:
 
1 θ3 −θ2
C21 =  −θ3 1 θ1 


θ2 −θ1 1
= 1 − θ× (24)

where
θ T = [θ1 θ2 θ3 ]
It is easy to show that the form of the rotation matrix for infinitesimal rotations (‘small
angle approximation’) does not depend on the order in which the rotations are performed.
For example, show that the same result is obtained for a 2-1-3 Euler sequence.
Euler Parameters
Euler’s Theorem. The most general motion of a rigid body with one point fixed is a
rotation about an axis through that point.
Let us denote, the axis by a = [a1 a2 a3 ]T and assume that it is a unit vector:

aT a = a21 + a22 + a23 ≡ 1 (25)

Kinematics of Spacecraft Motion 11

The angle of rotation is ϕ. We state, without proof, that the rotation matrix in this case is
given by
C21 = cos ϕ1 + (1 − cos ϕ)aaT − sin ϕa× (26)
It does not matter in which frame a is expressed because

C21 a = a (27)

The combination of variables

 
a sin(ϕ/2)
ϕ ϕ  1
η = cos , ε = a sin =  a2 sin(ϕ/2)  (28)

2 2
a3 sin(ϕ/2)

is particularly useful. The four parameters {η, ε} are called the Euler parameters associated
with a rotation. They are not independent because they satisfy the constraint

η 2 + ε21 + ε22 + ε23 = 1

The rotation matrix (26) can be expressed in terms of the Euler parameters as

C21 = (η 2 − εT ε)1 + 2εεT − 2ηε×

1 − 2(ε22 + ε23 ) 2(ε1 ε2 + ε3 η) 2(ε1 ε3 − ε2 η)
 

=  2(ε2 ε1 − ε3 η) 1 − 2(ε23 + ε21 ) 2(ε2 ε3 + ε1 η)  (29)

 

2(ε3 ε1 + ε2 η) 2(ε3 ε2 − ε1 η) 1 − 2(ε21 + ε22 )

Euler parameters are used most often in actual space applications. There are no singularities
associated with them and the calculation of the rotation matrix does not involve trigonomet-
ric functions (compare (29) with (22) and (23)) which is a significant numerical advantage.
The only drawback is the use of four parameters instead of three as is the case with Euler
angles.

1.4 Attitude Solution Given Angular Velocity

In the last section, we showed that the orientation of one frame F 2 with respect to another
F 1 could be parameterized in different ways. In other words, the rotation matrix could be
written as a function of Euler angles or Euler parameters. However, in most applications
the attitude changes with time and we must solve a different problem. Usually, the angular
velocity of one frame with respect to another, ω 21 , is known and one must calculate the
attitude from this information.
An example of this type of problem is Eq. (20),

Ċ21 = −ω ×
21 C21

If the initial condition of C21 is known and ω 21 (t) is given, then the equation can be integrated
to yield the attitude history in the form C21 (t). If C21 (t) is known than one can calculate
the corresponding Euler angles or Euler parameters.
Kinematics of Spacecraft Motion 12

Euler Angles
Consider the 3-2-1 Euler angle sequence and its associated rotation matrix (23). If this is
substituted into (21), one can show (eventually) that
  
1 0 − sin θ2 θ̇1
ω 21 =  0 cos θ1 sin θ1 cos θ2   θ̇2  (30)
  

0 − sin θ1 cos θ1 cos θ2 θ̇3

| {z }
= S(θ1 , θ2 )θ̇

which gives the angular velocity in terms of the Euler angles and the Euler rates, {θ̇1 , θ̇2 , θ̇3 }.
By inverting the matrix S, we arrive at a system of differential equations which can be
integrated to yield the Euler angles assuming ω 21 is known:

θ̇ = S−1 (θ1 , θ2 )ω 21
 
1 sin θ1 tan θ2 cos θ1 tan θ2
 0
=  cos θ1 − sin θ1   ω 21 (31)
0 sin θ1 sec θ2 cos θ1 sec θ2

Infinitesimal Rotations
The angular velocity for a sequence of small rotations can be calculated using (21) and (24):

ω× T
21 = −Ċ21 C21
×
= θ̇ (1 + θ × )
× ×
= θ̇ + θ̇ θ ×

If we assume that the angular rates θ̇i are small and neglect products of small angles and
×
rates, then the above reduces to ω ×
21 = θ̇ or

ω 21 = θ̇ (32)

Therefore, Z t
θ(t) = θ(0) + ω 21 (τ ) dτ
0
For small angles and rates, the components of angular velocity can be integrated directly to
yield the (small) attitude angles. In general, however, one must integrate (20) for C21 or an
equation of the form of (31) for the Euler angles.