AIAA 2016-0086 AIAA SciTech 4-8 January 2016, San Diego, California, USA AIAA Guidance, Navigation, and Control Conference

Spacecraft Attitude Stabilization using Magnetorquers with Separation between Measurement and Actuation Fabio Celani∗

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

Sapienza University of Rome, 00138 Rome, Italy In actual implementations of magnetic control laws for spacecraft attitude stabilization, the time in which Earth magnetic field is measured must be separated from the time in which magnetic dipole moment is generated. In this work we present a magnetic attitude stabilization law that takes into account of the latter separation. A proportional-derivative feedback is presented, and a design method for determining its parameters is obtained. The proposed method is validated through a case study.

Nomenclature Bb, Bi C incl I J mcoils n q = [qvT q4 ]T R T ∆ ∆a ∆m ω Ω

geomagnetic field in Fb and Fi coordinates respectively, T attitude matrix of Fb w.r.t. Fi orbit’s inclination, rad identity matrix spacecraft inertia matrix, kg m2 vector of magnetic dipole moments generated by three coils, A m2 orbital rate, rad/sec quaternion representation of attitude of Fb w.r.t. Fi orbital radius, m external torque, N m = ∆m + ∆a , s actuation time, s measurement time, s angular rate of Fb w.r.t. Fi in Fb coordinates, rad/s Right Ascension of the Ascending Node (RAAN), rad

I.

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 moment generated by those coils and the Earth’s magnetic field; in fact, the interaction with the Earth’s field generates a torque that attempts to align the total magnetic 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’s magnetic field. As a result, usually a different type of actuator accompanies 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 [1, Chapter 7]). Lately, attitude stabilization using only magnetorquers has been considered as a feasible option especially for low-cost micro-satellites and for satellites with a failure in the main attitude control system. In such scenario many control laws have been designed, and a survey of various approaches adopted can be found in ref.;2 in particular, Lyapunov-based design has been adopted in ref.’s.3–7 ∗ Assistant

Professor, School of Aerospace Engineering, Via Salaria 851.

1 of 15 American Institute of Aeronautics and Astronautics Copyright © 2015 by Fabio Celani. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

All control laws presented in the cited works are continuous-time; thus, in principle they require a continuous-time measurement of the geomagnetic field and a continuous-time generation of coils’ currents. However, in practical implementations the time in which Earth’s magnetic field is measured (measurement time) has to be separated from the time in which coils’ currents are generated (actuation time). The latter separation is necessary because currents flowing in the spacecraft’s magnetic coils generate a local magnetic field that impairs an accurate measurement of the geomagnetic field; as a result, when the Earth’s magnetic field is being measured no currents must flow in the coils, and consequently no magnetic actuation is possible; on the other hand, when magnetic actuation is active, it is not possible to obtain accurate measurements of the geomagnetic field. Consequently, in practical applications a periodic switch between measurement time and actuation time is implemented. Magnetic control laws compatible with the latter constraint could be obtained by discretizing continuos-time control algorithms; however, that approach could lead to a degradation of the performances of the control system. Then, in order to overcome such issue, it is important to design control laws that take into account of the separation constraint directly during the design phase. In the present work, a proportional-derivative feedback that satisfies the previous constraint is designed.

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

A.

Notations

For x ∈ Rn , kxk denotes the Euclidean norm of x. For a ∈ R3 , a× represents the skew symmetric matrix   0 −a3 a2   a× =  a3 (1) 0 −a1  −a2 a1 0 so that for b ∈ R3 , multiplication a× b is equal to the cross product a × b.

II.

Modeling and control algorithm

In order to describe the attitude dynamics of an Earth-orbiting rigid spacecraft, and in order to represent the geomagnetic field, it is useful to introduce the following reference frames. 1. Geocentric Inertial Frame Fi . A commonly used inertial frame for Earth orbits is the Geocentric Inertial Frame, whose origin is in the Earth’s center, its xi axis is the vernal equinox direction, its zi axis coincides with the Earth’s axis of rotation and points northward, and its yi axis completes an orthogonal right-handed frame (see [1, Section 2.6.1]). 2. Spacecraft body frame Fb . The origin of this right-handed orthogonal frame attached to the spacecraft, coincides with the satellite’s center of mass; its axes are chosen so that the (inertial) pointing objective is having Fb aligned with Fi . Since the pointing objective consists in aligning Fb to Fi , the focus will be on the relative kinematics and dynamics of the satellite with respect to the inertial frame. Let q = [q1 q2 q3 q4 ]T = [qvT q4 ]T with kqk = 1 be the unit quaternion representing rotation of Fb with respect to Fi ; then, the corresponding attitude matrix is given by C(q) = (q42 − qvT qv )I + 2qv qvT − 2q4 qv× (2) (see [8, Section 5.4]). Let 1 W (q) = 2

"

q4 I + qv× −qvT

# (3)

Then the relative attitude kinematics is given by q˙ = W (q)ω

(4)

where ω ∈ R3 is the angular rate of Fb with respect to Fi resolved in Fb (see [8, Section 5.5.3]). The attitude dynamics in body frame can be expressed by J ω˙ = −ω × Jω + T 2 of 15 American Institute of Aeronautics and Astronautics

(5)

where J ∈ R3×3 is the spacecraft inertia matrix, and T is the external torque expressed in Fb (see [8, Section 6.4]). The spacecraft is equipped with three magnetic coils aligned with the Fb axes which generate the magnetic attitude control torque Tcoils = mcoils × B b = −B b× mcoils (6) where mcoils ∈ R3 is the vector of magnetic dipole moments for the three coils, and B b is the geomagnetic field at spacecraft expressed in body frame Fb . Moreover, let Tdist ∈ R3 be the sum of disturbance torques acting on the spacecraft such as gravity gradient torque and torque due to magnetic dipole residual. Let B i be the geomagnetic field at spacecraft expressed in inertial frame Fi . Note that B i varies with time at least because of the spacecraft’s motion along the orbit. Then B b (q, t) = C(q)B i (t)

(7)

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

which shows explicitly the dependence of B b on both q and t. Grouping together equations (4) (5) (6) the following nonlinear time-varying system is obtained q˙ J ω˙

= W (q)ω = −ω × Jω − B b (q, t)× mcoils + Tdist

(8)

in which mcoils is the control input, and Tdist is a disturbance input which will be neglected for control design. In order to analyze and design control algorithms, it is important to characterize the time-dependence of B b (q, t) which is the same as characterizing the time-dependence of B i (t). Assume that the orbit is circular of radius R; then, adopting the so called dipole model of the geomagnetic field (see [9, Appendix H]) we obtain µm ˆ i )T rˆi (t))ˆ ri (t) − m ˆ i] (9) B i (t) = 3 [3((m R In equation (9), µm = 7.746 1015 Wb m is the total dipole strength (see ref.10 ), ri (t) is the spacecraft’s position vector resolved in Fi , and rˆi (t) is the vector of the direction cosines of ri (t); finally m ˆ i is the vector of direction cosines of the Earth’s magnetic dipole moment expressed in Fi . Here we set m ˆ i = [0 0 − 1]T ◦ which corresponds to having the dipole’s coelevation angle equal to 180 . Even if a more accurate value for that angle would be 170◦ (see ref.10 ), here we approximate it to 180◦ since this will substantially simplify future symbolic computations. Equation (9) shows that in order to characterize the time dependence of B i (t) it is necessary to determine an expression for ri (t) which is the spacecraft’s position vector resolved in Fi . Define a coordinate system xp , yp in the orbital’s plane whose origin is at Earth’s center; then, the position of satellite’s center of mass is clearly given by xp (t) = R cos(nt + φ0 ) (10) y p (t) = R sin(nt + φ0 ) where n is the orbital rate, and φ0 an initial phase. Then, coordinates of the satellite in inertial frame Fi can be easily obtained from (10) using an appropriate rotation matrix which depends on the orbit’s inclination incl and on the value Ω of the Right Ascension of the Ascending Node (RAAN) (see [1, Section 2.6.2]). Plugging into (9) the equations of the latter coordinates, an explicit expression for B i (t) can be obtained; it can be easily checked that B i (t) is periodic with period 2π/n. Consequently system (8) is a periodic nonlinear system with period 2π/n. As stated before, the control objective is driving the spacecraft so that Fb is aligned with Fi . From (2) it follows that C(q) = I for q = [qvT q4 ]T = ±¯ q where q¯ = [0 0 0 1]T . Thus, the objective is designing control strategies for mcoils so that qv → 0 and ω → 0. Ref.7 presents a stabilizing proportional-derivative control law which was obtained as modification of those in ref.’s5, 6 and is given by mcoils (t) = (B b (q(t), t)× )T (2 kp qv (t) + kd ω(t))

(11)

It is shown that picking kp > 0, kd > 0, and choosing  > 0 and small enough, local exponential stability of (q, ω) = (¯ q , 0) is achieved.

3 of 15 American Institute of Aeronautics and Astronautics

Feedback control law (11) seems implementable in practice since B b , qv , and ω can be measured using magnetometers, attitude sensors, and attitude rate sensors respectively. However, as explained in the introduction, the time in which B b is measured should be separated from the time in which magnetic control action is applied; in fact, magnetic torque is obtained by generating currents flowing through magnetic coils, and those currents create a local magnetic field which impairs an accurate measurement of B b . As a result, in order to take into account of the previous constraint, we should consider only magnetic feedback laws of the following type; in a first interval of length ∆m , B b is measured (along with q and ω) and mcoils is set to 0; then, in a successive interval of length ∆a , mcoils is generated and held constant because the measure of B b cannot be updated during that interval; the operations of measurement and actuation are repeated periodically. Thus, the considered scenario is represented by the time diagram in Fig. 1

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

0

...

actuation

measure ∆m

∆m + ∆ a

...

measure j(∆m + ∆a )

actuation

j(∆m + ∆a ) + ∆m (j + 1)(∆m + ∆a )

t

Figure 1. Time diagram (j=0,1,2,. . . ).

Note that the value of ∆m is determined mainly by the following physical fact; in order to measure B b , no current must flow in the magnetic coils; then, when the current in the magnetorquers is switched off, it is necessary to wait for the dissipation of existing currents through the coils before obtaining an accurate measurement of B b ; usually the latter dissipation time is much larger than the time necessary to read and process data from magnetometers, attitude sensors and attitude rate sensors; as a consequence, ∆m is mostly determined by the dissipation phenomenon just described. Ref.11 considers the case in which ∆m = 0 which corresponds to having ∆m negligible with respect to ∆a ; here, we study situations in which that does not occur and consequently ∆m > 0. In order to simplify forthcoming expressions let us define ∆ , ∆m + ∆a . A proportional-derivative feedback law that fulfills the separation requirement between measurement and actuation can be obtained by modifying (11) as follows

mcoils (t) =

    0

j∆ ≤ t < j∆ + ∆m

(B b (q(j∆ + ∆m ), j∆ + ∆m )× )T    2 ( kp qv (j∆ + ∆m ) + kd ω(j∆ + ∆m ))

j∆ + ∆m ≤ t < (j + 1)∆ j = 0, 1, 2, . . . (12)

III.

Feedback Design

In this section we will focus on feedback law (12) and will derive conditions on parameters kp , kd , , and intervals ∆m and ∆a which ensure that equilibrium (q, ω) = (¯ q , 0) is locally exponentially stable for closed-loop system (8) (12). In order to derive those conditions, it suffices considering the restriction of closed-loop system (8) (12) to the open set S3+ × R3 where S3+ = {q ∈ R4 | kqk = 1, q4 > 0}

(13)

Since on the latter set q4 = (1 − qvT qv )1/2 then the considered restriction is given by the interconnection of the following systems q˙v (t) = Wv (qv (t))ω(t) J ω(t) ˙ = −ω(t)× Jω(t) − [Cv (qv (t))B i (t)]× mcoils (t)

4 of 15 American Institute of Aeronautics and Astronautics

(14)

mcoils (t) =

    0

j∆ ≤ t < j∆ + ∆m

(Cv (qv (j∆ + ∆m ))B i (j∆ + ∆m )× )T    2 ( kp qv (j∆ + ∆m ) + kd ω(j∆ + ∆m ))

j∆ + ∆m ≤ t < (j + 1)∆ j = 0, 1, 2, . . . (15)

where

i
1/2 1h 1 − qvT qv I + qv× 2

1/2 × T Cv (qv ) = 1 − 2qv qv I + 2qv qvT − 2 1 − qvT qv qv

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

Wv (qv ) =

and where (7) has been used. It is simple to show that if (qv , ω) = (0, 0) is a locally exponentially stable equilibrium for (14) (15), then (q, ω) = (¯ q , 0) is a locally exponentially stable equilibrium for (8) (12). It will be shown next that the closed-loop system (14) (15) can be expressed as hybrid system of the type considered in ref.12 For that purpose define the following time sequence  k   ∆ k = 0, 2, 4, . . .  2 (16) tk =    k − 1 ∆ + ∆m k = 1, 3, 5, . . . 2 which represent the instants that appear in Fig. 1. Introduce the following discrete-time function p(k) =

1 (1 − (−1)k ) k = 0, 1, 2, . . . 2

(17)

Note that p(k) is equal to 0 for k even and to 1 for k odd. Then, time sequence (16) can be expressed by the following alternative form which will be useful for the sequel   k k−1 tk = (1 − p(k)) ∆ + p(k) ∆ + ∆m k = 0, 1, 2, . . . (18) 2 2 By using the previous expressions, it follows that equation (15) can be written in the following compact form mcoils (t) = p(k)(Cv (qv (tk )B i (tk )× )T (2 kp qv (tk ) + kd ω(tk ))

tk ≤ t < tk+1

k = 0, 1, 2, . . .

(19)

Note then that the closed-loop system (14) (19) is a nonlinear time-varying hybrid system of the type considered in ref.12 Next, consider the linear approximation of (14) (19) about the equilibrium (qv , ω) = (0, 0) as defined in ref.12 The latter approximation is given by q˙v (t) = 21 ω(t) J ω(t) ˙ = −B i (t)× mcoils (tk )

tk ≤ t < tk+1

mcoils (tk ) = p(k)(B i (tk )× )T (2 kp qv (tk ) + kd ω(tk )) i

(20)

k = 0, 1, 2, . . .

(21)

12

Note that since B is bounded (see (9)), then Theorem II.1 in ref. applies to the nonlinear time-varying hybrid system (14) (19); consequently, (qv , ω) = (0, 0) is a locally exponentially stable equilibrium for (14) (19) if and only if it is an exponentially stable equilibrium for (20) (21). Next consider the continuous-time system (20) and sample its state [qv (t)T ω(t)T ]T at t = tk . Setting ∗ qv (k) , qv (tk ), ω ∗ (k) , ω(tk ), and m∗coils (k) , mcoils (tk ), the following discrete-time system is obtained qv∗ (k + 1)

=

ω ∗ (k + 1)

=

1 qv∗ (k) + a(k, ∆m , ∆a )ω ∗ (k) − J −1 G1 (k, ∆m , ∆a )m∗coils (k) 2 ω ∗ (k) − J −1 G2 (k, ∆m , ∆a )m∗coils (k)

5 of 15 American Institute of Aeronautics and Astronautics

k = 0, 1, 2 . . .

(22)

where a(k, ∆m , ∆a ) , G1 (k, ∆m , ∆a ) , G2 (k, ∆m , ∆a ) ,

tk+1 − tk Z tk+1 1 (tk+1 − τ )B i (τ )× dτ 2 tk Z tk+1 B i (τ )× dτ

(23) (24) (25)

tk

Note that the explicit dependence of a, G1 , and G2 from k, ∆m , and ∆a can be obtained by replacing tk+1 and tk with the corresponding expressions that can be derived using (18). Defining B i∗ (k) , B i (tk ), and adopting the previously introduced notation, state-feedback (21) reads as

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

m∗coils (k) = p(k)(B i∗ (k)× )T (2 kp qv∗ (k) + kd ω ∗ (k))

k = 0, 1, 2, . . .

(26)

By a simple adaptation of Proposition 7 in ref.,13 it follows that if the linear time-varying discrete-time system (22) (26) is exponentially stable, then the linear time-varying hybrid system (20) (21) is exponentially stable, too. Based on the previous discussion, our objective has become studying stability of linear time-varying discrete-time system (22) (26). For that purpose it is useful to perform the following change of variables z1 (k) = qv∗ (k), z2 (k) = ω ∗ (k)/ obtaining z1 (k + 1)

=

z2 (k + 1)

=

1 z1 (k) +  a(k, ∆m , ∆a )z2 (k) − 2 p(k)J −1 G1 (k, ∆m , ∆a )(B i∗ (k)× )T (kp z1 (k) + kd z2 (k)) 2 z2 (k) − p(k)J −1 G2 (k, ∆m , ∆a )(B i∗ (k)× )T (kp z1 (k) + kd z2 (k)) (27)

Let M (k, ∆m , ∆a ) , p(k)G2 (k, ∆m , ∆a )(B i∗ (k)× )T

(28)

and consider the following averages aav (s, ∆m , ∆a ) ,

Mav (s, ∆m , ∆a ) ,

1 N →∞ N lim

1 N →∞ N lim

k=s+N X

a(k, ∆m , ∆a )

k=s+1 k=s+N X

M (k, ∆m , ∆a )

k=s+1

with s = 0, 1, 2, . . .. Both averages have been computed symbolically using MathematicaTM . It turns out that aav is independent of s and is equal to aav (∆m , ∆a ) =

∆m + ∆ a 2

Symbolic computations show also that Mav is independent of both s and ∆m . Its long expression is here omitted to save space. From equation (25), it is easy to see that for k odd G2 (k, ∆m , 0) = 0; thus, taking into account that for k even p(k) = 0, from (28) we obtain that M (k, ∆m , 0) = 0 . Then, it occurs that Mav (∆a ) vanishes at ∆a = 0, and consequently it can be expressed as Mav (∆a ) = Lav (∆a )∆a where Z Lav (∆a ) , 0

1

(29)

∂Mav (u∆a )du ∂∆a

The expression of Lav has also been computed symbolically using MathematicaTM , but its long expression is here omitted to save space. It turns out that L0av , lim∆a →0 Lav (∆a ) is symmetric, and that matrix L0av is positive definite if and only if the orbit is not equatorial, that is if its inclination satisfies incl 6= 0. Thus, we make the following hypothesis. 6 of 15 American Institute of Aeronautics and Astronautics

Assumption 1. The spacecraft’s orbit satisfies condition L0av > 0. Remark 1. In the special case of polar orbit (incl = 90◦ ) with zero RAAN (Ω = 0), the expression of Lav simplifies as follows   9 sin(2n∆a ) 18 sin(n∆a )2 4+ 0 −   n∆a n∆a    9 sin(2n∆a ) µ2    (30) 0 2 2 + Lav (∆a ) = 0   6   32R n∆a  18 sin(n∆ )2  9 sin(2n∆a ) a 0 n∆a n∆a The general fact that with incl 6= 0, matrix L0av , lim∆a →0 Lav (∆a ) is symmetric and positive definite can be easily checked in this special case.

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

At this point, consider the average system of (27) as defined in ref.14 which is given by z1 (k + 1)

=

z2 (k + 1)

=

∆m + ∆ a z2 (k) 4 z2 (k) − ∆a J −1 Lav (∆a )(kp z1 (k) + kd z2 (k)) z1 (k) + 

(31)

It can be verified that all assumptions of [14, Theorem 2.2.2] are fulfilled by system (27). Thus, the following proposition holds true. Proposition 1. If the average system (31) is exponentially stable for all 0 <  ≤ 0 , then there exists 0 < 1 ≤ 0 , such that system (27) is exponentially stable for all 0 <  ≤ 1 . Rewrite (31) in the following compact form z(k + 1) = z(k) + ∆a A(∆m , ∆a )z(k) with A(∆m , ∆a ) = A1 (∆a ) + where



∆m A2 ∆a 1 4I

0

A1 (∆a ) =  −kp J

−1

Lav (∆a ) −kd J

and



0

1 4I

0

0

A2 = 

(32)

−1

(33)  

(34)

Lav (∆a )

 

(35)

For matrix A1 the following lemma holds true. Lemma 1. Assume that kp > 0, kd > 0 and Assumption 1 is satisfied. Then there exists ∆∗a > 0 such that A1 (∆a ) is a Hurwitz matrix for all 0 < ∆a < ∆∗a . Proof. First, it will be shown that A01 , lim∆a →0 A1 (∆a ) is a Hurwitz matrix. For that purpose it is useful to introduce the continuous-time system w˙ = A01 w which in expanded form reads as follows w˙ 1 J w˙ 2

1 w2 4 = −L0av (kp w1 + kd w2 ) =

(36)

Introduce the Lyapunov function for (36) 1 V1 (w1 , w2 ) = 2kp w1T L0av w1 + w2T Jw2 2 then V˙ 1 (w1 , w2 ) = −kd w2T L0av w2 By La Salle’s invariance principle [15, Theorem 4.4], it follows that (36) is exponentially stable, and thus A01 is Hurwitz. Then, by continuity of A1 (∆a ) with respect to ∆a , the thesis follows immediately. 7 of 15 American Institute of Aeronautics and Astronautics

Remark 2. In the general case, the value of ∆∗a can be determined numerically; in fact, once we fix the orbital parameters (R, incl, and Ω which appear in Lav (∆a )), the inertia matrix J, and the controller’s gains kp and kd , it suffices studying the behavior with ∆a of the maximum of the real parts of the eigenvalues of A(∆a ). Moreover, an analytical expression of ∆∗a can be determined in the very special situation of polar orbit (incl = 90◦ ) and isoinertial spacecraft (J = J0 I). First note that in general, because of the symmetry of the scenario, attitude stability properties are independent of the value Ω of RAAN; as a result, there is no loss of generality in setting Ω = 0. Moreover, in the special case under consideration of polar orbit, the simplified expression of Lav (∆a ) reported in (30) applies. Then, it is possible to obtain an analytical expression of ∆∗a computing the characteristic polynomial of (34) having set J = J0 I, and applying Routh’s criterion. Standard but lengthy symbolic computations show that in the considered particular setting, matrix A(∆a ) is Hurwitz for all 0 < ∆a < 0.11/n. Then, 0.11/n, which corresponds to about one sixtieth of the orbital period, can be used as a first initial guess for ∆∗a even in more general situations. For average system (32) the following stability result holds true.

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

Theorem 1. Assume that Assumption 1 is satisfied. First fix kp > 0 and kd > 0, then fix 0 < ∆a < ∆∗a . Consider symmetric matrix P1 > 0 such that P1 A1 (∆a ) + A1 (∆a )T P1 = −I and let ∆∗m ,

∆a kP1 A2 + AT 2 P1 k

(37)

(38)

If 0 < ∆m < ∆∗m , then, there exists 0 > 0, such that system (32) is exponentially stable for all 0 <  ≤ 0 . Proof. First it will be shown that picking kp > 0, kd > 0, ∆a and ∆m as stated in the hypothesis, it follows that matrix A(∆m , ∆a ) (see (33)) is Hurwitz. Note that the choice of kp > 0, kd > 0 and ∆a guarantees that A1 (∆a ) is Hurwitz; then, consider the continuos time system   ∆m A2 w w˙ = A(∆m , ∆a )w = A1 (∆a ) + ∆a with the candidate Lyapunov function V2 (w) = wT P1 w. Since P1 is the solution of (37) we obtain    
∆m ∆m P1 A 2 + A T P w ≤ − 1 − V˙ 2 (w) = wT P1 A1 (∆a ) + A1 (∆a )T P1 + kwk2 1 2 ∆a ∆∗m Thus, if 0 < ∆m < ∆∗m , it follows that V˙ 2 is negative definite and consequently A(∆m , ∆a ) is Hurwitz. Then, consider symmetric matrix P > 0 such that P A(∆m , ∆a ) + A(∆m , ∆a )T P = −I

(39)

and let V3 (z) = z T P z be a candidate Lyapunov function for discrete-time system (32). Then, the following holds ∆V2 (z) , (z + ∆a A(∆m , ∆a )z)T P (z + ∆a A(∆m , ∆a )z) − z T P z ≤ ∆a (−1 + ∆a kA(∆m , ∆a )T P A(∆m , ∆a )k)kzk2 As a result, setting 0 =

1 2∆a kA(∆m , ∆a )T P A(∆m , ∆a )k

(40)

we obtain that system (32) is exponentially stable for all 0 <  ≤ 0 . Based on the previous theorem, fix kp > 0, kd > 0, ∆a and ∆m as indicated in the theorem, and consider the corresponding value of 0 given by (40); by Proposition 1 it follows that there exists 0 < 1 ≤ 0 such that system (27) is exponentially stable for all 0 <  ≤ 1 . In conclusion, from Theorem 1, Proposition 1, and the preceding discussion, the following main proposition is obtained.

8 of 15 American Institute of Aeronautics and Astronautics

Proposition 2. Consider the magnetically actuated spacecraft described by (8) and apply state-feedback (12) with kp > 0 and kd > 0. Then, under Assumption 1, there exists ∆∗a > 0 such that fixed first 0 < ∆a < ∆∗a and then 0 < ∆m < ∆∗m with ∆∗m given by (38), there exists 1 ≤ 0 , with 0 given by (40), such that for all 0 <  ≤ 1 equilibrium (q, ω) = (¯ q , 0) is locally exponentially stable for (8) (12).

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

Remark 3. Usually in practical situations time ∆m is not a design parameter for the control engineer; in fact, as explained before, ∆m is basically determined by the time that coils take to dissipate existing currents once the latter are switched off; as a result, ∆m is a parameter intrinsic to the controlled plant and not a controller’s parameter. As a consequence, it might seem that the previous proposition is not useful for control design since it seems to imply that ∆m can be chosen by control designer. Actually, the previous proposition can be used to do control design as follows. First fix gains kp > 0 and kd > 0; then, determine ∆∗a by using Lemma 1; for each 0 < ∆a < ∆∗a determine the corresponding ∆∗m (see Theorem 1). By plotting ∆∗m versus ∆a , we can graphically determine the set S of values of ∆a for which it occurs that ∆∗m ≥ ∆m . Thus, fix ∆a ∈ S, and finally choose  ≤ 0 sufficiently small so that stability is achieved. In the next section, the described designed procedure will be illustrated through a case study. Remark 4. The obtained stability results hold even if saturation on magnetic dipole moments is taken into account by replacing control (12) with

mcoils (t) =

 0   

j∆ ≤ t < j∆ + ∆m 

mcoils max sat mcoils1 max (B b (q(j∆ + ∆m ), j∆ + ∆m )× )T  
 2 ( kp qv (j∆ + ∆m ) + kd ω(j∆ + ∆m ))

j∆ + ∆m ≤ t < (j + 1)∆ j = 0, 1, 2, . . . (41)

where mcoils max is the saturation limit on each dipole moment, and sat : R3 → R3 is the standard saturation function defined as follows; given x ∈ R3 , the i-th component of sat(x) is equal to xi if |xi | ≤ 1, otherwise it is equal to either 1 or -1 depending on the sign of xi . The previous results still hold because saturation does not modify linearization (21).

IV.

Case Study

We consider the same case study presented in ref.7 in which the spacecraft’s inertia matrix is given by J = diag[27 17 25] kg m2 ; the inclination of the orbit is given by incl = 87◦ , and the orbit’s altitude is 450 km; the corresponding orbital period is about 5600 s. Without loss of generality the value Ω of RAAN is set equal to 0, whereas the initial phase φ0 (see (10)) has been randomly selected equal to φ0 = 0.94 rad. Consider an initial state characterized by attitude equal to to the target attitude q(0) = q¯, and by the following high initial angular rate ω(0) = [0.02 0.02 − 0.03]T rad/s

(42)

due for example to an impact with an object. Moreover, consider that that the time needed for the spacecraft to perform all required measurements is equal to ∆m = 3 s. Disturbance torques have been considered including gravity gradient torque and a residual magnetic moment equal to m0 = [0.5 0.5 0.5]T A m2 In ref.7 the continuous-time feedback (11) has been designed setting kp = 2 1011 , kd = 3 1011 , and  = 10−3 . Here we design feedback (12) by keeping kp = 2 1011 , kd = 3 1011 and by choosing parameters ∆a and  as follows. First, studying numerically the behavior with ∆a of the maximum of the real parts of the eigenvalues of A1 (∆a ) (see (34)), the plot in Fig. 2 is obtained; from the latter plot we determine the value ∆∗a = 1501 s for which it occurs that A1 (∆a ) is Hurwitz for all 0 < ∆a < ∆∗a . Then, for each 0 < ∆a < ∆∗a we determine the value ∆∗m given by (38) obtaining the plot in Fig. 3. Recall that for each ∆a , magnitude ∆∗m represents the maximum value for ∆m for which stability can be achieved. Then, since in our case study ∆m is fixed and equal to 3 s, the set S of values of ∆a for which stability is achievable is determined considering those values for which ∆∗m is greater than 3 s. In order to obtain set S it is useful 9 of 15 American Institute of Aeronautics and Astronautics

0.25

0.2

max(real(eig(A 1 (∆ a ))))

0.15

0.1

0.05

0

-0.05

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

-0.1

-0.15

-0.2 0

200

400

600

800

1000

1200

1400

1600

∆ a (sec)

Figure 2. Maximum of the real parts of the eigenvalues of A1 (∆a ) .

to look at the plots in Fig. 4 which present portions of the plot in Fig. 3. From Fig.’s 3 and 4, we see that stability can be achieved for all 3.49 s < ∆a < ∆∗a = 1501 s. Then, we set ∆a = 7 s to which there corresponds the value 0 = 1.1 10−3 (see (40)). By Proposition 2, it occurs that setting  ≤ 1.1 10−3 and small enough, equilibrium (q, ω) = (¯ q , 0) is locally exponentially stable for (8) (12). Proceeding by trial and error, we fix  = 10−3 . Corresponding simulation results are reported in Fig. 5. Note that when disturbance torques are not included, the spacecraft acquire the desired attitude within 5 orbits. In addition, simulations with disturbance torques show that the designed attitude control system is able to reject the considered disturbances to a significant extent. In Fig. 6 the time behavior of mcoils restricted to the time interval [0 30] s is plotted in order to display clearly the separation between measurement and actuation. The capability of the the proposed state-feedback to handle saturation in the magnetic actuators has been tested, too. This has been done by implementing control law (41) with mcoils max = 10 A m2 . The corresponding time behaviors are plotted in Fig. 7.

V.

Conclusion

In this work a magnetic control law for spacecraft attitude stabilizations has been presented; the control law possesses the important feature of separating measurement time from actuation time. A design method for the parameters of the controller has been obtained and it has been applied to a case study. Further research in the field could consist in designing a different feedback law which still satisfies the separation constraint, but that does not require measures of attitude rate i.e. a so called output feedback.

Acknowledgments The author acknowledges Prof. A. Nascetti for fruitful discussions.

References 1 Sidi,

M. J., Spacecraft dynamics and control, Cambridge University press, 1997. E. and Lovera, M., “Magnetic spacecraft attitude control: A survey and some new results,” Control Engineering Practice, Vol. 13, No. 3, 2005, pp. 357–371. 2 Silani,

10 of 15 American Institute of Aeronautics and Astronautics

1000 ∆ *m ∆ m=3 sec

900 800 700

sec

600 500 400 300

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

200 100 0 0

500

1000

1500

∆ a (sec)

Figure 3. Behavior of ∆∗m with respect to ∆a .

3 Arduini, C. and Baiocco, P., “Active Magnetic Damping Attitude Control for Gravity Gradient Stabilized Spacecraft,” Journal of Guidance, Control, and Dynamics, Vol. 20, No. 1, 1997, pp. 117–122. 4 Wi´ sniewski, R. and Blanke, M., “Fully magnetic attitude control for spacecraft subject to gravity gradient,” Automatica, Vol. 35, No. 7, 1999, pp. 1201–1214. 5 Lovera, M. and Astolfi, A., “Spacecraft attitude control using magnetic actuators,” Automatica, Vol. 40, No. 8, 2004, pp. 1405–1414. 6 Lovera, M. and Astolfi, A., “Global magnetic attitude control of inertially pointing spacecraft,” Journal of Guidance, Control, and Dynamics, Vol. 28, No. 5, 2005, pp. 1065–1067. 7 Celani, F., “Robust three-axis attitude stabilization for inertial pointing spacecraft using magnetorquers,” Acta Astronautica, Vol. 107, 2015, pp. 87–96. 8 Wie, B., Space vehicle dynamics and control, American institute of aeronautics and astronautics, 2008. 9 Wertz, J. R., editor, Spacecraft attitude determination and control, Kluwer Academic, 1978. 10 Rodriguez-Vazquez, A. L., Martin-Prats, M. A., and Bernelli-Zazzera, F., “Full magnetic satellite attitude control using ASRE method,” 1st IAA Conference on Dynamics and Control of Space Systems, 2012. 11 Celani, F., “Spacecraft attitude stabilization with piecewise-constant magnetic dipole moment,” 7th International Conference on Recent Advances in Space Technologies, 2015, pp. 169–174. 12 Mancilla Aguilar, J. L. and Garcia, R. A., “Stability analysis of a certain class of time-varying hybrid dynamical systems,” Latin American Applied Research, Vol. 33, No. 4, 2003, pp. 463–468. 13 Iglesias, P. A., “Input-output stability of sampled-data linear time-varying systems,” IEEE Transactions on Automatic Control, Vol. 40, No. 9, 1995, pp. 1646–1650. 14 Bai, E.-W., Fu, L.-C., and Sastry, S. S., “Averaging analysis for discrete time and sampled data adaptive systems.” IEEE Transactions on Circuits and Systems, Vol. 35, No. 2, 1988, pp. 137–148. 15 Khalil, H. K., Nonlinear systems, Prentice Hall, 2002.

11 of 15 American Institute of Aeronautics and Astronautics

60

4

∆ *m

∆ *m

∆ m=3 sec

∆ m=3 sec

50 3.5

3

40

2.5

sec

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

4.5

30 2

1.5

20

1 10 0.5

0 0

1

2

3

4

5

∆ a (sec)

0 1498

1499

1500

∆ a (sec)

Figure 4. Portions of the plot in Fig. 3.

12 of 15 American Institute of Aeronautics and Astronautics

1501

q1

1 0 -1 0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

q2

1 0 -1

q3

1 0 -1

q4

1 0

ωx (rad/s)

0.04 0.02 0 -0.02

0.03

ωy (rad/s)

0.02 0.01 0 -0.01

0.01

ωz (rad/s)

0 -0.01 -0.02 -0.03

100

(A m 2 )

m coils x

0 -100 -200 -300

600

(A m 2 )

m coils y

400 200 0 -200

(A m 2 )

200

m coils z

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

-1

0 -200 -400

orbit

Figure 5. Simulations without (solid blue lines) and with (red dashed lines) disturbance torques.

13 of 15 American Institute of Aeronautics and Astronautics

(A m 2 )

m coils x

-100 -200 -300 0

5

10

15

20

25

30

0

5

10

15

20

25

30

0

5

10

15

20

25

30

(A m 2 )

m coils y

600 400 200 0

200

(A m 2 )

100

m coils z

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

0

0 -100 -200

time (sec)

Figure 6. Time behavior of mcoils restricted to the time interval [0 30] s. Simulations without (solid blue lines) and with (red dashed lines) disturbance torques.

14 of 15 American Institute of Aeronautics and Astronautics

q1

1 0 -1 0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

q2

1 0 -1

q3

1 0 -1

q4

1 0

0.04

ωx (rad/s)

0.02 0 -0.02 -0.04

0.03

ωy (rad/s)

0.02 0.01 0 -0.01

0.04

ωz (rad/s)

0.02 0 -0.02 -0.04

10

(A m 2 )

m coils x

5 0 -5 -10

10

(A m 2 )

m coils y

5 0 -5 -10

10

(A m 2 )

5

m coils z

Downloaded by Fabio Celani on January 5, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-0086

-1

0 -5 -10

orbit

Figure 7. Controller subject to saturation. Simulations without (solid blue lines) and with (red dashed lines) disturbance torques.

15 of 15 American Institute of Aeronautics and Astronautics