Optimal Cones Intersection Technique - LAIRS [PDF]

Optimal Cones Intersection Technique

Daniele Mortari ∗ , Puneet Singla ∗∗ Department of Aerospace Engineering, Texas A&M University, College Station, TX 77843-3141.

Abstract This paper presents an extension to the classical Cones Intersection Technique, which was, and is, exhaustively used to estimate the spin axis direction of spin stabilized spacecraft. The resulting Optimal Cones Intersection Technique is derived from a demonstrated co-planarity condition and does not present ambiguities or singularities. The proposed method allows an optimal estimation of any onboard direction in a closed-form by fully complying with Wahba’s optimality criterion of spacecraft attitude. Numerical tests show that the Optimal Cones Intersection Technique is faster than ESOQ-2, the presently fastest optimal attitude estimator, for two vector observation case. Key words: Spacecraft Attitude Estimation PACS: 07.87.+v, 95.40.+s

∗ Associate Professor, Corresponding Author, E-mail: [email protected] ∗∗Post Doctoral Research Associate, E-mail: [email protected] URLs: http://aero.tamu.edu/people/mortari/ (Daniele Mortari), http://people.tamu.edu/∼puneet/ (Puneet Singla).

Preprint submitted to ACTA Astronautica

27 April 2006

1

Introduction

Spacecraft attitude determination is the process of estimating the orientation of a spacecraft by taking on-board observations of other celestial bodies or reference points, with respect to a Body Reference Frame (BRF). Leonard Euler in 1775 has shown [1] that the configuration of a rigid body can be fully defined by locating a cartesian set of coordinate fixed in rigid body (the BRF) relative to the some inertial coordinate axes. Three parameters are needed to define the origin of the BRF and another three parameters are needed to specify the orientation of the BRF with respect to an external inertial frame. That is why, generally two or more co-ordinate systems are defined for the attitude determination process: one on the vehicle body, the BRF, and the second can be, for instance, the Earth-Centered Inertial (ECI) reference frame [2]. The spacecraft attitude estimates must be calculated quickly and continuously during the entire operational life of a mission as it is an important task of spacecraft navigation, dynamics, and control. Generally, the attitude estimation is accomplished in a precise manner using the n unit vectors bi observed in BRF and the corresponding directions ri evaluated in ECI. The projection of the BRF observations, bi onto the ECI is given by an orthogonal matrix (C) called the Direction Cosine Matrix (DCM), or the orientation matrix, or the attitude matrix. So the attitude estimation problem corresponds to finding the elements of unknown attitude matrix, C, given the knowledge of the same set of directions in BRF and ECI. Furthermore, the BRF observations are constructed from the measurements available from the attitude sensors and therefore, are susceptible to sensor noise. In order to achieve high 2

precision attitude determination, these sensor errors, which tend to introduce error in the attitude, must be accounted for in the final design. Generally, all these sensor errors are modeled by zero mean Gaussian white noise of known standard deviation, σi . Many algorithms have been presented in literature to optimally estimate the spacecraft attitude from vector observations. Most of the single-point algorithms are described in references [3] and fully comply with the Wahba’s [4] attitude optimality problem 1 , i.e., to find the orthogonal matrix C, with determinant det(C) = +1, that minimizes the loss function L(C) =

n 1X ai | bi − Cri |2 2 i=1

(1)

where, ai are non-negative weights, bi and ri are reference unit vectors evaluated in BRF and ECI, respectively. Generally, non-negative weights, ai , are normalized to unity but if the weights are chosen to be the inverse of measured error variance, i.e. ai = σi−2 , then it has been demonstrated [7] that the Wahba Problem becomes a Maximum Likelihood Problem. It can be easily proved [8] that if the attitude sensors are ideal (no measurement noise), then for n ≥ 2 the attitude determination problem has a unique solution and, if the solution to Wahba problem exists, then this unique solution is also the solution of Wahba problem. But in real life, the attitude estimation problem is not so simple, and the vector observations, generally body vector measurements, are subject to sensor errors. Due to errors in vec1

The word optimal attitude is usually associated with an attitude satisfying

Wahba’s definition of attitude optimality. To our knowledge, only the Euler-q [5] and the OLAE [6] algorithms are based on different optimality definitions.

3

tor measurements, we may not even guarantee the existence of the solution. In this case, we can determine the first approximation of the spacecraft attitude using some ideal deterministic methods like CONES [9,10], TRIAD [11], and EULER-2 [12]. For many applications, the estimation of the three-axis spacecraft attitude is unnecessary as the knowledge of only one direction is required. The historical example of this case can be a spinning spacecraft. Most of the times, it is not really necessary to know the three-axis attitude of a spin stabilized spacecraft, but only the knowledge of the orientation of the spin axis in the ECI frame is required. Therefore, a mathematical tool, capable of a direct estimation of the pointing direction r defined in the ECI corresponding to a given spacecraft direction b defined in the BRF, would be very useful especially if the evaluation is fast and optimally accomplished. However, the necessity of having such a mathematical tool is not a problem restricted to spin stabilized spacecraft, but it pertains to almost all satellites as most of the upper stages, which take care of the last delicate phase of orbital injection, are spin stabilized. It should be noticed that even though the spin axis direction has been estimated, the orientation of the BRF is still undefined as the orientation of the other two axes is still unknown by a rotation angle. This paper provides the solution to the above requirements by presenting the Optimal Cones Intersection Technique (OCIT) algorithm [13], which optimally estimates any on board direction when n = 2 observed vectors are available. The OCIT mathematical model represents the extension of the ideal deterministic Cones Intersection Technique [9] (well known and widely used for more than thirty years for evaluating the spin axis orientation of spinning spacecraft) to the optimal solution (satisfying the Wahba’s optimality crite4

rion) and without the sign ambiguity characterizing the classic non-optimal solution [10]. The structure of this paper is as follows: first an overview of the classical ideal Cones Intersection Technique is presented, followed by the Optimal Cones Intersection Technique to take care of Wahba’s optimality criteria and sign ambiguity resolution. Finally, the proposed algorithm is validated by numerical speed tests comparison with ESOQ-2 [14], the presently fastest optimal attitude estimator.

2

The Cones Intersection Technique

From an historical point of view, the Cones Intersection Technique was devised by Carl Grubin and presented in Ref. [9]. Since then, this technique, which is summarized in the next section, has been largely applied to estimate the spin axis direction of spin-stabilized spacecraft.

2.1 Spin z-axis estimation

Ideal sensors are characterized by the conditions σi = 0,

and

C ri = bi

(2)

where i = 1, · · · , n. Using the orthogonality property of C it follows that riT rj = bTi bj

(3) 5

where (i, j = 1, · · · , n). From this equation, it can be inferred that the internal angle structure of the observed directions bi is identical to the internal angle structure of the reference directions ri . This is certainly false in the real case when the observations are made using sensors with limited accuracy (σi > 0). However, it is still important to obtain the first estimates of the spacecraft attitude using algorithms which solve for the ideal case. Based on this property and when n = 2, it is possible to evaluate the spin axis direction as the intersections of two cones [9]. To this end, let us consider a spin stabilized spacecraft spinning about the z body axis (our unknown direction), and has two attitude sensors (e.g., sun and horizon sensors) on-board which provide, at time t, the two directions, b1 and b2 , respectively. Further, let r1 and r2 be their counterparts as defined in the ECI. The problem is, therefore, mathematically described by following three relationships:

z T r1 = eT3 b1 = b1 (3)

z T r2 = eT3 b2 = b2 (3)

and

zTz = 1

(4)

where, e3 = {0 0 1}T denotes the z direction expressed in the BRF. Provided that r1 and r2 are not parallel, the z-axis can always be expressed as

z = az r1 + bz r2 + cz r1 × r2

(5)

It should be noticed that the first two terms in the above expression denote the component of z direction in the plane identified by r1 and r2 and the second term denotes the out-of-plane component of z direction. Furthermore, substituting Eq. (5) into Eqs. (4) and carrying out the necessary algebra, the 6

following values for az , bz , and cz can be easily found     b1 (3) − b2 (3) r1T r2    a = z   1 − (r1T r2 )2      T

b2 (3) − b1 (3) r1 r2

bz =    1 − (r1T r2 )2    s     1 − a2z − b2z − 2az bz r1T r2   c = ±  z T 2

(6)

1 − (r1 r2 )

This solution provides (as the intersections of the two cones) two possible directions for the spin axis z because Eqs. (4) are the equations associated with the problem of finding the intersections of two cones, the first with axis r1 and aperture b1 (3) = cos ϑ1 , and the second with axis r2 and aperture b2 (3) = cos ϑ2 . These directions are symmetric with respect to the plane identified by the r1 and r2 directions. These two solutions are represented by the sign ambiguity for cz , which is a problem to be solved. From an historical point of view the cz sign ambiguity problem was solved by looking at the subsequent solutions of the problem. Since the spin-axis drifts slowly, one of the two solutions is almost constant while the other one changes. Therefore, by keeping track of the constant solution we can discriminate which of the two solutions actually represents the spin axis.

2.2 How to Solve for Sign Ambiguity

In this section, a procedure to choose the correct solution (cz sign) is presented. Pre-multiplying Eq. (5) by (r1 ×r2 )T the following relationship can be obtained (r1 × r2 )T z = cz (r1 × r2 )T (r1 × r2 ) = cz | r1 × r2 |2 7

(7)

which implies that sign{cz } = sign{(r1 × r2 )T z}

(8)

Now, the directions r1 , r2 , and z, identify a spatial structure whose internal angles (in the ideal case) are independent from the chosen reference frame. This implies that Eq. (7) can also be written in the BRF, as follows (r1 × r2 )T z = (b1 × b2 )T e3 = b1 (1) b2 (2) − b1 (2) b2 (1)

(9)

Hence, Eqs. (7), (8), and(9), provide us with the complete expression for cz v u u 1 − a2 − b2 − 2az bz r T r2 1 z z cz = sign{b1 (1) b2 (2) − b1 (2) b2 (1)} t T 2

1 − (r1 r2 )

(10)

a solution with no sign ambiguity for cz . By equating Eq. (7) with Eq. (9) we can directly compute cz cz =

b1 (1) b2 (2) − b1 (2) b2 (1) b1 (1) b2 (2) − b1 (2) b2 (1) = | r1 × r2 |2 1 − (r1T r2 )2

(11)

with no sign ambiguity. Finally, the closed form solution for spin axis z can be written as z=

b2 (3) − b1 (3) r1T r2 b1 (3) − b2 (3) r1T r2 r + r2 + 1 1 − (r1T r2 )2 1 − (r1T r2 )2

(12)

b1 (1) b2 (2) − b1 (2) b2 (1) + r1 × r2 1 − (r1T r2 )2 Note that this method cannot be applied when the vector cross product between r1 and r2 vanishes, which occurs only when r1T r2 = ±1, since r1 and r2 are unit-vectors. The Cones Intersections Technique, however, is not responsible for this problem. In fact, the condition r1T r2 = ±1 does not allow to solve 8

the attitude problem by this and by any other method, because under this condition the problem becomes intrinsically unresolvable. In the other words, in the attitude determination problems the condition r1T r2 = ±1 implies the knowledge of only one direction in space and, therefore, there is no information on a direction orthogonal to r1 (or r2 ), which implies that the condition r1T r2 = ±1 coincides with that of having only one observation (n = 1).

2.3

x-axis estimation

The method can be applied for any body axis determination, as for instance for the x-axis

x = ax r1 + bx r2 + cx r1 × r2

(13)

By analogous procedure, the following three conditions      T T    x r1 = e1 b1 = b1 (1)      

(14)

x r2 = e1 b2 = b2 (1)              xT (r1 × r2 ) = eT 1 (b1 × b2 ) T

T

allow us to obtain the closed form expression for the x-axis b2 (1) − b1 (1) r1T r2 b1 (1) − b2 (1) r1T r2 r1 + r2 + x= 1 − (r1T r2 )2 1 − (r1T r2 )2 +

b1 (2)b2 (3) − b1 (3)b2 (2) r1 × r2 1 − (r1T r2 )2 9

(15)

2.4 Three axis attitude determination

Once z and x have been evaluated, then the direction cosine matrix can be written as

C = [ x z × x z ]T

(16)

where

z × x = [(az cx − ax cz ) r1T r2 + (bz cx − bx cz )] r1 + +[(ax cz − az cx ) − (bz cx − bx cz ) r1T r2 ] r2 +

(17)

+(az bx − ax bz ) r1 × r2

2.5 Estimation of a general on-board direction

Note that the method can be applied not only to determine one of the spacecraft coordinate axes, but any on-board direction. For this case, let r be the unknown direction in ECI, and b be the same direction (known) in BRF. In this case the equations are

rT r1 = bT b1

rT r2 = bT b2

and

rT (r1 × r2 ) = bT (b1 × b2 )

(18)

where we look for a solution having the form

r = ar r1 + br r2 + cr r1 × r2

(19) 10

which implies the general solution r=

bT (b1 − b2 r1T r2 ) bT (b2 − b1 r1T r2 ) bT (b1 × b2 ) r + r + r1 × r2 1 2 1 − (r1T r2 )2 1 − (r1T r2 )2 1 − (r1T r2 )2

(20)

It is interesting to see that, by replacing the condition rT r = 1 (which leads to a solution ambiguity!) by the condition rT (r1 × r2 ) = bT (b1 × b2 ), not only the ambiguity is resolved, but also the constraint for r to be a unit-vector is guaranteed. This is easily demonstrated by writing the three conditions of Eq. (18) into matrix form rT [ r1

r2

r1 × r2 ] = bT [ b1

b2

b1 × b2 ]

(21)

now, the solution of the TRIAD [11] method C = [ b1

b2

b1 × b2 ] [ r1

r2

r1 × r2 ]−1

(22)

allows us to re-write Eq. (21) as rT = bT C, which implies r = CTb

(23)

and, consequently, implies that if b is a unit-vector, then also r is a unit vector, because orthogonal transformations preserve lengths of vectors and angles between vectors.

3

Optimal Cones Intersection Technique

In the ideal case, that is, when no errors affect the observed and the reference directions, there exists a rotation matrix C such that the n relationships given in Eq. (2) hold. In the real case, unfortunately, these relationships do not 11

hold anymore because the spatial distribution of the bi differs from the spatial distribution of the ri . Hence, the problem to compute the orthogonal matrix C, that best satisfies Eq. (1), arises. Grace Wahba, in 1965 2 , posed this problem as an optimal problem for attitude: to find the direction cosine matrix C that minimizes the cost function L(C) given in Eq. (1). In the case of n independent observations, it is easy to prove that L(C) can be written as L(C) =

n X i=1

αi −

n X

αi bTi Cri = λ0 − G(C) = min

(24)

i=1

from which it is clear that L(C) is minimized when G(C) is maximized. The gain function G(C) can be written as G(C) =

n X

αi bTi Cri = · · · = tr[CB T ] = max

(25)

i=1

where, B =

n X

αi ri bTi is known as the attitude profile matrix. The relative

i=1

weights αi can be chosen normalized, λ0 = λ0 =

n X 1 2 i=1 σi

n X

αi = 1, or un-normalized,

i=1 −2 = σtot , where the σi2 represent the sensor error variances.

For all the existing attitude estimation algorithms, accuracy and speed are two very important requirements. As for the attitude accuracy, it can be improved by increasing the precision of the observed directions (with better hardware and better data processing), and by increasing the precision in the reference directions (with more accurate codes). With these given accuracies, then there is no further room to improve the attitude accuracy since the Wahba optimality 2

At that time graduate student at Stanford, she proposed this problem while on

a summer job with IBM.

12

criterion yields for C an accuracy upper limit. Therefore, particular attention should be given to the development of faster attitude estimation techniques. A fast solution for the OCIT is obtained thanks to a coplanarity condition (among four unit-vectors) that is demonstrated in the next subsection.

3.1 Co-Planarity Condition

As already specified, in the real case when n = 2, the condition r1T r2 = bT1 b2 does not hold anymore. Therefore, it is not possible to devise an attitude matrix C that satisfies both Cr1 = b1 and Cr2 = b2 . However, associated with the optimal attitude C, we can introduce two new direction vectors, ¯b1 and ¯b2 as follows C r1 = ¯b1 6= b1 ,

and

C r2 = ¯b2 6= b2

(26)

which implies that r1T r2 = ¯bT1 ¯b2 , while the loss function, G(C) can be written as G(¯b1 , ¯b2 ) =

2 X i=1

αi bTi C ri =

2 X

αi bTi ¯bi =

i=1

2 X

αi cos ϑi = max

(27)

i=1

This implies that, if we are able to compute the directions ¯b1 and ¯b2 , then the attitude estimation problem would become as easy as in the ideal case because we could apply the same solving formulations devised for the ideal case, where b1 is replaced by ¯b1 , and b2 is replaced by ¯b2 . The problem, therefore, is recasted into the problem of evaluating the directions ¯b1 and ¯b2 . To this end, it is useful to demonstrate that, when n = 2 only, the plane identified by ¯b1 and ¯b2 coincides with that identified by b1 and b2 . This implies a coplanarity 13

condition among the four directions b1 , b2 , ¯b1 , and ¯b2 . To this goal, the Wahba’s cost function given by Eq. (27), which is a function of two unknowns ¯b1 , and ¯b2 , is subjected to three constraints ¯bT¯b1 = 1, 1

¯bT¯b2 = 1, 2

¯bT¯b2 = rT r2 = cos ϑr . 1 1

and

(28)

The problem to maximize G(¯b1 , ¯b2 ) subject to these constraints can be solved by using the Lagrangian multiplier technique. The augmented cost function G∗ (¯b1 , ¯b2 ) is given by the following equation

G∗ (¯b1 , ¯b2 ) = α1 bT1 ¯b1 + α2 bT2 ¯b2 − λ12 (¯bT1 ¯b2 − cos ϑr )+

(29)

−λ1 (¯bT1 ¯b1 − 1) − λ2 (¯bT2 ¯b2 − 1)

Now, the stationary conditions, α1 b1 = λ12 ¯b2 + 2λ1 ¯b1

∂G∗ (¯b1 , ¯b2 ) ∂G∗ (¯b1 , ¯b2 ) = 0 and = 0, yield ∂¯b1 ∂¯b2

and

α2 b2 = λ12 ¯b1 + 2λ2 ¯b2

(30)

which implies that the four unit-vectors, b1 , b2 , ¯b1 , and ¯b2 , are linearly correlated. As a consequence of this property, they must be co-planar (see Fig. 1).

Based on the co-planarity condition it is possible to evaluate the ¯b1 and ¯b2 directions. In fact, the optimality condition, as expressed by Eq. (27), implies that G(ϑ1 , ϑ2 ) = α1 cos ϑ1 + α2 cos ϑ2 = max 14

(31)

Fig. 1. Coplanarity condition geometry

Now, since ϑ2 = ϑb − ϑr − ϑ1 , Eq. (31) can be rewritten as G(ϑ1 ) = α1 cos ϑ1 + α2 cos(ϑb − ϑr − ϑ1 ) = max

therefore, stationary condition

(32)

dG(ϑ1 ) = 0 yields to the following expression dϑ1

α1 sin ϑ1 = α2 sin ϑ2

(33)

Using the fact that ϑ2 = ϑb − ϑr − ϑ1 , where cos ϑr = r1T r2 = ¯bT1 ¯b2 and cos ϑb = bT1 b2 , then Eq. (33) can be solved for the following expression for ϑ1     α2 sin(ϑb − ϑr )    sin ϑ1 = q   2 2  α1 + α2 + 2 α1 α2 cos(ϑb − ϑr )    α1 + α2 cos(ϑb − ϑr )    cos ϑ1 = q    α2 + α2 + 2 α α cos(ϑ − ϑ ) 1

2

1

2

b

(34)

r

This allows us to evaluate ϑ2 = ϑb − ϑr − ϑ1 and, then, the solution becomes     b2 sin ϑ1 + b1 sin(ϑr + ϑ2 )   ¯   b1 =

sin ϑb

   b2 sin(ϑr + ϑ1 ) + b1 sin ϑ2    b2 = ¯

sin ϑb

15

(35)

The directions ¯b1 and ¯b2 can also be computedby rigid rotation. In fact, the unit vector ¯b2 can be obtained by the rigid rotation ¯b2 = R(nb , ϑr ) ¯b1

(36)

where R(nb , ϑr ) = I cos ϑr + (1 − cos ϑr ) nb nTb + n ˜ b sin ϑr

is the matrix performing rigid rotation about the axis nb =

(37) b1 × b2 by an | b1 × b2 |

angle ϑr . Equation (36) allows us to write the Wahba’s cost function as G(¯b1 ) = α1 bT1 ¯b1 + α2 bT2 R(nb , ϑr ) ¯b1 = y T¯b1 = max

(38)

In order to maximize G(¯b1 ), the unit-vector ¯b1 must be parallel to y. Therefore T ¯b1 = y = α1 b1 + α2 R (nb , ϑr ) b2 |y| | α1 b1 + α2 RT (nb , ϑr ) b2 |

(39)

and then Eq. (36) can be used to compute ¯b2 . Once ¯b1 and ¯b2 are evaluated, then using either Eq. (35) or Eqs. (36,39), it is possible to carry out a fast and optimal evaluation of any spacecraft axis by using the formulation developed for the ideal case, described by Eq. (20), by replacing the observations, b1 and b2 , by the directions ¯b1 and ¯b2 , respectively. Note that, the capability to evaluate ¯b1 and ¯b2 , allows an immediate extension of all the ideal attitude determination techniques that use only n = 2 observed directions, to fully comply with Wahba’s optimality definition. In particular, it is possible to optimize EULER-2 and TRIAD algorithms [12]. 16

3.2 Optimal Cones Intersection Technique for n > 2

In this section, we present the extension of OCIT to when more than 2 observations are available. It has been demonstrated [15] that, for n = 2 case, the maximum attitude error 3 εij , which could affect the attitude estimation, can be well estimated from the knowledge of the sensor precision [ σi , σj ] and the referenced vector pair ri and rj , without any knowledge of the spacecraft attitude. A simplified expression for this error is given by following expression v u u sin2 (kσi ) + sin2 (kσj ) + 2 sin(kσi ) sin(kσj ) r T r2 1 sin εij = t ≈ εij T 2

1 − (r1 r2 )

(40)

where k can be 3 or 4 as kσi denote the maximum sensor error possible. It should be noticed that the angle kσi represents the precision of the observed direction bi and the angle εij represents the precision of the optimal direction εij rij , estimated by CONES using the ith and j th observed directions. Hence, k εij 2 4 denotes the scalar variance for the OCIT algorithm i.e. = E(δφ ) . Now, k the idea is to evaluate the optimal pointing direction, r, as a weighted sum of all the possible optimal directions rij obtained by CONES using all possible observed vector pairs [ bi , bj ]. r=

n−1 X

n X

i=1 j=i+1

3

ξij rij =

Nc X

ξ` r`

(41)

`=1

The maximum attitude error is defined as the maximum error experienced by

an onboard direction due to the fact that the attitude estimation (C) does not coincide with the true attitude (T ). This error coincides with the principal angle of the corrective attitude matrix ∆ = T C T . · ¸ T 4 δφ = cos−1 trace(T C ) − 1 2

17

n(n − 1) , which is the number of all 2 the possible different combinations of n objects taken two by two. The relative

where the index `, ranges from 1 to Nc =

weights ξ` , appearing in Eq. (41), are computed from the maximum attitude errors εij (as the αi are computed from the kσi ) as follows 1 ξ` = µ ¶2 ε` ε¯ k

where

ε¯ =

Nc X k2 ij

ε2ij

(42)

The optimal solution, ropt , is then obtained by normalizing r

ropt =

r |r|

(43)

It is easy to check that ropt direction is optimal when n = 2 because it coincides with the CONES solution, while, for n > 2, this solution minimizes the loss function

L(r) =

Nc Nc X 1X ξ` ϑ2` ' 1 − ξ` r`T r 2 `=1 `=1

(44)

The loss function, L(r) can be interpreted as the energy function for a set of spherical springs of stiffness ξ` linking the solution r with the estimated values ϑ2 r` . Further, r`T r = cos ϑ` and ϑ` is so small that cos ϑ` ≈ 1 − ` . 2 Now, the loss function L(r) can be minimized by means of the Lagrangian multiplier technique. To this end, the augmented cost function, which includes the constraint for v to be a unit vector, can be written as 1 L∗ (r) = L(r) − λ(rT r − 1) 2

(45)

18

where λ is the Lagrangian multiplier. The stationarity condition implies Nc dL∗ (r) X = ξ` r` − λr = 0 dr `=1

(46)

which, without the need of computing the value of the Lagrangian multiplier, demonstrates that the solution r is the direction of the r vector provided by n(n − 1) Eq. (41). Obviously, the number Nc = of all possible combinations of 2 observed vector pairs becomes rapidly large as n increases. Therefore, when the maximization of the achieved accuracy is not required, the choice to estimate the pointing direction v by using a small subset of m < n estimations r` , is suggested. The m solutions r` are, obviously, those providing the minimum εij values.

4

Numerical Tests

The attitude estimation algorithm, OCIT, described in this paper is numerically tested for both accuracy and speed. This analysis is performed using MATLAB [16] on a 1.7 GHz Pentium PC. For the analysis presented below, the vector observation data set consisting of the true observed direction b, the inertial frame directions r and the true attitude matrix T , is generated randomly. The attitude sensor errors are modeled by zero mean white Gaussian noise process with known variance, (σi ). The measured directions bi are then obtained as bi = R(hi , kσi ) T ri , where R(hi , kσi ) performs the rigid rotation about the random axis hi (orthogonal to ri ) by the random angle kσi . For the accuracy test the number of vector observations, n, is varied from 2 to 10 and the attitude sensor accuracy, σ, is varied from 50 µrad to 0.05 rad 19

and thus covering the practical accuracy values for various existing attitude sensors. For one particular value of n and σ the attitude accuracy is computed using the following formula "

εmax = cos

−1

trace(T C T ) − 1 2

#

(47)

where, C is the estimated attitude matrix and T is the corresponding true attitude matrix. Figure 2 shows the plot of mean attitude accuracy result for 10, 000 Monte Carlo simulations. The solid line in plot represents the attitude accuracy for Wahba compliant ESOQ-2 and the solid dots “•” denote the attitude accuracy obtained by the OCIT algorithm as presented in this paper. The empty circular dots “◦” denote the minimum attitude accuracy obtained by considering two vector observations at a time. As expected, for n = 2 case, there is no differences between the attitude accuracy for the ESOQ-2 and OCIT algorithms as both of them fully complies with the Wahba optimality definition. But for n > 2 case, the attitude accuracy obtained by the OCIT is bit degraded from the ESOQ-2 accuracy as in this case the OCIT attitude is the weighted average of the Wahba complaint attitude obtained by considering two vector observations at a time. But it should be noticed that for all practical purposes this difference in accuracies is negligible. It is also clear from Figure 2 that the attitude accuracy obtained using the best pair of vector observations (i.e. for which εij is minimum) is better than the ESOQ-2 attitude accuracy when all the available vector observations are used at once. This result supports our earlier claim that for n > 2 case the maximum attitude accuracy can be obtained by using the m (< n) set of vector observations for which εij is minimum. 20

−1

10

5.0e−002 −2

10

Mean Attitude Accuracy

5.0e−003 −3

10

5.0e−004

−4

10

5.0e−005 −5

10

ESOQ2 OCIT Min Pairwise OCIT Accuracy

−6

10

2

3

4

5 6 7 8 Number of Vector Observations

9

10

Fig. 2. Attitude Accuracy Plot

0.2507 ms

ESOQ2

0.0960 ms

TRIAD

OCIT

0.08 ms

0

0.05

0.1

0.15 Time (ms)

0.2

0.25

0.3

Fig. 3. Time Comparison for OCIT

For computational speed test, the n = 2 case is considered and the OCIT algorithm is compared with the ESOQ-2 (presently fastest attitude estimation algorithm) and the TRIAD algorithms using the MATLAB function, cputime. 21

For speed test, the TRIAD algorithm is also modified by using the procedure listed in section 3.1 to make it fully complying with Wahba optimality condition. Figure 3 shows the CPU time bar chart for the ESOQ-2, TRIAD and OCIT algorithms for n = 2 case and 10, 000 Monte-Carlo simulations. From this Figure, it is clear that the OCIT and the TRIAD algorithms are more than 50% faster than the ESOQ-2 even for computing 3-axis attitude. It should be also noticed that this speed gain result is independent from the sensor precision (σi ). However, we should mention that the CPU time requirement for the ESOQ-2 generally do not changes with n but for the OCIT and the TRIAD increases with O(n2 ).

5

Conclusions

The Wahba optimal condition compliant extension to the classical Cones Intersection Technique is presented to estimate any generalized on-board direction or 3-axis spacecraft attitude. The resulting Optimal Cones Intersection Technique (OCIT) has a mathematical model without any sign ambiguities or singularities, and allows an optimal estimation of any onboard direction in a closed-form by fully complying with Wahba’s optimality criterion of spacecraft attitude. The advantage of the OCIT in terms of accuracy and computational speed, over the ESOQ-2, are demonstrated with the help of Monte Carlo simulations. The simulation results presented in this paper form the basis for considering the OCIT as one of the most suitable algorithm when a fast and optimal estimation of a pointing direction is required. 22

References

[1] L. Euler, Formulae Generales pro Translatione Quacunque Corporum Rigidorum, Novi Commentarii Academiae Scientiarum Petropolitanae 20 (1776) 189–207. [2] D. A. Vallado, Fundamentals of Astrodynamics and Applications, Vol. 2, McGraw–Hill, New York, 2001. [3] F. L. Markley, D. Mortari, Quaternion Attitude Estimation using Vector Observations, Journal of the Astronautical Sciences 48 (2/3) (2000) 359–380, special Issue: The Richard H. Battin Astrodynamics Symposium. [4] G. Wahba, Problem 65–1: A Least Squares Estimate of Spacecraft Attitude, SIAM Review 7 (3) (1965) 409. [5] D. Mortari, Euler–q Algorithm for Attitude Determination from Vector Observations, Journal of Guidance, Control, and Dynamics 21 (2) (1998) 328– 334. [6] D. Mortari, F. L. Markley, J. L. Junkins, Optimal Linear Attitude Estimator, in: Advances in the Astronautical Sciences, Vol. 105, Pt. I, 2000, pp. 465–478. [7] M. D. Shuster, Maximum Likelihood Estimate of Spacecraft Attitude, Journal of the Astronautical Sciences 37 (1) (1989) 79–88. [8] P. Singla, D. Mortari, J. L. Junkins, How to Avoid Singularity for Euler Angle Set?, Paper AAS 04-190 of the 2004 Space Flight Mechanics Meeting Conference, Maui, Hawaii (February 9–13, 2004). [9] C. Grubin, Simple Algorithm for Intersecting two Conical Surfaces, Journal of Spacecraft 14 (1977) 251–252. [10] J. R. Wertz, Spacecraft Attitude Determination and Control, D. Reidel Publishing Company, Boston, Ma, 1978.

23

[11] H. Black, A Passive System for Determining the Attitude of a Satellite, AIAA Journal 2 (1964) 1350–1351. [12] D. Mortari, EULER–2 and EULER–n Algorithms for Attitude Determination from Vector Observations, Space Technology 16 (5/6) (1996) 317–321. [13] D. Mortari, Optimal Cones Intersection Techniques for Attitude Pointing Error Evaluation, in: Advances in the Astronautical Sciences, Vol. 97, Pt. I, 1997, pp. 949–962. [14] D. Mortari, Second Estimator of the Optimal Quaternion, Journal of Guidance, Control, and Dynamics 23 (5) (2000) 885–888. [15] D. Mortari, The Attitude Error Estimator, in: International Conference on Dynamics and Control of Systems and Structures in Space 2002, King College, Cambridge, England, 2002. [16] Matlab reference guide, The MATH WORKs Inc. (October 1992).

24