ARTÍCULO ORIGINAL

Design and simulation of an attitude determination system based on the Extended Kalman Filter for Cube-Sat Colombia I

Diseño y Simulación de un Sistema de Determinación de Actitud basado en el Filtro Extendido de Kalman para el Cubesat Colombia I

Cesar A. Castellanos^*, Lilia Edith Aparicio Pico

Research Group GITEM, Engineering Faculty, Universidad Distrital FJC. Carrera 7 # 40-53. CP 110231. Bogotá DC., Colombia.

^*Autor de correspondencia: teléfono: + 57 + 1 + 814 87 97, correo electrónico: divertech07@gmail.com (C. Castellanos)

(Recibido el 16 de febrero de 2012. Aceptado el 14 de enero de 2014)

Abstract

The purpose of this document is to describe the design and simulation of an attitude determination system (ADS) for the CubeSat Colombia I picosatellite, based on an extended Kalman filter EKF [1,2]. In this development the propagation of the State of the system has been implemented through a Runge-Kutta procedure instead of traditional State Transition matrix method, improving predictions from the EKF filter without making direct measurements of angular velocity, it was obtained a reduced prediction time by incorporating a deterministic algorithm (the TRIAD algorithm) within the EKF. The results validated through simulation are those shown in this article.

Keywords: Attitude determination system, extended Kalman filter, non deterministic and deterministic algorithm, picosatellite Colombia I

Resumen

El propósito de este documento es describir el diseño y simulación de un Sistema de Determinación de Actitud (ADS) para el picosatélite Cubesat Colombia I, basado en el Filtro Extendido de Kalman (EKF) [1, 2] . En este desarrollo la propagación del estado del sistema ha sido implementada a través de un procedimiento Runge-Kutta a cambio del método tradicional por matriz de transición de estado, mejorando las predicciones del filtro EKF, sin hacer mediciones directas de velocidad angular. Se obtuvo una reducción del tiempo de predicción incorporando un algoritmo determinístico (El algoritmo TRIAD) dentro del EKF. Los resultados mostrados en este artículo se validaron mediante simulación.

Palabras clave: Sistema de determinación de actitud, filtro extendido de Kalman, algoritmo determinístico y no determinístico, picosatélite Colombia I

Introduction

The System for determination of a Satellite's Attitude (ADS, by the English acronym for Attitude Determination System) establishes the Cube-Sat's orientation at the point of its terrestrial orbit, with respect to different reference frames in particular one centered on Earth. This attitude will be used by another functional block, the ACS (Attitude Control System) [3, 4] to modify this attitude with the ultimate aim of pointing communication antennas towards the Earth, or any other pointing requirement. To achieve this goal, the Cube-Sat makes use of signals from sunlight and magnetic field sensors, mathematical models for satellite Orbit , terrestrial magnetic field , and position of the Sun by establishing a Solar Vector, and a model of the plant (the picosatellite). This information is combined optimally in a so-called extended Kalman filter algorithm to obtain a good estimate of the attitude of the picosatellite. These concepts are expanded in the following sections .

Necessary elements for attitude determination

Attitude Definition and representation

All attitude can be conceived as a rotation of a reference frame [5-7] with respect to another, some angle around an axis. The angle and rotation axis are required to establish the attitude. A common representation of attitudes is through quaternions. A Quaternion is an arrangement of four real numbers, three of which are multiplied by the imaginary numbers i, j and k, as shown in equation (1).

The first three elements are called the vector part of the Quaternion and the remaining part is known as a scalar of the quaternion. Another way to represent a Quaternion is as in equation (2) a 4x1 matrix:

The rotation axis and angle information are incorporated into the Quaternion with the formula shown below in equation (3):

Where n represents a unit vector in the direction of the axis of rotation in the frame of reference and is the angle rotated in radians or degrees.

Reference frames [4-7]

The orientation of the picosatellite is related to reference frames associated with the earth, the following reference frames and corresponding coordinate systems may be defined: The reference frame linked to the Cube-Sat is called the BODY. The origin of this frame is the picosatellite's center of mass and its axes are parallel to the faces of the picosatellite and mutually perpendicular. The reference frame associated with the Translational motion of the Earth considered as inertial, referred as ECI (by the English acronym Earth-Centered Inertial frame). The origin of this system is the center of mass of the Earth and the axes are oriented so: the x axis always points towards a direction in space known as the Aries point where the Sun's movement through the ecliptic crosses the equator plane from south to north. The z axis points towards the geographical North of the Earth and coincides with its rotation axis. The y axis completes the coordinate system according to right hand rule. This system is not linked to the Earth's rotation, but shares its travel around the Sun.

Measurement and reference vectors

To determine the attitude at a given point along the orbit, measurements are taken in the BODY coordinate system and compared to vectors given in the ECI frame. Measurements are given as two unit vectors [6-9]. The first of these is generated by a three-axis magnetometer located in the picosatellite and indicates the direction of the Earth's magnetic field. The second is an array of photo-cells located in each of the 6 faces of the picosatellite (it is a cube) and provide information about the position of the Sun.

The magnetic vector measurement is compared with a model of Earth's magnetic field (so-called World Magnetic Model, WMM) that specifies the magnetic vector in a fixed point along the orbit where the measurement is being made.

The solar vector obtained from measurements of light intensity provided by photocells inserted on 6 solar panels is compared with a mathematical model of the Sun's position from the satellite, that includes the Albedo effect (light reflected from earth's surface) [10,11].

Model of the plant and the system state vector

With the information provided by the four above-mentioned vectors (two measurement vectors and two model reference vectors) we can estimate the attitude with some deterministic and/or non deterministic algorithms, which require a mathematical model of the body. The modeling represents the picosatellite as a rigid body that can rotate around some axis about its center of mass. Under these conditions the kinematics and dynamics of the rigid body satellite are described by the following equations (4, 4a, 4b, 4c) [3, 7, 9]:

Kinematics:

With

and,

Dynamics:

The Satellite Model is then Non linear and time variant

With:

w_b= angular velocity vector of the body with respect to reference frame eci. Represented in body.

q = Quaternion rotation of body with respect to frame of reference eci.

I = Tensor of inertia of the picosatellite on principal axes (similar to body axes)

τ_b = Torque applied to the picosatellite (in body)

[w_bx]= angular velocity vector cross product operator expressed as a matrix. [w_bx] v = wxv, with v a vector 3 x 1.

The information that is appropriate to incorporate in a State vector (x) of this system is the attitude of the body (represented as a Quaternion) and the angular velocity of the body (represented as a 3-element vector), so x dimension is 7x1.

As the picosatellite according to the equations (4) is a time variant nonlinear continuous system the State equation (5) that determines its evolution over time is:

Where u(t) is the input to the system (torque), w(t) represents the process noise because of the uncertainties in the mathematical model of the plant. Noise has the following characteristics, equation (6):

With E representing w(t)'s expected value, and covariance in time as in equation (7)

In other words, is a white noise, Gaussian with zero mean, with Q(t) the covariance of the noise and δ(t) the Dirac delta function. Similarly, the system here treated, has a discreet and non-linear measurements model as expressed in the equation (8):

The added vector v_k, is due to sensor's noise and features the following characteristic expressed by equations (9) and (10):

Where R(i) is the covariance of the measurement noise, and δ_ij, the Kronecker delta, which is set to 1 if i = j and 0 otherwise. Equations (9) and (10) indicate a white noise.

An extended Kalman filter is used for this type of system. To perform the prediction of the evolution ofthe State x from a time t_k-1 until a time t_k, the equation (4) is evaluated without adding noise w(t). This evaluation is done usually in the EKF through the State transition matrix [1-3, 7] (ф) both for the prediction of the future state x, and the covariance matrix of State error P. This state transition matrix is obtained linearizing and discretizing the non-linear function f. The function f is linearized with respect to the State Vector x around the operation point along the path on the State Space equation (11):

And the transition matrix as in equation (12)

The prediction of the State variation is then given by equation (13):

With matrix values as in equation (14)

Being:

x_k = true state

= state estimate at instant k

Equally, Δz_k is the measurement's innovation as in (15)

With:

z_k = measurement; = measurement estimate

h() = measurement model as a function of state estimate

If measurement model h is linearized and discretized, then measurement error is equation (16):

With H_k = linearized and discretized h.

Figure 1 shows the Extended Kalman Filter algorithm used in this development.

To preset the proper functioning of the ADS a software simulation using MATLAB was implemented which initially generated good predictions around 90 cycles through the EKF. Due to the need for prediction times as short as possible, it was incorporated a deterministic algorithm called TRIAD [5,12] which reduced the initial prediction about 20 cycles, ie, close to one fifth. This algorithm is explained below.

Algorithm TRIAD [13]

In block 2 of the simulation software [8] was implemented a small algorithm called TRIAD [5,8] that uses the information provided by two reference vectors and two vectors of measurements to give an estimate of the orientation of a vehicle. The TRIAD algorithm is as follows: If there are two vectors of reference given in the ECI frame and called _ECI and _ECI which correspond in our application to the magnetic vector, according to the WMM, and the solar vector deducted according to the satellite's orbit, then one could transform these vectors to the components of the ideal measurement vectors and seen in the body framework are called _c and _c, through a rotation matrix equation (17):

However, measurements have noise, deviations, and generally it is not possible to know previously that rotation matrix. What makes the algorithm is to build two Ortho normal vector triads and through them represent rotations to an intermediate reference frame; rotations to (or from) the ECI reference framework as well as towards (or from) the body framework, from the intermediate framework, and by composition of rotations, is obtained. To achieve this, algorithm TRIAD assumes that one of the vectors of measurement is an exact transformation of the vector in the ECI frame. The first basis vector of the reference framework built with a triad is the following in equations (18) and (19):

The second component of the intermediate reference frame is built as a unit Ortho normal vector to both measurements as in (20),

The third component of the triad is found as an Orthonormal vector to the first two components as in (21):

The second triad is obtained similarly in equation (22) using the magnetic field and solar vector models _ECI, and _ECI:

The rotation matrix from intermediate reference frame to the body framework is therefore (23):

The rotation matrix from intermediate reference frame to the ECI framework is therefore (24):

Now, to make the transformation from ECI reference frame to body reference frame, is calculated the following matrix given by equation (25):

Simulation

With the purpose of assessing the effectiveness of the filter to correctly obtain the attitude and angular velocity of a body, a simulation was implemented using MATLAB. The measurements are simulated adding Gaussian noise to a trajectory known in the State Space that must be deducted by the EKF of the ADS. For the simulation is considered that the actual values in the ECI reference frame of magnetic field and solar vector do not change during the time of simulation which is less than 1 second, accordingly these vectors only change by variations in attitudes and/or Gaussian noise in each cycle or step of the EKF. After implementing the EKF, deviations were found in the prediction of the angular velocities obtained through the State transition matrix ɸ, and opted for a deterministic procedure to obtain the propagation of the State x, specifically following the algorithm of Runge- Kutta of order 4 [5,13], using the dynamic equations of the satellite and without measuring angular velocity in the system. This procedure512 corrected the predictions in the part of the State vector corresponding to the angular velocities, and obtained deviations ≤ 2.5%.

From different simulations carried out, only the simplest with zero angular velocity are shown. The results are satisfactory, the same for simulations for nonzero angular velocity and angular acceleration (constant Torque). Suppose there is the situation shown in Figure 2, where the coordinate system associated with the body is rotated 45 degrees (meaning counterclockwise rotation) around the z-axis with respect to the coordinate reference system (ECI). According to the discussion in the section of Quaternions, the initial attitude of this system is represented by the following Quaternion given by equation (26):

Magnetic reference vector , Blue color in figure 2 is located on the positive semi-major ECI axis x, and the reference solar vector , in Orange on the positive semi-major ECI y axis . Both are unit vectors.

According to the situation shown, vectors measured by sensors mounted in the body of the picosatellite are on the frame of reference BODY in the following position or representation given by equations (27) and (28):

and

To these measured vectors we added Gaussian noise, with a covariance given by the matrix R. The angular velocity between the two coordinate systems is zero, so the angular velocity vector is (29):

Simulation graphs on figure 3, show the behavior of the extracted angle from the state vector , produced by the filter for each step k:

As can be seen, the angle provided by the filter ranges 45 degrees. In this particular simulation, the average was 45,0032 degrees with a standard deviation of 0,0151 degrees.

Averages and deviations shown below are taken from step 20, and during that period the filter converges to a minimum in the covariance matrix values.

Similarly in Figure 4 it has been plotted, from EKF's estimated State Vector x, the estimated rotation axis vector components. Again, the result is consistent with the theoretical or true conditions.

On the following table 1 related average and standard deviations are given:

The results shown in Figure 5 and Table 2 were obtained for the angular velocity:

Indeed, after 20 steps filter converges to the values expected with a minimum deviation. Angular component satatistics may be seen in table 2.

Conclusion

The extended Kalman filter was used as the central algorithm for Attitude Determination of picosatellite Colombia I. It has been generated an attitude determination system simulation implemented in MATLAB which gives successful prediction of attitude, minimizing the calculation time by incorporating a known deterministic method, as is the TRIAD algorithm, and using Runge-Kutta's method to propagate the State vector without measuring angular velocity, using the dynamic equations of the satellite managed to correct the predictions ofangular velocity set out by the ''normal'' transition matrix procedure. However the origin of this deviation must be investigated in future applications. The state estimation (which is the output of the ADS system ) constitutes input to the ACS system which has been conceived to make use of a Sliding methodology with probable mixing with fuzzy concepts [14-21]. The ADS itself may incorporate an EKF improvement by the use of UKF (Unscented Kalman Filtering) , which filter better in general but having a higher computational cost [22, 23].

Acknowledgements

This work extensively uses the undergraduate electronic engineering thesis from the student Fredy Ruiz Ch., at UDFJC whom we thank for his very important contributions in this development. This thesis was directed by Eng. César A. Castellanos G. coauthor of this article. We make a special recognition to the research group in Control Engineering of Aalborg University for having provided information and tips that have made possible the development of the system. Similarly to the Universidad Distrital FJC Faculty of engineering, the Center for Research and Scientific Development CIDC for economic support and processing of various elements that make part of this system. To Dr. Rafael Wisniewski by his valuable visit to Colombia, the members of the Group GITEM Eng. Leonardo Plazas, Eng. José Paternina, undergraduate students who have collaborated us: Oscar Urrego, John E. Goyeneche, Robinson Cadena, Alejandro Beltrán, Nestor Pulido and others who have participated to a lesser degree but always valuable and relevant.

References

1. R. Brown, P. Hwang. Introduction to random signals and applied Kalman filtering. 3^rd ed. Ed. John Wiley and Sons. USA. 1999. pp. 190-242, 335-377.         [ Links ]

2. G. Welch, G. Bishop. An Introduction to the Kalman Filter, SIGGRAPH 2001 course 8. In Computer Graphics. Annual Conference on Computer Graphics & Interactive Techniques. Los Angeles, USA. 2001. pp. 23-33.         [ Links ]

3. R. Wisniewski. Satellite Attitude Control Using Only Electromagnetic Actuation. Ph.d. Thesis. Department of Control Engineering. University of Aalborg. Aalborg, Denmark. 1997. pp.71-85         [ Links ]

4. R.Wertz. Spacecraft Attitude Determination and Control. 4^th ed. Ed. Springer. Boston, USA. 2002. pp. 485-502         [ Links ]

5. Group 05 gr 833.Attitude Determination for AAUSAT- II. Department of Control Engineering, Aalborg University. Aalborg, Denmark. 2005.pp.55-70.         [ Links ]

6. K. Krogh, E. Schreder. Attitude Determination for AAU Cubesat. Control Engineering Department, University of Aalborg. Aalborg, Denmark. 2002. pp. 42-57.         [ Links ]

7. T. Bak. Spacecraft Attitude Determination - A magnetometer approach. Ph.d. Thesis. Control Engineering Department, University of Aalborg. Aalborg, Denmark. 1999. pp. 23-65.         [ Links ]

8. J. Maybeck, S. Peter. Stochastic models, estimation and control. University of Southern California. Ed. Academic Press. Ohio, USA. Vol 141. 1979. pp. 342-358.         [ Links ]

9. F. Ruiz. Design of an Attitude Determination System for Picosatellite Colombia I. Thesis. group GITEM. Universidad Distrital FJC. Bogotá DC., Colombia. 2010. pp.25-32.         [ Links ]

10. D. Bhanderi. Spacecraft Attitude Determination with Earth Albedo Corrected Sun Sensor Measurements. Ph.d. Thesis. Aalborg University. Aalborg, Denmark. 2005. pp. 52-105.         [ Links ]

11. A. Pontus. ''Attitude Estimation from Magnetometer and Earth-Albedo-Corrected Coarse Sun Sensor Measurements''. Acta Astronautica. Vol. 56. 2005. pp. 115-126.         [ Links ]

12. H. Black. ''A passive system for determining the attitude of a satellite''. American Institute of Aeronautics and Astronautics Journal. Vol. 2. 1964. pp. 1350-1351.         [ Links ]

13. Ogata, Katsuhiko. Modern, Control Engineering. 2^nd ed. Ed. Prentice-Hall. Minnesota, USA.1993. pp. 334-337.         [ Links ]

14. T. Arif. A Satellite Attitude Control Algorithm. 2010. 2^nd International Conference on Computer Research and Development. Huala Lumpur, Malaysia. 2010. pp. 666-670.         [ Links ]

15. G. Sun, W. Huo. Adaptive Fuzzy Predictive Control of Satellite Attitude Based on Hierarchical Fuzzy Systems. International Conference on Intelligent System Design and Engineering Application. Changsha, China. 2010. pp. 208-211.         [ Links ]

16. Y. Zhou, W. Huo. Quaternion-Based Direct Adaptive Fuzzy Predictive Control for Attitude Tracking of Satellites. IEEE International Conference on Intelligent Computing and Intelligent Systems. No. 2. Shanghai, China. 2009. pp. 510-516.         [ Links ]

17. C. Cheng, S. Shu. ''Application of Fuzzy Controllers for Spacecraft Attitude Control''. IEEE Transactions on Aerospace and Electronic Systems. Vol. 45. 2009. pp. 761-765.         [ Links ]

18. B. Rani, C. Gomathy, B. Sowmya, R. Narmadha. ''Attitude Determination and Control System of Sathyabamasat.'' Recent Advances in Space Technology Services and Climate Change. Proceedings of a meeting held 13-15 November. Chennai, India. 2010. pp. 326-329.         [ Links ]

19. C. Pukdeboon. Dynamic Output Feedback Sliding Mode Control for Spacecraft Attitude Manouvers. ECTI-CON 2011. 2011. pp. 545-548.         [ Links ]

20. C. Pukdeboon, A. Zinober. Optimal Sliding Mode Controllers for Attitude Tracking of Spacecraft. 18^th IEEE International Conference on Control Applications. St.Petersburg, Russia. 2009. pp. 1708-1713.         [ Links ]

21. C. Li, Y. Wang, L. Xu, Z. Zhang. Spacecraft Attitude Stabilization Using Optimal Sliding Mode Control. 3^rd International Symposium on Systems and Control in Aeronautics and Astronautics. Harbin, China. 2010. pp. 1085-1089.         [ Links ]

22. H. Shou, J. Sheu, J.Wang. Micro-Satellite Detumbling Mode Attitude Determination and Control: UKF Approach. 8^th IEEE International Conference on Control and Automation. Xiamen, China. 2010. pp. 673-678.         [ Links ]

23. N. Chaturvedi, A. Sanyal, N. McClamroch.''Rigid- Body Attitude Control Using Rotation Matrices for Continuous, Singularity-Free Control Laws.'' IEEE Control Systems Magazine. Vol 31. 2011. pp. 30-51.         [ Links ]