Robust Static Attitude Determination via Robust Optimization

Robust Static Attitude Determination via Robust Optimization

Proceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011 Robust Static A...

446KB Sizes 0 Downloads 111 Views

Proceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011

Robust Static Attitude Determination via Robust Optimization Shakil Ahmed ∗ Eric C. Kerrigan ∗∗ ∗

Department of Electrical and Electronic Engineering, Imperial College London, UK (e-mail: [email protected]). ∗∗ Department of Aeronautics and Department of Electrical and Electronic Engineering, Imperial College London, UK (e-mail: [email protected]).

Abstract: We address the problem of robust attitude determination using a static approach. In contrast to the dynamic approach, a static approach does not depend on the system dynamics. This approach only requires measurements of some vectors, such as the earth magnetic field, sun vector, etc in two different coordinate frames. These vectors, obtained from some sensor or a mathematical model, may not be accurate due to sensor errors and modeling inaccuracies. We consider all such errors as infinity-norm bounded uncertainties and our main focus is to obtain an attitude estimate, which is least sensitive to such uncertainties. We formulate a robust optimization (RO) problem with a quadratic cost and nonlinear constraints and propose a solution using a quaternion approach with an affine uncertainty parameterization. We transform the RO problem into a suboptimal minimization problem, which is non-convex but can be solved with a good initial guess using a nonlinear optimization solver. The results show a significant advantage of the proposed approach for the worst case uncertainties. Keywords: Transformation matrices, attitude algorithms, modeling errors, least square estimation, robust estimation 1. INTRODUCTION The main focus of this work is on robust attitude determination in the presence of uncertainties using a static approach, which is based on vector information. One main advantage of this approach is that it always gives an attitude estimate independent of the dynamic system nonlinearities, making it especially useful for highly nonlinear systems, such as a tumbling satellite, where dynamic approaches based on linear or nonlinear filtering suffer from divergence issues due to lack of good a priori state estimate (Crassidis, Markley, and Cheng, 2007). An added benefit of static estimation in such cases is that it can provide a good state initialization for dynamic filters, reducing the likelihood of divergence. To compute the attitude of an object i.e. its orientation in space with respect to some known reference, two coordinate frames are needed. One, which is fixed to the body of the object, is called the body frame, while the second is called the reference frame. Selection of the reference frame normally depends upon the control system requirements. Formally, attitude of an object is defined as a coordinate transformation that transforms reference coordinates into the body coordinates (Sidi, 2000). This transformation is obtained through a proper orthogonal transformation matrix C ∈ R3×3 , also known as a direction cosine or attitude matrix. Orthogonality of the transformation matrix dictates the constraint C T C = I, while the proper transformation matrix is defined by having the determinant equal to +1 to preserve the orientation in a rotation, thus imposing the constraint det(C) = +1. 978-3-902661-93-7/11/$20.00 © 2011 IFAC

The vector quantities normally used in static attitude determination are the earth magnetic field (EMF), sun vector, star vector, etc. Information of these vectors is required both in the body and reference frames in order to determine the attitude or transformation matrix. Normally the vectors in the body frame are measured by some sensor installed on the body of the object, while the same vector information in the reference frame is obtained from some mathematical model. With only one pair of information, i.e. a vector quantity in the body and reference frame, determination of the attitude matrix remains an underdetermined problem with infinite solutions. Adding a second pair makes this problem overdetermined (Shuster and Oh, 1981). A commonly used approach for such problems is to find an optimal solution that minimizes the weighted least square cost, first proposed by Wahba (1965) for satellite applications, given as: n

min C

1X 2 wi kbi − Cri k2 2 i=1

(1)

T

subject to C C = I, det(C) = 1, where bi ∈ R3 represents ith measurement in the body frame and i = 1, .., n, n being the total number of sensors, ri ∈ R3 is the corresponding vector in the reference frame obtained from some model, wi ∈ R are non-negative weights while k·k2 represents vector 2-norm. Many efficient solutions of this problem for satellite applications can be found in the literature, such as Keat (1977), Shuster and Oh (1981), Markley (1988), Mortari (1997) and Mortari (2000).

5807

10.3182/20110828-6-IT-1002.02001

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

A general approach used in solving the minimization problem (1) is to convert it into an equivalent maximization problem. For this, consider the cost function used in (1) for minimization and expand it as n n n X 1X 1X 2 wi bTi Cri . wi kbi − Cri k2 = wi {bTi bi +riT ri }− 2 i=1 2 i=1 i=1 Neglecting the constant term as it will not change argument of the optimization problem, an equivalent maximization problem can be written as n X max wi bTi Cri C (2) i=1

subject to C T C = I, det(C) = 1. To solve this maximization problem the most popular approach is based on Davenport’s q-method (Keat, 1977; Shuster and Oh, 1981). Two important steps of the qmethod are given now, which will be used for solving the problem addressed in this work. Step 1: Find an equivalent formulation of (2) in terms of a quaternion. This new formulation, first reported in Keat (1977), states that the maximization of (2) is equivalent to the following problem (see Appendix A for derivation): max

q T Kq

subject to

q T q = 1,

q

(3)

where q ∈ R4 represents a 4-parameter quaternion, while K ∈ R4×4 is a symmetric, indefinite and traceless matrix given as  T  B + B − tr(B)I z K := , (4) zT tr(B) Pn 3 T := where B ∈ R3×3 := i=1 wi bi ri and z ∈ R Pn i=1 wi (bi × ri ), while × represents cross product and tr(·) represents trace operator.

Step 2: Using the new formulation, the non-concave maximization problem can be converted into an eigenvalue problem. For this we first add the constraint q T q = 1 using a Lagrange multiplier λ in (3) that gives f (q, λ) = q T Kq − λ(q T q − 1),

The use of a quaternion in the new formulation imparts many benefits. It not only reduces the number of optimization variables from nine to four, but also avoids the constraint det(C) = +1 of (1) being inherent in the definition of quaternion. However, the main benefit obtained is the transformation of the optimization problem into an eigenvalue problem. This classical approach, based on (1) does not directly address the issue of robustness in the measured and model vectors. Mostly, sensitivity analysis have been presented for different algorithms with analytical expressions of the maximum error covariance under stochastic variation in the measurement vector, such as Shuster and Oh (1981), Shuster (2006), Markley (1988), Markley (2008). However, modeling errors in the reference vectors are generally not considered. These errors could be significant. In the case of the earth magnetic field, which is a commonly used sensor especially in satellite applications, errors between sophisticated models and the actual field can be significant (McLean et al., December, 2004; Roithmayr, 2004). The use of simple models, such as low order IGRF (Roithmayr, 2004), which are normally preferred due to less computational cost, results in a less accurate earth magnetic field vector in the reference frame, leading to errors in the attitude estimate. Attitude inaccuracy is further increased due to sensor errors, which are mainly due to noise and installation issues. The magnetic field sensing in the post launch tumbling phase is also a big source of uncertainty. We consider all such errors as ∞norm bounded uncertainties in the input vectors. The main contribution of this paper is consideration of the measurement and modeling uncertainties in the problem formulation for static attitude determination. We formulate a robust optimization (RO) problem for the uncertain inputs. The formulated min-max RO problem is relaxed and transformed to a suboptimal minimization problem using a suitable affine uncertainty parameterization. The resulting non-convex minimization problem with nonlinear cost and constraints is solved using a nonlinear optimization solver with a good initial guess.

(5)

To obtain a stationary point, we solve ∂f / ∂q = 0 and ∂f /∂λ = 0 and obtain the following expression, that has the same formulation as the standard eigenvalue problem Kq = λq. (6) Four eigenvectors of the symmetric matrix K are possible solutions of (6); however, the eigenvector corresponding to the maximum eigenvalue is the solution that will solve (3) (Shuster and Oh, 1981),(Markley and Mortari, 2000), i.e. Kqopt = λmax (K)qopt , (7) where qopt is the solution to (3) and λmax (K) is the maximum eigenvalue of the matrix K. Most of the work in static satellite attitude determination is based on this result and many efficient algorithms have been proposed, such as QUEST (Shuster and Oh, 1981), ESOQ1 (Mortari, 1997), ESOQ2 (Mortari, 2000). A survey paper by Markley and Mortari (2000) provides a general description of many of these algorithms.

2. ROBUST PROBLEM DESCRIPTION To formulate a robust attitude determination problem, we will represent an uncertain measurement vector in the body frame with ¯bi ∈ Bi ⊆ R3 and an uncertain reference vector with r¯i ∈ Ri ⊆ R3 , i = 1, . . . , n, where Bi and Ri are bounded uncertainty sets. To formulate the problem of determining the best uncertainty immunized transformation matrix for attitude determination, we will use the weighted least square approach as used in (1) for the nominal problem. We define the robust problem as n

2 1X min max wi ¯bi − C r¯i 2 C ¯ 2 bi ∈ Bi , r¯i ∈ Ri , i=1 (8) i = 1, . . . , n

subject to C T C = I, det(C) = 1. In order to take advantage of using a quaternion to simplify the optimization problem, as a first step, we

5808

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

reformulate (8) introducing the quaternion q using the same approach used to derive (3). We get ( n ) 1X T T T ¯ min max wi (¯bi ¯bi + r¯i r¯i ) − q Kq q ¯ 2 i=1 bi ∈ Bi , r¯i ∈ Ri , i = 1, . . . , n

rameterization. For this, we normalize βil and ρil with γbi and γri , respectively. The new normalized perturbation vector is δ := [δ1 , δ2 , δ3 ], such that kδk∞ ≤ 1. Using this normalization, we can represent the above mentioned sets Bi and Ri as

q T q = 1,

subject to

Bi = {¯bi = bi +

(9) where   ¯ +B ¯ T − tr(B)I ¯ B z¯ ¯ K := ¯ , z¯T tr(B) n n X X ¯ := wi (¯bi × r¯i ). wi¯bi r¯iT , z¯ := B

Ri = {¯ ri := ri + (10)

¯ ∈ R4×4 is also a In this formulation, similar to (4), K ¯ symmetric indefinite matrix, while B ∈ R3×3 and z¯ ∈ R3 depend on uncertain input vectors. In contrast to (2), the first term in (9) is a function of optimization variables and cannot be neglected in the optimization problem formulation. 3. UNCERTAINTY MODELING To develop a tractable method of solving (9), we define an affine uncertainty parameterization of the uncertainty sets Bi and Ri . A general description of such an affine parameterization is given now (Ben-Tal, Ghaoui, and Nemirovski, 2009). Let βi := [βi1 , βi2 , βi3 ]T and ρi := [ρi1 , ρi2 , ρi3 ]T are vectors of perturbation variables used for uncertainty parameterization, then Bi := {¯bi = bi + Ri := {¯ ri = ri +

3 X l=1

l=1

ρil rˆl : kρi k∞ ≤ γri },

δl˜bil : kδk∞ ≤ 1},

3 X l=1

δl r˜il : kδk∞ ≤ 1},

(12)

4. SOLUTION OF THE ROBUST PROBLEM To solve (9), we transform the robust min-max problem into a minimization problem. In this regard, firstly, we will ¯ similar derive an expression for the uncertain matrix K, to the K matrix of (3), using the presented uncertainty model. This will provide a basis for the subsequent results that describe the simplified but a suboptimal formulation obtained for the robust problem (9). Using the uncertainty model described in Section 3, matrix ¯ defined in (10), can be written as K ¯ =K+ K

Although, each perturbation variable in the body and reference frame has different bound, we can use normalization to obtain a single perturbation variable for uncertainty pa-

(13)



# z ˜ l ˜ l := , K ˜l ) z˜lT tr(B  " # ˜s + B ˜ sT − tr(B ˜ s )I B z˜ls s l l l ˜ Kl := , ˜ls ) z˜lsT tr(B  " # c cT c c ˜lk ˜lk ˜lk B + B − tr( B )I z ˜ c lk ˜ := K . lk cT c ˜lk z˜lk tr(B )

(11)

As we considered all vectors in the body and reference frame to be unit vectors, we choose ˆbl and rˆl as ˆb1 = rˆ1 = T T T [1 0 0] , ˆb2 = rˆ2 = [0 1 0] , ˆb3 = rˆ3 = [0 0 1] . With this choice, γbi and γri represent the maximum possible change in each value of bi and ri as a percentage of the input vector magnitude.

c ˜ l + δl2 K ˜ ls + δl δk K ˜ jk (δl K ),

where the superscripts s and c represent matrices with square and cross terms in δ,k = l+1 if k ≤ 3 or k = l+1−n ˜ l, K ˜ s, K ˜ c ∈ R4×4 are given as: if k > 3. The matrices K l lk "

where γbi and γri represent uncertainty bounds for the input vectors in the body frame and reference frame, respectively. This type of uncertainty is called an interval uncertainty and corresponding perturbation set represents a box (Ben-Tal et al., 2009). The interval uncertainty can be considered a suitable representation for a mix of different type of errors, such as sensor noise, having a stochastic interpretation, biases and misalignments, which are constant and modeling inaccuracies, etc. However, bounds for each measurement or model vector should be carefully chosen, as unnecessary large values may result in large residual between the non-robust and the robust solution for nominal cases.

3 X l=1

βilˆbl : kβi k∞ ≤ γbi },

3 X

l=1

where ˜bil := γbiˆbl and r˜il := γri rˆl .

i=1

i=1

3 X

˜l + B ˜ T − tr(B ˜l )I B l

(14)

The other terms used in these expressions are givens as: ˜l := B

n X

T ˜ s := wi (bi r˜il + ˜bil riT ), B l

c ˜lk B :=

T , wi˜bil r˜il

i=1

i=1

n X

n X

T T wi (˜bil r˜ik + ˜bik r˜il ),

i=1

z˜l :=

n X i=1

z˜lk :=

n X i=1

wi {(bi × r˜il ) + (˜bil × ri )}, z˜ls := wi {(˜bil × r˜ik ) +

n X i=1

n X i=1

wi {(˜bik × r˜il ).

wi (˜bil × r˜il ), (15)

˜ l depends on measurements, It can be observed that K ˜ s and K ˜ c depend on uncertainty modeling and but K l lk remain fixed for different measurement vectors, provided the uncertainty description remains same. We now present two results using (13), that will be used to reformulate (9) in terms of the proposed uncertainty model.

5809

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

Lemma 1. The formulation given in (9) is equivalent to T

min

T

T

{−q Kq+ max (p δ + δ Qδ)}

q

kδk∞ ≤1

T

subject to

(16)

q q = 1,

where p ∈ R3 and Q ∈ R3×3 depend on q and are defined as   ˜ 1q c˜1 − q T K ˜ 2 q , p := c˜2 − q T K T ˜ c˜3 − q K3 q   1 T ˜c 1 T ˜c s T ˜s q K q − q K q c ˜ − q K q − 12 13  1  1 2 2   1 1 T c s T s T c .  ˜ 12 q ˜2q ˜ 23 Q :=  − q K c˜2 − q K − q K q 2   21 1 T ˜c s T ˜s ˜c q − q K q c ˜ − q K q − qT K 13 23 3 3 2 2 In this equation, c˜l , c˜sl are constants that depend on the uncertain vectors and are given as: n

c˜l :=

1X T ri ), wi (bTi ˜bil + ˜bTil bi + riT r˜il + r˜il 2 i=1 n

c˜sl

1X T wi (˜bTil ˜bil + r˜il r˜il ). := 2 i=1

(17)

Proof. We start with the cost given in (9). Expanding the first term, we get n n 1X 1X wi (¯bTi ¯bi + r¯iT r¯i ) = wi (bTi bi + riT ri ) 2 i=1 2 i=1 +

3 X

(δl c˜l + δl2 c˜sl )

(18)

l=1

Using (18) and (13) and joining similar terms in δl we can write the required expression.  The robust problem (16) is always feasible because the solution of the nominal problem (3) always exists (Shuster and Oh, 1981). However, finding the optimal maximum value of pT δ + δ T Qδ in (16) is difficult as the term is nonconcave in δ (Q is positive semidefinite). Based on this fact, we determine an upper bound on the maximum of (pT δ + δ T Qδ) over δ. The result is given in the following lemma. Lemma 2. An upper bound on the maximization term appearing in (16) is max (pT δ + δ T Qδ) ≤ kpk1 + 3λmax (Q) 0 ≤ kδk ≤1

(19)



Proof. The left hand side in (19) can be written as max (pT δ + δ T Qδ) ≤ max pT δ + max δ T Qδ (20) kδk ≤1 kδk ≤1

kδk∞ ≤1





Using the H¨older dual norm (Boyd and Vandenberghe, 2004), the first term on the right hand side in (20) is given as max pT δ = kpk1 . (21) kδk ≤1 ∞

For the second term appearing in (20), since Q is a symmetric matrix, we can write the maximum eigenvalue of Q as (Boyd and Vandenberghe, 2004) λmax (Q) = sup δ T Qδ . kδk ≤1 2

(22)

Hence, we first replace the ∞-norm in the second term on the right hand √ side of (20) with the 2-norm using the inequality kδk2 ≤ 3 kδk∞ for δ ∈ R3 . We can write max δ T Qδ ≤

kδk∞ ≤1



max√ δ kδk2 ≤ 3

T

Qδ (23)

3λmax (Q),

Using (21) and (23), we can write (19).



Lemma 2 gives an upper bound on the maximum value of the quadratic term pT δ+δ T Qδ and can be used to simplify (16). Theorem 1. The min-max problem (16) can be approximated with the following minimization problem min

q,u1 ,u2 ,u3

−q T Kq + u1 + u2 + u3 + 3λmax (Q)

q T q = 1, −uj ≤ pj ≤ uj , j = 1, 2, 3, where pj , j = 1, 2, 3 are elements of vector p. subject to

(24)

Proof. In formulation (16) given in Lemma 1, replacing the max term with the upper bound given in Lemma 2 and using the P fact that a set of 2n + 1 linear inequalities −uj ≤ xj ≤ uj , j uj ≤ 1, j = 1, . . . , 3 represent the nonlinear P inequality j |xj | ≤ 1 (Ben-Tal, Ghaoui, and Nemirovski, 2009, definition 1.3.1), we can write (24).  The argument q of the approximate robust problem (24) approaches the argument of the nominal problem (3) when the bound on the uncertainty approaches zero. This becomes evident if, instead of taking kδk∞ ≤ 1, we assume a bound α on each element of the normalized vector δ. In that case we can write (19) as 0≤

max (pT δ + δ T Qδ) ≤ α kpk1 + 3α2 λmax (Q). (25)

kδk∞ ≤α

If α tends to zero, the right hand side of (25) will also tend to zero, making the robust problem equivalent to (3). Since the cost and constraints in (24) are nonlinear, we solve it using a nonlinear optimization solver with a good initial guess. This initial guess for q is obtained from the non-robust form of the problem using the q-Method or QUEST (Shuster and Oh, 1981), etc. 5. SIMULATION RESULTS We consider the problem of attitude estimation for a low cost CubeSat (Heidt et al., 2001), a pico class of satellite having three units each of dimension 10 × 10 × 10 cm3 and weight less than 1 kg. CubeSats are mainly used for scientific missions and currently most of the satellites are being developed in universities. Being low cost, interest in such satellites is increasing for a wider class of applications; however, there are many challenges from an estimation and control point of view to successfully achieve a given objective. The CubeSat under consideration is assumed to be moving in a circular orbit having a radius of 650 km from the center of the earth. We are considering that only two measurements are available in the satellite, namely the earth magnetic field and the sun vector. For the earth magnetic field, two magnetometers are installed, one inside the satellite, which is mainly used in the post-

5810

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

launch phase when the satellite is recovering from launch disturbances, while the second is installed on an extended boom, which is deployed once the satellite has de-tumbled and achieved an equilibrium. The sun vector is sensed by a pair of sun sensors installed on the satellite body. Both of these measurements are in the body frame. We used the orbit frame as the reference frame. For the reference earth magnetic field vector the first order IGRF model (Roithmayr, 2004) is used, while the reference sun vector is obtained using a simplified sun model based on sun ephemeris (Wertz, 1978).

non-robust robust 0.02

0.015

J

0.01

0.005

0 1

1

0.5

0.5

0

0

−0.5

−0.5 −1

δ2

−1 δ1

Fig. 1. Robust and non-robust design comparison for given set of uncertain vectors at one instant. non-robust

φ

50

robust

actual

0 −50 0

500

1000

1500

0

500

1000

1500

0

500

1000

1500

θ

50

0

−50 100

ψ

Sensor measurements are not accurate due to many factors, such as sensor noise, misalignments due to installation errors, biases, etc. Especially in the post-launch tumbling phase measurement errors further increase due to the use of an internal magnetometer installed on-board the satellite, which suffers from interaction with the magnetic field generated by surrounding electronics. Similarly, the reference vectors for both measurements are also not exact since the mathematical models used to obtain measurement vectors in the reference frame are based on loworder approximations. For simulation purpose, we consider that the bound on uncertainty in the measurements and reference vectors i.e. γbi and γri , respectively to be 20% of the input magnitude. We take equal weights, such that P n i=1 wi = 1. We need to solve the optimization problem formulated in (24). Since this formulation includes nonlinear objective function and constraints, we used the nonlinear optimization solver (fmincon) of MATLAB. Two plots are presented to see the effectiveness of the proposed methodology. Firstly, the problem was solved for a given uncertain data set for one instant. To obtain the uncertain data set, we added a uniformally distributed random error in the range of ±γbi and ±γri in each measured and model vectors obtained from a simulation without error. The vectors with added uncertainty are given as:

0

−100

T

b1 = [−0.590 −0.367 0.719] ,

Time (sec)

T

r1 = [−0.583 −0.344 0.736] , T

b2 = [−0.721 0.021 0.693] , T

r2 = [−0.495 −0.0004 0.869] .

Figures 1 presents a comparison between the non-robust and the robust approach for one time instant. For this comparison, a non-robust solution was obtained satisfying (3), while a robust solution was found satisfying (24). Then, two components of the normalized uncertainty vector δ, i.e. δ1 and δ2 are varied, while keeping δ3 constant P for simulan tion purpose, and a variation in the cost J := 12 i=1 k¯bi − 2 C r¯i k2 is calculated for both solutions. This cost is plotted along the z-axis for both solutions, while x- and y-axis represent δ1 and δ2 . It can be observed from both plots that the cost J obtained using the nominal solution has large variation when δ1 and δ2 are close to set bounds, while the J obtained using the robust solution shows the best immunization against uncertainty. However, the cost obtained for the robust solution has an offset when compared with the cost value obtained using the nominal solution without error. Figure 2 presents a performance comparison of the robust and non-robust approaches in the presence of uncertain-

Fig. 2. Comparison of attitude angles obtained using nonrobust and robust algorithms. The dotted line shows the data without errors while the other two cases include errors within the uncertainty bound ±γbi or ±γri .

ties, using in-orbit data obtained from nonlinear closedloop simulation for the satellite. The error free data obtained from simulation was corrupted by adding uniformly distributed random errors in the range of ±γbi and ±γri in corresponding vectors. The simulation was run with a sample time of 1 second with initial roll, pitch and yaw body rates of 0.05, 0.15 and 0 deg/s and roll, pitch and yaw angles of 30, 0, 0 deg, respectively. It can be observed that due to the uncertainties in the input information, the non-robust approach gives large errors in the estimated attitude angles, while the robust approach gives much better performance, limiting the maximum attitude error to a smaller band. 6. CONCLUSIONS We presented a formulation for robust attitude determination with ∞-norm bounded uncertainty in the input vectors. The robust optimization problem, using a weighted

5811

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

least square approach, has a nonlinear cost and constraints. The formulated min-max optimization problem was transformed into a suboptimal minimization problem due to the use of an upper bound and solved using a nonlinear optimization solver with a good initial guess obtained from the nominal solution. The results showed that the robust approach has significant benefit over the nominal approach in the situations when inputs are uncertain with bounded uncertainty. The benefit is generally maximum in the worst case scenarios. However, the robust solution for the case without error has a larger residual than the nonrobust solution, which is expected in the robust optimal solutions. However, the main issue is increased computational burden in solving the optimization problem. Further work will focus on finding a computationally efficient solution of this problem.

Shuster, M.D. (2006). The generalized wahba problem. Journal of Astronautical Sciences, 54(2), 245–259. Sidi, M.J. (2000). Spacecraft Dynamics and Control. Cambridge University Press, UK, 2nd edition. Wahba, G. (1965). A least square estimate of spacecraft attitude. SIAM Review, 7(3), 409. Wertz, J.R. (1978). Spacecraft Attitude Determination and Control. D. Reidel Publishing Company, Holland. Appendix A. DAVENPORT TRANSFORMATION To derive the Davenport transformation, consider the cost function given in (2) n X

wi bTi Cri =

The first author is thankful to the Higher Education Commission of Pakistan for supporting his PhD studies. REFERENCES Ben-Tal, A., Ghaoui, L.E., and Nemirovski, A. (2009). Robust Optimization. Princeton University Press, USA. Boyd, S. and Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press. Crassidis, J.L., Markley, F.L., and Cheng, Y. (2007). Survey of nonlinear attitude estimation methods. Journal of Guidance, Control and Dynamics, 30(1), 12–28. Heidt, H., Suari, J.P., Moore, A.S., Nakasuka, S., and Twiggs, R.J. (2001). Cubesat: A new generation of pico-satellite for education and industry low cost space experimentation. In 15th Annual/USU Conference on Small Satellite. SSC01-VIIIb-5. Keat, J.E. (1977). Analysis of least-squares attitude determination routine DOAOP. Technical Report NASA/CSC/TM-77/6034. Markley, F.L. (1988). Attitude determination using vector observations and the singular value decomposition. Journal of Astronatutical Sciences, 36(3), 245–258. Markley, F.L. (2008). Optimal attitude matrix from two vector measurements. Journal of Guidance, Control and Dynamics, 31(3), 765–768. Markley, F.L. and Mortari, D. (2000). New developments in quaternion estimation from vector observation. In Advances in the Astronautical Sciences, volume 106, p 373–393. McLean, S., Macmillan, S., Maus, S., Lesur, V., Thomson, A., and Dater, D. (December, 2004). The US/UK world magnetic model for 2005-2010. Technical Report NOAA, NESDIS/NGDC-1. Mortari, D. (1997). ESOQ: A closed form solution to the Wahba problem. Journal of the Astronautical Sciences, 45(2), 195–204. Mortari, D. (2000). Second estimator of the optimal quaternion. Journal of Guidance, Control and Dynamics, 23(5), 885–888. Roithmayr, C.M. (2004). Contributions of spherical harmonics to magnetic and gravitational fields. Technical Report NASA/TM-2004-213007. Shuster, M.D. and Oh, S. (1981). Three-axis attitude determination from vector observation. Journal of guidance and Control, 4(1), 70–77.

tr(wi bTi Cri ).

(A.1)

i=1

i=1

ACKNOWLEDGEMENTS

n X

Now we use two properties of the trace. Firstly trace is P invariant P under cyclic permutations, and secondly, tr ( i Ai ) = i tr(Ai )∀A ∈ Rn×n . Using these properties we can write n X wi bTi Cri = tr(CB T ), (A.2) i=1

Pn where B := i=1 wi bi riT . Now we represent C using the  T T quaternion q ∈ R4 := qT q4 , where q := [q1 q2 q3 ] , written as (Shuster and Oh, 1981)

C = (q42 − qT q)I + 2qqT + 2q4 Q. (A.3) In the above equation Q := −q×, where × represents vector cross product, is given as # 0 q3 −q2 Q = −q3 0 q1 . q2 −q1 0 "

(A.4)

Substituting (A.3) in (A.2), we get n X i=1

wi bTi Cri = (q42 − qT q)tr(B T ) + 2tr(qqT B T ) + 2q4 tr(QB T ).

(A.5)

The second term on the right hand side of the above equation can be written as 2tr(qqT B T ) = qT (B T + B)q.

(A.6)

The last term can be written as 2q4 tr(QB T ) = q4 (qT z + z T q), (A.7) Pn where z := i=1 wi (bi × ri ). Substituting (A.6) and (A.7) in (A.5), we get n X i=1

wi bTi Cri



= q

T

    B T + B − tr(B)I z q q4 zT tr(B) q4

= q T Kq,

(A.8)

4×4

which is the required form of (A.1), where K ∈ R symmetric, traceless matrix given as   T B + B − tr(B)I z . K := zT tr(B)

5812

is a

(A.9)