GranaSAT-I: Kinematic and Dynamic Model for Attitude Determination and Control System

GranaSAT-I: Kinematic and Dynamic Model for Attitude Determination and Control System Luis Marchante Arjona Trabajo Fin de Grado en Física Tutor: Andrés María Roldán Aranda Departamento de Electrónica y Tecnología de Computadores Universidad de Granada 1 Abstract This Bachelor Thesis introduces a Theoretical Description and Simulink implementation of dierent ADCS (Attitude Determination and Control System), a fundamental feature of spacecrafts. As it makes possible to control the orientation and stabilize the satellite within its orbit. Which is as complex as fascinating. Mathematical Formalism and Physical Model are presented in order to fully understand this phenomenon and build it to make it work in Simulink, creating a simulation of a Digital Dynamical System that would be present at the on board computer of the satellite GranaSAT-I. In this context, the Aerospace GranaSAT Group at UGR is building the GranaSAT-I. Este Trabajo de Final de Grado presenta una descripción teórica y su implementación Simulink de diferentes ADCS (sistema de control y determinación de la orientación), una característica fundamental de las naves espaciales. Ya que hace posible controlar la orientación y estabilizar el satélite dentro de su órbita. Algo que resulta tan complejo como fascinante. El formalismo matemático y el modelo físico se presentan para comprender completamente este fenómeno y construir un modelo para hacerlo funcionar en Simulink, creando una simulación de un sistema digital dinámico que estará presente en el ordenador a bordo del satélite GranaSAT-I. En este contexto, el Grupo Aeroespacial GranaSAT de la UGR está construyendo el GranaSAT-I. 2 Agradecimientos Este trabajo, como tal, supone cerrar un capítulo de mi formación académica, pero también personal. Física es una carrera que potencia aspectos tales como el espíritu crítico, el razonamiento, creatividad, abstracción... además de fraternidad o trabajar en grupo, pero va más allá de todo eso, en esencia, hace que saque lo mejor y lo peor de uno, haciéndole a uno más humano. Respecto a este proyecto, quiero reconocer la atenta supervisión y consejos del Profesor Andrés Roldán que ha sabido con mucho esfuerzo, crear la semilla de lo que un día llegará a materializarse como un satélite. Es absolutamente necesario recordar el apoyo de mis padres, incondicional, sincero y afable. Este trabajo está completamente dedicado a mi Madre y a mi Padre. Con pundonor. Construir un satélite no es trivial, y en particular un ADCS. Es necesario tener ambición y visión, pero sobretodo paciencia y atención en los detalles. Todo surge como una idea, que se puede explicar con ecuaciones físicas y modelos que describen la realidad hasta cierto punto, pero precisamente, poder modicar una pequeña parte de esta realidad para que haga exactamente lo que uno quiera, es complejo. Es decir, plasmar esas ideas en algo tan real que se puede construir y trascender los límites clásicos de la existencia humana, es algo sencillamente fascinante... Per Aspera, Ad Astra. 3 CONTENTS CONTENTS Contents 1 Introduction 6 1.1 Attitude control for GranaSAT-I . . . . . . . . . . . . . . . . . . 6 1.2 ImportanceofADCS......................... 6 1.3 Structure of this work . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Mathematical Formalism 8 2.1 ReferenceFrames........................... 8 2.2 Vector Transformations . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Euler angles and Quaternions . . . . . . . . . . . . . . . . . . . . 10 2.4 Lyapunov stability . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4.1 Positive denite function . . . . . . . . . . . . . . . . . . 11 2.4.2 Stability of autonomous systems . . . . . . . . . . . . . . 11 2.4.3 Uniform global stability . . . . . . . . . . . . . . . . . . . 12 2.4.4 Energy of the satellite . . . . . . . . . . . . . . . . . . . . 12 2.4.5 Lyapunov Function . . . . . . . . . . . . . . . . . . . . . . 13 3 Satellite Model 14 3.1 Orbital Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2 Dynamics ............................... 16 3.3 Kinematics .............................. 18 3.4 Environmental and Actuators Torques . . . . . . . . . . . . . . . 19 3.4.1 Environmental Torques . . . . . . . . . . . . . . . . . . . 19 3.4.2 Actuators Torques . . . . . . . . . . . . . . . . . . . . . . 22 4 Simulation Results 24 4.1 Angular Velocity Feedback Controller . . . . . . . . . . . . . . . 25 4.2 Attitude Feedback Controller . . . . . . . . . . . . . . . . . . . . 27 4.3 BdotController ........................... 29 5 Conclusions 31 5.1 Energy................................. 32 5.2 MagneticMoment .......................... 33 5.3 Stability................................ 34 6 Future Work and Recommendations 36 4 CONTENTS CONTENTS NOTATION v 3x1 vector vT vector transpose vb vector expressed in b (body) frame Rb a 3x3 rotation matrix from frame a to b RT matrix transpose qb a 4x1 quaternion from frame a to b ωb ib angular velocity from body to inertia expressed in body frame ωb ob angular velocity from body to orbit expressed in body frame ωo io angular velocity from orbit to inertia expressed in orbit frame ωb io angular velocity from orbit to inertia expressed in body frame I 3x3 inertia tensor I3x3 3x3 identity matrix S(v) 3x3 skew-symmetric form of a vector 5 1 INTRODUCTION 1 Introduction 1.1 Attitude control for GranaSAT-I The attitude of a spacecraft is its orientation in space. [1] The attitude of the satellite tends to vary under the action of couples, which may be external, due to radiation pressure or atmospheric drag on solar panels, or internal, due to mechanical motion of the instrument motors. A stabilization system is thus required to maintain the satellite in the right position relative to the local orbital frame. Figure 1: GranaSAT-I pointing ground station. Potential implementation of a LED matrix and a camera in the cubesat. 1.2 Importance of ADCS The attitude determination and control system (ADCS) plays an indispensable part in satellite on-orbit operation which could greatly aect the satellite's performance. When the satellite is dropped by its launcher, e.g. ISS, it certainly will not be stable but it will twist upon itself. The rst mission of this system is thus to stabilize the satellite. 6 1.3 Structure of this work 1 INTRODUCTION All in all, ADCS allows us performing precise measurements or observing towards Earth's surface in orbit, which requires the satellite to be stable, as well as for receiving or transmitting telemetry data or collecting sun's energy. 1.3 Structure of this work For a proper understanding of the movement and orientation of the spacecraft, with GranaSAT-I characteristics (4), is natural to dene, rstly, a Reference Frame in which to study Kinematic (equations that express how position changes for a given velocity) and Dynamic (equations describing how velocity changes for a given force) properties of the satellite. Throughout this work there will be dened ve dierent reference frames, being of special importance Orbital frame and Body Frame , B (2.1), and from where most of equations are dened. In the Satellite Model section, using Euler's Equations for rigid body dynamics in B for computing Attitude Dynamics results to be useful for small maneuvers, but has singularities, therefore Quaternions are introduced. (2.3). Then, environmental torques are introduced, important to notice that aerodynamic pressure/drag and solar radiation pressure are not modeled in Simulink. In Mathematical Model section, environmental and actuators torques are modeled, thus obtaining the equations necessary for Attitude Determination and Stabilization through Lyapunov's Stability (2.4). Finally, using Simulink (6) environment to perform various simulations in order to test three dierent Attitude Control Algorithms (4.1), yields consistent results that lead to Conclusions (5) about energy, applied torque and stabilization. Additionally, Recommendations (6) for future work are presented. Appendix (6) is not necessary to read in order to understand this work. 7 2 MATHEMATICAL FORMALISM 2 Mathematical Formalism 2.1 Reference Frames This section introduces dierent reference frames for representing satellite position and attitude. Essentially, we have: Inertial frame ( ECI ): The Earth-centered inertial frame is an inertial frame for terrestrial navigation, which means that it is a non-accelerated reference frame in which Newton's laws of motion apply. The origin of the frame is located at the center of the Earth. Here sub-index i is used. Orbit frame ( O ): The orbit frame has it's origin in the satellites center of mass. The xo -axis points in the normal direction of the orbital plane, while the zo -axis points to the Earth center and the yo -axis completes the right-hand orthogonal system. It is a non-inertial reference frame. Body frame ( B ): The body-xed reference frame is a moving coordinate frame which is xed on the vessel. The axes are locked in the satellite, xb -axis is forward, zb -axis is downwards and the yb -axis completes the right-hand orthogonal system. The origin is situated at the center of mass. This coordinate system rotates with the satellite, xb , yb and zb coincide with principal axes of the moment of inertia tensor. Figure 2: Inertial frame (I), Orbit frame (O) and Body frame (B) representations. 8 2.2 Vector Transformations 2 MATHEMATICAL FORMALISM 2.2 Vector Transformations This section contains the main principles when transforming between dierent reference frames, and introduces unit quaternions and the inertia matrix. Rotation Matrix In terms of generic angles between each axes from one reference frame to the other, θ , the rotation matrix may be written: Ra0 a=  cos(θx0x)cos(θx0y)cos(θx0z) cos(θy0x)cos(θy0y)cos(θy0z) cos(θz0x)cos(θz0y)cos(θz0z)   (1) Cross product operator The vector cross product is dened by [22] : λ×a=S(λ)a (2) where S is dened as: S(λ) = −S(λ)T=  0−λ3λ2 λ30−λ1 −λ2λ10  ;λ=  λ1 λ2 λ3   (3) The super-index T means the transpose of the matrix. The angle-axis parameterization of the rotation matrix, Rλ,θ , corresponding to a rotation θ about the λ -axis: Rλ,θ =I3×3+sin(θ)S(λ) + (1 −cos(θ))S2(λ) (4) is an useful parameterization. Rotation matrix dierential equation The time derivative of the rotation matrix between a frame a and a frame b is: ˙ Ra b=Ra bS(ωa ba) (5) ωa ba is the angular velocity of a with respect to b expressed in a . 9 3.2 Dynamics 3 SATELLITE MODEL Keplerian Elements The parameters discussed above dene the orbit and the position of the satellite on the orbit. These parameters constitute the six orbital elements, also known as the Keplerian elements. They are generally organized in the following order: a, e, i, Ω, ω, M The parameter a has dimensions of length, whilst the ve others ( e and the four angles) are dimensionless. We should ask why there are 6 parameters. Here are two equivalent reasons:  Three points dene the position of a solid in space. Once the point O is xed, 6 parameters (2 times 3 position coordinates) dene the two other points.  The position of a point (3 position coordinates) and its velocity (3 velocity components) at a given time can provide the initial conditions required to integrate the equations of motion, thereby dening the position of a point on its trajectory. Since its axes are specied (OP orients the ellipse), the ellipse is characterized by two parameters:  Length of the semi-major axis and eccentricity, respectively: a and e Using standard astronomical notation, the three Euler angles are :  Ω , called here the right ascension of the ascending node or the longitude of the ascending node: Ω = (Ox,ON) .  i , called here the inclination. This is the dihedral angle i = (E,P) between the equatorial and orbital planes, i = (Oz,OZ) .  ω , called here the argument of the perigee: ω = (ON,OP o ) . 3.2 Dynamics A satellite can be regarded as an ideal rigid body. The dynamic model of the satellite is derived using a Newton-Euler formulation, where the angular momentum change is related to applied torque. Euler's second law states that the rate of change of angular momentum L about a point that is xed in an inertial reference frame (often the mass center of the 16 3.2 Dynamics 3 SATELLITE MODEL body), is equal to the sum of the external moments of force (torques) acting on that body about that point: d~ L dt =~τ (24) Angular momentum can be expressed as: ~ Lo=Io~ωo (25) Therefore, the rate of change of angular momentum vector is broken up in two parts. Change in angular momentum due to angular acceleration I˙ ~ω and change in angular momentum due to inertia tensor rotation (considering a generic rotation matrix R=R(t) ) which using rotation matrix dierential equation (5), ˙ R=RS(ω) , S(ω)R=ω×R , results in : ~ω ×(I~ω) . ~τ =d~ L dt =d(I~ω) dt =I˙ ~ω +d(RIRT) dt ~ω =I˙ ~ω + ( ˙ RIRT+RI ˙ RT)~ω (26) Using anticommutativity of cross product a×b=−(b×a) and a property of skew-symmetric operator ( S=−ST ): I˙ ~ω + ( ˙ RIRT+RI ˙ RT)~ω =I˙ ~ω +S(~ω)RIRT~ω +RIRTST(~ω)~ω = (27) =I˙ ~ω +~ω ×(I~ω) + I~ω ×~ω =I˙ ~ω +~ω ×(I~ω) (28) Having: ~τ =I˙ ~ω +~ω ×(I~ω) (29) The satellite model is: I˙ ~ωb ib +˙ ~ωb ib ×(I˙ ~ωb ib) = ~τb (30) where I is the moment of inertia, ~ωb ib is the angular velocity from body to inertia decomposed in body frame and ~τb are the torques acting on the satellite also decomposed in body frame (such as gravitational torque). Another way of representing the previous equation is : I˙ ~ωb ib +S(~ωb ib)I~ωb ib =~τb (31) The angular velocity of the satellite relative to the inertial frame is expressed in the body frame according to: ~ωb ib =~ωb io +~ωb ob =Rb o~ωo io +~ωb ob (32) Where ~ωb io angular velocity from orbit to inertia decomposed in body frame. 17 3.3 Kinematics 3 SATELLITE MODEL 3.3 Kinematics The kinematics describes the satellite's orientation in space and is derived by integration of the angular velocity. First, is presented the most general case of a rotation matrix. A rotation matrix may also be referred to as a direction cosine matrix, because the elements of this matrix are the cosines of the unsigned angles between both reference frames. Let a be a reference frame with unitary vectors ~ex, ~ey, ~ez and a0 other reference frame with unitary vectors ~ex 0, ~ey 0, ~ez 0 . Then orientation from a to a0 is completely determined by transformation matrix Ra0 a , this is the direction cosine matrix (DCM): Ra0 a=  ~ex 0·~ex~ex 0~ ·ey~ex 0~ ·ez ~ey 0·~ex~ey 0·~ey~ey 0~ ·ez ~ez 0~ ·ex~ez 0·~ey~ez 0·~ez   (33) In terms of generic angles between each axes from one reference frame to the other, θ , the rotation matrix may be written: Ra0 a=  cos(θx0x)cos(θx0y)cos(θx0z) cos(θy0x)cos(θy0y)cos(θy0z) cos(θz0x)cos(θz0y)cos(θz0z)   (34) Euler's Theorem: Any two independent orthonormal coordinate frames may be related by a minimum sequence of rotations (less than four) about coordinate axes, where no two successive rotations may be about the same axis. Then, it is possible to bring a rigid body into an arbitrary orientation by performing three successive rotations. The composition of three rotations, one over each (x, y, z) axes, also suppose θ=θ(t) : R(t) =   cos(θ3)sin(θ3) 0 −sin(θ3)cos(θ3) 0 0 0 1    cos(θ2) 0 −sin(θ2) 0 1 0 sin(θ2) 0 cos(θ2)    1 0 0 0cos(θ1)sin(θ1) 0−sin(θ1)cos(θ1)   (35) where θ1 is a rotation over x -axis (roll usually represented by φ ), θ2 is a rotation over y -axis (pitch, θ ) and θ3 is a rotation over z -axis (yaw, ψ ). Finally, kinematics in this matrix form yields: dR(t) dt =R(t)  0−ω3ω2 ω30−ω1 −ω2ω10   (36) 18 3.4 Environmental and Actuators Torques 3 SATELLITE MODEL Euler angles representation of kinematics:   ˙ θ1 ˙ θ2 ˙ θ3  =1 cos(θ2)  cos(θ3)sin(θ1)sin(θ2)cos(θ1)sin(θ2) 0cos(θ1)cos(θ2)−sin(θ1)cos(θ2) 0sin(θ1)cos(θ1)    ω1 ω2 ω3   (37) Euler angles representation is useful for small maneuvers, but has singularities at cos(θ2) . This is the reason we use quaternions, applying a conversion from Euler angles to quaternions:     ˙q1 ˙q2 ˙q3 ˙q4     =1 2     0ω3−ω2ω1 −ω30ω1ω2 ω2−ω10ω3 −ω1−ω2−ω30         q1 q2 q3 q4     (38) Simplied: ˙q=˙η ˙ε=1 2−εT ηI3×3+S(ε)ωb ob (39) To nd the rotation velocity for the body frame relative to the orbit frame: ~ωb ob =~ωb ib −~ωocb 1=~ωb ib −Rb o~ωo io (40) Where c are columns in Rb o= (cb 1cb 2cb 3) . 3.4 Environmental and Actuators Torques 3.4.1 Environmental Torques In order to design the attitude control and prediction system, environmental disturbance torques acting on the spacecraft shall be modeled suciently. The torques must be modeled as a function of time, the spacecraft's position and attitude so that they can be integrated to Euler's equations and any other mathematical models. The dominant sources of environmental disturbance torques on the spacecraft attitude are the solar radiation pressure, aerodynamic drag and Earth's gravitational and magnetic elds. 19 3.4 Environmental and Actuators Torques 3 SATELLITE MODEL Figure 5: Environmental disturbance torques as a function of altitude. Diagrams like this are strongly dependent on the mass and geometry of the spacecraft, although proves to be useful for having an idea of the typical order of magnitude of these torques. The solar radiation pressure is eective on attitude of the satellite for altitudes higher than 1000 km. The gravity gradient disturbance are most signicant below 1000 km. Aerodynamic perturbations are most eective below 500 km and negligible over 1000 km altitudes. Gravity Gradient There are many mathematical models for gravity gradient torque. The most common one can be derived (careful derivation is explained in [1] p. 530-533) assuming homogeneous mass distribution of the Earth, the gravity gradient is: ~τgrav =3µ R3 o ~ue×(I~ue) (41) where µ =3.986 ·1014 m 3· s −2 is the Earth`s gravitational coecient, Ro is the distance from Earth's center (m), I is the inertia tensor and nally, ~ue is the unit vector towards nadir, i.e., downward-facing viewing geometry, usually pointing the center of Earth (its opposite is the zenith). For example, a spacecraft in a low earth orbit (LEO: altitude between 200 km and 1500 km) has 3µ R3 o ∼ =4·10−6 s −2 with the moments of inertia of the Space 20 3.4 Environmental and Actuators Torques 3 SATELLITE MODEL Shuttle are on the order of 106 kg · m −2 , so the gravitational torques on this large vehicle are on the order of 1 N · m. In body frame: ~τgrav = 3ω2 o~cb 3×(I~cb 3) ; ω2 o=µ R3 o (42) where cb 3 is, again, the third column of the rotation matrix, Rb o , which transforms zb into zo , using quaternions: c3=  2(ε1ε3−ηε2) 2(ε2ε3+ηε1) 1−2(ε2 1+ε2 2)   (43) Thus: ~τb grav = 3ω2 o  2(Iz−Iy)(ε2ε3+ηε1)(1 −2(ε2 1+ε2 2)) 2(Ix−Iz)(ε1ε3−ηε2)(1 −2(ε2 1+ε2 2)) 2(Iy−Ix)(ε1ε3−ηε2)(ε2ε3+ηε1)   (44) Solar Radiation Pressure Radiation pressure is the pressure exerted upon any surface exposed to electromagnetic radiation. For example, if the eects of the sun's radiation pressure on the spacecraft of the Viking program had been ignored, the spacecraft would have missed Mars orbit by about 15,000 kilometers being the average distance between Earth and Mars about 225 million km. The intensity of the solar radiation varies over time, this makes the determination of its energy and frequencies dicult. Most analysis uses the solar radiation constant SF [2] : SF = 1353 W m2 (45) the force of solar pressure per unit area is then given by: pSR =SF c= 4.51 ·10−6N · m−2 (46) The expression of the force due to solar radiation on the satellite then becomes: FSR =−pSR ·CR·A·r (47) 21 3.4 Environmental and Actuators Torques 3 SATELLITE MODEL where A is the exposed area to the Sun. (  is the symbol for the sun,  for Earth) The reectivity, CR , indicates how the satellite reects incoming radiation, and its value is between 0.0 and 2.0. Because CR is time variant and the constant change in orientation of the object to the sun [6] . The resulting magnitude of the torque can be expressed as: τSR =FSR(cpSR −cg) (48) where cpSR is the center of solar radiation pressure and cg the center of gravity. For more information on this topic, see [22, 4]. Aerodynamic Pressure Satellites orbiting the Earth at low altitude will be inuenced by the air density. This disturbance is most eective on satellites orbiting below 400-500 km. This may reduce the velocity of the satellite, and the result will be lower altitude for the satellite. The torque is written [12]: ~τaero =1 2ρV 2CdAinc(~uv×((cp−cg)~uv)) (49) where ρ is the atmospheric density, Cd is the drag coecient, Ainc is the area perpendicular to uv , which is the unit vector in velocity, V , direction. Faero =1 2ρmV2CdAinc (50) this force is known as lift force and can be easily obtained from basic uid dynamics. Magnetic Disturbance This torque is resulted from the interaction of Earth's magnetic eld and spacecraft's residual magnetic eld. If ~m is the sum of all magnetic moments in the satellite, the torque acting on the satellite [3] : ~τ =~m ×~ B (51) where ~ B is Earth's magnetic eld vector can be described using IGRF or Dipole Model, see [3, 4]. ~m is caused by satellite-generated current loops, permanent magnets or induced magnets which should be computed. 3.4.2 Actuators Torques Reaction wheels, momentum wheels, or magnetic torquers are devices used for the changing satellite's angular momentum. They are simply used on spacecraft for several aims: to add stability against disturbance torques, to absorb cyclic torques, and to transfer momentum to the satellite body for slewing maneuvers. 22 3.4 Environmental and Actuators Torques 3 SATELLITE MODEL Magnetic Torquers Of special interest, as GranaSAT-I will use this type of torquer. Torque produced by the magnetic torquers in body frame is : ~τb m=~mb×~ Bb (52) ~mb is the magnetic dipole moment generated by the torquer, ~ Bb is the local geomagnetic eld vector, relative to the satellite. It is interesting to mention that these kind of magnetic actuators can only create a torque within a plane (perpendicular to the local magnetic eld), which may represent a limitation as for ADCS with these torquers. Magnetic dipole moment is given by : ~mb=  NxixAx NyiyAy NzizAz  =  mb x mb y mb z   (53) where, Nk is the number of windings in the torquer, Ak is the span area of the coil, and ik the torquerer current. Using skew-symmetric operator : ~τb m=S(~mb)~ Bb=  Bb zmb y−Bb ymb z Bb xmb z−Bb zmb x Bb ymb x−Bb xmb y   (54) 23 4 SIMULATION RESULTS 4 Simulation Results GranaSAT-I is going to be a nanosatellite (10x10x10 cm cube) with approximately 1 kg mass , supposed to carry a camera as the payload, which should take pictures of the city of Granada. It will also have a very powerful LED matrix that must be seen from Earth surface and follow a LEO (Low Earth Orbit), that is typically consider from 200 km to 1500 km. As well as having an (diagonal) inertia matrix: Ixx = 0.0018 , Iyy = 0.0017 , Izz = 0.0015 (m 2· kg) which is given by a Solidworks 3D model prototype. This simulated satellite follows a LEO orbit, typical in cubesats, at h= 750 km with very small eccentricity e= 0.005 as well as an inclination of i= 17.2 º . Magnetorquers physical properties in each axis are,      Coil Resistances 50 Ω Coil Areas 88 cm2 Intesity Limits 100 mA Inputs of the simulations, are essentially the quaternion q , and ~ωb ib , angular velocity from body to inertia decomposed in body frame. Various initial values haven been tested, as we would not know them until the satellite in orbit starts measuring, they are chosen randomly, those represented in the following gures are the result of,      Orientation :Quaternion q = [−0.001 0.957 0.0928 −0.275] Orientation :Euler Angles [ϕ θ ψ] = [176.688 º −31.769 º −12.020 º ] Angular V elocity ~ωb ib = [4 2 1] rad · s−1 The dierent controllers are tested on the Simulink model of this project. Henceforth, a comparative and detailed study as well as a latter discussion is given in this section. All graphs represent time in orbits, in this scenario, Torbit = 5940 s. 24 4.1 Angular Velocity Feedback Controller 4 SIMULATION RESULTS 4.1 Angular Velocity Feedback Controller The angular velocity feedback controller or Wisniewski controller [6], ~mb=H~ωb ob ×~ Bb (55) with a gain H= 500 . 0.1 0.2 0.3 0.4 0.5 0.6 -1.5 -1 -0.5 0 0.5 1 1.5 x 10-3 t (orbits) m b [Am 2 ] Applied magnetic moment by Angular Velocity Feedback Controller mx my mz 0.1 0.2 0.3 0.4 0.5 0.6 0 1 2 3 4 5 6 7 8 x 10-3 t (orbits) Energy [J] Energy consumed by Angular Velocity Feedback Controller Figure 6: Angular Velocity Feedback Controller. Topmost, applied torque. Down, energy consumed. 25 5.1 Energy 5 CONCLUSIONS 5.1 Energy Energetic eciency is of capital importance in a spacecraft. These controllers consume generally, in the order of some mJ, not B dot controller, that presents a higher energy consumption, about one order of magnitude higher. It is a result that makes sense looking a the magnetic moment generated by the coils. Additionally, the angular velocity and attitude controller present similar results, as their mathematical expression are analogous. Nevertheless, the most ecient controller turns out to be the attitude feedback controller, but as it will be explained is not as reliable for stabilization than the others. 0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 10 12 14 16 18 time (orbits) Energy [J] Total Energy Figure 12: Total Energy (expressed in equation 18) of the satellite in orbit dissipates over time. 32 5.2 Magnetic Moment 5 CONCLUSIONS 5.2 Magnetic Moment The applied magnetic moment by the controller makes possible to obtain the desired satellite's orientation. At rst, torque is applied with higher intensity, as an initial impulse, is interesting to notice how it is particularly higher for the B dot. As well, as how it is applied in dierent axis to reach stabilization. Magnetic moment is applied until the satellite is stable and with much less intensity when is stabilized, this allows the angular velocity to be aligned with Earth's local magnetic eld. Typical values shows it ranges from [-1.5 ·10−3 to 1.5 ·10−3 ] A · m 2 . 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -4 -3 -2 -1 0 1 2 3 4 5 x 10-5 time (orbits) Magnetic Field [T] Earth's Magnetic Field from Body Frame Bb x Bb y Bb z Figure 13: Earth Magnetic Field components as seen from body frame. Before stabilization and the mentioned alignment of the angular velocity, the satellite measures higher variations due to its own rotation. 33 5.3 Stability 5 CONCLUSIONS 5.3 Stability A fundamental requirement for controllers is for them to be able to reduce angular speed and align body and orbital frames, that means Rb o=I3x3 . According to equation (32), ~ωb ib =~ωb io +~ωb ob =Rb o~ωo io +~ωb ob Therefore, ~ωb ib =~ωb io +~ωb ob =~ωo io +~ωb ob with Rb o deduced in [6], Rb o= 2   1 2ε3−ε2 −ε31 2ε1 ε2−ε11 2   For stability to be studied, quaternions, or equivalently, Euler angles are shown. Euler angles, to reach stability, need to tend to [0 0 0], simultaneously, quaternions to q= [±1 0 0 0] . The three controllers give an acceptable result, angular speed is reduce to less than 0.001 rad/s. But total alignment is not fully reached. Not only that, but the attitude feedback controller tends to be unstable and after convergence to q= [±1 0 0 0] , starts to deviate slowly with small angular velocity. Therefore, this last controller is not recommended, in spite of its power eciency. After an orbit, and a slight dierence, the most stable controller is the Wisniewski controller as it gives q= [−0.9915 0.1258 0.0122 −0.03015] while B dot presents q= [−0.9823 0.1810 0.0182 −0.0521] . 34 5.3 Stability 5 CONCLUSIONS 0 0.1 0.2 0.25 -6 -4 -2 0 2 4 6 time (orbits)  b ob [rad/s] Angular Velocity x y z Figure 14: Convergence of components of ~ωb ob , angular velocity of the body frame with respect to orbit frame (expressed in body frame), to ~ωb ob = [0 0 0] rad/s. 35 6 FUTURE WORK AND RECOMMENDATIONS 6 Future Work and Recommendations This work is possible due to the maintained eort of the Aerospace GranaSAT Group carried out by my colleagues and, to some extent, to my previous collaborations to GranaSAT-I project, as in designing a 3D model in Solidworks of a prototype rst small rocket (it measured around 50 cm high), or such as simple launching tests or study of groundtrack control and nally this ADCS. This is a rst approximation of a attitude control system that will be some day controlling a real spacecraft. Continuing this work means having in consideration all the results given and improvements that can be added. A more accurate model can be created, for example, adding aerodynamic drag to the Simulink model, an easy implementable task. Or changing World Magnetic Model to IGRF. A more precise propagator, could be a good idea, although through references can be found not to be of paramount importance (providing a decent propagator). One practical idea would be to implement a function to read TLE (two-line element set) automatically so the propagator would not need to be introduced keplerian elements and more data, of course modications to the Simulink model shall be added in that scenario, for example regarding aerodynamic drag, that stronlgy depends on altitude. More complex enhancements, are the complete migration of the simulinlk code to C (for microcontroller) or HDVL (for FPGA), as the program is itself a feedback loop system and that can represent some complications, none the less, to encourge people that will follow this path, a migration to C of the propagator is given in the appendix. This, however, might be more doable, with less diculty than implementing a linearized system with a linear quadratic estimator, like Kalman lter. Again, in the appendix, the linearized model is written. 36 6 FUTURE WORK AND RECOMMENDATIONS ADCS Budget Estimation Although a nanosatellite budget can range from $2.000 to $20.000, more affordable options are taking place due to cheaper high performance electronics, making possible to design a low cost ADCS. Therefore, based on [27, 25, 16, 17, 18], a rough estimation on building this attitude control system is given: Magnetorquers Gyroscope Coils Microcontroller GPS Total High estimate $9,256.95 1 $71.87 3 $4,976.88 5 $112.20 7 $9,495 9 $23,912.90 Low estimate $22.44 2 $55.67 4 $12.45 6 $9.46 8 $99.00 10 $199.02 1 ISIS Magnetorquer Board 2 Honeywell HMC 1052L 3 ADXRS453BEYZ Analog Devices 4 ADXRS453BRGZ Analog Devices 5 CubeWheel Small Cube Sat Shop 6 MilliBird50 Didel 7 TE0887-03M Trenz Electronic GmbH 8 PIC24FJ256GA110 9 OEM4-G2L NovAtel 10 Venus838FLPx-L Navspark Further information of an overall cost a this type of mission can be found in [27] NASA Cost Symposium AMES. 37 6 FUTURE WORK AND RECOMMENDATIONS Gantt Chart of this Project 38 REFERENCES REFERENCES References [1] James R. Wertz. Spacecraft Attitude Determination and Control. Kluwer Academic, 2002. [2] Fundamentals of Spacecraft Attitude Determination and Control. F. Landis Markley, John L. Crassidis. Springer, 2014. [3] Orbital Mechanics for Engineering Students. Howard D. Curtis. Elsevier, 2005. [4] Satellites Orbits and Missions. Michel Capderou. Springer, 2005. [5] Introduction to Simulink ® with Engineering Applications. Orchard Publications. 2006. [6] Satellite Attitude Control Using Only Electromagnetic Actuation Ph.D. Thesis Rafaª Wisniewski 1996. [7] Rotaciones, secuencia aeroespacial y cuaterniones. Una revisión de las relaciones fundamentales. Eduardo Serrano, Ricardo Oscar Sirne y Guillermo La Mura. [8] AstroNet-II: The Astrodynamics Network. Attitude dynamics and control. University of Roma Tor Vergata. James Biggs, 2013. [9] Mathematical Model of the Solar Radiation Force and Torques Acting on the Components of a Spacecraft. R. M. Georgevic. NASA, 1971. [10] Attitude Determination and Control Systems Scott R. Starin, NASA Goddard Space Flight Center John Eterno, Southwest Research Institute, 2011. [11] Attitude control for the Norwegian student satellite nCube. Master Thesis. Eli Jerpseth Øverby, 2004. [12] LEO Satellites: Attitude Determination and Control Components: Some Linear Attitude Control Techniques Graduate School of Natural and Applied Sciences Of Middle East Technical University. Ceren Kaplan, 2006. [13] Controls Algorithm For A Satellite Using Earth's Magnetic Field: Orbit Maneuvers And Attitude Positioning Karthik Ganesh University of Central Florida, 2007. [14] NASA Cost Symposium AMES Cost Model Ames Micro/Nanosatellites Cost Model, 2015. [15] Satellite Attitude Determination with Low-Cost Sensors by John C. Springmann, 2013. [16] Hardware Review of an On Board Controller for a Cubesat Norwegian University of Science and Technology Trondheim Magne Alver Normann, 2015. 39 REFERENCES REFERENCES [17] Review of the hardware description of the camera module prototype for NTNU Test Satellite (NUTS) By Jon Kalevi Oltedal Fall Semester, 2015. [18] Attitude Determination and Control System design of KufaSat. Mohammed Chessab Mahdi and Abdal-Razak Shehab, 2014. [19] Useful to visulize Euler to Quaternions: http://quaternions.online [20] http://www.kostasalexis.com/frame-rotations-and-representations.html [21] Comunicaciones por Satélite Curso 2009/10 Comunicaciones por Satélite. © Ramón Martínez, Miguel Calvo UPM. Parámetros orbitales Formato TLE de NORAD. [22] Spacecraft Attitude Dynamics and Control Christopher D. Hall January 12, 2003 [23] Astronautica y Vehículos Espaciales. Rafael Vazquez Valenzuela Departamento de Ingeniera Aeroespacial Escuela Superior de Ingenieros, Universidad de Sevilla. 2012. [24] Satellite Dynamics and Control in a Quaternion Formulation - Lecture note for course 31365 Spacecraft Dynamics and Control at DTU. 2010. [25] Pose Estimation using Linearized Rotations and Quaternion Algebra. Timothy Barfoota, James R. Forbesa, Paul T. Furgalea .University of Toronto Institute for Aerospace Studies. [26] Basic Navigational Mathematics, Reference Frames and the Earth's Geometry. Springer-Verlag Berlin Heidelberg, 2013. [27] Modeling and Simulation for Automatic Control Olav Egeland and Jan Tommy Gravdahl. Norwegian University of Science and Technology Trondheim, Norway. 2002. 40 REFERENCES REFERENCES APPENDIX Linearization of the Satellite Model (for future implementation) What Is Linearization? (Brief Explanation) Linearization is a linear approximation of a nonlinear system that is valid in a small region around the operating point, using a Taylor series. For example, suppose that the nonlinear function is y=x2 . Linearizing this nonlinear function about the operating point x= 1 , y= 1 results in a linear function y= 2x−1 . Near the operating point, y= 2x−1 is a good approximation to y=x2 . Away from the operating point, the approximation is poor. The actual region of validity depends on the nonlinear model. Extending the concept of linearization to dynamic systems, we can write continuoustime nonlinear dierential equations in this form : ˙x(t) = f(x(t), u(t), t) y(t) = g(x(t), u(t), t) (58) In these equations, x(t) represents the system states, u(t) represents the inputs to the system, and y(t) represents the outputs of the system. A linearized model of this system is valid in a small region around the operating point we linearize. To represent the linearized model, dene new variables centered about the operating point : δx(t) = x(t) = xo δu(t) = u(t)−uo δy(t) = y(t) = yo (59) The linearized model in terms of δx(t) , δu(t) , and δy(t) is valid when the values of these variables are small : δ˙x(t) = Aδx(t) + Bδu(t) δy(t) = Cδx(t) + Dδu(t) (60) Kinematics The kinematics of the satellite is already described : 41 REFERENCES REFERENCES Omega_o es la velocidad angular orbital del satélite, expresada en rad/s. R_B_O es la matriz de rotación de Orbit a Body frame. Torque Gravitatorio 1 R_B_O 3*(1.243*10^-3)^2 3*omega_o^2 A B C Cross Product C = AxB 3x3 Cross Product 1 tau_g Select Columns Seleccionar C3 1 -0 0 Display -0 0 0 Display1 fcn I Tensor de Inercia [m^2*kg] Matrix Multiply Matrix Multiply 0.1043 0 0 Display2 tau_grav 48 REFERENCES REFERENCES Da vector campo magnetico en la orbit frame, justo lo que necesito, funciona OK. Los valores del campo magnetico son los esperados comparando con la bibliografía. Aquí se calcula la derivada del campo magnetico terrestre en la body frame Cuidado con unidades, he metido essto para que sean Teslas y no nT. Corrección de la altura que da el propagador, ya que toma como referencia el centro de la Tierra, no la altura desde la superficie terrestre. 1 R_B_O 3 B_o Matrix Multiply Matrix Multiply du/dt Derivative1 Scope B_o Matrix Multiply Matrix Multiply1 2 w_B_OB A B C Cross Product C = AxB 3x3 Cross Product Add Scope B_b 1 B_b 2 B_bdot h (m) μ (deg) l (deg) Magnetic Field (nT) World Magnetic Model 2010 3 LLA_pos |u| Abs -KECI to h 1e-9 Gain1 B_b To Workspace 49 REFERENCES REFERENCES Momento magnetico a aplicar Velocidad angular de Body frame respecto de la Orbit frame expresada en Body frame Campo magnético terrestre en body frame. En el bloque de control se encuentra el algoritmo que calcula cuánto momento magnetico hay que aplicar a los magnetorquers que lleva el satélite para que este se alinee localmente con el campo magnetico de la Tierra. En este caso, se muestra el Wisniew ski controller. 1 w_B_OB 2 B_b 1 m A B C Cross Product C = AxB 3x3 Cross Product 500 Gain Matrix H Scope m -0.002669 -0.06222 -0.0009568 Display m 2.096e-05 -1.028e-06 8.336e-06 Display w_B_OB -2.561e-28 -0.09128 5.935 Display B_b m To Workspace1 50 REFERENCES REFERENCES Momento magnetico a aplicar Velocidad angular de Body frame respecto de la Orbit frame expresada en Body frame Campo magnético terrestre en body frame. En el bloque de control se encuentra el algoritmo que calcula cuánto momento magnetico hay que aplicar a los magnetorquers que lleva el satélite para que este se alinee localmente con el campo magnetico de la Tierra. En este caso, se muestra el Wisniew ski controller. 1 w_B_OB 2 B_b 1 m A B C Cross Product C = AxB 3x3 Cross Product 500 Gain Matrix H Scope m -0.002669 -0.06222 -0.0009568 Display m 2.096e-05 -1.028e-06 8.336e-06 Display w_B_OB -2.561e-28 -0.09128 5.935 Display B_b m To Workspace1 51 REFERENCES REFERENCES B DOT CONTROLLER 2 w_B_OB 1 B_b 1 m A B C Cross Product C = AxB 3x3 Cross Product Dot Product 1.25e-8 Gain Scope1 Scope2 Divide m To Workspace1 52 REFERENCES REFERENCES ATTITUDE FEEDBACK CONTROLLER 1 w_B_OB 2 B_b 1 m Scope m 3 q w_B_OB B_b Epsilon m Attitude Control q Epsilon Demux m To Workspace5 53 REFERENCES REFERENCES B DOT CONTROLLER 2 w_B_OB 1 B_b 1 m A B C Cross Product C = AxB 3x3 Cross Product Dot Product 1.25e-8 Gain Scope1 Scope2 Divide m To Workspace1 54 REFERENCES REFERENCES ATTITUDE FEEDBACK CONTROLLER 1 w_B_OB 2 B_b 1 m A B C Cross Product C = AxB 3x3 Cross Product 500 Gain Scope1 Scope3 3 Epsilon A B C Cross Product C = AxB 3x3 Cross Product1 0.9 Alpha 55 REFERENCES REFERENCES derivada de velocidad angular de la body frame respecto a la inertial frame velocidad angular de la body frame respecto a la inertial frame Torque total aplicado sobre el satelite Este bloque resuelve la ecuación diferencial dada en Eq. 30, que modela la dinámica del satelite Se inicializa a una velocidad angular aleatoria, ya que en la realidad no es un parámetro que desconocemos inicialmente. 1 tau Add A B C Cross Product C = AxB w x (Iw) 1 w_B_IB_dot 2 w_B_IB 1 s xo Integrator fcn u v w y Tensor de Inercia Invertido fcn u v w y Tensor de Inercia fcn y W_B_IB Inicial 56 REFERENCES REFERENCES velocidad angular de la body frame respecto a la inertial frame derivada quaternion componentes del quaternion. Primero, se introduce el quaternion inicial, se normaliza, va a "Orbit to Body" donde se obtiene R_B_O y la velocidad angular desdeada para resolver la Eq. 58 que nos da la derivada del quaternion, es decir, cómo cambia la orientación del satélite. Seguidamente, se descompone el quaternion en su parte escalar (eta) y vectorial (e1,e2,e3) en el demux. Se calcula las tres operaciones necesarias descritad en Eq. 58 quaternion Matriz de rotacion de Orbit a Body velocidad angular de Body frame respecto de la Orbit frame expresada en Body frame Los quateriones siempre tienen que estar normalizados a la unidad. 1 w_B_IB w_B_IB q w_B_OB R_B_O Orbit to body 3 R_B_O 4 w_B_OB Dot Product -0.5 Gain4 1 q 2 q_dot qnormal(q) Quaternion Normalize 1 s xo Integrator y Quaternion Inicial 57