Gain Selection for Attitude Stabilization of Earth-Pointing Spacecraft Using Magnetorquers

This paper considers a feedback control law that achieves attitude stabilization for Earth-pointing spacecraft using only magnetorquers as torque actuators. The control law is proportional derivative (PD)-like with matrix gains, and it guarantees asymptotic stability. The PD matrix gains are determined through the numerical solution of a periodic linear quadric regulator problem. A case study shows the effectiveness of the considered control law, and specifically of the gain selection method, in a simplified simulation scenario.


Introduction
Spacecraft attitude control can be obtained by adopting several actuation mechanisms. Among them electromagnetic actuators are widely used for generation of attitude control torques on satellites flying in low Earth orbits. They consist of planar current-driven coils rigidly placed on the spacecraft typically along three orthogonal axes, and they operate on the basis of the interaction between the magnetic dipole moment generated by those coils and the Earth magnetic field. In fact, the interaction with the Earth field generates a torque that attempts to align the total magnetic dipole moment in the direction of the field. Magnetic actuators, also known as magnetorquers, have the important limitation that control torque is constrained to belong to the plane orthogonal to the Earth magnetic field. As a result, usually actuators of a different type accompany magnetorquers to provide full three-axis control. Moreover, magnetorquers are often used for angular momentum dumping when reaction or momentum-bias wheels are employed for accurate attitude control (see [11,Chapter 7]). Lately, attitude stabilization using only magnetorquers has been considered as a feasible option especially for low-cost micro-and nano-satellites and for satellites with a failure in the main attitude control system. As a result, several papers have been dedicated to the design of stabilizing control laws in such a setting. In some of those works stability is achieved in the case of inertial pointing (see, e.g. [2,4,5,7]), whereas in other papers it is achieved for Earth (or nadir) pointing spacecraft (see, e.g. [6,8,12,16,17]). In all of the cited papers PD-like control algorithms are employed. In [2,4,5,16] the corresponding scalar control gains are determined substantially on a trialand-error basis. In [6,7,12] the scalar gains are selected so to maximize the degree of stability of the closed-loop system. In [8,17] matrix gains are employed and are determined by solving linear quadratic regulator problems.
In the present work, attitude stabilization of an Earthpointing spacecraft using only magnetorquers is considered, and a stabilizing PD-like control law is applied. Matrix gains are employed since they provide more degrees of freedom in the control action with respect to scalar gains. The matrix gain selection is performed through the use of a result on periodic linear quadratic regulators presented in [13]. The proposed approach is more promising than those in [8,17].

Coordinate Frames
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: The Earth-pointing objective is having F b aligned with F o . The notation used throughout the paper is presented in Table 1.

Attitude Kinematics
Since in this work we deal with an Earth-pointing problem, the focus will be on the kinematics of the satellite with respect to the orbital frame. Let with ‖q‖ = 1 , be the quaternion representing rotation of F b w.r.t. F o . Then, the corresponding attitude matrix that transforms vectors of coordinates in F o into vectors of coordinates in F b is given by (see [15,Sect. 5.4]) where I 3×3 is the 3 × 3 identity matrix, and the notation v The attitude kinematics is given by (see [15,Sect. 5.5

.3])
Clearly Table 1 Notation 4 Vector part and scalar part of q C(q) Attitude matrix corresponding to q (see Eq.

Attitude Dynamics
The attitude dynamics can be expressed in body frame as gg is the gravity gradient torque, T b coils is the external torque induced by magnetic coils, and T dist is a disturbance torque.
The gravity gradient torque in body coordinates is given by (see [ bi . For that purpose note that from (4) and (5) it , which corresponds to characterizing the timedependence of b i (t) and R oi (t) . By adopting the so-called dipole model of the geomagnetic field (see [14, Appendix H]), and letting R orbit denote the radius of the circular orbit, obtain: In equation (15), m is the total dipole strength, 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 magnetic dipole expressed in F i , which is set equal to where m is the dipole coelevation, e = 360.99 deg/day is the Earth average rotation rate, and 0 is the right ascension of the dipole at time t = 0 . We use m = 7.746 10 15 Wb m and m = 170.0 deg as reported in [9]. To characterize the time dependence of b i (t) in (15), one needs to determine an expression for r i (t) which is the spacecraft position vector resolved in F i . Define a coordinate system x p , y p in the orbital plane, whose 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 where 0 is the argument of the spacecraft at time t = 0 . Let incl be the orbit inclination and let denote the Right Ascension of the Ascending Node (RAAN) of the orbit (see [11,Sect. 2.6.2]). Then, the coordinates of the satellite center of mass in the inertial frame can be obtained as follows where is the rotation matrix corresponding to a rotation around the x-axis of magnitude and is the rotation matrix corresponding to a rotation around the z-axis of magnitude (see [11,Sect. 2.6.2]). Combining Eqs (15) -(20), the expression of b i (t) can be easily obtained. Moreover, it holds that Then, by using Eq. (14), 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. In addition, if coelevation of the Earth magnetic dipole is approximated to m = 180 deg, which corresponds to having the Earth magnetic dipole aligned with the Earth rotation axis, then expression of b o (t) simplifies as follows Thus, in the latter simplified scenario b o (t) is periodic with period T orbit = 2 ∕n. (17) x p (t) = R orbit cos(nt + 0 ) y p (t) = R orbit sin(nt + 0 ),

Control Law
The control objective is driving the spacecraft so that F b is aligned with F o . From (1) it follows that C(q) = I 3×3 for q = [q v q 4 ] T = ±q where q = [0 0 0 1] T . Thus, the control objective can be equivalently stated as having that q v → 0 and b bo → 0. In this paper the following control law is considered with K p and K d being 3 × 3 matrices. It is a PD-like control law in which the cross product with b b guarantees that m coils is orthogonal to b b . This is useful since a component of m coils parallel to b b does not provide any contribution to T coils (see Eq. (9)).
Note that employing matrix gains K p and K d provides more degrees of freedom to the control action with respect to adopting scalar gains as in [2,[4][5][6][7]. Those additional degrees can be beneficial to the stabilizing action considering that, in the general case, the dynamics about each single body axis are different and are coupled. On the other hand, selecting appropriate values for two 3 × 3 matrix gains is clearly more challenging than doing it for two scalar gains.

Gain Selection
The goal of this section is presenting a systematic method for determining effective values for the 3 × 3 matrix gains K p and K d .
Neglect disturbance torque by setting T dist = 0 , and rewrite control law in Eq. (22) in the following equivalent forms The resulting closed-loop system is given by Eqs. On the latter set the following holds Then, the restriction is given by the reduced order system obtained considering the following equation along with the equations obtained from Eqs. (13) and (23) replacing C(q) with C v (q v ) which is given by Considering that b o changes with time (see Eq. (14)), the linear approximation of the latter restriction around (q v , b bo ) = (0, 0) can be represented as the following linear time-varying system (see also [8,17]) subject to state-feedback For selecting matrix gains K p and K d (or equivalently matrix gain K), set the coelevation angle of the magnetic dipole that generates the geomagnetic field to m = 180 deg in which case b o (t) is periodic with period T orbit = 2 ∕n (see Eq. (21)). As a consequence, system in Eq. (28) is periodic with period T orbit , too. Then, matrix gain K is chosen so that the closedloop system given by Eqs. (28) (29) is asymptotically stable and control (29) minimizes the following quadratic performance index 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 initial state x 0 which is assumed to be a random variable having zero mean and known covariance X 0 = E{x 0 x T 0 } [13]. Note that parameter X 0 can be used by the designer that has some a priori knowledge on the stochastic distribution of the initial state x 0 . If this type of information is not available, usually one sets X 0 = 2 0 I 6×6 in which 0 > 0 is a design parameter. Clearly, matrices Q and R represent design parameters, too.

Remark 1
The periodic linear quadratic regulator problem considered in this work compares favorably to the linear quadratic regulator problem formulated in [8]. In fact, in [8] the time horizon is limited to only one orbital period whereas the effective maneuvering time can easily involve several orbital periods. In addition, the cost function adopted in [8] includes a penalty of the state at time t = T orbit which does not seem to have a clear physical meaning. Moreover, compared to the linear quadratic method presented in [17], the proposed approach has the advantage of adopting constant matrix gains K p and K d rather than time-varying ones which simplifies the practical implementation of the control law.

Solution to the Periodic Linear Quadratic Regulator Problem
Reference [13] shows that a matrix K that solves the periodic linear quadratic regulator problem just formulated, can be determined by solving an equivalent simpler optimization problem which can be formulated as follows. Consider the dynamic matrix of the closed-loop system in Eqs. (28) (29) which is equal to Ā (t) = A − B(t)K and denote by Ā(t, ) the corresponding transition matrix which is the solution of the following matrix differential equation Moreover, denote by Ā the monodromy matrix associated to Ā (t) which is given by Ā =Ā(T orbit , 0) . Then, it can be shown that K is an optimal point of the following optimization problem where P 0 is the solution of the following discrete-time Lyapunov equation in which W is obtained by solving the following periodic Lyapunov differential equation with initial condition Z(0) = 0 6×6 and setting W = Z(T orbit ) . Matrix Q appearing in the previous differential equation is given by Q = Q + K T RK . The considered optimization problem is well defined if all eigenvalues of Ā are within the open unit disk in the complex plane which is equivalent to requiring that closed-loop system in Eqs. (28) (29) is asymptotically stable [1]. In fact, the latter condition guarantees that the solution P 0 of Eq. (30) exists and is unique. As a result, from a numerical point of view it is necessary to start the optimization procedure with an initial guess for K that renders the closed-loop system in Eqs. (28) (29) asymptotically stable. 1

Case Study
For simulation consider Tigrisat, a 3U Cubesat (see Fig. 1) designed and manufactured at the School of Aerospace Engineering of Sapienza University of Rome. Its principal moments of inertia are as follows 2 The satellite follows a circular orbit having incl = 97 deg and altitude h = 629 km which corresponds to the following orbital period T orbit = 5832 sec [12]. The Right Ascension of the Ascending Node is equal to = 68.5 deg, whereas the argument 0 of the spacecraft at time t = 0 (see Eq. (17)) has been randomly selected and set equal to 0 = 1.60 rad. The approach described in Sect. 3 has been adopted for determining the matrix gains K p and K d . The following design parameters have been selected for that purpose Q = X 0 = I 6×6 , R = I 3×3 . The optimization procedure has been initialized by choosing K p = 300 I 3×3 and K d = 1.8 ⋅ 10 4 I 3×3 . The initializing values have been chosen on a trial-and-error basis and are such that the characteristic multipliers of Ā (t) have all negative real part. The optimization procedure determines the following optimal gains It can be verified that all characteristic exponents of closedloop system in Eqs. (28) (29) have negative real part. Thus, the latter linear system is asymptotically stable. As a result, (q, b bo ) = (q, 0) is a (locally) exponentially stable equilibrium of closed-loop system in Eqs. (3) (13) (23).
Consider an initial condition characterized by attitude equal to the target attitude q(0) =q , and by the following initial angular velocity b bo (0) = [10 −3 10 −3 10 −3 ] T rad/s due for example to an impact with an object.
Two simulation are performed for the closed-loop system. In the first one nominal conditions are considered by employing the nominal value J nom = diag[J x J y J z ] for the spacecraft inertia matrix (see Eq. (31)), and by setting the dipole coelevation as m = 180 • , and the disturbance torque as T dist = 0 . In the second simulation the following perturbations are introduced. For the spacecraft inertia matrix the following perturbed value is used The dipole coelevation is set as m = 170 • (see [9]), and the disturbance torque is set as T dist = m 0 × b b which represents the torque due to a residual dipole. A realistic magnitude for the residual dipole of a 3U cubesat is 3 ⋅ 10 −4 A m 2 . Moreover, since the residual dipole is typically normal to the electronic boards, it is mostly aligned along the longitudinal axis. Note that the residual dipole torque is usually the largest among the disturbance torques acting on a nanosatellite flying in a low Earth orbit. Thus, for simplicity other disturbance torques such as those due to aerodynamic drag and solar radiation pressure, are not included in the simulation.
Plots of the Euler angles are included even if redundant since they are useful to quantify the error in steady-state. In nominal conditions convergence to the desired attitude is achieved within 5 orbital periods. Simulations with perturbed conditions show that the proposed attitude control system is able to compensate for the considered perturbations to a significant extent. As a matter of fact, in steadystate the maximum roll, pitch, and yaw error are, respectively, about 2°, 4°, and 5°.
The time histories of m coils in Fig. 5 show that the required maximum magnitude for the magnetic dipole generated by a single coil is lower that 4 10 -3 A m 2 . Such a magnitude can be easily achieved by Tigrisat that can generate dipoles up to 0.22 A m 2 for m coils x and 0.696 A m 2 for m coils y and m coils z 3 The scope of the considered numerical study is limited. In fact, a more accurate study could consider a more realistic model for the geomagnetic field, a time-varying residual magnetic dipole, additional disturbance torques such as those due to aerodynamic drag and solar radiation pressure, and errors in attitude determination and in angular velocity sensors (see [3]). However, such more complete study goes beyond the scope of this work whose focus is mostly theoretical.

Conclusion
This work studies the problem of attitude stabilization of an Earth-pointing spacecraft using only magnetorquers. A PD-like control law with matrix gains is employed. The matrix gains selection is based on the numerical solution of a periodic linear quadratic problem. The periodic linear quadratic problem adopted in this paper compares favorably under several aspects to other linear quadratic regulator problems employed in existing works for gain selection. A case study illustrates the effectiveness of the proposed method considering a simplified simulation scenario. More accurate simulation environments are not considered since that would go beyond the scope of this work whose focus is mostly theoretical.