Hardware Testing of Attitude Control for the NTNU Test
Hardware Testing of Attitude Control for the
NTNU Test Satellite
Marius Fløttum Westgaard
December 2014
Project Thesis
Department of Engineering Cybernetics
Norwegian University of Science and Technology
Supervisor 1: Professor Jan Tommy Gravdahl
Supervisor 2: Roger Birkeland
i
Abstract
NUTS, the NTNU test satellite, is a student satellite that are being developed at NTNU. It is
part of the Norwegian student satellite program that are being run by the Norwegian center for
space related education. The satellite are going to be carrying a camera that will be used to take
pictures of the earth. To achieve this the satellite needs to have a attitude determination and
control system, ADCS. This project are going to look at the attitude control system.
A prototype for the ADCS have been built earlier, and several control algorithms have been
tested in Matlab. Here the focus is on hardware testing of a control algorithm on the ADCS
prototype. To do this the prototype is being suspended in a set of Helmholtz coils that creates a
strong magnetic field for the satellite to interact with.
Because of the natural limitations of the test setup two tests have been performed. It is found
that the controller is able to start the satellite up from no speed and rotate it in the correct direction, and that it is able to slow the satellite down. These are positive results that indicates that
the control algorithm will work when the satellite is launched.
ii
Acknowledgement
I would like to thank Henrik Rudi Haave and Antoine Francois Xavier Pignède who also worked
on the NTNU test satellite for all the help they have provided during this semester.
Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
i
Acknowledgement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ii
1 Introduction
1
1.1 The NUTS Project . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.3 Contributions of this Project . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2 Background Theory
6
2.1 Reference Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.1.1 Earth-Centered Inertial Frame . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.1.2 Earth-Centered Earth-Fixed Frame . . . . . . . . . . . . . . . . . . . . . . . .
7
2.1.3 Orbit Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.1.4 Body-Fixed Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2 Attitude Representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2.1 The Rotation Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.2.2 Euler Angels and Quaternions . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
3 Satellite Theory
12
3.1 Satellite Dynamics and Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
3.2 Magnetorquers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
3.3 Environmental Disturbances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
4 Control
16
4.1 Nonlinear Controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iii
17
CONTENTS
iv
4.2 Linear-Quadratic Regulator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
4.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
5 Simulation Results
22
5.1 Nonlinear Controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
5.1.1 Nadir Pointing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
5.1.2 Target Pointing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
6 Hardware Testing
28
6.1 Hardware Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
6.2 Software Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
6.3 Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
6.3.1 Velocity Reduction Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
6.3.2 Velocity Increase Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
7 Discussion and Conclusion
37
7.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
7.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
7.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
A Acronyms
41
B Low Pass Filter for Gyroscope Measurements
42
C Matlab Code
44
D C Code
57
List of Figures
1.1 Concept image showing NUTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2 ADCS prototype with satellite frame built in [14] . . . . . . . . . . . . . . . . . . . . .
5
2.1 The body fixed coordinate system. The magnetorquers will be mounted on the
positive sides . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
5.1 Attitude of the satellite during nadir pointing simulation . . . . . . . . . . . . . . . .
25
5.2 Angular velocity of the satellite during nadir pointing simulation . . . . . . . . . . .
25
5.3 Attitude of the satellite during target pointing simulation . . . . . . . . . . . . . . .
27
6.1 The test setup used for this project . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
6.2 Angular velocity of the satellite, with and without the controller turned on, around
the x-axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
6.3 Angular velocity of the satellite, with and without the controller turned on, around
the y-axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
6.4 Angular velocity of the satellite, with and without the controller turned on, around
the z-axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
B.1 Angular velocity data from the gyroscope and some estimates . . . . . . . . . . . .
43
v
List of Tables
3.1 Magnetorquers Specifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
4.1 Attitude Control System Requirements . . . . . . . . . . . . . . . . . . . . . . . . . .
17
5.1 Satellite Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
5.2 Initial Values for Nadir Pointing Simulation . . . . . . . . . . . . . . . . . . . . . . . .
24
5.3 Initial Values for Target Pointing Simulation . . . . . . . . . . . . . . . . . . . . . . .
26
5.4 Desired Values for Target Pointing Simulation . . . . . . . . . . . . . . . . . . . . . .
26
vi
Chapter 1
Introduction
This project thesis is a study of the attitude control system that is being developed for the NTNU
test satellite, NUTS. Much work have already been done, and the NUTS Project is now in a late
stage where the equipment and algorithms needs to be tested and final specifications have to
be determined.
1.1 The NUTS Project
The Norwegian University of Science and Technology test satellite(NUTS) aims to send a double CubeSat [17] to space by 2016. The satellite is built and designed by students at NTNU as
part of the Norwegian Student Satellite Program, ANSAT, which is run by NAROM (Norwegian
Centre for Space-related Education). Two other universities are also in the same program, the
Univeristy of Oslo are building a CubeSat called CubeSTAR and Narvik University College have
built one called HiNCube. There have also been built two CubeSats before in cooperation between these institutions, the NCUBE-1 and the NCUBE-2. Many of the ideas and methods used
in the NUTS project is built on what was learned from the previous CubeSat projects.
The main contributors to the project are final year master students and volunteers from several
different fields of study, including electronics, communications, physics, cybernetics, mechanics and more. Most of the components and parts are being designed and manufactured in-house
at NTNU. As a part of the national student satellite program, the NUTS CubeSat project was ini1
CHAPTER 1. INTRODUCTION
2
tiated in 2010. The project goal is to give hands-on experience to students within different fields
of satellite technology. This includes planning, specification, design, construction, launch and
operation of the satellite. The main parts of the satellite bus are:
• Power system (solar cells and batteries)
• Attitude determination and control system
• Mechanical system (structures and mechanisms)
• TT&C (telemetry, tracking and command)
• Thermal system
• On-board data handling
• Payload and experiment
The NUTS CubeSat is using the AVR32 UC3 MCU as main processor. There are two processors
on board, one for the main computer and one for the communications system. The system is
running the FreeRTOS operating system. The payload of the satellite are going to be a CMOS
camera, and it will also carry two radios, making it possible for anyone with an amateur radio
license to listen to it. The project has its own home page at [11].
In my project I am going to look at attitude control of the satellite. The complete system is
called a Attitude Determination and Control System, ADCS. This paper will look at and test a
control algorithm for the satellite on hardware. The testing will be done on a previously built
prototype of the ADCS that are placed in a satellite frame that are going to be suspended from
a wire in a magnetic field. The work will be built on what previous students have found mathematically and in simulations. There will also be done new simulations for this project to review
the controller and get familiar with it. The overall goal of the satellites ADCS is to point the camera towards the earth. The goal of the tests that are to be performed here is to see whether or
not the control system can change the satellites attitude, and see how it behaves. We will also
investigate whether or not the ADCS can be used to perform target pointing and not just nadir
CHAPTER 1. INTRODUCTION
3
pointing. Target pointing here means pointing the satellites camera at a specific target on Earth,
like Trondheim.
Figure 1.1: Concept image showing NUTS
1.2 Previous Work
Since we are now approaching the final stages of the NUTS project there have obviously been a
lot of work done previously. The ADSC have already been designed and a prototype have been
built, and different control methods, both for detumbling and pointing have been investigated.
It was early on decided that magnetorquers are going to be the actuators for the satellite, these
are small and power efficient, making them ideal for CubeSats. They are very popular, have been
used before and are also now being used by CubeSTAR [16] and HiNCube [12]. Much of the work
done for NUTS on the magnetorquers and control is built on [15].
The ADCS can be split into two parts, a attitude determination part and a control part. To de-
CHAPTER 1. INTRODUCTION
4
termine the attitude we are measuring the earths magnetic field, using a set of photo diodes to
determine the sun vector and using a gyroscope to find the satellites angular speed. Since we
are using measurements of the magnetic field to find the attitude we can not use the magnetorquers to control the satellite at the same time, thus control and attitude determination needs
to be two separate parts. The main focus for this project thesis is the control part of the ADCS.
The attitude determination is being tested in a parallel project, much of the work there is built
on [7]. The control system consists of two parts, detumbling and pointing. When the satellite is
released from the launch vehicle it will have a unknown spin, so the first thing that have to be
done is to stabilize the satellite with the detumbling algorithm. Once the satellite is stable we
can determine its attitude, and start pointing it in the desired direction. There have already been
tested a detumbling algorithm on the prototype [14], and there is also another parallel project
ongoing that are investigating other detumbling algorithms. Thus the focus for this project will
be on the pointing controller that are being used after the detumbling controller is finished.
The prototype for the ADCS was built in [14], there have also been built two Helmoltz coils [8] to
create a strong magnetic field to test the satellite prototype in. Unfortunately the environment
we can create is quite different from the one the satellite will actually operate in, but the tests
are still useful and necessary. The algorithms for pointing have all been designed, some in [15]
and then they where expanded upon in [18]. The newest work done was by [2] who introduced
another controller and ran simulations to compare the controllers. Based on the simulations
previously done a PD-like nonlinear controller seems to be the best choice for the pointing. For
comparison HiNCube uses a nonlinear controller [12](similar to the one recommended in [2])
while CubeSTAR uses a linear-qadratic regulator(LQR) [16]. [2] also investigated LQR but found
that the nonlinear controller would perform better with the expected disturbances.
1.3 Contributions of this Project
The main contribution of this project to the overall NUTS project is hardware implementation
and testing of the nonlinear controller first developed in [15]. It have been implemented on
hardware before, but never on a functioning prototype of the satellite as have been done here.
CHAPTER 1. INTRODUCTION
5
Figure 1.2: ADCS prototype with satellite frame built in [14]
This is the first time the controllers performance have been tested in a test setup.
This project also looks at simulations of the nonlinear controller with the power output and
magnetorquers that are intended for the final specification of the satellite. Hence the simulations done here are probably closer to how the satellite will actually behave then what previous
simulations have been.
Chapter 2
Background Theory
This chapter will cover reference frames and different ways to represent a satellites attitude. It
is assumed that the reader have a technical background since not all concepts are explained in
depth.
2.1 Reference Frames
A reference frame is a coordinate system that are being used to represent the position and orientation of a object. For satellite navigation it is clearly necessary to specify where the satellite is,
and how it is oriented (what its attitude is), but this have to be done with respect to some other
point in space, reference frames allows us to do this. Which frame one choose to use depends
on what one is trying to describe.
2.1.1 Earth-Centered Inertial Frame
The Earth-centered inertial (ECI) frame is a coordinate frame that have its origin at the center
of the Earth. Its axis are defined by the vectors {i } = (x i , y i , z i ). The z i axis is defined as a vector
directed from the center of the Earth towards the celestial north pole. The x i axis points in the
vernal equinox direction. This is a vector pointing from the center of Earth to the center of the
Sun on the vernal equinox. The vernal equinox is the time of year when the Sun crosses the
Earths equatorial plane going from south to north. This is usually on the 20th or 21st of March,
6
CHAPTER 2. BACKGROUND THEORY
7
the sun is then directly over the equator. I.e the center of the Sun lies in the same plane as the
Earths equator. The y i axis is defined by the right hand rule to complete the coordinate system.
2.1.2 Earth-Centered Earth-Fixed Frame
The Earth-centered Earth-fixed (ECEF) frame is like the ECI frame except that it is rotational. Its
axis are defined by the vectors {e} = (x e , y e , z e ). The z e axis is defined as a vector directed from
the center of the Earth towards the celestial north pole. The x i axis points in the direction of 0°
longitude and 0° latitude, thus the frame rotates with Earth around the z e axis with a angular
ad
speed of ωe = 7.2921 × 10−5 rsec
. Again the y e axis is defined by the right hand rule to complete
the coordinate system.
2.1.3 Orbit Frame
The orbit frame is a moving coordinate frame that have its origin at the center of the satellite. Its
axis are defined by the vectors {o} = (x o , y o , z o ). The z o axis is defined as a vector pointing in the
nadir direction (towards the center of the Earth). The x o axis is defined so that it points in the
orbit normal direction,e.g along the satellites linear velocity vector. Again the y o axis is defined
by the right hand rule to complete the coordinate system.
2.1.4 Body-Fixed Frame
The body-fixed frame is also a moving coordinate frame with its origin at the center of the satellite. Its axis are defined by the vectors {b} = (x b , y b , z b ). Unlike the orbit frame the body-fixed
frame is fixed to the satellite, as the name implies, it moves and rotates with the satellite. We will
define it so that z b is the axis of minimum inertia, y b is the axis of maximum inertia and the x b
axis is defined by the right hand rule to complete the coordinate system. x b is the roll axis, y b is
the pitch axis and z b is the yaw axis. The reasoning behind this choice of axis follows from the
stability analysis done in [15], the most important parts of the analysis is repeated in Chapter
4.3. This frame is also referred to as the body frame.
For the hardware test done on the prototype ADCS in this project we have defined the body
CHAPTER 2. BACKGROUND THEORY
8
frame so that the z b axis points from the center of the satellite trough the small coil. The other
axis are also defined so that they point from the center of the satellite and trough the coils, thus
you can find the x b axis and y b axis by using the right hand rule. The camera will be mounted
on the z face of the satellite.
Figure 2.1: The body fixed coordinate system. The magnetorquers will be mounted on the positive sides
2.2 Attitude Representation
The attitude of a satellite (or any craft) describes how the craft is oriented relative to some other
reference frame. In this project Euler angels and Quaternions are being used to represent the
satellites attitude and they will be defined here, as will the rotation matrix.
CHAPTER 2. BACKGROUND THEORY
9
2.2.1 The Rotation Matrix
The rotation matrix is used to rotate or transform a vector from one reference frame to another
or to rotate a vector within a reference frame. A rotation of the vector v from frame a to b is
denoted as v b = Rba v a . The rotation matrix is an element in SO(3) which is formally defined as:
n
SO(3) := R|R ∈ R3 ,
R is orthogonal and d et R = 1
o
(2.1)
A useful way to look at the rotation matrix is to describe it as a rotation β about a unit vector λ
as follows:
h
i
Rλ,β = I3×3 + sin(β)S(λ) + 1 − cos(β) S2 (λ)
(2.2)
Where S is a skew-symmetric matrix defined as:
−λ3 λ2
0
S(λ) = −S(λ)T =
0
−λ1
λ3
,
−λ2 λ1
0
λ1
λ=
λ2
λ3
(2.3)
For more on the rotation matrix see Chapter 6 in [4].
In this project the rotation matrix from the body frame b, to the orbit frame o, can be written as:
x b · x o
Rob =
y b · xo
zb · xo
xb · y o
yb · yo
zb · y o
x b · z o c 11 c 12 c 13
i
h
c
= c c c
=
y b · zo
c
c
1
2
3
21 22 23
zb · zo
c 31 c 32 c 33
(2.4)
The elements of the first matrix are the x, y and z components defining the satellites attitude
in the body and orbit frame, as defined in Chapter 2.1. The elements in the second matrix are
called direction cosines. When we are talking about the satellites attitude we are referring to
how the body frame is oriented relative to the orbit frame. This rotation matrix can thus be used
to describe the satellites attitude.
CHAPTER 2. BACKGROUND THEORY
10
2.2.2 Euler Angels and Quaternions
The most common way to describe the attitude of any system are by using the Euler angles, roll
(φ), pitch (θ) and yaw (ψ). As mentioned earlier, for this satellite we have that roll describes rotation around the x b axis, pitch describes rotation around the y b axis and yaw describes rotation
around the z b axis. The Euler angel vector is:
φ
Θ=
θ
ψ
(2.5)
One can with this notation easily rewrite the already mentioned rotation matrix into a set of
principal rotations defined as follows
0
0
1
Rx,φ =
0 cos(φ) − sin(φ)
0 sin(φ) cos(φ)
cos(θ) 0 sin(θ)
R y,θ =
0
1
0
− sin(θ) 0 cos(θ)
cos(ψ) − sin(ψ) 0
Rz,ψ =
sin(ψ) cos(ψ) 0
0
0
1
(2.6)
(2.7)
(2.8)
Rx,φ then is a rotation of φ around the x-axis, R y,θ is a rotation of θ around the y-axis and Rz,ψ is
a rotation of ψ around the z-axis. You can then write the full rotation matrix as
Rca = Rx,φ R y,θ Rz,ψ
(2.9)
This attitude representation is easy to understand, but you get a problem with singularities at
some angles, for example at a pitch angle of θ = ±90°. At this angle a singularity appears in the
rotation matrix, and it can not be used. For our satellite system, that can rotate more than 90° in
CHAPTER 2. BACKGROUND THEORY
11
all directions, this is then not a suitable way to describe the attitude in a controller.
The alternative used for this project are quaternions, they are harder to interpret directly but
have no singularities and are thus much better to use for a satellite. A quaternion q is a comh
iT
plex number with one real part η and three imaginary parts ² = ²1 , ²2 , ²3 . The set Q of unit
quaternions can be formally defined as
½
T
Q := q|q q = 1,
h
q = η, ²T
iT
,
¾
² ∈ R and η ∈ R
3
(2.10)
The real and imaginary parts of the unit quaternion are defined as
h
where λ = λ1 , λ2 , λ4
iT
µ ¶
β
η := cos
2
(2.11)
µ ¶
h
iT
β
² = ²1 , ²2 , ²3 = λ sin
2
(2.12)
is a unit vector that satisfies
²
λ = ±p
²T ²
if
p
²T ² 6= 0
(2.13)
This is the same λ and β as in Equation 2.2. For more on Euler angels and quaternions see
Chapter 6 in [4].
Chapter 3
Satellite Theory
In this chapter the satellite dynamics and kinematics will be stated. The theory behind the magnetic actuators and how they work will be explained and we will look at the environmental disturbances that can affect the satellite.
3.1 Satellite Dynamics and Kinematics
Assuming that we can treat the satellite as a rigid body its equation of motion for rotational
movement can be expressed as [15]:
˙ bib + ωbib × (Iωbib )
τb = Iω
(3.1)
Where ωbib is the angular velocity of the satellite in the body frame with respect to the ECI frame
given in the body frame. I is the satellites moment of inertia matrix around its mass center
and τb is a moment affecting the satellite, given in the body frame. The angular velocity of the
satellite ωbib can be expressed as:
ωbib = ωbob + ωbio = ωbob + Rbo ωoio
Where ωoio = [0 − ωo 0] and ωo =
q
GM
R o3
(3.2)
is the satellites orbital speed. G is the gravitational con-
stant, M is the mass of the earth and R o = R s + R e , where R e is the earths radius and R s is the
satellites altitude. Here ωbob is the angular velocity of the satellite in the body frame with respect
12
CHAPTER 3. SATELLITE THEORY
13
to the orbit frame, given in the body frame
The kinematic differential equation for quaternions are given as [3]:
1
η˙ = − ²T ωbob
2
1
²˙ = (ηI + S(²))ωbob
2
(3.3)
(3.4)
This Equation, along with Equation 3.4 and 3.2 is all that is needed to model the satellite and
how it will react to torque input.
3.2 Magnetorquers
The satellite uses magnetic actuators, called magnetorquers. These magnetorquers are coils
that have been winded up a number of times, when a current runs trough them they will create a
magnetic moment that interacts with the Earths magnetic field and creates torque. The satellite
will have three magnetorquers mounted such that you have one on the z side, one on the x
side and one on the y side. Running the current trough the coils creates a magnetic moment m
perpendicular to the coil that reacts with the magnetic field B and creates torque according to:
τb = mb × Bb
(3.5)
Since all the measurements and calculations are either being done in or transformed to the body
frame we use the body superscript. To find the direction of the magnetic moment one can look
at the current and use the right-hand rule. The magnetic field created by the coils will try to
align itself with the local geomagnetic field, this will create the torque, making the satellite turn.
The magnetic moment is given by:
N x i x A x
mb =
N y i y A y
Nz i z A z
(3.6)
CHAPTER 3. SATELLITE THEORY
14
Where N is the number of windings in the coil, A is the area the coil covers and i is electric
current. A problem that comes with using magnetorquers is that it gives us an underactuated
system if one of the axis are aligned with the local magnetic field. If that happens and we want
to turn the satellite around that axis the coils will provide no torque and the system is underactuated. This leads to the system being uncontrollable whenever one of the axis are parallel to
the local magnetic field.
For the test setup of the satellite the magetorquers have the specifications given in Table 3.1. For
more information on how the magnetorquers where built and how these values where found see
[14]. The table are included again here since some of these values are needed in the simulations.
Notice that the X- and Y-coils are identical, while the Z-coil is different, this is because we want
the actuators to create the same magnetic moment when we apply the same power to all three.
Since A z < (A x = A y ) we need to have either i z > (i x = i y ) and/or N z > (N x = N y ) to produce the
same magnetic moment, as can be seen from Equation 3.6.
Table 3.1: Magnetorquers Specifications
Parameter
X and Y-coil
Z-coil
Windings
155
238
R theoretical
26.7 Ω
19.2 Ω
R measured
27.2 Ω
18.5 Ω
f
590 Hz
580 Hz
L
8.53 mH 9.13 mH
τ
313 µs
494 µs
For the test setup the magnetic actuators interacts with the magnetic field generated by the
Helmholtz coils. The Helmholtz coils are just two large coils, they can produce a magnetic field
of about 1.3 milliTesla, which is much higher then the average field strength in a low Earth orbit.
This magnetic field should be strong enough to allow us to test the satellites pointing controller.
The Helmholtz coils are being driven by a 10 Volt power supply. Unfortunately the magnetometer on the satellite only have a operating range of ±0.8 milliTesla [14], thus we have to keep the
field generated by the Helmholtz coils below that magnitude.
CHAPTER 3. SATELLITE THEORY
15
3.3 Environmental Disturbances
Even tough the satellite are going to operate in a low Earth orbit there are several environmental
factors that are going to affect it and create disturbance torques. The most important disturbances are
• Gravity
• Aerodynamic torque
• Solar radiation torque
• Internal magnetic dipole moments and other electrical noise
The most important to model is gravity, since, as we are going to discuss in Chapter 4.3, the
controller is based on the satellite working with the torque generated by the gravity of Earth.
Simulations have been done with all these disturbances in [2]. Here we are only going to repeat
the simulations with gravity.
The gravitational pull of the earth will act as a torque on any satellite in orbit that are non symmetric in the form of a gravity gradient torque. For a unsymmetrical satellite the Earths gravity
will pull differently on different parts of the satellite, creating torque according to [6]:
τbg = 3ω2o S(c3 )Ic3
(3.7)
Where ωo is the satellites orbital speed, and c3 comes from the rotation matrix given in equation
2.4.
Chapter 4
Control
The payload of the satellite is going to be a camera, the objective being to take pictures of Earth.
Originally the payload was supposed to be a IR camera, but this have been changed due to cost
considerations. To be able to take these pictures the satellites camera needs to be pointed downwards towards Earth, pointing the camera straight down towards the center of Earth is also
called pointing it in the nadir direction, or nadir pointing. The originally planned IR camera
also required us to point the satellite in the nadir direction, so the stability and pointing demands that have been found in previous work is still valid. The satellite does not need to be
pointing straight down, small deviations would still allow us to take pictures of the Earth. The
requirements for the pointing accuracy of the ADCS given in [2] are repeated here in Table 4.1.
The B-dot detumbling controller will not be presented here, see [2] and [14] for more information on that. The controller presented here will be used after the B-dot controller have been
able to stop the satellites rotation (or at least almost stop it) to point the satellite towards Earth.
The controller are based on the work done in [15], [18] and [2]. Since the last work was done on
the controller in [2] the camera have as mentioned been changed from a IR-camera to a regular
camera. With the new camera we want to take pictures of Trondheim, we will therefore investigate if the controller can be modified to point the camera towards a given target instead of just
in the nadir direction.
16
CHAPTER 4. CONTROL
17
Table 4.1: Attitude Control System Requirements
Roll
Pitch
Yaw
Pointing Accuracy [deg]
±25
±25
−
−3
−3
Angulvar Velocity [rad/sec] ±1 · 10
±1 · 10
±1 · 10−3
4.1 Nonlinear Controller
To get the satellite to point in the nadir direction we need to align the body frame with the orbit
frame. Using quaternions to express the attitude of the satellite this orientation can be expressed
as q = [1 0 0 0]T . The quaternion are a measure of how the satellites body frame deviates from
the orbit frame.
If we want to point the satellite at a certain target we have to look at the quaternion error. This
is, as the name suggests, a measure of the difference between two quaternions. Say that we have
two quaternions q1 and q2 , the difference q0 is then given by q0 = q−1
1 · q2 . Define qd esi r ed and
qd et er mi ned as the attitude we want the satellite to have and the current attitude of the satellite
respectively, the error is then given by [19]:
qer r or = q−1
d et er mi ned · qd esi r ed
(4.1)
When this takes the form qer r or = [1 0 0 0]T the current attitude is aligned with the desired attitude. This is similar to aligning the body frame with the orbit frame.
The desired nonlinear controller was found in [15] and is given by:
τbd = −p² − d ωbob ,
p, d > 0
(4.2)
Again ωbob is the angular velocity of the satellite in the body frame with respect to the orbit frame,
given in the body frame. ² is the last three elements of the quaternion vector as defined in Equation 2.12. ² represents the satellites current attitude for nadir pointing, for target pointing ² are
the three last elements of the quaternion error vector.
CHAPTER 4. CONTROL
18
The magnetic torque τbd will with this equation be in the plane perpendicular to Bb , and in
reality τbd will not be in that plane very often. So we have to project τbd to the plane that is
perpendicular to Bb , this gives us the correct moment in the vector space. This procedure is
described in detail in [15], the resulting equation for the magnetic moment is:
mb = −
1
³
|Bb |2
p(Bb × ²) + d (Bb × ωbob )
´
(4.3)
The magnetic torque becomes:
τb = −
1
|Bb |2
³
´
p(Bb × ²) + d (Bb × ωbob ) × Bb
(4.4)
4.2 Linear-Quadratic Regulator
In [2] a Linear-Quadratic Regulator was also investigated. It has not been time to deal with it
properly in this project, but it is mentioned here since it have been done some investigation
into its properties. For an actual description of the controller see [2].
On of the steps of the optimal LQR controller is solving the following equation for P:
0 = AT P + PA − PBR−1 BT P + Q
(4.5)
Where Q and R are weighting matrices while A and B comes from the linear system x˙ = Ax + Bu.
This is the algebraic Riccati equation, and solving it is quite computationally heavy.
The linearized state matrices for the system, A and B, are sparse. So it has been investigated if
this could lead to a explicit solution fo the algebraic Riccati equation, unfortunately none have
been found. After doing some research into how this have been done on other CubeSats it have
been fond that most of them uses a constant gain LQR such as the one used for CubeStar [16].
If this controller is to be implemented and tested on hardware a constant gain LQR is probably
the way to go. This is something that can be investigated in a future project.
CHAPTER 4. CONTROL
19
4.3 Stability
We can analyze the stability properties of the satellite by deriving a Lyapunov function and applying Lyapunov theory. A Lyapunov function based on the satellites potential and kinetic energy was found in [15]:
¤ 1 2£
¤
3 £
1
2
2
2
2
V = (ωbob )T Iωbob + ω2o (I x − I z )c 13
+ (I y − I z )c 23
+ ωo (I y − I x )c 12
+ (I y − I z )c 32
2
2
2
(4.6)
The potential energy is dominated by the gravitational force caused by the Earth [2].
The time derivative of this function is found to be:
V˙ = (ωbob )T τb
(4.7)
h
i
The state vector used for the stability analysis is defined as x = (ωb )T c 13 c 23 c 12 c 32 . We want
ob
to use Lyapunov theory to prove that x = 0 is a uniformly stable equilibrium point. The theory behind this can be found in [9]. In short the requirements for stability is that V have to be
positive definite, that is V = 0 for x = 0, and V > 0 for x 6= 0. The time derivative, V˙ , have to be
negative semidefinite, that is V˙ ≤ 0, but here the equality is not constrained to only x = 0.
V is positive definite for I y > I x > I z , the function will only equal zero for x = 0 and be larger
than zero for any other values of x. Thus the satellite need to be built in such a way that it satisfies the inequality:
I y > Ix > Iz
(4.8)
If we provide no torque from the coils we will get the trivial solution of the Lyapunov derivate
V˙ = 0, which is enough to satisfy the Lyapunov conditions, making the equilibrium point x = 0
uniformly stable. This tough will correspond to more than one attitude, all the attitudes where
c 13 = c 23 = c 12 = c 32 = 0 will be a stable attitude. Since the rotation matrix is orthogonal and have
a determinant equal to one we can find which attitudes this corresponds to. We can easily find
that the satellite (without any torque produced by the magnetorquers) is stable for any attitude
where z b is parallel with z o and x b is parallel with x o . If I x = I y we will have stability for every
CHAPTER 4. CONTROL
20
attitude where z b is parallel with z o .
If the satellite is built according to Equation 4.8, the gravitational pull of the Earth would in
theory pull the satellite into one of these orientations. Thus the gravity gradient torque in Equation 3.7 is a stabilizing factor on our satellite.
Looking at the case where we do actually use the magnetorquers to generate torque the time
derivative of the Lyapunov function can be written as (by inserting 4.4 into 4.7):
V˙ = −
1
|Bb |2
p(Bb × ²)T (Bb × ωbob ) −
1
|Bb |2
d (Bb × ωbob )T (Bb × ωbob )
(4.9)
The second term will be less than zero for (Bb × ωbob ) 6= 0, unfortunately ² can be both positive
and negative, which makes it hard to say anything about the first term. If this term is positive
and dominates the second term we can not guarantee that V˙ will stay negative. Thus we can not
conclude that V˙ is negative semidefinite.
By looking at the desired control law, equation 4.2, and the satellite dynamics, equation 3.1, it is
possible to find a constraint on the gain p so that the term −p² dominates the angular velocity
and the gravitational moment [15]. This will be true if we chose p such that:
p > 8ω2o (I y − I z )
(4.10)
We want this term to dominate the other terms so that the satellite only will be stable for the
attitude we want, that is the body frame aligned with the orbit frame, and not any of the other
stable equilibriums found above. This however does not mean that we can prove uniform stability for the controller, in fact this have not been shown in any of the projects related to NUTS
so far.
Soglo, [15], goes on to linearize the system around ωbob = 0, ² = 0 and η = 1 and look at the
systems eigenvalues. We know that a linear system is stable if all the real parts of the eigenvalues are less than or equal to zero, and asymptotically stable if the eigenvalues are all less than
CHAPTER 4. CONTROL
21
zero. The eigenvalues have been calculated numerically in Matlab and found to be negative, but
the system parameters in [15] is very different from the ones that belong to the NUTS satellite.
The simulations that have been done on this system, in [18], [2] and [1] have all shown satisfactory performance of the controller. The analysis that have been done shows that the satellite
do have certain stable orientations, and steeps have been taken to try to make sure that the attitude will converge to the right one. We thus have reason to believe that the satellite with the
controller is globally asymptotically stable, but unfortunately no analytic proof exists.
Chapter 5
Simulation Results
In this chapter we will look at numerical simulations done in Matlab to test the performance of
the controller. The code used is based on that of [18] and [2]. To do the calculations from quaternions to Euler angels and back we are using the MSS toolbox. The toolbox have been developed
for control and simulation of marine craft, but the functions for dealing with quaternions, Euler
angels and rotation matrices are very useful for any application that are operating with these
parameters. For more on the MSS toolbox see the Marine Systems Simulator homepage [10].
We are using Matlabs built in od e45 function to simulate the system, which is given by the state
matrix:
h
x = (ωb )T
ib
η ²T
λ W
iT
(5.1)
This is the state matrix first used in [15] and the one that have been used throughout the entire
NUTS project.
ωbib is the angular velocity of the satellite in the body frame with respect to the ECI frame,
given in the body frame, η is the first element of the unit quaternion vector, ² are the three
last elements of the unit quaternion, λ is the latitude and W is the satellites power consumption(originally called p). The last state does not have any effect on the satellites orientation
or control, but are included to monitor the power usage. It can be related to watt-hours as
Wh =
W
.
3600
We are not actually using it to monitor power usage in this project, but the function
can be nice to have to compare the power consumption of different controllers. This have been
22
CHAPTER 5. SIMULATION RESULTS
23
done before, and might come in useful at a later stage so the code part have not been removed.
£
¤
We set the initial conditions for the satellite in ωbob and Euler angels φ θ ψ , which are the attitude variables of the satellite in the body frame with respect to the orbit frame, and transform
them to the variables in the state matrix using the MSS toolbox. The initial latitude, λ, and power
consumption, W , is set to zero in all the simulations. The reason we need to include the latitude,
λ is that it is needed in the IGRF model that are being used to calculate the magnetic field, for
more information see [18].
All the Matlab code can be found in Appendix C. The satellite parameters can be found in Table
5.1. Since the satellite have not been built yet all these variables are based on assumptions and
approximations. The inertia are calculated assuming uniform mass distribution, then changed
a bit to satisfy the stability requirement I y > I x > I z .
Table 5.1: Satellite Parameters
Parameter
Value
m [kg]
2.6
2
I x [kg m ]
0.0098
I y [kg m 2 ]
0.0108
2
I z [kg m ]
0.0043
Voltage [V]
3.3
Rs[km] (orbit height )
600
Orbit Period [min]
96.54
5.1 Nonlinear Controller
To test the controller we first need to find the gains. Let us first look at the requirement given by
Equation 4.10, with our parameters and assuming a orbit height of 600 km we get that p > 6.292·
10−8 . It is also stated in [15] that the value for p should be small. We can see from simulations
done in [2] and [18] that the gain d is usually a factor of 1000 larger than p. After doing new
simulations this assumption seems to be good.
CHAPTER 5. SIMULATION RESULTS
24
5.1.1 Nadir Pointing
After trying and testing lots of different values the best gains seem to be p = 7 · 10−8 and d =
7 · 10−5 . We will present results of just one simulation here, but all combinations of positive and
negative signs on the initial values have been tried, as well as different values. If you compare
the values used for p and d here to the ones used in [2] you can see that these are higher, this is
because it was discovered that the satellite did not reach the desired pointing direction for some
initial values with the gains used in [2]. The initial values used here are shown in Table 5.2.
Table 5.2: Initial Values for Nadir Pointing Simulation
ωbob = [0.05 − 0.02 0.045]
φ = 80◦
θ = −50◦
ψ = 130◦
The amount of time the controller needs to stabilize the satellite in the nadir pointing direction is heavily influenced by the initial parameters, in the simulations it have varied from 1.5
orbits and up to 5 orbits. Of course not all possible initial conditions have been tried and tested,
but enough to make us confident that the controller will function as we want. The result from
these initial values can be seen in Figure 5.1. Here it takes almost 3 orbits for the satellite to
reach the desired attitude.
As can be seen in Figure 5.2 the controller dampens out the angular velocities of the satellite
much faster than it reaches the desired attitude. This is as expected since the dampening term
d is a factor of 1000 larger than p.
Yaw (z) [degrees]
Pitch (y) [degrees]
Roll (x) [degrees]
CHAPTER 5. SIMULATION RESULTS
25
Euler Angeles
200
0
−200
0
0.5
1
1.5
Time [orbits]
2
2.5
3
0
0.5
1
1.5
Time [orbits]
2
2.5
3
0
0.5
1
1.5
Time [orbits]
2
2.5
3
100
0
−100
200
0
−200
Figure 5.1: Attitude of the satellite during nadir pointing simulation
zcomp [rad/sec]
ycomp [rad/sec]
xcomp [rad/sec]
Angular Velocities, ωB
OB
0.05
0
−0.05
0
0.5
1
1.5
Time [orbits]
2
2.5
3
0
0.5
1
1.5
Time [orbits]
2
2.5
3
0
0.5
1
1.5
Time [orbits]
2
2.5
3
0.05
0
−0.05
0.05
0
−0.05
Figure 5.2: Angular velocity of the satellite during nadir pointing simulation
5.1.2 Target Pointing
For the target pointing variant of the nonlinear controller we will use the same initial values,
they can be found in Table 5.3. The desired pointing direction can be found in Table 5.4.
CHAPTER 5. SIMULATION RESULTS
26
Table 5.3: Initial Values for Target Pointing Simulation
ωbob = [0.05 − 0.02 0.045]
φ = 80◦
θ = −50◦
ψ = 130◦
Table 5.4: Desired Values for Target Pointing Simulation
φ = 30◦
θ = 30◦
ψ = 0◦
Several other initial conditions and desired pointing directions have also been tried. There is
no point in rotating the satellite around the z-axis. What we want to do is change the satellites
attitude around the x- and y-axis to tilt it upwards and sidewards toward the horizon of the Earth
to be able to point the camera at a certain location on the Earths surface.
The result can be found in Figure 5.3. It can be seen that the satellite does not reach the desired attitude. This is as expected from the stability analysis done previously, where we found
that the satellite only had certain stable attitudes. We also know from simulations that the magnetorquers provide approximately the same amount of torque as the gravity forces working on
the satellite. When we are doing nadir pointing we are working with this gravity gradient torque,
but when we try to do target pointing we are working against it. Thus we can not create enough
torque to put the satellite into other stable attitudes than the ones found for the satellite when
it generates no magnetic moment.
Yaw (z) [degrees]
Pitch (y) [degrees]
Roll (x) [degrees]
CHAPTER 5. SIMULATION RESULTS
27
Euler Angeles
200
0
−200
0
2
4
6
8
10
Time [orbits]
12
14
16
0
2
4
6
8
10
Time [orbits]
12
14
16
0
2
4
6
8
10
Time [orbits]
12
14
16
100
0
−100
200
0
−200
Figure 5.3: Attitude of the satellite during target pointing simulation
Chapter 6
Hardware Testing
The purpose of this project is to test the nonlinear control algorithm on hardware using the
ADCS prototype developed in [14] and a set of Helmholtz coils to create a magnetic field. In this
chapter we go trough the tests that have been performed and the issues that came up during
testing. To record the data we use the Bluetooth device already connected [14] and LabView. The
LabView software have been expanded upon a bit from [14] to also include transfer of estimation
data [5].
Figure 6.1: The test setup used for this project
28
CHAPTER 6. HARDWARE TESTING
29
6.1 Hardware Issues
Originally the satellite was attached to the coils trough a magnetic bearing and string. This is
fine for testing detumbling and the attitude determination algorithms, but provides to much
friction to be able to properly test the control algorithm. If the satellite is hanging still it is not
able to overcome the friction in the magnetic bearing and start moving again by itself by using
the magnetorquers. Thus another alternative for suspending the satellite had to be found.
For some of the tests we ended up using a fishing string. The string was fastened to the Helmholtz
coils trough the magnetic bearing, and then the satellite prototype was attached to the other
end of the string. The fishing line then works as a torsion spring with a very low spring constant,
meaning that it can rotate around itself easily. This allowed us to test what happens when the
satellite is not rotating but the attitude is wrong.
Another problem is that the satellite can not create enough torque to overcome gravity, this
means that the satellite is not able to tilt itself upwards or downwards in any way, it can only
spin perpendicularly to the wire it is suspended from. This means that we can only test one axis
at a time.
6.2 Software Issues
To test the control algorithm properly we also need a functioning and reliable attitude determination system. This system will provide the controller with the attitude, in the form of a quaternion, and the angular velocity. A extended Kalman filter for attitude determination was tested in
a parallel project, unfortunately this estimator can be unreliable when the satellite spins or is to
far away from the desired attitude and it does not provide any estimate of the angular velocity,
just the attitude. Thus we need to get the angular velocity directly from the gyroscope. For more
information on the attitude determination system see [5].
A issue with the C code implementation itself is setting a output gain. The code is implemented
using floating point values, these then needs to be transformed to a integer value before being
CHAPTER 6. HARDWARE TESTING
30
sent to the coil output function. This transformation means multiplying the floating point value
with a number, hereafter called the output gain, to get a integer value. The coil output function
controls how long the available current will go trough the magnetorquers and in what direction
[14]. Originally the range of the integer value was from -2047 to 2047. To give a larger output
range this have been changed to the maximum value that a 16 bit integer can hold, e.g -32,677
to 32,767.
This was deemed to be an advantage after looking at the Matlab simulation outputs. The magnetic moments set up by the controller in the simulations varies a lot in size, ranging from the
max value 0.2708Am 2 to as low as 1 · 10−7 Am 2 . To be able to capture as much of this range as
possible it was decided to increase the integer value that the coil output function could take.
We first tried to tune the output gain such that setting the output to the max value of 32 767
would be done when the controller wanted to output the maximum magnetic moment we can
set up, 0.2708Am 2 . Thus the output gain can be found from:
g ai n =
i nt 16_t max
b
m max
=
32767
= 121001
0.2708
(6.1)
Unfortunately this gain is not high enough, since we would get a value of less than zero for many
of the moments that the controller outputs. This is a problem since we could end up getting a
value off less than zero when the satellite is standing still with a attitude that is different from
the desired one. Since the coil output function only takes integer values it would not set up any
magnetic moment, and the satellite would never get to the correct attitude.
6.3 Testing
Because of the issues mentioned above testing of the control algorithm is hard, and we have
therefore settled on two tests that are to be performed. As can be seen from Figure 5.2 the controller quickly dampens out the satellites initial spin (if any) and then continues slowly to the
desired attitude. Thus we need to check that the controller can slow the satellite down. Sec-
CHAPTER 6. HARDWARE TESTING
31
ondly we want to see if the controller can start from zero initial speed and set up a spin in the
correct direction to reach a desired attitude.
We found that the different tests had to be performed with different output gains. This comes
from the fact that the p and d gains in the controller are so different, Equation 4.3. When we are
testing velocity reduction we do not want the output gains to be saturated all the time. But to
get enough magnetic moment when testing velocity increase we need the output gain set quite
high, so that it can reach the saturation limit (32767) even when there are no speed.
We ended up using the following gain for the velocity increase testing:
g ai n = 12100100000000
(6.2)
For the velocity decrease testing we used:
g ai n = 121001000000
(6.3)
These gains are simply found by testing, we took the gain found in Equation 6.1 and added zeros
to it until the controller worked as we wanted it to.
Clearly we can not use different gains for different cases on the actual controller when the satellite is launched. But to get the satellite to start moving by itself in our lab set up we needed the
output gains to be close or equal to the maximum value. This is far from the case in the simulations where the satellite can start spinning from a very small magnetic moment. To avoid the
outputs being constantly saturated we had to use a lower gain for the velocity decrease testing.
A good value for the output gain needs to be found before the satellite is launched. However
finding a good value for this gain is not straight forward.
In [2] it was developed a low-pass filter for the magnetometer measurements, all of the tests
done here have been tried both with and without this filter. The final results are presented without the filter since it did not seem to have any effect on the performance. The magnetometer
CHAPTER 6. HARDWARE TESTING
32
data is really smooth as it is, so the filter did not really do much.
The C code can be found in Appendix D. Input-output testing of the C code against the Matlab code have been done for a wide variety of input combinations. The outputs from the Matlab
code and the C code have then been compared. The outputs from the two implementations
always matched, which is a strong indication that the code is implemented correctly in C (assuming the Matlab implementation is correct).
6.3.1 Velocity Reduction Testing
Here we are going to look at the controllers ability to reduce any spin the satellite might have.
From simulations we have seen that the controller tries to damp out the angular velocity initially, and after that move the satellite slowly to the desired attitude. This means that the controller works almost like a detumbling controller.
For this test we are using the state estimator even tough it does not necessarily provide the correct data, this is to differentiate the nonlinear controller from a detumbling controller. We could
also have given the controller a constant quaternion argument, but using the constantly changing attitude estimate from the attitude determination system will be closer to how the system
are really going to work even tough the estimate might be wrong.
We are using the magnetic bearing to suspend the satellite for this test. The friction in the bearing is not a problem for showing a reduction in angular velocity, and actually provides less error
sources here than the fishing string would have done. The fishing string would get wounded up
and start providing a torque in the opposite direction, whereas the magnetic bearing provides a
almost constant friction when spinning. This will affect the satellite in the same way with and
without the controller turned on.
The figures showing the results for the different axis can be seen below, Figure 6.2, Figure 6.3
ad
and Figure 6.4. The satellite was spun up to about 5 rsec
, which is the fastest angular speed the
gyroscope can measure. We have then looked at what happens when the speed starts dropping
CHAPTER 6. HARDWARE TESTING
33
ad
below 5 rsec
. We can clearly see that the angular velocity is reduced faster with the controller on,
than with the controller off. These results are similar to what was found in [13], a parallel project
where testing was done on different detumbling controllers. One of the detumbling controllers
tested in [13] was:
mb = −d (Bb × ωbob )
1
(6.4)
|Bb |2
Which is just the dampening term for the nonlinear controller that are being tested here.
Spin around the x−axis
5
Controller off
Controller on
4
rad/s
3
2
1
0
−1
0
5
10
15
20
25
Time(s)
30
35
40
45
Figure 6.2: Angular velocity of the satellite, with and without the controller turned on, around
the x-axis
CHAPTER 6. HARDWARE TESTING
34
Spin around the y−axis
5
Controller off
Controller on
4
rad/s
3
2
1
0
−1
0
10
20
30
40
50
60
70
Time(s)
Figure 6.3: Angular velocity of the satellite, with and without the controller turned on, around
the y-axis
Spin around the z−axis
5
Controller off
Controller on
4
rad/s
3
2
1
0
−1
0
5
10
15
20
Time(s)
25
30
35
40
Figure 6.4: Angular velocity of the satellite, with and without the controller turned on, around
the z-axis
CHAPTER 6. HARDWARE TESTING
35
6.3.2 Velocity Increase Testing
Here we are looking at whats going to happen when the satellite is hanging still and the orientation is wrong. We do this by setting up the input to the controller so that its attitude is 45°
wrong around the axis we want to test. The satellite should then start spinning around that axis
in the direction that provides the shortest rotation to reach the correct attitude. This means that
we are not using the estimator, we are only feeding the controller with a quaternion argument
corresponding to this attitude. We use the following arguments for the different axis:
qx−axi s = [0.9239 0.3827 0 0]T
(6.5)
q y−axi s = [0.9239 0 0.3827 0]T
(6.6)
qz−axi s = [0.9239 0 0 0.3827]T
(6.7)
This corresponds to a 45° error in the positive rotation direction.
We have also had to input false gyroscope data to the controller for this test. When the satellite is not rotating the gyroscope measurements should be zero, unfortunately this is not the
ad
case. The data that are output from the gyroscope are extremely noisy, varying from −0.06 rsec
ad
, this is a problem since the gain d is so much higher than the gain p. What happens
to 0.06 rsec
is that the coil output starts oscillating around zero so that we get no consistent torque in any
direction and the satellite does not begin to rotate.
A low pass filter have been investigated, see Appendix B, but it did not reduce the noise sufficiently. This is less of a problem when the satellite rotates, since the overall trend then will be
in the correct direction and the calculated magnetic moment will not just oscillate around zero.
Thus we are feeding the controller with a constant vector for angular velocity ωbib = [0 0 0]. We
are hoping that when a estimator that also provides estimates of the angular velocity is implemented these values will be less noisy and will work with the controller.
In this part we are using the fishing line to suspend the satellite in the Helmholtz coils. Since
CHAPTER 6. HARDWARE TESTING
36
this works as a spring we have to be careful that the movement we are seeing does not come
from any moment originating in the fishing line itself, therefore we let the satellite hang for a
while before turning on the controller.
Since the estimator is not reliable enough we unfortunately have no plots of this, but it have
been observed that the satellite starts to rotate in the correct direction. We could have plotted
the angular velocity directly from the gyroscope, but the speeds are so small that this is mostly
noise, making the plots meaningless. Videos of the satellite are included in the appendix in the
folder Tests/VelocityIncrease.
In these videos the satellite have been hanging untouched from the fishing line for a while so
that the velocity is greatly reduced. Unfortunately the line produces so little resistance that getting the satellite to be completely still was incredibly hard. Thus the satellite had a small periodic
rotation when we started the control algorithm of about ±10° around a equilibrium point. When
we started the controller the satellite started to spin multiple rounds around the axis we were
testing.
In all of these tests the positive side for each satellite points upwards, e.g for the x-axis test the
positive x-axis points upwards. Since the input to the satellite corresponds to the attitude being
45° wrong about that axis it should start spinning in the clockwise direction, which it does.
Chapter 7
Discussion and Conclusion
7.1 Discussion
The simulations done in Matlab are mostly just rehashing of previous work, and everything
worked as expected, it reached the desired attitude. We did discover some errors with some
of the previous work that have been done, but this have been corrected here. A issue that have
come up is that the ADCS system uses a special type of pulse width modulation to control the
magnetic moment [14]. This means that what we are controlling is not the amplitude of the
magnetic moment, as we do in the simulations, but how long the current trough the coils are set
to the maximum value. How this will affect the stability of the satellite is unknown and should
be investigated in a future project.
We have also looked at target pointing, that is pointing the satellites camera in some direction
other than the nadir direction, but still pointing it at Earth. This have been found to be impossible, simply because the satellite have certain stable attitudes and target pointing would mean
leaving these stable attitudes. This can not be done since the torque we are generating with the
magnetorquers are of the same order as the gravity gradient torque. When we are doing nadir
pointing we are working with this torque, when we try to do target pointing we are working
against it, and we simply do not have enough power to do that.
Due to the technical limitations of the Helmholtz coils and the test setup in general, and the
37
CHAPTER 7. DISCUSSION AND CONCLUSION
38
fact that we have no properly working attitude determination system yet the tests that have
been done here is the best we can do so far. Hopefully when the attitude determination system
is fully implemented the controller can be tested more thoroughly. This should be done before
the satellite is launched. The test that have been done here have been velocity increase and decrease tests. Both tests have yielded positive results.
Even with a functioning attitude determination system testing the control algorithm properly
will probably be difficult, if not impossible. More thorough testing needs to be done, but we
have no way of simulating the free fall conditions the satellite will experience in orbit. Since
it can not generate enough torque to overcome gravity this means that we can not test control
around all axes simultaneously. There exists air-cushions and different other devices to overcome this, but none that we have access to here.
To better test the satellite we need a new bearing to suspend it in the magnetic field generated by the Helmholtz coils. The fishing line works to a certain extent, but since it behaves like
a torsion feather it also introduces other error sources. A better bearing with less friction and a
good attitude determination system could allow for much more thorough testing of the attitude
control system.
7.2 Conclusion
The simulations showed that the controller was able to steer the satellite to the correct attitude.
Even tough not every possible initial value have been tried, the simulations done here and previously points to the satellite being globally stable with the controller. It should also be noted
that this is only the case when we model the satellite with the gravity gradient torque as the only
external disturbance. More disturbances have been investigated in [2].
We have also found that it is impossible to use the controller to do target pointing, this follows
CHAPTER 7. DISCUSSION AND CONCLUSION
39
almost directly from the stability analysis and have also been shown in simulations.
The hardware tests performed here have all yielded positive results. The prototype satellite behaves as expected with the controller, it damps out angular velocity, and it manages to start a
rotation in the correct direction around at least one axis when the satellite is standing still.
Several simplifications have been done when testing the hardware, this is mainly since we do
not have a fully functioning attitude determination system and because the bearing used to
suspend the satellite is not ideal. From the tests done here we have increased confidence that
the controller will work as expected when the satellite is launched. It is however necessary with
more testing before we can be sure that the controller will be able to control the satellite to the
desired attitude.
7.3 Future Work
New simulations where the actual usage of the magnetorquers are being mimicked needs to be
done to see how this affects the performance compared to the ideal usage of the magnetorquers
that have been looked at so far. This also needs to include the switching between the attitude
determination system and the attitude control system. It could also be a advantage to include
some noise on the gyroscope data in the simulations to see how this affects the controller.
The control algorithm should be tested on hardware together with a fully functional attitude
determination system. A framework for how these interact, and how the ADCS interacts with
the OBC and the rest of the satellite also needs to be implemented and tested. Some other way
to suspend the satellite between the Helmholtz coils that can provide less friction should also
be found to be able to do more thorough tests on the whole system.
The linear quadratic regulator can also be investigated further and tested on hardware. Simulations have shown that the nonlinear controller that have been tested here are better than
LQR, but testing and comparing both on hardware could be beneficial.
CHAPTER 7. DISCUSSION AND CONCLUSION
40
More work needs to be done on tuning of the controller gains and the output gain. This should
be based on the new simulations and hardware testing. Since we can not simulate the vacuum
of space or the free fall condition the satellite will experience down here on Earth finding good
values for the output can be challenging, but it have to be done before launch.
Appendix A
Acronyms
NTNU Norwegian University of Science and Technology
NUTS NTNU Test Satellite
ANSAT Norwegian Student Satellite Program
MCU Microcontroller Unit
CMOS Complementary metal–oxide–semiconductor
NAROM Norwegian Center for Space-Related Education
ADCS Attitude Determination and Control System
LQR Linear Quadratic Regulator
ECI Earth-Centered Inertial
ECEF Earth-Centered Earth-Fixed
IR Infrared
MSS Marine Systems Simulator
41
Appendix B
Low Pass Filter for Gyroscope Measurements
The gyroscope measurements where quickly found to be noisy. This is a problem for the controller since the gain for the dampening part is d = 7 · 10−5 while the gain on the proportional
part is only p = 7 · 10−8 . This means that a small variation in the measurement of ωbib will result
in errors on the calculated control gain. This is mostly a problem when the satellite is not moving, since the gyroscope measurements then will oscillate around zero, leading to our output
doing the same.
This is a problem that should be fixed by the attitude estimator. When the attitude determination system is fully functional the controller will receive the estimated angular velocity from
the system together with the attitude (quaternion), instead of getting it directly from the gyroscope. But since the estimator that have been developed so far does not provide any estimates
of ωbib or ωbob we need some other way to smooth the gyroscope data, or find some other way
around the problem. A simple low pass filter have been tested in Matlab by recording a data
series from the gyroscope when the ADCS prototype was hanging still.
A simple low pass filter for our angular velocities have the following form:
ˆ i = αωi + (1 − α)ω
ˆ i −1
ω
42
(B.1)
APPENDIX B. LOW PASS FILTER FOR GYROSCOPE MEASUREMENTS
43
ˆ i is the current estimate, ωi is the current measurement, and ω
ˆ i −1 is the previous estiWhere ω
mate. α is a design variable. We tried with several different values of α. The original data and
some estimates can be seen in Figure B.1. The low pass filter was also tested on hardware, but it
Angular Velocity (x)
x−axis
x−axis estimate, a = 0.6
x−axis estimate, a = 0.2
Angular Velocity (rad/sec)
0.015
0.01
0.005
0
−0.005
−0.01
0
5
10
15
20
25
Time (seconds)
30
35
Angular Velocity (y)
Angular Velocity (rad/sec)
45
y−axis
y−axis estimate, a = 0.6
y−axis estimate, a = 0.2
0.1
0.05
0
−0.05
−0.1
0
5
10
15
20
25
Time (seconds)
30
35
40
45
z−axis
z−axis estimate, a = 0.6
z−axis estimate, a = 0.2
Angular Velocity (z)
0.03
Angular Velocity (rad/sec)
40
0.02
0.01
0
−0.01
−0.02
0
5
10
15
20
25
Time (seconds)
30
35
40
45
Figure B.1: Angular velocity data from the gyroscope and some estimates
did not smooth the signal enough to make it usable directly in the control algorithm. This can
also be seen from the figure, where we see that we only managed to halve the amplitude of the
noise. This is still enough noise to interfere with the control algorithm. Thus other ways around
the problem had to be found.
It was also found that the gyroscope had a built in low pass filter, there where some options
for tunning it, but none that gave us the necessary results. If the estimate we end up getting
from the attitude estimator also turns out to be noisy around zero we need to look at other ways
of tunning the controller.
Appendix C
Matlab Code
%**************************************************************************
% File main.m
% Main script. Contains the following sections:
%
% - Simulate (integrate) system
%
% Written by Zdenko Tudor, 2011. Edited by Marius Westgaard, 2014
%**************************************************************************
clear all; close all; clc;
addpath('myIGRF');
% Initialize storage variables (For testing and plotting purposes only)
global VAR;
% Create parameters struct
P = parameters();
% Load neccessary parameters
orbitPeriod = P.orbitPeriod;
% For defining simulation length
scalePlot = P.scalePlot;
% For x-axes scaling when plotting
w_o = P.w_o;
44
APPENDIX C. MATLAB CODE
45
ix_max = P.ix_max;
iy_max = P.iy_max;
iz_max = P.iz_max;
Nx = P.Nx;
Ny = P.Ny;
Nz = P.Nz;
Ax = P.Ax;
Ay = P.Ay;
Az = P.Az;
m_B_x_max = Nx*Ax*ix_max;
m_B_y_max = Ny*Ay*iy_max;
m_B_z_max = Nz*Az*iz_max;
% Simulation parameters
numOfOrbits = 3;
% Simulation time in number of orbits
tSpan = [0,orbitPeriod*numOfOrbits];% Time span for ODE45
% Initial parameters
w_B_OB = [0.05; -0.02; 0.045];
% Initial rotational velocity
eulAng = [80;-50;130];
% Initial orientation
initLat = 0;
% Initial latitude
initJoule = 0;
% Always to be set to zero
% Transform initial parameters
qua = euler2q(eulAng(1)*pi/180, eulAng(2)*pi/180, eulAng(3)*pi/180);
S_eps = [ 0, -qua(4), qua(3);
qua(4), 0, -qua(2);
-qua(3), qua(2), 0];
R_O_B = eye(3) + 2*qua(1)*S_eps + 2*S_eps^2;
R_B_O = R_O_B';
w_B_IB = w_B_OB + R_B_O*[0; -w_o; 0];
% State vector x = [wx wy wz eta eps1 eps2 eps3 latitude totalJoule];
xInit = [w_B_IB;qua;initLat;initJoule];
APPENDIX C. MATLAB CODE
46
% ODE45 simulaton settings
options = odeset('MaxStep',1,'OutputFcn',@outFcn,'Refine',1);
%options = odeset('MaxStep',0.1,'OutputFcn',@outFcn,'RelTol',1e-5,'AbsTol',1e-5,'Refine',1)
% Simulate system
[tout,yout] = ode45(@(t,x)nonlinearSatellite(t,x,P),tSpan,xInit,options);
%**************************************************************************
% File outFcn.m
%
% status = outFcn(t,~,flag)
%
% Is run when ODE45 successfully completes a step.
% This function is used implicitly by ODE45, when ODE settings are
% properly set.
%
% Example:
% options = odeset('OutputFcn',@outFcn,'Refine',1);
% Refine sets how many steps should be created for each successfull step.
% With this set to 1 length of all vectors created by the OutputFcn will be
% the same as the length of ODE timevector TOUT.
%
% Written by Zdenko Tudor, 2011
%**************************************************************************
function status = outFcn(t,~,flag)
persistent ite;
global tmpVAR VAR;
switch flag
case 'init'
%initialize arrays etc.
ite = 1;
VAR.moment = zeros(10000,3);
VAR.torque = zeros(10000,3);
VAR.J = zeros(10000,1);
APPENDIX C. MATLAB CODE
VAR.t = zeros(10000,1);
status = 0;
case 'done'
VAR.moment(ite+1:end,:) = [];
VAR.torque(ite+1:end,:) = [];
VAR.J(ite+1:end,:) = [];
VAR.t(ite+1:end,:) = [];
status = 0;
otherwise
ite = ite+1;
VAR.moment(ite,:)=tmpVAR.moment';
VAR.torque(ite,:)=tmpVAR.torque';
VAR.W(ite,:)=tmpVAR.W;
VAR.t(ite,1)=t(end);
status = 0;
end
end
%**************************************************************************
% File parameters.m
%
% P = parameters()
%
% Creates a parameter struct P, which is used to pass the parameters
% between different functions.
%
% Written by Zdenko Tudor, 2011. Edited by Marius Westgaard, 2014
%**************************************************************************
function P = parameters()
% Coil data
% Battery voltage
47
APPENDIX C. MATLAB CODE
48
V = 3.3;
% Number of coil windings (turns)
Nx = 155;
Ny = 155;
Nz = 238;
% [m^2] Coil area
Ax = 0.08*0.18;
Ay = 0.08*0.18;
Az = 0.08^2;
% [Ohm] Coil resistance
Rx = 27.2;
Ry = 27.2;
Rz = 18.5;
% [A] Maximum current through coils
ix_max = V/Rx;
iy_max = V/Ry;
iz_max = V/Rz;
% Satellite dimensions
x = 0.1;
% X-axis length
y = 0.1;
% Y-axis length
z = 0.2;
% Z-axis length
% Other constants
m = 2.6;
% [kg] Satellite mass
Re = 6371.2e3;
% [m] Earth radius
Rs = 600e3;
% [m] Satellite altitude
Rc = Re+Rs;
% [m] Distance from earth center to satellite
G = 6.67428e-11;
% Earth gravitational constant
M = 5.972e24;
% Earth mass
APPENDIX C. MATLAB CODE
49
% Calculated constants
Ix = (m*0.9/12)*(y^2+z^2);
% X-axis inertia
Iy = (m/12)*(x^2+z^2);
% Y-axis inertia
Iz = (m/12)*(x^2+y^2);
% Z-axis inertia
I = diag([Ix Iy Iz]);
% Inertia matrix
w_o = sqrt(G*M/Rc^3);
% Satellite angular velocity relative Earth
orbitPeriod = (2*pi)/(w_o); % For defining simulation length
scalePlot = 1/orbitPeriod;
% For x-axes scaling when plotting
% Aerodynamic drag constants
jx = 0.03;
% X-axis length for centre of pressure
jy = 0.03;
% Y-axis length for centre of pressure
jz = 0.083;
% Z-axis length for centre of pressure
Jx = (m/12)*(jy^2+jz^2);
% X-axis inertia for centre of pressure
Jy = (m/12)*(jx^2+jz^2);
% Y-axis inertia for centre of pressure
Jz = (m/12)*(jx^2+jy^2);
% Z-axis inertia for centre of pressure
J = diag([Jx Jy Jz]);
% Inertia matrix for centre of pressure
A_drag = z*sqrt(x^2 + y^2); % Maximum area for calc disturbance drags
vel = 2*pi*Rc/orbitPeriod;
% Velocity
rho_a = 4.89e-13;
% Atmospheric density
c_p = [0.02 0.02 0.02]';
% Relative location of centre of pressure
c_p_x = [
0
-c_p(3)
c_p(3)
-c_p(2)
0
c_p(2)
-c_p(1)
c_p(1)
0 ];
% Store variables
P.w_o = w_o;
P.V = V; P.I = I; P.J = J;
P.Nx = Nx; P.Ny = Ny; P.Nz = Nz;
APPENDIX C. MATLAB CODE
P.Ax = Ax; P.Ay = Ay; P.Az = Az;
P.Rx = Rx; P.Ry = Ry; P.Rz = Rz;
P.ix_max = ix_max; P.iy_max = iy_max; P.iz_max = iz_max;
P.Re = Re; P.Rc = Rc;
P.orbitPeriod = orbitPeriod; P.scalePlot = scalePlot;
P.A_drag = A_drag; P.vel = vel;
P.rho_a = rho_a; P.c_p_x = c_p_x;
end
%**************************************************************************
% File nonlinearSatellite.m
%
% xDot = nonlinearSatellite(T,X,P)
%
% Simulates the earth-satellite system. To be used with an ordinary
% differential equation solver such as MATLAB's Runge Kutta(4,5) method,
% ODE45.
% The input:
% T - time vector (implicit)
% X - state vector, initial state has to be provided. The state vector is
% of the form X = [W^B_IB' Q' LAMBDA JOULE]. W^B_IB is a 3x1 vector
% containing the three rotational velocity components. Q is the quaternion
% vector of the form Q = [eta eps1 eps2 eps3]. LAMBDA is the latitude of
% the satellite and JOULE is only a testing paramemter which will return
% the total Joule consumption of the satellite during the simulation.
% Initial JOULE should always be set to zero.
% P - struct containg necessary parameters (w_o, I, Rc, coil data, voltage
% and current ratings).
%
% Example:
% [TOUT,YOUT] = ode45(@(t,x)nonlinearSatellite(t,x,P),tSpan,xInit,options)
% this produces:
% TOUT - time vector for the integration
% YOUT - Integrated xDot vector for each timestep size(YOUT)=length(TOUT)x9
%
50
APPENDIX C. MATLAB CODE
51
% Written by Zdenko Tudor, 2011. Edited by Marius Westgaard, 2014
%**************************************************************************
function xDot = nonlinearSatellite(t,x,P)
global R_B_O;
global tmpVAR;
% Load necessary parameters from P struct.
w_o = P.w_o;
I = P.I;
J = P.J;
rho_a = P.rho_a;
A_drag = P.A_drag;
c_p_x = P.c_p_x;
vel = P.vel;
% Normalization of quartenions
x(4:7) = x(4:7)./norm(x(4:7));
% State vector x = [];
w_B_IB = x(1:3);
% Angular velocity vector
eta = x(4);
% Euler parameter eta
eps = x(5:7);
% Euler parameter epsilon
lat = x(8);
% Lambda
q = x(4:7);
% Epsilon crossmatrix, skew-symmetric
S_eps = [ 0, -eps(3), eps(2);
eps(3), 0, -eps(1);
-eps(2), eps(1), 0];
% Rotation matrices
R_O_B = eye(3) + 2*eta*S_eps + 2*S_eps^2;
R_B_O = R_O_B';
% Angular velocities of frames relative each other
APPENDIX C. MATLAB CODE
52
w_O_IO = [0;-w_o;0];
w_B_OB = w_B_IB-R_B_O*w_O_IO;
% Geomagnetic field in Orbit frame and in Body frame
[B_O] = IGRF(10,10,lat,0,P);
% IGRF model
B_B = R_B_O*B_O;
% Chose desired attitude
desiredEulAng = [30;30;0]*pi/180;
q_desired = euler2q(desiredEulAng(1), desiredEulAng(2), desiredEulAng(3));
% Choose controller
tic
[tau_m, m_B, W] = pd_nadir_controller(P, w_B_IB, B_B, q);
%[tau_m, m_B, W] = pd_desired_controller(P, w_B_IB, B_B, q, q_desired);
%[tau_m, m_B, W] = lq_nadir_controller(P, w_B_IB, B_B, q);
toc
%-----Storing variables used for testing purposes only-----Start----tmpVAR.torque = tau_m;
tmpVAR.moment=m_B;
tmpVAR.W = W;
xDot(9,1) = W;
%-----Storing variables used for testing purposes only------End------
% Gravitational torque
c3 = R_B_O(:,3);
tau_g = 3*w_o^2*cross(c3,I*c3);
% Aerodynamic drag torque
Vr = R_B_O*[-1 0 0]';
Vr_x = [
Vr(3)
-Vr(2)
0
-Vr(3)
0
Vr(1)
Vr(2)
-Vr(1)
0 ];
tau_a = rho_a*vel*(vel*A_drag*c_p_x*Vr - (I + Vr_x*J)*w_B_OB);
APPENDIX C. MATLAB CODE
% Total disturbance torques
tau_dist = tau_g; % + tau_a;
% Angular acceleration
wDot_B_IB=I\(tau_m+tau_dist-cross(w_B_IB,I*w_B_IB));
% Eta
etaDot = -0.5*eps'*w_B_OB;
% Epsilon
epsDot = 0.5*(eta*eye(3)+S_eps)*w_B_OB;
% Variable assignment
xDot(1:3,1) = wDot_B_IB;
xDot(4,1) = etaDot;
%qDot(1);
xDot(5:7,1) = epsDot;
%q(2:4);
xDot(8,1) = w_o;
end
%**************************************************************************
% File pd_nadir_controller.m
%
% [tau_m, m_B, W] = pd_nadir_controller(P, w_B_IB, B_B, q)
%
% Calculates the torque, magetic moment and power usage generated
% by the controller. To be used with nonlinearSatellite.m
%
% Written by Marius Westgaard, 2014
%**************************************************************************
function [tau_m, m_B, W] = pd_nadir_controller(P, w_B_IB, B_B, q)
% Epsilon crossmatrix, skew-symmetric
S_eps = [ 0, -q(4), q(3);
q(4), 0, -q(2);
-q(3), q(2), 0];
53
APPENDIX C. MATLAB CODE
% Rotation matrices
R_O_B = eye(3) + 2*q(1)*S_eps + 2*S_eps^2;
R_B_O = R_O_B';
% Angular velocities of frames relative each other
w_o = P.w_o;
w_O_IO = [0;-w_o;0];
w_B_OB = w_B_IB-R_B_O*w_O_IO;
% Controller gains
p = 7e-8;
d = 7e-5;
% Moment set up by coils before scaling (saturating current)
m_B = (1/norm(B_B,2)^2)*(-p*cross(B_B,q(2:4))-d*cross(B_B,w_B_OB));
% Moment set up by coils after scaling
[m_B, W] = currentScaling(P,m_B);
% Torque set up by coils
tau_m = cross(m_B,B_B);
end
%**************************************************************************
% File pd_desired_controller.m
%
% [tau_m, m_B, W] = pd_desired_controller(P, w_B_IB, B_B, q, q_desired)
%
% Calculates the torque, magetic moment and power usage generated
% by the controller. To be used with nonlinearSatellite.m
%
% Written by Marius Westgaard, 2014
%**************************************************************************
54
APPENDIX C. MATLAB CODE
55
function [tau_m, m_B, W] = pd_desired_controller(P, w_B_IB, B_B, q, q_desired)
% Quaternion difference
q_diff = quatmultiply(q', quatconj(q_desired'))';
% Epsilon crossmatrix, skew-symmetric
S_eps = [ 0, -q(4), q(3);
q(4), 0, -q(2);
-q(3), q(2), 0];
% Rotation matrices
R_O_B = eye(3) + 2*q(1)*S_eps + 2*S_eps^2;
R_B_O = R_O_B';
% Angular velocities of frames relative each other
w_o = P.w_o;
w_O_IO = [0;-w_o;0];
w_B_OB = w_B_IB-R_B_O*w_O_IO;
% Controller gains
p = 7e-7;
d = 7e-4;
% Moment set up by coils before scaling (saturating current)
m_B = (1/norm(B_B,2)^2)*(-p*cross(B_B,q_diff(2:4))-d*cross(B_B,w_B_OB));
% Moment set up by coils after scaling
[m_B, W] = currentScaling(P,m_B);
% Torque set up by coils
tau_m = cross(m_B,B_B);
end
%**************************************************************************
% File currentScaling.m
%
APPENDIX C. MATLAB CODE
% [m_B, W] = currentScaling(P,m_B)
%
% Scales the magnetic moment m_B so that it does not use a higher current
% on each coil than what is available for the real system
%
% Written by Zdenko Tudor, 2011. Edited by Marius Westgaard, 2014
%**************************************************************************
function [m_B, W] = currentScaling(P,m_B)
% Load coil parameters from parameter struct
Nx = P.Nx; Ny = P.Ny; Nz = P.Nz;
Ax = P.Ax; Ay = P.Ay; Az = P.Az;
ix_max = P.ix_max; iy_max = P.iy_max; iz_max = P.iz_max;
V = P.V;
% Current scaling
% Creates ratio variable containing a coil's: max current/wanted current.
% The lowest value in ratio variable, if below 1 is the highest current
% violation, thus all the currents should be scaled by this factor.
ix = m_B(1)/(Nx*Ax); ratio(1)=ix_max/abs(ix);
iy = m_B(2)/(Ny*Ay); ratio(2)=iy_max/abs(iy);
iz = m_B(3)/(Nz*Az); ratio(3)=iz_max/abs(iz);
if min(ratio)<1
m_B = m_B*min(ratio);
end
% Power consumption (Joules)
Wx = abs(V*m_B(1)/(Nx*Ax));
Wy = abs(V*m_B(2)/(Ny*Ay));
Wz = abs(V*m_B(3)/(Nz*Az));
W = Wx+Wy+Wz;
end
56
Appendix D
C Code
/*
* adcs.c
*
* Created: 18.10.2013 14:16:37
*
Author: Oyvind Rein
* Edited: 15.11.2014
* By: Marius Westgaard and Henrik Rudi Haave
*/
#define F_CPU (32000000UL)
#include <avr/io.h>
#define MODE_IDLE 0
#define MODE_DETUMBLING 1
#define MODE_TUMBLING 2
#define MODE_ESTIMATE 3
#define MODE_CONTROL 4
#define RX_BUFFER_SIZE 200
#define PI 3.1415926535F
#include "coil_driver/coil_driver.h"
#include "init/clk_init.h"
57
APPENDIX D. C CODE
58
#include "bluetooth/bluetooth.h"
#include "sun_sensor/sun_sensor.h"
#include "magnetometer/magnetometer.h"
#include "gyroscope/gyroscope.h"
#include <util/delay.h>
#include <avr/interrupt.h>
#include "b_dot/b_dot.h"
#include "tumbling/tumbling.h"
#include "estimator/EKF.h"
#include "estimator/sun_vect.h"
#include "estimator/normlize.h"
#include "estimator/rtwtypes.h"
#include "pd_controller/pd_controller.h"
int mode = MODE_IDLE;
volatile int new_cycle;
volatile int MODE = MODE_IDLE;
char rx_buffer[RX_BUFFER_SIZE]; // Must be less that 256 since rx_buffer_length is uint8_t
uint8_t rx_message_length = 0;
uint8_t message_received = 0;
int16_t sensor_data[18];
int16_t output_vector[3];
// Testing Variables
float q[4];
int16_t w[3];
int compareStrings(char *a, char *b);
int init();
int main(void)
{
init();
APPENDIX D. C CODE
59
while(1){
while(!new_cycle);
new_cycle = 0;
poll(sensor_data);
add_to_log(sensor_data);
switch(mode){
case MODE_IDLE:
output_vector[0] = 0;
output_vector[1] = 0;
output_vector[2] = 0;
coil_output(output_vector);
break;
case MODE_DETUMBLING:
b_dot_update(&sensor_data[0], output_vector);
coil_output(output_vector);
break;
case MODE_TUMBLING:
b_dot_update(&sensor_data[0], output_vector);
output_vector[0] = (-1) * output_vector[0];
output_vector[1] = (-1) * output_vector[1];
output_vector[2] = (-1) * output_vector[2];
coil_output(output_vector);
break;
case MODE_ESTIMATE:
PORTB_OUTSET = 0x01;
// The coordinate system used is based on what was defined by rein for
// EKF_U is input
// EKF_Y is output
// memcpy(&sun_vect_step_U.In1[0], &output_vector[6], 6 * (sizeof(int16
for (uint8_t i = 0; i<6 ; i++){
sun_vect_U.In1[i] = sensor_data[i+6];
}
for (uint8_t i = 0; i<3 ; i++){
normlize_U.in[i] = sensor_data[i];
//EKF_U.z_n[i+3] = [m_x m_y m_z
APPENDIX D. C CODE
60
}
normlize_step();
sun_vect_step();// this builds and normalizes a sun vector [s_x s_y s_z
for (uint8_t i = 0; i<3 ; i++){
EKF_U.W[i] = ((real32_T)sensor_data[i+3]*0.00875F*PI)/(180.0F) ;
EKF_U.z_n[i] = sun_vect_Y.Out1[i];
//sun_vect_Y.Out1[] = [s_x s_y
EKF_U.z_n[i+3] = normlize_Y.Out[i];
}
EKF_U.r_m[0]=0.99F;//-0.36F;
EKF_U.r_m[1]=-0.12F;//-0.01F;
EKF_U.r_m[2]=0.11F;//0.93F;
EKF_U.r_s[0]= 0.0001F;
EKF_U.r_s[1]=-1.0F;
EKF_U.r_s[2]= 0.0001F;
EKF_U.h = 0.250F;
static const float Q[3] = {0.04F, 0.0000000F, 0.0000000F}; // i am cons
static const float R[2] = {0.002F, 0.002F,}; // [sun mag]
static const float q_0[4] = {1.0F , 0.0F, 0.0F, 0.0F};
memcpy(EKF_U.Q,Q,sizeof(Q));
memcpy(EKF_U.R,R,sizeof(R));
memcpy(EKF_U.q_0,q_0,sizeof(q_0));
EKF_step(); // estimation time is unknown, converting this to a fixed p
PORTB_OUTCLR = 0x01;
//EKF_Y.x[] // at the moment the plan is to memcopy the floating point
pd_controller_update(&sensor_data[0], &sensor_data[3], &EKF_Y.x[0], out
coil_output(output_vector);
break;
APPENDIX D. C CODE
61
case MODE_CONTROL:
q[0] = 0.9239;
q[1] = 0.0;
q[2] = 0.0;
q[3] = 0.3827;
w[0] = 0;
w[1] = 0;
w[2] = 0;
pd_controller_update(&sensor_data[0], &w, &q, output_vector);
coil_output(output_vector);
break;
default:
break;
}
read_message();
while(!coils_are_off());
_delay_ms(1);
if (new_cycle == 1){
printf("Timer is too fast!\n");
}
}
}
int init(){
ckl_init_32MHz();
_delay_ms(1000);
coil_driver_init();
bluetooth_init();
sun_sensor_init();
magnetometer_init();
gyro_init();
start_timer();
sei();
APPENDIX D. C CODE
62
printf("Starting\n");
b_dot_init();
pd_controller_init();
///////////Setting a port to find timing///////////
PORTB_DIRSET = 0x01;
///////////-----------------------------///////////
return 1;
}
int poll(int16_t sensor_data[18]){
magnetometer_poll(&sensor_data[0]);
gyro_poll(&sensor_data[3]);
sun_sensor_poll(&sensor_data[6]);
}
void start_timer(){
TCC1.CTRLA = TC_CLKSEL_DIV256_gc;
TCC1.INTCTRLA = TC_OVFINTLVL_LO_gc;
TCC1.PER = 500000; // PER = 12500, fCPU = 32M and DIV = 256 gives 10hz. 1/0.25s = 4hz (
PMIC.CTRL = PMIC_LOLVLEN_bm;
TCC1.INTFLAGS = 0x01;
}
int read_message(){
if (!message_received){
return 1;
}
if (compareStrings(rx_buffer, "i")){
mode = MODE_IDLE;
}
else if (compareStrings(rx_buffer, "d")){
mode = MODE_DETUMBLING;
}
else if (compareStrings(rx_buffer, "get_log")){
get_log();
APPENDIX D. C CODE
63
}
else if (compareStrings(rx_buffer, "clear_log")){
clear_log();
}
else if (compareStrings(rx_buffer, "get_data")){
for(int i = 0; i < 12; i++){
printf("%d,", sensor_data[i]);
}
int16_t temp_coil = output_vector[0];
output_vector[0] = output_vector[1];
output_vector[1] = temp_coil;
output_vector[0] *= -1;
output_vector[1] *= -1;
output_vector[2] *= -1;
printf("%d,%d,%d,", output_vector[0], output_vector[1], output_vector[2]);
printf("%ld,%ld,%ld,%ld,",EKF_Y.x[0],EKF_Y.x[1],EKF_Y.x[2],EKF_Y.x[3]);
printf("%ld,%ld,%ld,",EKF_Y.x[4],EKF_Y.x[5],EKF_Y.x[6]);
printf("%ld,%ld,%ld",EKF_Y.x[7],EKF_Y.x[8],EKF_Y.x[9]);
//printf("%ld,%ld,%ld,",EKF_Y.z_hat[3],EKF_Y.z_hat[4],EKF_Y.z_hat[5]);
//printf("%ld,%ld,%ld",EKF_Y.z_hat[0],EKF_Y.z_hat[1],EKF_Y.z_hat[2]);
printf("\n");
}
else if (compareStrings(rx_buffer, "get_EKF")){
printf("EKF %f, %f, %f, %f\n", EKF_Y.x[0], EKF_Y.x[1], EKF_Y.x[2], EKF_Y.x[3]);
printf("sun_vect %f, %f, %f\n", sun_vect_Y.Out1[0], sun_vect_Y.Out1[1], sun_vect_Y.
printf("mag_vect %f, %f, %f\n", normlize_Y.Out[0], normlize_Y.Out[1], normlize_Y.Ou
printf("gyro vect %f %f %f \n",EKF_U.W[0], EKF_U.W[1], EKF_U.W[2]);
printf("Q %f %f %f and R %f %f\n", EKF_U.Q[0], EKF_U.Q[1], EKF_U.Q[2], EKF_U.R[0],
}
else if (compareStrings(rx_buffer, "t")){
mode = MODE_TUMBLING;
}
else if (compareStrings(rx_buffer, "e")){
mode = MODE_ESTIMATE;
}
else if (compareStrings(rx_buffer, "c")){
APPENDIX D. C CODE
mode = MODE_CONTROL;
}
message_received = 0;
rx_message_length = 0;
return 1;
}
int compareStrings(char *a, char *b){
int x = 0;
while(1){
if(a[x] != b[x]){
return 0;
}else if(a[x] == '\0'){
return 1;
}
x++;
}
}
ISR(USARTF0_RXC_vect){
if (rx_message_length < RX_BUFFER_SIZE && !(message_received)){
char c =
USARTF0_DATA;
rx_buffer[rx_message_length] = c;
if (c == '\r' || c == '\n'){
message_received = 1;
rx_buffer[rx_message_length] = '\0';
}else{
rx_message_length++;
}
}
}
ISR(TCC1_OVF_vect){
new_cycle = 1;
}
64
APPENDIX D. C CODE
65
/*
* pd_controller.c
*
* Created: 15.10.2014 19:17:22
*
Author: Marius Westgaard
*/
#include <avr/io.h>
#include <math.h>
static const float Re = 6371200.0F;
// Re is the earth radius
static const float Rs = 0.0F;
// Rs is the satellite orbit height, set to 0 for t
static const float gyro_scaling = (0.00875/180.0)*3.1415926535; // Found in datasheet
static const float magnetometer_scaling = 1.0/230.0;
volatile float m_B[3];
volatile float w_B_OB[3];
volatile float w_o;
volatile float B_B_filt[3];
volatile float w_B_IB_filt[3];
float omega_previous[3];
float B_previous[3];
// Controller constants
#define p 0.00000007
// Gain found by simulations
#define d 0.00007
// Gain found by simulations
#define gain_x 121000100000.0F
// Gain found by experimenting
#define gain_y 121000100000.0F
// Gain found by experimenting
#define gain_z 121000100000.0F
// Gain found by experimenting
#define COIL_MAX 32767.0F
// LP filter constants
#define a_B 0.8
#define b_B 0.2
#define a_omega 0.3
#define b_omega 0.7
// Found in datasheet
APPENDIX D. C CODE
66
int pd_controller_init(){
for(int i = 0; i < 3; i++){
m_B[i] = 0;
w_B_OB[i] = 0;
B_B_filt[i] = 0;
w_B_IB_filt[i] = 0;
omega_previous[i] = 0;
B_previous[i] = 0;
}
w_o = sqrt((6.67428*pow(10,-11)*5.972*pow(10,24))/pow((Re+Rs),3));
w_o = 0;
// Orbit speed
// For testing
return 1;
}
int pd_controller_update(int16_t B_B[3], int16_t w_B_IB[3], float q[4], int16_t coil_output
b_filter(B_B, B_B_filt);
omega_filter(w_B_IB, w_B_IB_filt);
w_B_OB[0] = w_B_IB_filt[0]+(2*q[0]*q[3]+2*q[1]*q[2])*w_o;
w_B_OB[1] = w_B_IB_filt[1]+(1-2*q[1]*q[1]-2*q[3]*q[3])*w_o;
w_B_OB[2] = w_B_IB_filt[2]+(-2*q[0]*q[1]+2*q[2]*q[3])*w_o;
m_B[0] = (-p*(B_B_filt[1]*q[3] - B_B_filt[2]*q[2]) -d*(B_B_filt[1]*w_B_OB[2] - B_B_filt
m_B[1] = (-p*(B_B_filt[2]*q[1] - B_B_filt[0]*q[3]) -d*(B_B_filt[2]*w_B_OB[0] - B_B_filt
m_B[2] = (-p*(B_B_filt[0]*q[2] - B_B_filt[1]*q[1]) -d*(B_B_filt[0]*w_B_OB[1] - B_B_filt
m_B[0] = m_B[0]/(B_B_filt[0]*B_B_filt[0] + B_B_filt[1]*B_B_filt[1] + B_B_filt[2]*B_B_fi
m_B[1] = m_B[1]/(B_B_filt[0]*B_B_filt[0] + B_B_filt[1]*B_B_filt[1] + B_B_filt[2]*B_B_fi
m_B[2] = m_B[2]/(B_B_filt[0]*B_B_filt[0] + B_B_filt[1]*B_B_filt[1] + B_B_filt[2]*B_B_fi
// Make sure there are no overflow for the x gain
if (m_B[0] > (COIL_MAX/gain_x)) {
coil_output[0] = 32767;
}
else if (m_B[0] < -(COIL_MAX/gain_x)) {
APPENDIX D. C CODE
coil_output[0] = -32767;
}
else {
coil_output[0] = gain_x*m_B[0];
}
// Make sure there are no overflow for the y gain
if (m_B[1] > (COIL_MAX/gain_y)) {
coil_output[1] = 32767;
}
else if (m_B[1] < -(COIL_MAX/gain_y)) {
coil_output[1] = -32767;
}
else {
coil_output[1] = gain_y*m_B[1];
}
// Make sure there are no overflow for the z gain
if (m_B[2] > (COIL_MAX/gain_x)) {
coil_output[2] = 32767;
}
else if (m_B[2] < -(COIL_MAX/gain_z)) {
coil_output[2] = -32767;
}
else {
coil_output[2] = gain_z*m_B[2];
}
return 1;
}
void b_filter(int16_t B[3], float B_filt[3]) {
B_filt[0] = magnetometer_scaling*B[0];
B_filt[1] = magnetometer_scaling*B[1];
B_filt[2] = magnetometer_scaling*B[2];
67
APPENDIX D. C CODE
68
// Have been tested with and without the LP-filter
/*
B_filt[0] = a_B*magnetometer_scaling*B[0] + b_B*B_previous[0];
B_filt[1] = a_B*magnetometer_scaling*B[1] + b_B*B_previous[1];
B_filt[2] = a_B*magnetometer_scaling*B[2] + b_B*B_previous[2];
B_previous[0] = B_filt[0];
B_previous[1] = B_filt[1];
B_previous[2] = B_filt[2];
*/
}
void omega_filter(int16_t omega[3], float omega_filt[3]) {
omega_filt[0] = gyro_scaling*omega[0];
omega_filt[1] = gyro_scaling*omega[1];
omega_filt[2] = gyro_scaling*omega[2];
// Have been tested with and without the LP-filter
/*
omega_filt[0] = a_omega*gyro_scaling*omega[0] + b_omega*omega_previous[0];
omega_filt[1] = a_omega*gyro_scaling*omega[1] + b_omega*omega_previous[1];
omega_filt[2] = a_omega*gyro_scaling*omega[2] + b_omega*omega_previous[2];
omega_previous[0] = omega_filt[0];
omega_previous[1] = omega_filt[1];
omega_previous[2] = omega_filt[2];
*/
}
/*
* coil_driver.c
*
* Created: 28.02.2014 11:35:11
*
Author: Oyvind Rein
* Edited: 15.11.2014
*
*/
By: Marius Westgaard
APPENDIX D. C CODE
69
#include "coil_driver.h"
#include <avr/io.h>
#include <avr/interrupt.h>
#define COIL_MAX 32767
#define COIL_PER 22000 //11000
int coils_off;
int coil_driver_init(){
PORTC_DIR |= 0x07; //
Set pin CXSLEEP, CYSLEEP and CZSLEEP to output.
PORTD_DIR |= 0x0F; // Set pin CX_PHASE, CX_PWM, CY_PHASE and CY_PWM to output.
PORTE_DIR |= 0x0C; // Set pin CZ_PHASE and CZ_PWM to output
PORTE_DIR &= ~0x13; // Set pin CX_BR_FAULT, CY_BR_FAULT and CZ_BR_FAULT to input.
set_coil(COIL_X, COIL_OFF);
set_coil(COIL_Y, COIL_OFF);
set_coil(COIL_Z, COIL_OFF);
coils_off = 1;
PORTC_OUTSET |= 0x07; // Set pin CXSLEEP, CYSLEEP and CZSLEEP to high. Disables sleep m
return 1;
}
void coil_output(int16_t output_vector[3]){ // Values from -2048 to 2047. 12 bit accuracy
// Change it so that the axes are correct
int16_t temp_coil = output_vector[0];
output_vector[0] = output_vector[1];
output_vector[1] = temp_coil;
output_vector[0] *= -1;
output_vector[1] *= -1;
output_vector[2] *= -1;
coils_off = 0;
int16_t coil_x_output = output_vector[0];
int16_t coil_y_output = output_vector[1];
int16_t coil_z_output = output_vector[2];
if (coil_x_output < 0){
APPENDIX D. C CODE
70
coil_x_output = -1*coil_x_output;
set_coil(COIL_X, COIL_NEG);
}else{
set_coil(COIL_X, COIL_POS);
}
if (coil_y_output < 0){
coil_y_output = -1*coil_y_output;
set_coil(COIL_Y, COIL_NEG);
}else{
set_coil(COIL_Y, COIL_POS);
}
if (coil_z_output < 0){
coil_z_output = -1*coil_z_output;
set_coil(COIL_Z, COIL_NEG);
}else{
set_coil(COIL_Z, COIL_POS);
}
if (coil_x_output > COIL_MAX) {
coil_x_output = COIL_MAX;
}
if (coil_y_output > COIL_MAX) {
coil_y_output = COIL_MAX;
}
if (coil_z_output > COIL_MAX) {
coil_z_output = COIL_MAX;
}
TCC0.CTRLA = TC_CLKSEL_DIV256_gc;
TCC0.PER = COIL_PER; // PER = 500, fCPU = 32M and DIV = 64 gives 1000hz sampling
PMIC.CTRL = PMIC_LOLVLEN_bm;
TCC0.INTFLAGS = 0b01110001;
TCC0.CCA = (COIL_PER*(float)coil_x_output)/COIL_MAX;
TCC0.CCB = (COIL_PER*(float)coil_y_output)/COIL_MAX;
TCC0.CCC = (COIL_PER*(float)coil_z_output)/COIL_MAX;
TCC0.CTRLFSET = TC_CMD_RESTART_gc; // Restart timer.
TCC0.INTCTRLB = 0b00010101; // Set interrupt level to low for CCA, CCB and CCC
TCC0.INTCTRLA = TC_OVFINTLVL_LO_gc;
APPENDIX D. C CODE
}
ISR(TCC0_CCA_vect){
set_coil(COIL_X, COIL_OFF);
}
ISR(TCC0_CCB_vect){
set_coil(COIL_Y, COIL_OFF);
}
ISR(TCC0_CCC_vect){
set_coil(COIL_Z, COIL_OFF);
}
ISR(TCC0_OVF_vect){
set_coil(COIL_X, COIL_OFF);
set_coil(COIL_Y, COIL_OFF);
set_coil(COIL_Z, COIL_OFF);
TCC0.INTCTRLA = 0x00; // Turn off interrupts
TCC0.INTCTRLB = 0x00;
coils_off = 1;
}
int coils_are_off(){
return coils_off;
}
void set_coil(int coil, int dir){
if (coil == COIL_X){
if (dir == COIL_OFF){
PORTD_OUT &= ~CX_PWM;
}
else if (dir == COIL_POS){
PORTD_OUT &= ~CX_PHASE;
PORTD_OUTSET |= CX_PWM;
71
APPENDIX D. C CODE
}
else if (dir == COIL_NEG){
PORTD_OUTSET |= CX_PHASE;
PORTD_OUTSET |= CX_PWM;
}
}
else if (coil == COIL_Y){
if (dir == COIL_OFF){
PORTD_OUT &= ~CY_PWM;
}
else if (dir == COIL_POS){
PORTD_OUT &= ~CY_PHASE;
PORTD_OUTSET |= CY_PWM;
}
else if (dir == COIL_NEG){
PORTD_OUTSET |= CY_PHASE;
PORTD_OUTSET |= CY_PWM;
}
}
else if (coil == COIL_Z){
if (dir == COIL_OFF){
PORTE_OUT &= ~CZ_PWM;
}
else if (dir == COIL_POS){
PORTE_OUT &= ~CZ_PHASE;
PORTE_OUTSET |= CZ_PWM;
}
else if (dir == COIL_NEG){
PORTE_OUTSET |= CZ_PHASE;
PORTE_OUTSET |= CZ_PWM;
}
}
return;
}
72
Bibliography
[1] Alvenes, F. (2013). Attitude Controller-Observer Design for the NTNU Test Satellite. Master
thesis, NTNU.
[2] Bråthen, G. (2013). Design of Attitude Control System of a Double CubeSat. Master thesis,
NTNU.
[3] Fossen, T. I. (2011). Handbook of Marine Craft Hydrodynamics and Motion Control. Wiley.
[4] Gravdahl, J. T. and Egeland, O. (2002). Modeling and simulation for automatic control. Marine Cybernetics.
[5] Haave, H. R. (2014). Implementation of Attitude Estimation and a Look at the UWE CubeSat.
Project thesis, NTNU.
[6] Hughes, P. C. (1986). Spacecraft Attitude Dynamics. Wiley.
[7] Jenssen, K. L. r. and Yabar, K. H. (2011). Development , Implementation and Testing of Two
Attitude Estimation Methods for Cube Satellites. Master thesis, NTNU.
[8] Josephsen, S., Bergslid, T. S., Gulbrandsen, F., Hetland, A., and Sogge, N. (2012). Eksperter i
Team NUTS - studentsatellitt Modell for demonstrasjon av ADCS. Technical report, NTNU.
[9] Khalil, H. K. (2002). Nonlinear Systems, Third Edition. Prentice Hall.
[10] MSS. http://www.marinecontrol.org/.
[11] NUTS. http://nuts.cubesat.no/.
73
BIBLIOGRAPHY
74
[12] Oland, E., Aas, A., Steihaug, T. M., Mathisen, S. V., and Vedal, F. (2009). A Design Guide for
Attitude Determination and Control Systems for Picosatellites. Technical report, HiN.
[13] Pignède, A. (2014). Detumbling of the NTNU Test Satellite. Project thesis, NTNU.
[14] Rein, O. y. (2014). Developing an ADCS Prototype for NTNU Test Satellite. Master thesis,
NTNU.
[15] Soglo, P. K. r. (1994). 3-Aksestyring av Gravitasjonsstabilisert Satelitt Ved Bruk av Magnetspoler. Master thesis, NTH.
[16] Stray, F. (2010). Attitude Control of a Nano Satellite. Master thesis, UiO.
[17] The CubeSat Program (2009). Cubesat design specification. The CubeSat Program, California Polytechnic State . . . .
[18] Tudor, Z. (2011). Design and Implementation of Attitude Control for 3-axes Magnetic Coil
Stabilization of a Spacecraft. Master thesis, NTNU.
[19] Vik, B. r. (2014). Integrated Satellite and Inertial Navigation Systems. Department of Engineering Cybernetics, NTNU, Trondheim.
© Copyright 2024