IAA-CU-15-07-12 Optimal Feedback Gains for Spacecraft Attitude Stabilization Using Magnetorquers Fabio Celani∗, Renato Bruni†
Abstract Attitude stabilization of spacecraft using only magnetorqers as torque actuators can be achieved by a proportional-derivative (PD)like control algorithm. The gains of this algorithm are usually determined by a trial-and-error approach within the large search space of the possible values of the gains. Hence, only a small portion of this space can actually be explored. We propose here an innovative systematic approach for finding the gains: they should be computed as those that minimize the settling time of the attitude error. However, the relationship between gains and settling time cannot be explicitly given in analytical form, and this leads to several difficulties in specifying and in solving an explicit optimization model. Therefore, we use an optimization algorithm that guarantees, up to numerical precision, to find the global minimum of a generic function subject to simple bounds. This algorithm is able to work without the need for analytically writing such a function: it only needs to be able to compute it in a number of points. The proposed algorithm is employed in a case study where it converges at the numerical optimum in very reasonable times; thus, it can be used for solving similar stabilization problems.
1
Introduction
Spacecraft attitude control can be obtained by adopting several actuation mechanisms. Among them, electromagnetic actuators are widely used for ∗ School of Aerospace Engineering, Sapienza University of Rome, Italy,
[email protected] † Department of Computer Control and Management Engineering Antonio Ruberti, Sapienza University of Rome, Italy,
[email protected]
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
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. They operate on the basis of the interaction between the magnetic dipole moment generated by those coils and the Earth’s magnetic field: the interaction with the Earth’s 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 a plane orthogonal to the Earth’s magnetic field. Therefore, a different type of actuator often supports magnetorquers to provide full three-axis control. In any case, magnetorquers are frequently 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 [2]; in particular, Lyapunov-based design has been adopted in [3, 4, 5]. The latter works present theoretical investigations showing that, if certain feedback gains are chosen positive and a scaling gain is picked positive and sufficiently small, then attitude stabilization is achieved. However, they do not contain precise guidelines for determining appropriate values for those parameters. In practice, a trial-and-error approach is usually pursued in order to find numerical values for the feedback parameters that guarantee satisfactory performances of the resulting closed-loop system. Since the feedback parameters may range over large intervals, the corresponding search space is very large. As a consequence, the described approach still suffers from two main limitations: (a) it is often very time-consuming; (b) it is not systematic. This means that, even if a satisfactory solution is obtained, it is not know whether better solutions could be obtained by protracting the search, and not even how much the solution could be improved. This could be particularly inconvenient in practical applications. In this work we propose an innovative systematic approach for determining the gains of the proportional-derivative (PD)-like feedback presented in [5]: they should be computed as those that minimize the settling time of the attitude error. However, the relationship between gains and settling time cannot be explicitly given in analytical form, and this leads to several difficulties in specifying and in solving an explicit optimization model. Therefore, we use an optimization algorithm that guarantees, up
IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
to numerical precision, to find the global minimum of a generic function subject to simple bounds. This algorithm is able to work without the need for analytically writing such a function: it only needs to be able to compute it in a number of points. This method is a modification of the standard Lipschitzian approach that eliminates the need to specify a Lipschitz constant. This is possible by carrying out simultaneous searches using all possible constants from zero to infinity. The work is organized as follows: Section 2 defines the spacecraft model and the control algorithm; Section 3 focuses on the determination of the feedback gains by using the mentioned optimization approach; Section 4 reports a case study of determination of the gains using the proposed optimization method.
2
Spacecraft model and control algorithm
In the rest of the paper the following notation is adopted. Symbol I represents the identity matrix. For a ∈ R3 , a× represents the skew symmetric matrix 0 −a3 a2 0 −a1 a× = a3 (1) −a2 a1 0 so that, for b = [b1 b2 b3 ]T ∈ R3 , multiplication a× b is equal to the cross product a × b. 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, Chapter 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 IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
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 [6, Section 5.4]). To obtain a more compact notation, define an expression W (q) as follows 1 q4 I + qv× . (3) W (q) = −qvT 2 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 [6, Section 5.5.3]). The attitude dynamics in body frame can be expressed by J ω˙ = −ω × Jω + T (5) where J ∈ R3×3 is the spacecraft inertia matrix, and T is the control torque expressed in Fb (see [6, Section 6.4]). The spacecraft is equipped with three magnetic coils aligned with the Fb axes which generate the magnetic attitude control torque T = 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 (see [6, Section 12.17]). 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)
which shows explicitly the dependence of B b on both q and t. Grouping together equations (4) (5) (6) the following nonlinear timevarying system is obtained q˙ = W (q)ω J ω˙ = −ω × Jω − B b (q, t)× mcoils
IAA-CU-15-07-12
(8)
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
in which mcoils is the control input. 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 [7, Appendix H]) we obtain B i (t) =
µm [3((m ˆ i (t))T rˆi (t))ˆ ri (t) − m ˆ i (t)] R3
(9)
In equation (9), µm is the total dipole strength, 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 (t) is the vector of the direction cosines of the Earth’s magnetic dipole expressed in Fi which is set equal to sin(θm ) cos(ωe t + α0 ) (10) m ˆ i (t) = sin(θm ) sin(ωe t + α0 ) cos(θm ) where θm is the dipole’s coelevation, ωe = 360.99 deg/day is the Earth’s average rotation rate, and α0 is the right ascension of the dipole at time t = 0. Clearly, Earth’s rotation has been taken into account in equation (10). The following numerical values have been reported in [8] µm = 7.746 1015 Wb m and θm = 170.0◦ . 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 given by xp (t) = R cos(nt + φ0 ) y p (t) = R sin(nt + φ0 )
(11)
where n is the orbital rate, and φ0 an initial phase. The coordinates of the satellite in inertial frame Fi can be easily obtained from (11) 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. As stated before, the control objective is driving the spacecraft so that Fb is aligned with Fi . From equation (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 IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
control strategies for mcoils so that qv → 0 and ω → 0. Reference [5] presents the following stabilizing proportional-derivative (PD)-like control law which was obtained as modification of those in [3] and [4] mcoils = −B b (q, t) × (2 kp qv + kd ω) (12) In order to state the stabilizing properties of the above feedback, it is useful to introduce the following matrix Γi (t) = B i (t)T B i (t)I − B i (t)B i (t)T
(13)
along with its average Γiav
1 = lim T →∞ T
Z
T
Γi (τ )dτ
(14)
0
Then, the following stabilization result has been obtained in [5]. Theorem 1 Consider the magnetically actuated spacecraft described by (8). Apply the proportional-derivative (PD)-like control law (12) with kp > 0 and kd > 0. If the spacecraft’s orbit satisfies Γiav > 0, then there exists ∗ > 0 such that for any 0 < < ∗ , equilibrium (q, ω) = (¯ q , 0) is locally exponentially stable for (8) (12). Remark 2 In [5] it is shown that condition Γiav > 0 is not fulfilled only by low inclination orbits.
3
Determination of the feedback gains by optimization
The previous result does not give specific indications on how to choose the feedback gains kp and kd , and the scaling factor . It only states that kp and kd have to be positive, and has to be positive and smaller than an upper bound ∗ > 0 whose value is very hard to compute in almost all practical cases. As a matter of fact, in practical implementations, appropriate values for kp , kd and are mostly determined by a trial-anderror approach, that basically works as follows. A significant initial state for the spacecraft is considered, some first-guess positive values are chosen for kp , kd and , and simulations are run. Now, in case the spacecraft’s attitude does not converge to the desired one, the value of is lowered till convergence is achieved. If the corresponding time behaviors are not
IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
satisfactory in terms, for example, of overshoot, speed of convergence, and control inputs’ amplitudes, the values of kp and kd are modified according to the following guidelines used in proportional-derivative control: kp is raised in order to increase the speed of convergence at the expense of larger amplitudes of the control inputs; kd is raised in case the response displays an excessive overshoot. After having modified kp and kd , it might be necessary and/or appropriate to assign a different value to . However, the described trial-and-error approach suffers from several limitations. First, this approach is quite time-consuming. Each simulation may require short time, but running a very large number of simulations can often become very demanding. Second, and more important, it is not systematic. This means that, when a satisfactory value of the gains is finally obtained, it is not know whether protracting the search could lead to new values of the gains that provide an overall better performance of the closed-loop system or not. Moreover, in case one chooses to protract the search, it is not known how long that should be protracted, and not even the amount of the possible improvements that this additional work could produce. In any case, unless performing an exhaustive search for all the possible values of the gains, it is alway possible to neglect values of the gains providing an overall better performance. On the other hand, an exhaustive search is almost always impossible to perform, because the search space of the possible values of the gains is too large to do this in reasonable time; the latter aspect will be illustrated for the case study in Section 4. Consequently, we propose in this work the following new approach for determining the feedback gains. First, we introduce the settling-time ts for quaternion q. Since the desired attitude is reached when qv = 0, then we define the settling time tsi for component qi i = 1, 2, 3 as the time needed for |qi | to become and stay smaller than 0.05. Hence, we define the settling-time for quaternion q as ts , max tsi i=1,2,3
in other words, the settling time ts is determined by the slowest component of q. Furthermore, rather then expressing feedback (12) using three parameters kp , kd , and , we rewrite the same feedback in terms of two gains Kp > 0 and Kd > 0 as follows mcoils = −B b (q, t) × (Kp qv + Kd ω)
(15)
Now, the idea is setting the spacecraft’s initial state equal to a particularly IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
significant value, and determining Kp and Kd so that convergence to the desired attitude is achieved as fast as possible. This objective can be formulated by means of the following optimization problem min ts Kp > 0 Kd > 0
(16)
Observe that the solution to problem (16) would probably lead to large values for Kp and Kd since, roughly speaking, the larger Kp and Kd the faster the corresponding closed-loop system. However, using large values of the gains imply that the required time behavior of the magnetic dipole moment is likely to exceeds the maximum value that can be physically generated by magnetorquers. Thus, in order to overcome the latter issue, the optimization problem (16) is solved by considering the following saturated version of feedback (12) 1 ? b × T mcoils = mcoils sat (B (q, t) ) (Kp qv + Kd ω) (17) m?coils in which m?coils is the saturation limit on each magnetic 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 . Note that, as stated in [5, Remark 5], local exponential stability of equilibrium (q, ω) = (¯ q , 0) is guaranteed for the closed-loop system even when saturated feedback (17) is employed. After having defined the optimization problem, we need to solve it numerically. However, even though ts obviously depend on Kp and Kd , it is practically impossible to express this relationship in analytic from. Consequently, many optimization approaches that are able to find the minimum of specific (and easy) classes of functions would not be applicable in this case. Hence, we renounce to give explicit analytic form to such a relationship, and use an optimization routine based on the algorithm described in [9] to solve problem (16). This algorithm is able to work without the need for analytically writing the function under optimization. On the contrary, it needs to be able to compute the function in a number of points, for example by means of a simulation. In addition, the algorithm requires that the feasible set is a bounded rectangle (or a bounded hyperrectangle in the general case). Thus, for the attitude concp trol problem under consideration, it is neccessary that upper bounds K IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
cd for gains Kp and Kd are determined. In the case study in Section and K cp and K cd is illustrated. 4, a possible approach for finding upper bounds K Thus, the optimization problem (16) becomes min ts cp 0 < Kp ≤ K cd 0 < Kd ≤ K
(18)
A description of the optimization routine that will be used to find a solution of (18) is as follows. A major approach in the optimization of a generic function, that may even be not differentiable, is the use of Lipschitzian methods (see, e.g., [10]). In this case, an obstacle may be the need of specifying the so called Lipschitz constant, that is a bound on the rate of change of the function under optimization. However, by opportunely dividing the feasible region into sample rectangles and by selecting them in a suitably order, the search may be conducted without knowing the Lipschitz constant. This corresponds in practice to carrying out simultaneous searches using all possible values of the Lipschitz constant from zero to infinity. In more detail, the algorithm works as follows. The feasible region starts as a single rectangle that is internally normalized to a unit rectangle. The algorithm partitions the rectangle into a collection of smaller rectangles and evaluates the objective function at their central points. The choice of the central point is essential for reducing the number of function evaluations. Potentially optimal rectangles are identified, and only these rectangles are passed to the next iteration of the algorithm for further partitioning and investigation. The algorithm stops when the size of the rectangles becomes too small, or when it reaches the maximum number of iterations. This deterministic algorithm will converge to the global optimum of the function if the sampling is dense enough, even if the search process may require a large amount of function evaluations. Therefore, given enough iterations and up to numerical precision of the computing hardware, we are able to guarantee to find the optimum values of Kp and Kd .
4
Case study
We consider the same case study presented in [5] in which the spacecraft’s inertia matrix is given by J = diag[27, 17, 25] kg m2 ; the saturation level for each magnetic dipole moment is equal to m?coils = 10 A m2 . The IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
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 phases α0 (see (10)) and φ0 (see (11)) have been randomly selected and set equal to α0 = 4.54 rad and φ0 = 0.94 rad. In [5] it has been checked numerically that for the considered orbit, condition Γiav > 0 is fulfilled. 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
(19)
The latter initial state is significant since it can occur when the spacecraft possesses the desired attitude with no angular momentum, but an impact with an object (e.g. a piece of debris or a tiny meteorite) happens and consequently the spacecraft’s angular rate changes istantaneously. In [5], gains Kp and Kd were determined by trial-and-error and set equal to Kp = 2 × 105 and Kd = 3 × 108 . The corresponding time behaviors are reported in Fig. 1, and the settling time turns out to be equal to ts = 4.0179 orbit. Based on the previous values of Kp and Kd , the upper cp = 108 and K cd = 109 . bounds which appear in (18) are set equal to K cp and K cd , it is practically impossible to Note that, with such values of K find a significant solution of (18) by conducting an exhaustive search. In fact, a huge number of couples of real values should be explored (theoretically infinite). Even by limiting the exhaustive search to the integer values, the number of combinations is of the order of 1017 , and even by running 103 simulations per second, the total time would be 1014 sec that is about 3 million years. On the other hand, by applying the optimization method proposed in Section 3, in only 569 secs the following optimal values are obtained for the gains Kp = 8.6274864002 × 104 , Kd = 4.6478557565395 × 107 . The corresponding time evolutions are represented in Fig. 2, and the obtained optimal settling time is equal to ts = 1.0334 orbit. Thus, with respect to the trial-and-error approach, a significant improvement in convergence speed to the target attitude is achieved. Note that, even if the values of the gains determined by trial-and-error are not very far from the optimal ones, the reduction of the settling time is quite large. To achieve a comparable reduction of the settling time by tuning the gains by trial-and-error would have been significantly more time consuming, and would have not assured success.
IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
q1
1 0 -1 0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
q2
1 0 -1
q3
1 0 -1
q4
1 0 -1
ωx
(rad/s)
0.05
0
-0.05 0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
ωy
(rad/s)
0.05
0
-0.05
ωz
(rad/s)
0.05
0
-0.05
(A m2 )
m coils x
10
0
-10 0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
(A m2 )
m coils y
10
0
-10
(A m2 )
m coils z
10
0
-10
orbit
Figure 1: Evolutions with feedback gains determined by trial-and-error.
IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
q1
1 0 -1 0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
q2
1 0 -1
q3
1 0 -1
q4
1 0 -1
ωx
(rad/s)
0.05
0
-0.05 0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
ωy
(rad/s)
0.05
0
-0.05
ωz
(rad/s)
0.05
0
-0.05
(A m2 )
m coils x
10
0
-10 0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
(A m2 )
m coils y
10
0
-10
(A m2 )
m coils z
10
0
-10
orbit
Figure 2: Evolutions with optimal feedback gains.
IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
5
Conclusions
A PD-like feedback control that achieves attitude stabilization of a spacecraft using only magnetorquers is considered in this work, and a systematic method for determining the PD gains is presented. The method consists in setting the gains as those that minimize the settling time. Even though this optimization problem cannot be expressed in analytical form, an optimization routine providing the numerical solution was developed. It is shown through a case-study that the feedback gains determined by the proposed method lead to a closed-loop system that is substantially faster than by using gains determined through a trial-and-error approach. A saturation on the magnetic dipole moment is included in the spacecraft model in order to take into account of the limits on the currents flowing in the magnetorquers. We remark that, in the proposed procedure, the gains are optimized with respect to a fixed initial state set equal to a significant value. As a result, optimality of the gains cannot be guaranteed for different initial states, even though the reported gains should maintain near-optimality for initial states that are similar to the considered ones. It will be subject of future research to propose a more refined method for finding the gains, so that some type of optimality is ensured for all initial states belonging to a significant fixed set.
References [1] M. J. Sidi. Spacecraft dynamics and control. Cambridge University Press, 1997. [2] E. Silani and M. Lovera. Magnetic spacecraft attitude control: A survey and some new results. Control Engineering Practice, 13(3):357– 371, 2005. [3] M. Lovera and A. Astolfi. Spacecraft attitude control using magnetic actuators. Automatica, 40(8):1405–1414, 2004. [4] M. Lovera and A. Astolfi. Global magnetic attitude control of inertially pointing spacecraft. Journal of Guidance, Control, and Dynamics, 28(5):1065–1067, 2005. [5] F. Celani. Robust three-axis attitude stabilization for inertial pointing spacecraft using magnetorquers. Acta Astronautica, 107:87–96, 2015. IAA-CU-15-07-12
3rd IAA Conference On University Satellite Missions And Cubesat Workshop
[6] B. Wie. Space vehicle dynamics and control. American institute of aeronautics and astronautics, 2008. [7] J. R. Wertz, editor. Spacecraft attitude determination and control. Kluwer Academic, 1978. [8] A. L. Rodriguez-Vazquez, M. A. Martin-Prats, and F. BernelliZazzera. Full magnetic satellite attitude control using ASRE method. In 1st IAA Conference on Dynamics and Control of Space Systems, 2012. [9] B. E. Stuckman. D. R. Jones, C. D. Perttunen. Lipschitzian optimization without the Lipschitz constant. Journal of Optimization Theory and Application, 79(1):157–181, 1993. [10] J. D. Pinter. Global optimization in action - continuous and Lipschitz optimization: Algorithms, implementations and applications. Kluwer Academic Publishers, Dordrecht, 1996.
IAA-CU-15-07-12
