Next Article in Journal
Real-Time Precise Orbit Determination for LEO between Kinematic and Reduced-Dynamic with Ambiguity Resolution
Next Article in Special Issue
Attitude Stabilization of a Satellite Having Only Electromagnetic Actuation Using Oscillating Controls
Previous Article in Journal
Modal Analysis of Conceptual Microsatellite Design Employing Perforated Structural Components for Mass Reduction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Quaternion versus Rotation Matrix Feedback for Spacecraft Attitude Stabilization Using Magnetorquers †

School of Aerospace Engineering, Sapienza University of Rome, via Salaria 851, 00138 Roma, Italy
This paper is an extended version of my paper published in Fifth IAA Conference on University Satellite Missions and Cubesat Workshop 2020, Rome, Italy, 28–31 January 2020.
Aerospace 2022, 9(1), 24; https://doi.org/10.3390/aerospace9010024
Submission received: 8 December 2021 / Revised: 27 December 2021 / Accepted: 29 December 2021 / Published: 4 January 2022
(This article belongs to the Special Issue Spacecraft Attitude Control Using Magnetic Actuators)

Abstract

:
The purpose of this paper is to compare performances between stabilization algorithms of quaternion plus attitude rate feedback and rotation matrix plus attitude rate feedback for an Earth-pointing spacecraft with magnetorquers as the only torque actuators. From a mathematical point of view, an important difference between the two stabilizing laws is that only quaternion feedback can exhibit an undesired behavior known as the unwinding phenomenon. A numerical case study is considered, and two Monte Carlo campaigns are carried out: one in nominal conditions and one in perturbed conditions. It turns out that quaternion feedback compares more favorably in terms of the speed of convergence in both campaigns, and it requires less energy in perturbed conditions.

1. Introduction

Magnetic actuators, also known as magnetorquers, are widely used as torque actuators for attitude control of spacecraft in low Earth orbits. The interest in such devices is due to the following reasons: (i) they are simple, reliable, and low cost; (ii) they need only renewable electrical power to be operated; and (iii) magnetorquers save system weight with respect to any other class of torque actuators. On the other hand, magnetorquers have the following important limitation: (i) the control torque is constrained to belong to the plane orthogonal to the Earth’s magnetic field; (ii) the maximum torque they can generate is substantially smaller than for other types of torque actuators. Due to these limitations, different types of actuators often accompany magnetorquers to provide full three-axis control, and a considerable amount of work has been dedicated to the design of magnetic control laws in the latter setting (see [1] and references therein). However, it is possible to stabilize the attitude of spacecraft in low Earth orbit by using only magnetorquers if the orbit inclination is not too low. Using such actuators convergence to the nominal attitude occurs more slowly compared to other torque actuators (see [2]).
Simple control algorithms that achieves attitude stabilization using only magnetorquers, are based on quaternion and attitude rate feedback. It has been shown that such control laws achieves stabilization for both inertial pointing spacecraft [2,3,4,5] and Earth (or nadir) pointing spacecraft [6,7,8,9,10]. However, using quaternion and attitude rate feedback can result in undesired unwinding phenomenon for the spacecraft. Such phenomenona originate from the fact that one of the two opposite quaternions that represent the nominal spacecraft attitude is unstable. Thus, if, during attitude evolution, the current quaternion becomes very close to the latter unstable quaternion, then the spacecraft attitude rather than converging to the nominal one is most likely pushed away from it [11]. Clearly, the unwinding phenomenon results in both slower convergence and higher energy consumption in performing attitude stabilization.
An alternative attitude stabilizing law that does not generate such an undesired phenomenon is obtained by replacing quaternion feedback with an appropriate rotation matrix feedback [11]. Since the desired attitude is represented by only one rotation matrix that is made asymptotically stable, clearly no unwinding phenomenon can occur.
The goal of this paper is to compare performances of quaternion and attitude rate feedback with those of rotation matrix and attitude rate feedback. The comparison is performed numerically by using the spacecraft model presented in [12] in terms of both speed of convergence and energy consumption.
The rest of the paper is organized as follows. In the next section, the models of the spacecraft and of the geomagnetic field are presented. In Section 3, the two stabilizing laws under examination are discussed. The numerical comparison of the control laws is performed in Section 4. Concluding remarks end the paper. Two appendices present the methods adopted for selecting the control gains and for compensating for residual magnetic dipole moments.

2. Models

2.1. Coordinate Frames

In order to describe the attitude dynamics of a rigid Earth-pointing spacecraft in a circular orbit and to represent the geomagnetic field, it is useful to introduce the following reference frames:
  • Geocentric Inertial Frame F i . A commonly used inertial frame for Earth orbits is the Geocentric Inertial Frame, for which its origin is in the Earth’s center, its x i axis is the vernal equinox direction, the z i axis coincides with the Earth axis of rotation and points northward, and the y i axis completes the right-handed frame (see ([13], Section 2.6.1)).
  • Orbital frame F o . The origin of this frame, also called Local Horizontal Local Vertical (LVLH), is in the satellites center of mass. The x o axis is along the spacecraft velocity. The z o axis points to the Earth’s center, and the y o axis is perpendicular to the orbit plane and completes the right-handed frame.
  • Spacecraft body frame F b . The origin of this frame coincides with the satellite center of mass. Its axes are principal axes of inertia of the spacecraft.
The nominal Earth-pointing attitude corresponds to F b aligned with F o .
The notation used throughout the paper is presented in Table 1.

2.2. Attitude Kinematics

Since we deal with an Earth-pointing spacecraft in this work, the focus will be on the kinematics of the spacecraft with respect to the orbital frame. The attitude kinematics equations depend on the representation adopted for the attitude of F b with respect to F o .
If the attitude is represented by the following quaternion:
q = [ q 1 q 2 q 3 q 4 ] T = [ q v T q 4 ] T
with q = 1 , then the attitude kinematics is given by the following (see Section 5.5.3 of [14]):
q ˙ v = 1 2 ( q v × + q 4 I 3 × 3 ) ω b o b q ˙ 4 = 1 2 q v T ω b o b
where I 3 × 3 is the 3 × 3 identity matrix, and the notation v × for v = [ v 1 v 2 v 3 ] T represents the skew symmetric matrix:
v × = 0 v 3 v 2 v 3 0 v 1 v 2 v 1 0
so that for w = [ w 1 w 2 w 3 ] T , it occurs that v × w = v × w .
On the other hand, if the attitude is represented by rotation matrix R b o , which transforms vectors of coordinates in F o into vectors of coordinates in F b , then the the attitude kinematics is given by the following (see ([15], Section 1.4.1)).
R ˙ b o = ( ω b o b ) × R b o
The relation between q and R b o is given by the following.
R b o = ( q 4 2 q v T q v ) I 3 + 2 q v q v T 2 q 4 q v ×
Clearly, the following is obtained.
ω b o b = ω b i b ω o i b
Moreover, since a circular orbit is considered, then ω o i o = [ 0 n 0 ] T where n is the constant orbital rate, and, consequently, the following is obtained.
ω o i b = R b o ω o i o

2.3. Attitude Dynamics and Geomagnetic Field

Attitude dynamics can be expressed in a body frame as follows:
J ω ˙ b i b = ω b i b × ( J ω b i b ) + T g g b + T c o i l s b + T d i s t b
where J = diag [ J x , J y , J z ] is the spacecraft inertia matrix, T g g b , T c o i l s b , and T d i s t b are the body coordinates of the gravity gradient torque, the control torque induced by magnetic coils, and the disturbance torque, respectively.
The gravity gradient torque in body coordinates is given by the following (see ([14], Section 6.10)):
T g g b = 3 n 2 z ^ o b × J z ^ o b
where z ^ o b denotes the unit vector corresponding to the z o -axis of F o resolved in the body frame, which can be expressed as follows:
z ^ o b = R b o e 3
where the following is the case.
e 3 = 0 0 1
The spacecraft is equipped with three magnetic coils aligned with the F b -axes, which generate the following magnetic control torque:
T c o i l s b = m c o i l s × b b = ( b b ) × m c o i l s
where m c o i l s is the vector obtained by stacking the magnetic moments of the three coils, and b b is the geomagnetic field at the spacecraft expressed in F b . Clearly, the relation between b b and b o is given by the following.
b b = R b o b o
For what follows, it useful to rewrite Equation (7) in terms of ω b o b instead of ω b i b . For that purpose, note that from Equations (5) and (6) it follows that ω b i b = ω b o b + R b o ω o i o . Thus, since ω o i o is constant, the following holds.
ω ˙ b i b = ω ˙ b o b + R ˙ b o ω o i o
By using Equation (3), we obtain the following.
R ˙ b o ω o i o = ( ω b o b ) × R b o ω o i o = ( ω b o b ) × ω o i b = ω o i b × ω b o b = ( R b o ω o i o ) × ω b o b
Thus, from Equations (7)–(13), we obtain the following.
J ω ˙ b o b = J ( R b o ω o i o ) × ( ω b o b ) × J + ( J R b o ω o i o ) × ( R b o ω o i o ) × J ω b o b ( R b o ω o i o ) × J R b o ω o i o + 3 n 2 R b o e 3 × J R b o e 3 ( R b o b o ) × m c o i l s + T d i s t b
Let b i be the geomagnetic field at spacecraft expressed in inertial frame F i and let R o i be the rotation matrix that transforms vectors of coordinates in F i into vectors of coordinates in F o . Note that b i and R o i vary with time at least because of the spacecraft’s motion along the orbit. Then, the following is the case:
b o ( t ) = R o i ( t ) b i ( t )
which shows explicitly the dependence of b o on t. In order to study Equation (14), it is important to characterize the time-dependence of b o ( t ) , which corresponds to characterizing the time-dependence of b i ( t ) and R o i ( t ) . By adopting the so-called inclined dipole model of the geomagnetic field (see ([16], Appendix H)) and letting R o r b i t denote the radius of the circular orbit, we obtain the following.
b i ( t ) = μ i d R o r b i t 3 [ 3 ( ( m ^ i ( t ) ) T r ^ i ( t ) ) r ^ i ( t ) m ^ i ( t ) ] .
In Equation (16), μ i d is the total strength of the inclined dipole, r i ( t ) is the spacecraft position vector resolved in F i , and r ^ i ( t ) is the vector of the direction cosines of r i ( t ) . The components of vector m ^ i ( t ) are the direction cosines of the Earth’s magnetic dipole expressed in F i , which is set equal to the following:
m ^ i ( t ) = sin ( θ i d ) cos ( ω e t + α 0 ) sin ( θ i d ) sin ( ω e t + α 0 ) cos ( θ i d ) ,
where θ i d is the coelevation of the inclined dipole, ω e is the Earth average rotation rate, and α 0 is the right ascension of the dipole at time t = 0 .
In order to characterize the time dependence of b i ( t ) in (16), one needs to determine an expression for r i ( t ) , which is the spacecraft’s position vector resolved in F i . Define a coordinate system x p , y p in the orbital plane, for which its origin is at the center of the Earth and with the x p axis coinciding with the line of nodes. Then, the position of satellite centre of mass is given by the following:
x p ( t ) = R o r b i t cos ( n t + ϕ 0 ) y p ( t ) = R o r b i t sin ( n t + ϕ 0 ) ,
where ϕ 0 is the argument of the spacecraft at time t = 0 . Let i n c l be the orbit inclination and let Ω dentote the right ascension of the Ascending Node (RAAN) of the orbit (see ([13], Section 2.6.2)). Then, the coordinates of the satellite center of mass in the inertial frame can be obtained as follows:
r i ( t ) = R z ( Ω ) R x ( i n c l ) x p ( t ) y p ( t ) 0
where the following:
R x ( φ ) = 1 0 0 0 cos ( φ ) sin ( φ ) 0 sin ( φ ) cos ( φ )
is the rotation matrix corresponding to a rotation around the x-axis of magnitude φ and the following:
R z ( ψ ) = cos ( ψ ) sin ( ψ ) 0 sin ( ψ ) cos ( ψ ) 0 0 0 1
is the rotation matrix corresponding to a rotation around the z-axis of magnitude ψ (see ([13], Section 2.6.2)).
By combining Equations (16)–(21), the expression of b i ( t ) can be easily obtained. Moreover, the following holds.
R o i ( t ) = 0 1 0 0 0 1 1 0 0 R z ( n t + ϕ 0 ) R x ( i n c l ) R z ( Ω )
Thus, by using Equation (15), an explicit expression for b o ( t ) can be derived. It is easy to see that b o ( t ) can be expressed as a sum of sinusoidal functions of t having different frequencies since sinusoidal functions having angular frequencies n and ω e appear in the previous equations.
A simpler model of the geomagnetic field is the axial dipole model in which the Earth’s magnetic dipole is aligned with the Earth’s rotation axis (see [17]). Thus, the axial dipole model is obtained by setting m ^ i = [ 0 0 1 ] T in Equation (16) and by replacing μ i d with μ a d . Using such a model, the expression of b o ( t ) is simplified as follows.
b o ( t ) = μ a d R o r b i t 3 sin ( i n c l ) cos ( n t + ϕ 0 ) cos ( i n c l ) 2 sin ( i n c l ) sin ( n t + ϕ 0 )
The latter equation shows that the adoption of a simpler model results in a sinusoidal b o ( t ) with period T o r b i t = 2 π / n .
The spacecraft and orbit data employed in the following numerical study are obtained from [12]. The geomagnetic field data are obtained from [17]. Both data are reported in Table 2.

2.4. Disturbance Torques

The most significant disturbance torques acting on a spacecraft in low Earth orbit are modeled as follows (see [12]). The residual magnetic torque in body coordinates is given by the following:
T r m b = m r m b × b b
in which m r m b is the body coordinate of the residual magnetic dipole moment due to onboard electrical components. The aerodynamic torque is modeled as follows: T a b = r a p b × F a b where r a p b is the body coordinate of the vector from the center of mass to the center of pressure, and F a b is the body coordinate of the aerodynamic force acting on the spacecraft. A simplified model for F a b is considered by setting F a b = 1 2 C D S D ρ v b v b in which C D is the drag coefficient, S D is the area of the spacecraft cross section, ρ is the atmospheric density at orbit altitude, and v b is the body coordinate of spacecraft velocity with respect to the air, which is approximated along with spacecraft velocity. The solar radiation pressure torque is modeled as T s r p b = r s p b × F s r p b , where r s p b is the body coordinate of the vector from the center of mass to the center of solar pressure, and F s r p b is the body coordinate of the force due to solar radiation pressure, which is modeled as F s r p b = ( Φ s / c ) ( 1 + q s ) A s s ^ b . In the last equation, Φ s is the solar flux density constant, c is the speed of light, q s is the reflectance factor, A s is the area of the sunlit surface, which is assumed constant in the worst case scenario, and s ^ b is the body coordinate of the unit vector from the spacecraft to the Sun. Note that s ^ b = R b o R o i s ^ i , where s ^ i are the coordinates in the inertial frame F i of the same unit vector.
The disturbance torques data are obtained from [12] and reported in Table 3.

3. Attitude Stabilizing Laws

In this section, two attitude stabilizing laws are considered, and their performances are compared in Section 4.

3.1. Quaternion and Attitude Rate Feedback

Consider the following quaternion and attitude rate feedback:
m c o i l s = m m a x sat 1 m m a x b b × ( K p q v + K d ω b o b )
where sat is the standard saturation function, and m m a x is the maximum control dipole moment for which its numerical value is reported in Table 2. The numerical values of the matrix gains K p and K d are determined by solving the linear quadratic optimal control problem described in Appendix A, which results in the values for K p and K d reported in Equation (A3).
Consider the closed-loop system given by Equations (1), (14), and (24) and set T d i s t b = 0 and b o , as in Equation (22). For the latter system the following properties hold. Equilibrium ( q , ω b o b ) = ( q ¯ , 0 ) is asymptotically stable where q ¯ = [ 0 0 0 1 ] T . In fact, it can be verified that the linear approximation of the latter periodic system about ( q , ω b o b ) = ( q ¯ , 0 ) is given by Equations (A1) and (A2) and has all characteristic exponents with negative real parts. On the other hand, equilibrium ( q , ω b o b ) = ( q ¯ , 0 ) is unstable since the linear approximation about ( q , ω b o b ) = ( q ¯ , 0 ) has three positive characteristic exponents [18].
It is well known that both opposite quaternions q ¯ and q ¯ correspond to the nominal attitude having F b aligned with F o . Thus, the fact that ( q , ω b o b ) = ( q ¯ , 0 ) is unstable can result in the so-called unwinding phenomenon. In fact, if, during the spacecraft motion, it occurs that ( q , ω b o b ) becomes very close to ( q ¯ , 0 ) , then the spacecraft attitude rather than converging to the nominal one is most likely pushed away from it [11]. Such a behavior typically results in both slower convergence and higher energy consumption in performing attitude maneuvers. Consequently, it represents a weak point of the proposed control law.

3.2. Rotation Matrix and Attitude Rate Feedback

Consider the following rotation matrix and attitude rate feedback [11].
m c o i l s = m m a x sat 1 m m a x b b × K p 4 i = 1 3 ( e i × R b o T e i ) + K d ω b o b
In the previous equation, the numerical values of K p and K d are chosen as in Equation (A3), and e 1 , e 2 , and e 3 are the columns of I 3 × 3 .
Consider the closed-loop system given by Equations (1), (14), and (25) and set T d i s t b = 0 and b o as in Equation (22). Clearly, the obtained system is periodic with period T o r b i t too. Moreover, for the latter system, equilibrium ( R b o , ω b o b ) = ( I 3 × 3 , 0 ) is asymptotically stable. In fact, it can be verified that the linear approximation of the latter periodic system about ( R b o , ω b o b ) = ( I 3 × 3 , 0 ) has all characteristic exponents with a negative real part. Actually, the linear approximation of the closed-loop system about ( R b o , ω b o b ) = ( I 3 × 3 , 0 ) is identical to the linear approximation about ( q , ω b o b ) = ( q ¯ , 0 ) considered in the previous section. This is a consequence of the fact that both the linear approximation of the stabilizing law in Equation (24) about ( q , ω b o b ) = ( q ¯ , 0 ) and the linear approximation of the stabilizing law in Equation (25) about ( R b o , ω b o b ) = ( I 3 × 3 , 0 ) can be expressed as follows:
m c o i l s = b b × 1 2 K p ζ + K d ω b o b
where ζ = [ ϕ θ ψ ] T are the 3-2-1 Euler angles representing the attitude of F b with respect to F o (see [11]). Thus, the two control laws locally are identical. As a result, given an initial attitude sufficiently close to the nominal one and an initial angular velocity sufficiently small, if the initial attitude is represented by the quaternion closer to q ¯ , then the corresponding evolutions of the spacecraft attitude obtained by applying either control laws in Equation (24) or the one in Equation (25) will be very close to each other.
Clearly, the control law in Equation (25) does not generate the unwinding phenomenon. In fact, the desired attitude is achieved only when ( R b o , ω b o b ) = ( I 3 × 3 , 0 ) , which is an asymptotically stable equilibrium. As a result, it is interesting to compare the performances of the stabilizing laws in Equations (24) and (25) to observe how much the unwinding phenomenon affects the performances of the control law in Equation (24).

4. Quaternion versus Rotation Matrix Feedback

The purpose of this section is to compare, numerically, the performances of the quaternion and attitude rate feedback in Equation (24) with those of the rotation matrix and attitude rate feedback in Equation (25). The comparison is carried out by using all the numerical values indicated in the previous sections. The focus of the comparison is on both speed of convergence to the nominal attitude and energy required to perform the control action. The results of two simulation cases are presented. In the first one, nominal conditions are considered with the main aim of studying the effects of the unwinding phenomenon on the speed of convergence and on energy consumption. In the second case, perturbations due to disturbance torques and to the adoption of the tilted dipole model are introduced so that a comparison in a more realistic scenario is obtained.

4.1. Nominal Case

In this subsection, the comparison is executed by neglecting disturbance torque T d i s t b and considering the aligned dipole model for the geomagnetic field. In these conditions, the closed-loop systems are periodic with period T o r b i t , and the stability results stated in Section 4 hold. Thus, it is of interest to determine in such a scenario the effect of the unwinding phenomenon on the performances. Two performance indexes are considered. The first one is the settling time t s for the principal angle, which is defined as the time needed for the principal angle between F b and F o to become and stay smaller or equal to 1 deg. Clearly, the settling time measures the speed of convergence relative to the nominal attitude. The second performance index is given by the following:
E c o i l s = 0 t f i n m c o i l s ( t ) 2 d t
which provides an indirect and rough measure of the energy required to perform the attitude control action. Time t f i n in the previous equation is set to 30 T o r b i t so that the spacecraft is able to achieve the nominal attitude.
A Monte Carlo campaign has been run with the intent of obtaining some useful statistical information on the performances of the two control laws. Specifically, 100 simulations have been run. In each run, the initial rotation matrix R b o ( 0 ) has been selected randomly. The initial quaternion q ( 0 ) has been obtained by converting R b o ( 0 ) into a quaternion and selecting the quaternion with a non-negative scalar component q 4 ( 0 ) . The latter choice should reduce the possibility for the winding phenomenon to occur. The initial angular velocity has been selected randomly within the following set ω b o b ( 0 ) 20 deg/s. The argument ϕ 0 of the spacecraft at time t = 0 (see Equation (18)) has been randomly selected too.
The results of the campaign are reported in Table 4.
The first row (Q) refers to quaternion feedback, and the second row (RM) refers to rotation matrix feedback. Symbol t ¯ s denotes the mean value of t s , and E ¯ c o i l s denotes the mean value of E c o i l s . Symbol t s ( RM ) t s ( Q ) denotes the number of runs (expressed in percentage on the total number) in which t s ( RM ) t s ( Q ) . Similarly, E c o i l s ( RM ) E c o i l s ( Q ) denotes the number of runs (expressed in percentage on the total number) in which E c o i l s ( RM ) E c o i l s ( Q ) .
The statistics show that, in terms of settling time, quaternion feedback performs better than rotation matrix feedback. In terms of energy consumption, the two control laws are practically equivalent even if in 65% of the simulation runs rotation matrix feedback uses less energy than quaternion feedback.
Figure 1 displays the time histories of the principal angle for a single run obtained by using quaternion feedback and rotation matrix feedback.
It turns out that, for this specific run, the settling time obtained with quaternion feedback is equal to 11 orbital periods, whereas by using rotation matrix feedback, the settling time is equal to 14 orbital periods. The value of E c o i l s is equal to 2.50 104 A2 m4 s for both control laws.

4.2. Perturbed Case

In this subsection, the performance comparison is carried out by considering a more realistic scenario in which the disturbance torques described in Section 2.4 act on the spacecraft and in which the inclined dipole model is adopted. Moreover, in this case the residual magnetic torque T r m b is approximately compensated by modifying the stabilizing laws in Equations (24) and (25) as follows.
m c o i l s = m m a x sat 1 m m a x b b × ( K p q v + K d ω b o b ) m ^ r m m c o i l s = m m a x sat 1 m m a x b b × K p 4 i = 1 3 ( e i × R b o T e i ) + K d ω b o b m ^ r m b
The term m ^ r m b represents an estimation of the residual dipole moment m r m b obtained by means of a Kalman filter, as reported in Appendix B.
In such a scenario, it is not appropriate to choose the settling time for the principal angle as the performance index that measures the speed of convergence. In fact, because of the presence of T d i s t b , the principal angle will not converge to 0. It will oscillate about some residual value. Thus, it more appropriate to select the Integral Time Absolute Error (ITAE) [19] as index for the speed of convergence for the principal angle, which is defined as follows.
ITAE = 0 t f i n t principal angle ( t ) d t
In ITAE, the principal angle is weighted by time so that a high value of the principal angle at later times is penalized more than that at earlier times. In this perturbed case, t f i n is set equal to 40 T o r b i t so that the spacecraft’s attitude can reach a steady-state condition. On the other hand, index in Equation (27) is still used to measure the energy required to apply the control action.
A new Monte Carlo campaign including 100 runs has been executed. The initial rotation matrix R b o ( 0 ) , initial quaternion q ( 0 ) , initial angular velocity ω b o b ( 0 ) 20 deg/s, and initial argument of the spacecraft ϕ 0 have been all randomly selected, as explained in the previous subsection.
The results of the campaign are reported in Table 5.
In the table, ITAE ¯ denotes the mean value of ITAE, and ITAE ( RM ) ITAE ( Q ) denotes the number of runs (expressed in percentage on the total number) in which ITAE ( RM ) ITAE ( Q ) .
The statistics show that in terms of both ITAE and energy consumption, quaternion feedback performs better than rotation matrix feedback.
Figure 2 displays the time histories of the principal angle for a single run obtained by using quaternion feedback and rotation matrix feedback.
It turns out that for this specific run, the ITAE obtained with quaternion feedback is equal to 3.247 1011 deg sec, whereas the ITAE is equal to 3.255 1011 deg sec when using rotation matrix feedback. Moreover, the values of E c o i l s are, respectively, equal to 2.017 104 A2 m4 s and 2.062 104 A2 m4 s.

5. Conclusions

In this paper, a numerical comparison between quaternion plus attitude rate feedback and rotation matrix plus attitude rate feedback has been performed for a case study of an Earth-pointing spacecraft equipped with magnetorquers only as torque actuators. Monte Carlo campaigns show that in both nominal and perturbed conditions quaternion feedback achieves a faster convergence to the nominal attitude. Moreover, in terms of energy consumption, the two stabilizing laws are comparable in nominal conditions. However, quaternion feedback requires less energy in perturbed conditions. Note that, by using quaternion feedback, the unwinding phenomenon can occur, and its occurrence typically results in slower convergence and higher energy consumption. Thus, the fact that, in this study, quaternion feedback performs better both in terms of speed of convergence and of energy consumption testifies that the unwinding phenomenon most likely has occurred in rare simulation runs of Monte Carlo campaigns that have been executed.

Funding

This research received no external funding.

Conflicts of Interest

The author declare no conflict of interest.

Appendix A. Gain Selection

The goal of this appendix to present a method for determining effective values for the 3 × 3 matrix gains K p and K d . The method was obtained in [20] (see also [7]).
Let us neglect disturbance torque by setting T d i s t b = 0 , and consider the linear approximation of the closed-loop system given by Equations (1) and (14) and either (24) or (25). The linear approximation with either control laws is given by the following linear time-varying system:
x ˙ = A x + B ( t ) u
which is subject to the state-feedback.
u = K x
In the above equations, x = [ q v T ( ω b o b ) T ] T .
A = 0 3 × 3 1 2 I 3 × 3 A 21 A 22 B ( t ) = 0 3 × 3 J 1 ( ( b o ( t ) ) × ) 2 K = [ K p K d ]
A 21 = 8 n 2 σ x 0 0 0 6 n 2 σ y 0 0 0 2 n 2 σ z A 22 = 0 0 n 1 σ x 0 0 0 n ( 1 + σ z ) 0 0
σ x = J y J z J x σ y = J z J x J y σ z = J x J y J z
For selecting matrix gains K p and K d (or equivalently matrix gain K), consider the aligned dipole model of the geomagnetic field for which case b o ( t ) is periodic with period T o r b i t = 2 π / n (see Equation (22)). Then, matrix gain K is chosen so that the closed-loop system given by Equations (A1) and (A2) is asymptotically stable and control (A2) minimizes the following quadratic performance index:
E 0 x T ( t ) Q x ( t ) + u ( t ) T R u ( t )
where Q is a 6 × 6 positive semidefinite matrix, R is a 3 × 3 positive definite matrix, and the expectation operator E { · } is taken over the set of possible values for the initial state x 0 , which is assumed to be a random variable having zero mean and known covariance X 0 = E { x 0 x 0 T } . Table A1 reports the values selected for those matrices as well as the initial guesses K p 0 K d 0 for K p and K d , respectively.
Table A1. Parameters for gain selection.
Table A1. Parameters for gain selection.
QR X 0 K p 0 K d 0
I 6 × 6 I 6 × 6 I 6 × 6 7 10 3 I 3 × 3 9 10 6 I 3 × 3
Thus, the approach presented in [7,20] is employed for finding a solution to the above linear quadratic problem, and the following gains are obtained.
K p = 10 3 6.9970 0.0003 0.0031 0.0001 7.0000 0.0006 0.0037 0.0003 6.9880 K d = 9 × 10 6 0 0 0 9 × 10 6 0 0 0 9 × 10 6

Appendix B. Compensation for the Residual Magnetic Torque

This appendix presents the method employed to compensate for the residual magnetic torque in Equation (23). The approach is based on obtaining an estimate m ^ r m b of m r m b and plugging it into the stabilizing laws in Equation (28).
In order to estimate m r m b , we employ the method presented in [21], which uses a Kalman filter. Consider the following model of attitude dynamics:
ω ˙ b i b = J 1 ω b i b × ( J ω b i b ) + m r m b × b b + m c o i l s × b b
which is obtained from Equations (7), (10), and (23) and in which the gravity gradient torque T g g b , aerodynamic torque T a b , and solar radiation pressure torque T s r p b are all neglected. In the above equation, m c o i l s is given by either of the stabilizing laws in Equation (28). Moreover, m r m b is modeled as constant. Thus, the following is the case.
m ˙ r m b = 0
Note that when the spacecraft attitude is in nominal conditions, ( R b o , ω b o b ) = ( I 3 × 3 , 0 ) and, consequently, ω b i b = ω o i o (see Equations (5) and (6)). We introduce Δ ω = ω b i b ω o i o and determine the equation for Δ ω ˙ , and we linearize it about Δ ω = 0 , obtaining the following.
Δ ω ˙ = J 1 ( J ω o i o ) × ( ω o i o ) × J Δ ω J 1 ( b b ) × m r m b J 1 ( b b ) × m c o i l s + w 1
Note that in the linearized equation, the zero mean process noise w 1 has been added to take into account all the approximations that were introduced. Similarly, Equation (A5) is modified as follows:
m ˙ r m b = w 2
where the zero mean process noise w 2 has been included to take into account possible time variations of the residual magnetic dipole moment.
Consider discrete times t k = k Δ t k = 0 , 1 , , where Δ t is the sampling interval, and integrate Equations (A6) and (A7) over interval [ t k t k + 1 ] . If Δ t is sufficiently small b b ( t ) , m r m b ( t ) , m c o i l s ( t ) , w 1 ( t ) , and w 2 ( t ) can be approximately considered constant over the integration interval and are equal to b b ( t k ) , m r m b ( t k ) , m c o i l s ( t k ) , w 1 ( t k ) , and w 2 ( t k ) , respectively. Thus, the following equation is obtained:
Δ ω ( t k + 1 ) m r m b ( t k + 1 ) = e A 1 Δ t A 2 J 1 ( b b ( t k ) ) × 0 3 × 3 I 3 × 3 Δ ω ( t k ) m r m b ( t k ) + A 2 J 1 ( b b ( t k ) ) × 0 3 × 3 m c o i l s ( t k ) + A 2 0 3 × 3 0 3 × 3 Δ t I 3 × 3 w 1 ( t k ) w 2 ( t k )
where the following is the case.
A 1 = J 1 ( J ω o i o ) × ( ω o i o ) × J A 2 = 0 Δ t e A 1 ξ d ξ
Moreover, since ω b i b can be measured, consider the following observation equation:
Δ ω o b s ( t k ) = I 3 × 3 0 3 × 3 Δ ω ( t k ) m r m b ( t k ) + v ( t k )
where v ( t k ) is a zero mean measurement noise. Thus, the estimation of m r m b ( t k ) denoted by m ^ r m b ( t k ) is obtained by designing a Kalman filter to estimate the state of Equation (A8) with observation Equation (A9) (see [22]).
The following data are employed to design the Kalman filter (see [12]). The sampling time is selected to Δ t = 0.1 s. The estimates are initialized to Δ ω ^ ( 0 ) = ω b i b ( 0 ) ω o i o , m ^ r m b ( 0 ) = 0 3 × 1 . The corresponding covariance matrices are set to P 10 = diag ( 10 9 10 9 10 9 ) rad 2 / s 2 and P 20 = diag ( 10 5 10 5 10 5 ) A 2 m 4 . The covariance matrices of w 1 , w 2 , and v are chosen as Q 1 = 10 13 I 3 × 3 rad 2 / s 4 , Q 2 = 10 13 I 3 × 3 A 2 m 4 / s 2 , and R = 10 8 I 3 × 3 rad 2 / s 2 respectively.
In order to evaluate the performance of the proposed estimation method, the time histories of of m ^ r m b obtained in the Monte Carlo campaign for the perturbed case using rotation matrix feedback are displayed in Figure A1. It is convenient to recall that, in the campaign, the residual magnetic dipole moment is constant and equal to m r m b = [ 0.15 0.12 0.1 ] T (see Table 3). Similar time behaviors are obtained in the Monte Carlo campaign for the perturbed case using quaternion feedback and are not reported to save space. An important aspect shown in Figure A1 is that m ^ r m b achieves steady state after few orbital periods.
Figure A1. Time histories of m ^ r m b obtained in the Monte Carlo campaign for the perturbed case using rotation matrix feedback.
Figure A1. Time histories of m ^ r m b obtained in the Monte Carlo campaign for the perturbed case using rotation matrix feedback.
Aerospace 09 00024 g0a1
In addition, Table A2 reports the statistics for m ^ r m b ( t f ) obtained from the Monte Carlo campaign for the perturbed case using either quaternion feedback or rotation matrix feedback.
Table A2. Statistics on m ^ r m b ( t f ) from the Monte Carlo campaign for the perturbed case using either quaternion feedback (Q) or rotation matrix feedback (RM). In the campaign, the residual magnetic dipole moment is constant and equal to m r m b = [ 0.15 0.12 0.1 ] T (see Table 3).
Table A2. Statistics on m ^ r m b ( t f ) from the Monte Carlo campaign for the perturbed case using either quaternion feedback (Q) or rotation matrix feedback (RM). In the campaign, the residual magnetic dipole moment is constant and equal to m r m b = [ 0.15 0.12 0.1 ] T (see Table 3).
Qmean value (A m2)standard deviation (A m2)
m ^ r m x b ( t f ) 0.140.75
m ^ r m y b ( t f ) −0.0870.83
m ^ r m z b ( t f ) −0.110.78
RMmean value (A m2)standard deviation (A m2)
m ^ r m x b ( t f ) 0.180.83
m ^ r m y b ( t f ) −0.0631.00
m ^ r m z b ( t f ) −0.130.88
The large standard deviations show that the estimation of m r m b is not accurate. Consequently, the compensation of the residual magnetic torque is not accurate either. Nevertheless, the proposed compensation is still highly beneficial. In fact, by running several simulations in perturbed conditions without operating any compensation of m r m b , it can be observed that no convergence to the nominal Earth-pointing attitude is achieved, as displayed in Figure A2, which is obtained from one of those simulations.
Figure A2. Time histories of principal angle using quaternion feedback (Q, blue solid line) and using rotation matrix feedback (RM, red dashed line) obtained in a simulation in perturbed conditions without compensation of the residual magnetic torque.
Figure A2. Time histories of principal angle using quaternion feedback (Q, blue solid line) and using rotation matrix feedback (RM, red dashed line) obtained in a simulation in perturbed conditions without compensation of the residual magnetic torque.
Aerospace 09 00024 g0a2

References

  1. Ovchinnikov, M.Y.; Roldugin, D. A survey on active magnetic attitude control algorithms for small satellites. Prog. Aerosp. Sci. 2019, 109, 100546. [Google Scholar] [CrossRef]
  2. Celani, F. Robust three-axis attitude stabilization for inertial pointing spacecraft using magnetorquers. Acta Astronaut. 2015, 107, 87–96. [Google Scholar] [CrossRef] [Green Version]
  3. Lovera, M.; Astolfi, A. Spacecraft attitude control using magnetic actuators. Automatica 2004, 40, 1405–1414. [Google Scholar] [CrossRef]
  4. Lovera, M.; Astolfi, A. Global magnetic attitude control of inertially pointing spacecraft. J. Guid. Control. Dyn. 2005, 28, 1065–1067. [Google Scholar] [CrossRef]
  5. Ovchinnikov, M.; Roldugin, D.; Penkov, V. Three-axis active magnetic attitude control asymptotical study. Acta Astronaut. 2015, 110, 279–286. [Google Scholar] [CrossRef]
  6. Wiśniewski, R.; Blanke, M. Fully magnetic attitude control for spacecraft subject to gravity gradient. Automatica 1999, 35, 1201–1214. [Google Scholar] [CrossRef]
  7. Celani, F. Gain Selection for Attitude Stabilization of Earth-Pointing Spacecraft Using Magnetorquers. Aerotec. Missili Spaz. 2021, 100, 15–24. [Google Scholar] [CrossRef]
  8. Ovchinnikov, M.; Roldugin, D.; Ivanov, D.; Penkov, V. Choosing control parameters for three axis magnetic stabilization in orbital frame. Acta Astronaut. 2015, 116, 74–77. [Google Scholar] [CrossRef]
  9. Psiaki, M.L. Magnetic Torquer Attitude Control via Asymptotic Periodic Linear Quadratic Regulation. J. Guid. Control. Dyn. 2001, 24, 386–394. [Google Scholar] [CrossRef]
  10. Yang, Y. An efficient algorithm for periodic Riccati equation with periodically time-varying input matrix. Automatica 2017, 78, 103–109. [Google Scholar] [CrossRef] [Green Version]
  11. Chaturvedi, N.A.; Sanyal, A.K.; McClamroch, N.H. Rigid-body attitude control. IEEE Control Syst. Mag. 2011, 31, 30–51. [Google Scholar]
  12. Avanzini, G.; de Angelis, E.L.; Giulietti, F. Two-timescale magnetic attitude control of Low-Earth-Orbit spacecraft. Aerosp. Sci. Technol. 2021, 116, 106884. [Google Scholar] [CrossRef]
  13. Sidi, M.J. Spacecraft Dynamics and Control; Cambridge University Press: Cambridge, MA, USA, 1997. [Google Scholar]
  14. Wie, B. Space Vehicle Dynamics and Control; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2008. [Google Scholar]
  15. de Ruiter, A.H.J.; Damaren, C.J.; Forbes, J.R. Spacecraft Dynamics and Control: An Introduction; Wiley: Chichester, UK, 2013. [Google Scholar]
  16. Wertz, J.R. (Ed.) Spacecraft Attitude Determination and Control; Kluwer Academic: Dordrecht, Holland, 1978. [Google Scholar]
  17. Alken, P.; Thébault, E.; Beggan, C.D.; Amit, H.; Aubert, J.; Baerenzung, J.; Bondar, T.; Brown, W.; Califf, S.; Chambodut, A.; et al. International geomagnetic reference field: The thirteenth generation. Earth Planets Space 2021, 73, 1–25. [Google Scholar] [CrossRef]
  18. Yakubovich, V.A.; Starzhinskii, V.M. Linear Differential Equations with Periodic Coefficients; J. Wiley & Sons: Jerusalem, Israel, 1975. [Google Scholar]
  19. Graham, D.; Lathrop, R.C. The synthesis of “optimum” transient response: Criteria and standard forms. Trans. Am. Inst. Electr. Eng. Part II Appl. Ind. 1953, 72, 273–288. [Google Scholar] [CrossRef] [Green Version]
  20. Viganò, L.; Bergamasco, M.; Lovera, M.; Varga, A. Optimal periodic output feedback control: A continuous-time approach and a case study. Int. J. Control 2010, 83, 897–914. [Google Scholar] [CrossRef]
  21. Inamori, T.; Sako, N.; Nakasuka, S. Magnetic dipole moment estimation and compensation for an accurate attitude control in nano-satellite missions. Acta Astronaut. 2011, 68, 2038–2046. [Google Scholar] [CrossRef]
  22. Anderson, B.D.; Moore, J.B. Optimal Filtering; Prentice-Hall: Englewood Cliffs, NJ, USA, 1979. [Google Scholar]
Figure 1. Time histories of principal angle using quaternion feedback (Q, blue solid line) and using rotation matrix feedback (RM, red dashed line) (single run of the Monte Carlo campaign for the nominal case).
Figure 1. Time histories of principal angle using quaternion feedback (Q, blue solid line) and using rotation matrix feedback (RM, red dashed line) (single run of the Monte Carlo campaign for the nominal case).
Aerospace 09 00024 g001
Figure 2. Time histories of principal angle using quaternion feedback (Q, blue solid line) and using rotation matrix feedback (RM, red dashed line) (single run of the Monte Carlo campaign for the perturbed case).
Figure 2. Time histories of principal angle using quaternion feedback (Q, blue solid line) and using rotation matrix feedback (RM, red dashed line) (single run of the Monte Carlo campaign for the perturbed case).
Aerospace 09 00024 g002
Table 1. Notation.
Table 1. Notation.
v i , v o , v b vector of components of Eucledian vector v with respect to frames F i , F o , and F b
respectively
v × skew-symmetric matrix corresponding to v = [ v 1 v 2 v 3 ] T (see Equation (2))
ω b o angular velocity of F b w.r.t. F o
ω b i angular velocity of F b w.r.t. F i
ω o i angular velocity of F o w.r.t. F i
b geomagnetic field at spacecraft
e 3 = [ 0 0 1 ] T
qquaternion representing rotation of F b w.r.t. F o
q v , q 4 vector part and scalar part of q
q ¯ = [ 0 0 0 1 ] T
r vector from the Earth center to the satellite center of mass
z ^ o unit Eucledian vector corresponding to the z o -axis of F o
I n × n n × n identity matrix
0 m × n m × n zero matrix
Table 2. Spacecraft, orbit, and geomagnetic field data.
Table 2. Spacecraft, orbit, and geomagnetic field data.
J x = J z (kg m2) J y (kg m2) R o r b i t (km) i n c l (deg) Ω (deg) μ i d (Wb m) θ i d (deg)
1.4162.08617021981377.71 1015171
ω e (deg/day) μ a d (Wb m) m m a x (A m2)
360.997.60 10153.5
Table 3. Disturbance torques data.
Table 3. Disturbance torques data.
m r m b (A m2) r a p b = r s p b (m) C D S D (m2) ρ (kg/m3) Φ s (W/m2)
[ 0.15 0.12 0.1 ] T [ 0.0082 0.003 0.0492 ] T 2.20.226.39 10−131367
q s s ^ i (m) A s (m2)
0.8 [ 0.578 0.578 0.578 ] T 0.33
Table 4. Statistics of the Monte Carlo campaign for the nominal case.
Table 4. Statistics of the Monte Carlo campaign for the nominal case.
control law t ¯ s (orbital period) E ¯ c o i l s (A2 m4 s)
Q13.86.20 104
RM15.7 (+14%)6.19 104 (−0.1%)
t s ( RM ) t s ( Q ) E c o i l s ( RM ) E c o i l s ( Q )
47%65%
Table 5. Statistics of the Monte Carlo campaign for the perturbed case.
Table 5. Statistics of the Monte Carlo campaign for the perturbed case.
control law ITAE ¯ (deg sec) E ¯ c o i l s (A2 m4 sec)
Q1.3 10121.21 106
RM2.1 1012 (+69%)1.36 10 6 (+12%)
ITAE ( RM ) ITAE ( Q ) E c o i l s ( RM ) E c o i l s ( Q )
30%32%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Celani, F. Quaternion versus Rotation Matrix Feedback for Spacecraft Attitude Stabilization Using Magnetorquers. Aerospace 2022, 9, 24. https://doi.org/10.3390/aerospace9010024

AMA Style

Celani F. Quaternion versus Rotation Matrix Feedback for Spacecraft Attitude Stabilization Using Magnetorquers. Aerospace. 2022; 9(1):24. https://doi.org/10.3390/aerospace9010024

Chicago/Turabian Style

Celani, Fabio. 2022. "Quaternion versus Rotation Matrix Feedback for Spacecraft Attitude Stabilization Using Magnetorquers" Aerospace 9, no. 1: 24. https://doi.org/10.3390/aerospace9010024

APA Style

Celani, F. (2022). Quaternion versus Rotation Matrix Feedback for Spacecraft Attitude Stabilization Using Magnetorquers. Aerospace, 9(1), 24. https://doi.org/10.3390/aerospace9010024

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop