A comparative study of spacecraft attitude determination and estimation algorithms (a cost–benefit approach)

Aerospace Science and Technology 26 (2013) 211–215

Contents lists available at SciVerse ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

A comparative study of spacecraft attitude determination and estimation algorithms (a cost–benefit approach) Tamer Mekky Ahmed Habib ∗ Egyptian National Authority for Remote Sensing and Space Science, Cairo, Egypt

a r t i c l e

i n f o

Article history: Received 28 September 2011 Received in revised form 8 April 2012 Accepted 19 April 2012 Available online 4 May 2012 Keywords: Three-axis Attitude Determination Estimation Magnetometer Gyro Disturbance Extended Kalman filter Cost–benefit

a b s t r a c t The primary goal of this work is to investigate various performance indices of spacecraft attitude determination and estimation algorithms to perform a cost–benefit analysis. This analysis could help the designer to choose whether to use an attitude determination algorithm or an attitude estimation algorithm. The choice depends mainly on the available computational resources and the required attitude determination accuracy. These performance indices include accuracy, and execution time. In order to determine these performance indices, a test-case spacecraft is utilized. The test-case spacecraft is subject to aerodynamics forces and moments, solar radiation pressure forces and moments, earth’s oblateness till the fourth order (i.e. J 4 ), residual magnetic dipole moments, and gravity gradient moments. The available spacecraft attitude sensors are magnetometer and gyro. © 2012 Elsevier Masson SAS. All rights reserved.

1. Introduction Magnetometer and gyro are commonly used sensors for spacecraft attitude estimation and attitude determination. Attitude determination as described in [15], is the process of extracting attitude parameters from measurements of at least two kinds of physical quantities. Sometimes this process is referred to as deterministic method or point-to-point approach [8]. Attitude estimation is defined also in [15], as the process of extracting attitude parameters from measurements of a single physical quantity or more. During presentation of the work given in [14], a question arose, why to use attitude estimation algorithms while it is possible to use attitude determination algorithms. Attitude determination algorithms (such as TRIAD, Q-method, etc.) have a lot of advantages such as low computational efforts, and independence of spacecraft attitude motion models. On the other hand, they couldn’t deal with sensor measurements if the number of independent sensors was less than two sensors. Attitude estimation algorithms (such as Kalman filter, and its derivatives) have the advantage of being able to deal with even measurements of a single sensor. These algorithms usually make use of modeling the spacecraft attitude motion. Conversely, attitude estimation algorithms are usually characterized by enormous computational effort. Dur-

*

Tel.: +20 1224207654. E-mail address: [email protected].

1270-9638/$ – see front matter © 2012 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ast.2012.04.005

ing some spacecraft operational modes such as the stand-by mode, only measurements of the magnetometer are available. If this is the case, an attitude estimation algorithm must be used. In the high accuracy spacecraft operational modes such as the imaging mode, multi-sensor measurements are available. In this situation, either attitude estimation algorithms or attitude determination algorithms could be used and the need arise to determine which kind of algorithm that should be used. The designer of the software related to the attitude determination and control subsystem (ADCS) usually sticks to attitude estimation algorithms and disregards totally the attitude determination algorithms in spite of their advantages. Such situation is encountered during the design of the first Egyptian scientific spacecraft which was launched in 2007, EGYPTSAT-1. The literature doesn’t provide an explanation to the answer of the questions: when should we use attitude determination algorithms and when should we use attitude estimation algorithms? Refs. [8,12,7], discussed both of attitude determination and estimation algorithms but with no answer to such questions. Refs. [13] and [4] compare only attitude estimation algorithms and don’t make use of attitude determination methods at all. Ref. [5] uses small angles approximation, and hence isn’t valid for large initial attitude estimation error. Ref. [6] suffers from large tracking error standard deviations that reach up to 11 degrees in the yaw angle. In Ref. [11] the attitude estimates converges after a long period typically, one orbit, while as the presented extended Kalman filter converges after about 0.25 orbits. Ref. [3] uses the unscented Kalman filter which

212

T.M.A. Habib / Aerospace Science and Technology 26 (2013) 211–215

is characterized by enormous computational load more than the extended Kalman filter. The contribution of the research at hand comes from addressing a problem that hasn’t been addressed before based on quantitative measures. The current attitude estimation and determination algorithms reported in the literature don’t answer the question: what is the cost (in terms of computational load) of a certain attitude determination algorithm with respect to attitude estimation algorithms and what is the gain achieved? And how to help the software designer of the ADCS to determine an affordable (in terms of computational load versus accuracy) method to extract spacecraft attitude angles based on magnetometer and gyro measurements. The main objective of this research is to provide a detailed cost– benefit analysis for attitude estimation and attitude determination algorithms with no prior information and no small angle approximations taking into account many disturbances such as aerodynamics forces and moments, solar radiation pressure forces and moments, earth’s oblateness till the fourth order (i.e. J 4 ), residual magnetic dipole moments, and gravity gradient moments. This analysis is applied to large initial attitude estimation errors (typically about 180◦ for the yaw angle, 175◦ for the roll and 85◦ for the pitch angle) so as to assure the convergence of the estimation algorithms from lost in space initial conditions. The sensors utilized in the research at hand are the same like those of [14], but with exclusion of the GPS receiver and the star sensor. 2. Attitude determination algorithms Several attitude determination algorithms exist in the literature such as the algebraic method (known also as TRIAD) [15], the modified algebraic method [12], the optimized TRIAD algorithm [1], and the Q-method [2]. The general attitude determination problem is defined as: given n abstract unit vectors resolved in a body and reference Cartesian coordinate system (in the body coordinate system these vectors are denoted by b i , i = 1, 2, 3, . . . , n these vectors are resolved into a reference axes as r i , i = 1, 2, 3, . . . , n); find an orthogonal 3 × 3 transformation matrix T R → B that minimizes the cost function

J1 = c



ai |b i − T R → B r i |2

(1)

where c is a constant, ai is the weight assigned to the ith pair. Each attitude determination algorithm solves this problem in a different way. The following well known transformations are needed. The first transformation is given by

B B = T I →B B I

(2)

where B B and B I are the body and inertial components of the earth’s magnetic field defined in the body frame of reference and the inertial frame of reference respectively. T I → B is the transformation matrix from the inertial to body axes. The second transformation is given by

T I → B B˚ I = B˚ B + ω × B B

B˚ (•) j +1 − B˚ (•) j

T

(4)

where  T is the time step between successive measurements, and j is a counter representing the time step number. Thus, the vector pairs are [9]

( B I , B B ) and ( B˚ I , B˚ B + ω × B B )

The algorithm starts by computing the following vectors

E 1 = b1 /|b1 | E2 =

(5)

So, if we multiplied B I by T I → B we get B B , and similarly if we multiplied B˚ I by T I → B we get B˚ B + ω × B B . Consequently, the attitude matrix, T I → B , could be constructed from these pairs.

(6)

E 1 × b2

(7)

| E 1 × b2 | E3 = E1 × E2

(8)

and their corresponding vectors in the reference axes

G 1 = r1 /|r1 | G 1 × r2 G2 = |G 1 × r2 |

(9) (10)

G3 = G1 × G2

(11)

The direction cosine matrix that transforms a vector from the reference frame to the body frame is given by:

T R → B = E 1 .G 1T + E 2 .G 2T + E 3 .G 3T

(12)

T

where (•) denotes the transpose. Eqs. (6) to (12) are applicable only if b1 and b2 are not parallel, and so are r1 and r2 , i.e.,

|b1 .b2 | < 1,

|r1 .r2 | < 1

(13)

The body matrix M b , and the reference matrix M r , could now be formed as,

Mb = [ E 1

E2

E3 ] ,

Mr = [ G 1

G2

G3 ]

(14)

Thus, the transformation matrix T R → B , is related to the body and the reference matrices by the relation

T R →B Mr = Mb

(15)

Consequently,

T R → B = M b M r−1

(16)

Since the matrix M r is orthogonal, we could write Eq. (16) as

T R → B = M b M rT

(17)

Thus, the only condition that must hold is that the matrix M r , has an inverse, which is guaranteed by the condition given in the inequality (13). The inequality, (13), assures that the rows of the M r matrix are linearly independent. 2.2. Q-method In this method, the cost function given in Eq. (1) is applied as [2]

(3)

where ω is the inertial angular velocity of the spacecraft body, and B˚ (•) denotes the first derivative of B (•) w.r.t. time computed from

B˚ (•) j +1 =

2.1. TRIAD algorithm

J2 =

1 2

ai |b i − T R → B r i |2

(18)

In Eq. (18), the need arise to solve for the quaternion vector, (q), rather than the attitude matrix, T R → B . Thus Eq. (18) is rewritten as

J2 =

2 1  ai b i − T R → B (q)r i  2

(19)

We are seeking for a unit length the quaternion that minimizes the cost function given in Eq. (19). The required quaternion vector is the eigenvector of a certain matrix, K , corresponding to the largest eigenvalue. Construction of the K matrix starts by defining σ , B, S, and z, as follows [2]

T.M.A. Habib / Aerospace Science and Technology 26 (2013) 211–215

σ=

n 

ai b iT r i

(20)

ai b i r iT

(21)

i =1 n 

B=

z=

n  i =1



K=

S − σ I3

z

zT

σ

Xˆ k− : A priori state estimate at a time step k. Ak :

The state transition matrix. The measurement matrix.

Xˆ k :

A posteriori state estimate at a time step k.

Pk:

A posteriori estimate error covariance at a time step k.

(23)

Q k:

The discrete process noise covariance.

Rk :

The discrete measurement noise covariance.

(24)

zk :

The measurement vector provided by measurement devices.

zˆ k :

The estimated measurement vector.

(22)

ai b i × r i

with

Hk :

i =1

S = B + BT

213



with I 3 , defined as a third order identity matrix. The eigenvalues of the K matrix are computed and the eigenvector corresponding to the largest eigenvalue is thought to be the quaternion vector. A property of the quaternion vector is that

The state transition matrix is given by

q = −q

with  T defined as the sampling time interval and

(25)

Thus, during computation of the required quaternion, one could obtain the negative quaternion. This could be explained by the noting the relation between the quaternion and the Euler axis of rotation [15]. Any general rotations about three axes could be viewed as a single rotation by an angle, Φ , about the axis of the resultant rotation (i.e., Euler axis, eˆ = [e 1 e 2 e 3 ] T ). The relation between the quaternion vector components and the Euler axis and angle is given as

q1 = e 1 sin

Φ

q2 = e 2 sin

Φ

q3 = e 3 sin

Φ

q4 = cos

(26)

2

(27)

2

(28)

2

Φ

(29)

2

From the last four equations it is evident that the rotation about the Euler axis, eˆ , by an angle, Φ , is equivalent to the rotation about the axis, −ˆe , by an angle, −Φ . Thus, if the scalar part of the corresponding quaternion resulting from the Q-method is negative, the required quaternion vector sign should be changed. 3. Attitude estimation algorithm (the extended Kalman filter) Attitude estimation algorithms such as the family of Kalman filter are used to minimize a certain cost function other than that specified in Eq. (1). The cost function of the Kalman filter family is given by





J 2 = E ( Xˆ k − X k ) T ( Xˆ k − X k )

(30)

where E is the expectation, X k is the true plant state vector, and Xˆ k , is the estimated state vector. The structure of the extended Kalman filter given in [14] is

Xˆ k− = f ( Xˆ k−1 ) −

Pk =

(31)

A k−1 P k−1 A kT−1

 T

K k = P k− H k

+ Q k −1 − 1 H k P k− H kT + R k

(32) (33)

Xˆ k = X k + K k [ zk − zˆ k ]

ˆ−

−

(34) T

P k = [I − Kk Hk ] P k [I − Kk Hk ] +

K k R k K kT

(35)













A k−1 Xˆ k+−1 = I + F k−1 Xˆ k+−1  T





F k−1 Xˆ k+−1 = ∂ f /∂ X |( Xˆ + ) k −1



(36)

(37)

The measurement matrix H k is related to the measurement vector, h, by the relation

Hk =

 ∂ h( X )  ∂ X  X = Xˆ

(38)

The same model of [14] for the simulated and test-case spacecraft is utilized with some modifications. The same notations of [14] are used. The state vector is reduced to be X = [q1 q2 q3 q4 ωx ω y ωz ]T . The measurement vector is cut to be

h = [ b xb

b yb

b zb

ωx ω y ω z ]

(39)

The test-case spacecraft is subject to many disturbances such as solar radiation pressure forces and moments, aerodynamics forces and moments, earth’s oblateness till the fourth order (i.e. J 4 ), gravity gradient moments, and residual magnetic dipole moments. For further details the reader should return to [14]. 4. Test-case parameters, and results The initial orbital parameters of the test-case spacecraft and its simulated counterpart are selected to be identical to those found in [14]. These parameters are based on design simulations of the EGYPTSAT-1 spacecraft which was launched in 2007. The initial attitude angles are chosen very large to indicate lost in space conditions, just to assure the ability of the algorithms to converge despite of the absence of prior attitude knowledge. The very large initial attitude angles assure that no small angles approximations exist in the methodologies presented. The initial orbital parameters of the test-case spacecraft are a (semimajor axis) = 7 139 200 m, e (orbit eccentricity) = 0, i (orbit inclination) = 101.085◦ , Ω (right ascension of ascending node) = 339.5◦ , ω A (argument of perigee) = 69◦ , υ (true anomaly) = 6◦ , φ (roll angle) = −165◦ , ψ (yaw angle) = 170◦ , and θ (pitch angle) = 85◦ . The estimated spacecraft parameters are initialized with a (semi-major axis) = 7 039 200 m, e (orbit eccentricity) = 0, i (orbit inclination) = 98.85◦ , Ω (right ascension of ascending node) = 337.5◦ , ω A = 69◦ , and υ = 0◦ , φ (roll angle) = 0◦ , ψ (yaw angle) = 0◦ , and θ (pitch angle) = 0◦ . Epoch time: 1/4/2013 0h:0m:0s. Time step ( T ) = 4 seconds. The measurement, and process, noise covariance matrices are those given in [14] after excluding those related to omitted measurements and states. The same also applies for the initial estimation error covariance matrices.

214

T.M.A. Habib / Aerospace Science and Technology 26 (2013) 211–215

Fig. 1. Error of the attitude determination and estimation algorithms.

Table 1 Comparison between accuracy and execution time among different algorithms.

TRIAD algorithm Q-method Extended Kalman filter

Yaw angle error standard deviation (deg.)

Pitch angle error standard deviation (deg.)

Roll angle error standard deviation (deg.)

Average execution time (s)

1.2517523 1.2517516 0.9241491

0.8943686 0.8943650 0.8725242

5.7003591 5.7703597 1.5811642

0.0009 0.0010 0.5303

The estimation error of attitude angles is shown in Fig. 1. The convergence criterion is chosen according to the real design process of EGYPTSAT-1 spacecraft. The most critical phases of the spacecraft life time were the detumbling and attitude acquisition phases respectively, and the attitude determination (and/or estimation) algorithm must be able to deal with these phases which are considered much more like lost in space conditions. An attitude estimation (or determination) algorithm is considered to be successful (i.e. convergent) if it could deal with no prior knowledge of the attitude angles and its error falls between ±7 degrees. As seen in Fig. 1, the estimation error is converging (according to the converging criterion described) within about 10 seconds despite large initial attitude estimation error. The attitude determination error of the TRIAD and the Q-method is shown also in the same figures. We should note that the error for the TRIAD and the Q-method is superimposed. If the convergence criterion is changed the convergence time of the estimator will be affected. In addition, a better tuning of the extended Kalman filter should be done throughout real coded genetic algorithm or the simplex algorithm to obtain the optimum values of either, the attitude estimation parameters, or the parameters of the attitude determination algorithms, such a methodology is given in [10]. Table 1 shows a comparison between accuracy and the average execution time among different algorithms. The simulations of the test-case spacecraft lasted for 1500 s approximately. The standard deviation of the errors of the attitude angles is computed for the time period 870–1500 s. The average execution time along the 1500 s time span is given in the fourth column of Table 1. In this table, the errors of the Q-method are very close to those of the TRIAD algorithm. So the number of significant digits is chosen just to indicate the first number that differs between the superimposed results. The TRIAD algorithm has the least execution time with the worst accuracy compared to the extended Kalman filter. The computational time of the Q-method is more than that of the TRIAD and less than that of the extended Kalman filter with no significant enhancement of the spacecraft attitude error standard deviation. The main advantage of the Q-method is its ability to

deal with measurements coming from more than two sensors. The extended Kalman filter shows a very good attitude error standard deviation (actually four times better than the TRIAD algorithm); on the other hand it has a huge computational load (approximately six times the TRIAD algorithm execution time). The high accuracy resulting from attitude estimation algorithms is due to the usage of the satellite angular motion model (dynamics and kinematics of the attitude motion). While as attitude determination methods don’t use these models at all. Using these models makes attitude estimates more accurate and more time consuming than attitude determination algorithms. The practical significance for achieving better accuracy for the attitude estimation methods described is the necessity of achieving high accuracy attitude angles knowledge required in the spacecraft imaging mode or before execution of orbital maneuver. 5. Conclusion Attitude estimation algorithms are characterized by high accuracy and high computational load with respect to attitude determination algorithms. This is due to its usage of spacecraft attitude motion models in addition to the computational load required by the rest of the extended Kalman filter structure. Attitude determination algorithms on the other hand don’t use these models at all. This behavior increases the computational speed considerably on the cost of accuracy. In addition, simplicity of the attitude determination algorithms makes them favorable for satellite on board software. The questions that haven’t been addressed before in the literature: when to use attitude determination methods and when to use attitude determination methods is easily answered. The answer depends mainly on the available computational resources and the required attitude determination accuracy. As shown in the testcase results, attitude determination algorithms such as the TRIAD algorithms and the Q-method are characterized by low computation effort. On the other hand, they provide low accuracies of spacecraft attitude angles. Extended Kalman filter make use of modeling the spacecraft to increase such accuracies (actually four

T.M.A. Habib / Aerospace Science and Technology 26 (2013) 211–215

times better than the TRIAD algorithm). Conversely, it needs a huge computational load (approximately six times the TRIAD algorithm execution time). The computational time of the Q-method is more than that of the TRIAD with no significant enhancement of the spacecraft attitude error standard deviation. The main advantage of the Q-method is its ability to deal with measurements coming from more than two sensors. Thus the cost and benefit of each attitude knowledge methodology is now clear for the software designer of the attitude determination and control subsystem. This could be of high crucial importance if the available computational resources are low. References [1] Itzhack Y. Bar-Itzhack, Optimized TRIAD algorithm for attitude determination, Journal of Guidance Control and Dynamics 20 (1) (1996) 208–211. [2] Itzhack Y. Bar-Itzhack, New method for extracting the quaternion from a rotation matrix, Journal of Guidance Control and Dynamics 23 (6) (2000) 1085– 1087. [3] J.L. Crassidis, F.L. Markley, Unscented filtering for spacecraft attitude estimation, Journal of Guidance, Control and Dynamics 26 (4) (2003) 536–542, http://www.acsu.buffalo.edu/~johnc/uf_att.pdf. [4] O. Diaz, Analysis and comparison of extended and unscented Kalman filtering methods for spacecraft attitude determination, MSc thesis, Naval Postgraduate School, 2010. [5] M. Francois, P.K. Parimal, M.K. Psiaki, Three axis attitude determination via Kalman filtering of magnetometer data, in: Flight Mechanics/Estimation Theory Symposium, Goddard Space Flight Centre, 1988, pp. 344–367.

215

[6] D. Gebre-Egziabher, G.H. Elkaim, J.D. Powell, B.W. Parkinson, A gyro-free quaternion-based attitude determination system suitable for implementation using low cost sensors, in: Position Location and Navigation Symposium, IEEE, 2000. [7] L. Kuylen, F. Simsky, Attitude determination methods used in the PolaRx2@ multi-antenna GPS receiver, ION GNSS, 2005. [8] S. Marques, Small satellites attitude determination methods, MSc thesis, Instituto Superior Técnico, Universidade Técnica de Lisboa, 2000. [9] G.A. Natanson, M.S. Challa, J. Deutschmann, D.F. Baker, Magnetometer-only attitude and rate determination for a gyroless spacecraft, in: Proceedings of the Third International Symposium on Space Mission Operations and Ground Data Systems, in: NASA Conference Publications, vol. 3281, NASA, Greenbelt, MD, November 1994, pp. 791–798. [10] P. Powell, Automated tuning of an extended Kalman filter using the downhill simplex algorithm, Journal of Guidance, Control and Dynamics 25 (5) (2002) 901–908. [11] M.L. Psiaki, Attitude-determination filtering via extended quaternion estimation, Journal of Guidance, Control and Dynamics 23 (2) (2000) 206–214. [12] M. Tamer, The global positioning system application to satellite position and attitude determination, MSc thesis, Aerospace Department, Cairo University, 2003. [13] M. Tamer, New algorithms of nonlinear spacecraft attitude control via attitude, angular velocity, and orbit estimation based on the Earth’s magnetic field, PhD thesis, Cairo University, 2009. [14] M. Tamer, Fast converging with high accuracy estimates of satellite attitude and orbit based on magnetometer augmented with gyro, star sensor and GPS via extended Kalman filter, The Egyptian Journal of Remote Sensing and Space Sciences (2011), http://dx.doi.org/10.1016/j.ejrs.2011.06.002. [15] J.R. Wertz, Spacecraft Attitude Determination and Control, D. Reidel Publishing Company, 1997.