- Email: [email protected]
H∞ attitude control of magnetically actuated satellites*
Proceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011
H∞ attitude control of magnetically actuated satellites Andrea Maria Zanchettin ∗ Marco Lovera ∗ ∗
Dipartimento di Elettronica e Informazione, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133, Milano, Italy (e-mail: {zanchettin, lovera}@elet.polimi.it)
Abstract: The problem of attitude regulation for magnetically actuated spacecraft is considered and a randomized approach for the solution of the robust periodic state-feedback control problem is proposed. The algorithm is applied to the tuning of "projection-based" attitude controllers. The performance of the designed controllers is illustrated in a detailed simulation study. Keywords: Guidance, navigation and control of vehicles; Autonomous systems. 1. INTRODUCTION Magnetic torquers are an effective and reliable technology for the attitude control of small satellites in low Earth orbit. The principle of operation of such actuators is the interaction between the magnetic dipole generated by the torquers and the magnetic field of the Earth. The main difficulty in the design of attitude control laws based on magnetic torquers is that the torques they generate are instantaneously constrained to lie in the plane orthogonal to the local direction of the geomagnetic field vector. This implies that the attitude control problem is formulated over a time-varying model. Furthermore, the open-loop attitude dynamics is usually marginally stable (i.e., undamped) or unstable and the time-variability of the dynamics introduced by the magnetic actuators is characterised by a very long period. In recent years this control problem has been studied extensively and two main lines of work can be identified in the literature (at least as far as linear regulation near a nominal attitude is concerned). Methods based on averaged models have been proposed in, e.g., Stickler and Alfriend [1976], Hablani [1995]: the time-varying model of the magnetically actuated spacecraft is replaced with an approximate time-invariant one obtained using averaging techniques, thus leading to a much simplified design problem. The main drawbacks of this approach are that the designer has to verify a posteriori the achieved stability degree and performance level and that averaging implies limitations in closed-loop performance. Methods based on periodic models, on the other hand, exploit the quasi-periodic variability of the geomagnetic field. Most of the recent work on this problem has focused on the use of optimal and robust periodic control theory for the design of state and output feedback regulators (see, e.g., Psiaki [2001], Silani and Lovera [2005], Yan et al. [2005], Wood et al. [2006]). While periodic control design methods have the advantage of guaranteeing closed-loop stability a priori, they lead to the synthesis of time-periodic regulators, which are difficult to implement and operate in practice for a number of reasons (see, e.g., Silani and Lovera [2005] for a detailed discussion). In view of this, we will focus on the problem of tuning so-called "projection based" controllers for magnetic attitude regulation (first proposed in Martel et al. [1988b]), a ⋆ This work has been supported by the MIUR project "‘New Methods for Identification and Adaptive Control for Industrial Systems"’.
978-3-902661-93-7/11/$20.00 © 2011 IFAC
controller structure widely used in the practice of magnetic attitude control. Previous work on this problem has relied on optimisation-based tuning of fixed structure LQ-like controllers (see Pulecchi et al. [2010], Vigano et al. [2010]. Optimizationbased design is a very active research area aiming at the development of computationally efficient methods for the solution of challenging control issues (e.g. static output feedback) in terms of constrained optimization problems. Among others, randomized methods, see, e.g., Tempo et al. [2005], Zhigljavsky and Zilinskas [2008] and non-smooth optimization techniques, see Boyd and Barratt [1991] and Lewis [2007], appear today as the most promising paths to the development of viable design methods applicable to practical engineering problems. In the light of the above discussion, the aim of this paper is to present the results obtained in the design of control laws parameterised via constant gains with guaranteed nominal stability and H∞ performance, using an optimisation-based method. 2. BACKGROUND ON LINEAR PERIODIC SYSTEMS In this Section some basic definitions associated with LTP systems are provided (see, e.g., Bittanti and Colaneri [2009] for an in-depth treatment of this topic) and a framework extending frequency domain analysis to continuous-time LTP systems is reviewed. Consider the LTP system x˙ = A (t) x (1) where matrix A (·) is T -periodic, namely A (t + T ) = A (t). Then, given the initial condition x (t = τ ) = x (τ ) at time τ the system (1) has the solution (2) x (t) = φA (t, τ ) x (τ ) where φA (·, ·) is the so-called transition matrix (compare with φA (·, ·) = eA(t−τ ) for LTI systems). The transition matrix can be computed solving the ODE system φ˙A (t, τ ) = A (t) φA (t, τ ) , (3) φA (τ , τ ) = I. Consider ψA (τ ) = φA (τ + T, τ ), the so called monodromy matrix, then the following similarity holds:
ψA (τ1 ) = φA (τ1 , τ2 ) ψA (τ2 ) φA (τ1 , τ2 )−1 (4) from which it follows that the eigenvalues of ψA (·) are timeinvariant quantities. Such eigenvalues are called characteristic
8479
10.3182/20110828-6-IT-1002.02751
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
• Earth Centered Inertial reference axes (ECI). The Earth’s center is the origin of these axes. The positive X-axis points in the vernal equinox direction. The Z-axis points in the direction of the North Pole. The Y-axis completes the right-handed orthogonal triad. • Orbital Axes (X0 , Y0 , Z0 ). The origin of these axes is at the satellite center of mass. The X-axis points to the Earth’s center; the Y-axis points in the direction of the orbital velocity vector. The Z-axis is normal to the satellite orbit plane. • Satellite Body Axes. The origin of these axes is at the satellite center of mass; in nominal Earth-pointing conditions the Xb (yaw), Yb (roll) and Zb (pitch) axes are aligned with the corresponding orbital axes.
multipliers. Again, unlike the LTI case, the time-varying eigenvalues of matrix A (·) have nothing to do with the stability of system (1). However, exponential stability can be checked using the characteristic multipliers. The following Theorem holds. Theorem 2.1. System (1) is uniformly exponentially stable if and only if its characteristic multipliers - i.e., the eigenvalues of ψA (τ ) - lie in the open unit disc. As for stable LTI systems, a frequency response operator can be defined for stable LTP systems. The frequency domain representation for LTI systems is related with the fact that sinusoidal inputs are mapped into sinusoidal outputs at the same frequency with possibly different amplitude and phase. This notion can be extended to LTP systems by making reference to a more general class of signals: the set of exponentiallymodulated periodic signals, i.e., h (t) = ∑ hk esk t (5) k∈Z
where sk = s + jkΩ, Ω = 2π /T , which can be arranged as follows T H = . . . hT−2 hT−1 hT0 hT1 hT2 . . . . (6) Consider the following T -periodic stable system x(t) ˙ = A (t) x(t) + B (t) u(t) (7) y(t) = C (t) x(t) + D (t) u(t). The key idea, see Wereley and Hall [1990], further developed and formalized in Zhou and Hagiwara [2002a], is to expand the complex Fourier series of the dynamics matrices A (t) =
∑ Am e jmΩt
(8)
The angular kinematics and dynamics of the spacecraft are T modeled using as state variables the quaternion q = q1 qTR = T [q1 q2 q3 q4 ] describing the attitude of the body axes with respect to the orbital axes, and the inertial angular velocity vector T ω = [ωx ωy ωz ] , with respect to the body axes. With respect to the selected variables, the nominal, Earth-pointing, equilibrium T corresponds to the attitude quaternion q¯ = [1 0 0 0] and to the T angular rate ω¯ = [0 0 −Ω0 ] , where Ω0 is the orbital angular rate. Define the state vector x = [δ qTR δ ω T ]T formed with small displacements of the vector part qR of the attitude quaternion q T from the nominal values q¯R = [0 0 0] and small deviations of the body rates from the nominal values ω¯ x = ω¯ y = 0, ω¯ z = −Ω0 . Then the attitude dynamics can be linearized and a local linear model for the attitude can be defined as in Pulecchi et al. [2010]
m∈Z
and similarly for B (t) ,C (t) and D (t) and perform a harmonic balance approach. It was shown that defining .. . . .. .. . . . . . . . A0 A−1 A−2 . . . A = . . . A1 A0 A−1 . . . (9) . . . A A A . . . 2 1 0 .. .. .. . . . . . .
and similarly B, C and D, the harmonic transfer function G (s) that describes the input-output relation is 1
G (s) = C [sI − (A − N )]−1 B + D (10) where N = blkdiag { jnΩI} , n ∈ Z. Note that the operator defined in (10) is a doubly-infinite dimensional one. On the other hand, it is no longer time-varying and this fact can be exploited to obtain input-output information related with system (7) without numerical integration.
x(t) ˙ = Ax(t) + BT [Tmag (t) + Td (t)]
where Tmag is the magnetic control torque vector and Td is any disturbance torque vector. Taking into account that Tmag can be written as Tmag (t) = m(t) × b(t) = S(b(t))m(t) = −b(t) × m(t) T
equation (11) can be equivalently written as
x(t) ˙ = Ax(t) + Bcm (t)m(t) + BT Td (t),
In the following denoted with G (s) = A(t),B(t),C(t),D(t) .
(14)
where
3. MODELING AND CONTROL OF MAGNETICALLY ACTUATED SATELLITES
1
(12)
where b = [bx by bz ] is the geomagnetic field vector (in body frame), m is the dipole vector of the magnetic torquers and " # 0 bz −by S(b) = −bz 0 bx , (13) by −bx 0
For practical purposes, however, truncation must be, of course, considered. Given a stable system, performance measures such as the H∞ or the H2 norm can be thus computed using a suitably truncated frequency response operator, see, e.g., Zhou and Hagiwara [2002a] or Zhou and Hagiwara [2002b] for details.
In order to represent the attitude motion of an Earth-pointing spacecraft on a circular orbit the following reference systems are adopted:
(11)
0 Ω0 0 A= 0 0 0 0 0 0 BT = I −1 xx 0 0
−Ω0 0 0.5 0 0 0 0 0 0.5 0 0 0 0.5 0 0 0 0 0 Wx 0 Wy 0 0 −6ky Ω20 0 0 +6kz Ω20 0 0 0 0 0 0 0 0 0 , Bcm (t) = BT S(b(t)) 0 0 −1 Iyy 0 −1 0 Izz
(15)
(16)
Ixx −Iyy Izz −Ixx Iyy , kz = Izz , Wx = −kx Ω0 + kwx Ω, J Wy = −ky Ω0 −kwy Ω, kwx = Ixx , kwy = IJyy . Here, Ω is the nominal
and kx =
8480
Iyy −Izz Ixx , ky
=
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
wheel speed 2 . A very common structure for magnetic attitude control laws, which goes back to classical papers such as Martel et al. [1988a], consists of control laws of the kind 1 m(t) = ST (b(t))Tid (t) (17) kb (t) k2 where b(t) is the measurement, at time t, of the geomagnetic field b and Tid (t) is an "ideal" control torque to be determined on the basis of a suitable (usually static) feedback of state or output variables, according to the specific attitude control architecture of the considered spacecraft. The above control law can be readily given a simple geometric interpretation. Recalling (12) we can express the torque generated by the magnetic coils as 1 S(b(t))ST (b(t))Tid (t) (18) Tmag (t) = kb(t)k2 which can be easily interpreted as the projection of vector Tid onto the plane orthogonal to the direction of the magnetic field vector b (hence the name of "projection-based" controllers). 4. H∞ CONTROLLER DESIGN 4.1 Formulation of the design problem Given the satellite attitude dynamics discussed so far and assuming the value of the geomagnetic field b as well as the state vector x are measurable, the attitude dynamics can be written as a LTP system as follows: x˙ (t) = Ax (t) + BT S (b (t)) m (t) + BT Td (t) 1 = Ax (t) + BT S (b (t)) ST (b (t)) u (t) + BT Td (t) kb (t)k2 | {z } B(t)
= Ax (t) + B (t) u (t) + BT Td (t) .
(19) Under the above assumptions, we can design a static state feedback controller Tid = u = Kx. The aim of the controller is to stabilize the system and possibly optimize some given performance criterion. Since the model is obtained through linearization and many effects, such as a magnetic residual dipole, see Corno and Lovera [2009], are not taken into account, the controller should be able to attenuate the effect of any torque disturbance acting on the attitude variables. To this end, we can define the following variable to be used as a performance measure: Z (20) x (t) = Cσ x (t) , σ > 0 z (t) = σK
which implies ρ ψA+B(t)K (0) < 1, where by ρ (X) we denote the spectral radius of the matrix X. In other terms, the tuning of the control law u = Kx can be formalized as follows: σ minkGzT k subject to ρ ψA+B(t)K (0) < 1. (22) d ∞ K
It should be noticed that, since the cost function in (21) is in general non-smooth, gradient-based descent algorithms such as the ones adopted in Pulecchi et al. [2010], Vigano et al. [2010] might fail in converging to its minimum. In order to circumvent this difficulty, a dedicated approach to the solution of this problem is proposed, which is outlined in the following subsection. 4.2 Controller tuning algorithm
A significant amount of work has been carried out to cope with non-smooth optimization, yielding either deterministic or stochastic methods. In this work, we adopted a stochastic method to solve this optimization problem. Many randomized method exists, however, but the key issue is that a control law is optimal if in its neighborhood a better control law cannot be found (or, equivalently, can be found with null probability). From this simple observation, it is possible to design a procedure to optimize a certain criterion. Assume an initial stabilising controller K (i) is available, e.g., tuned using common techniques for LTI systems applied to the linearized time-averaged model, see e.g. Stickler and Alfriend [1976], or as the output of the previous iteration step. To compute a new controller K (i+1) , the idea is to pick-up randomly sampled controllers in the neighborhood of K (i) and select among them the best one in terms of the cost function in (22). When it is no longer possible to compute better controllers, the algorithm ends with the optimal one. The actual algorithm is further detailed in the following. For the sake of brevity, we
define the cost function to be min σ imized as J (K) = GzTd . As for the optimization algorithm, ∞ the controller K is regarded as a column vector.
Tuning algorithm Set the parameter 0 ≪ rmax < 1 to check the convergence, the number N ≫ 1 of M ONTE C ARLO samples to be generated and the initial step size µ .
where Z = [I3 0] and σ has been considered in order to limit the control effort 3 . A natural criterion to achieve an acceptable disturbance rejection is the H∞ norm of the input-output operator from the disturbance Td to the performance measure z, namely:
σ
GzT ≡ sup kzk2 , (21) d ∞ Td ∈L2 kTd k2
σ = [A + B (t) K, B ,C , 0]. where GzT T σ d The controller gain K can be tuned in order to minimize the cost function in (21) subject, of course, to the stability constraint, 2
The assumption of a diagonal inertia matrix is made only for ease of presentation. 3 Notice that zT z = xT Z T Z + σ 2 K T K x = xT Z T Zx + σ 2 uT u.
8481
(1) given a stabilizing controller K (i) generate N controllers using the following M ONTE C ARLO generation, namely for j = 1, . . . , N !
∆K (i)
(i) j j (i)
K ← K + µ K η + (23)
∆K (i) where η ∼ N (018×1, I18 ) and ∆K (i) = K (i) − K (i−1) ; (2) compute the corresponding norms, namely c j ← J K j ; ˆ (3) select the candidate controller K j (24) jˆ ← argmin J K j j
(4) compute an approximate of the rejection ratio r(i+1) o n ˆ (25) r(i+1) ← prob J K j > J K (i) ˆ (5) if ρ ψA+B(t)Kˆ (0) < 1 then set K (i+1) ← K j and i ← i+ 1 else decrease the step µ and goto (1); (6) if r(i+1) > r(i) decrease the step µ , else if r(i+1) < r(i) increase it;
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
• the H∞ norm of the frequency response operator from the disturbance torque Td to the orientation error δ qR = Zx
kδ qR k2 c1 = Gδ qR Td ∞ ≡ sup (28) Td ∈L2 kTd k2
(7) if r(i+1) ≥ rmax then stop and return Kσ ← K (i+1) , else goto (1); Notice that ∆K (i) /∆K (i) can be regarded as an approximation of the steepest descent direction, see Figure 1. Thus, the candidate controllers K j ’s generated within the step (1) are randomly selected by possibly taking into account the knowledge of the steepest descend direction estimated in the previous step. For
• the H2 norm of the frequency response operator from the disturbance torque Td to the magnetic dipole vector m
kmk∞ c2 = GmTd 2 ≡ sup (29) kT d k∞ Td ∈L∞ • the spectral radius of the characteristic multipliers of A + B (t) Kσ ρ = max λi ψA+B(t)Kσ (0) (30) i
which are reported in Tab. 1 as functions of the selected values of σ . The design model in (19) has been then simulated in Input-output norms c1 (σ ) c2 (σ )
σ 60 120 240
5. SIMULATION STUDY The considered spacecraft is of the type described in Section 3 and operates in a near polar orbit (87 deg inclination) with an altitude of 450 km and a corresponding orbital period of T = 5614.8 s = 2π /Ω0. For the given orbit, the following nominal (periodic) model for the components of the geomagnetic field has been considered: # " 7 cos(Ω0t) + 48 sin(Ω0t) −6 23 cos(Ω0t) − 2 sin(Ω0t) T. (26) b(t) = 10 5 The numerical values for the parameters of the linearized model are as follows: • satellite inertia tensor (kg m2 ): I = diag (Ixx , Iyy , Izz ) = diag (35, 17, 25)
ρ (σ )
12072.45 11039.75 9348.98
0.00283 0.00345 0.00598
Table 1. Computed norms for the closed-loop system.
Fig. 1. One iteration of the tuning algorithm for two decision variables K1 and K2 . the sake of completeness, note that the computed controller is actually a random matrix. To check whether the algorithm converges to the best controller, the procedure has been implemented and run several times, for each given value of σ , ending with a standard deviation in the cost function of about ±0.7%. In other words, the randomized operation of the optimization algorithm is such that the minimality of the cost function is guaranteed within a very small tolerance.
239.19 283.96 399.32
closed-loop using the control law u = Kσ x. The following initial perturbation has been applied to the angular rate vector: T δ ω (0) = 1 × 10−3 2 × 10−3 0 rad/s (31)
while a residual dipole of 1 Am2 along each axis has been considered as an unmeasurable disturbance, i.e., the following disturbance in (19) has been taken into account T Td (t) = S (b (t)) [1 1 1] .
(32) q Time histories of the norm of the orientation error δ qTR δ qR = √ √ xT Z T Zx as well as norm of the control dipole mT m are shown during the first orbit, i.e. t ∈ [0, T ], and during the following four orbits, t ∈ [T, 5T ], in Fig. 2, and Fig. 3, respectively. Figure 2 shows that for decreasing σ , the transient related with a non-zero initial condition is faster and more damped. On the other hand, increasing σ the required control effort increases, as well. As for the effects of disturbances, it can be noticed 0.03
0.02
0.01
0 0
(27)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.4 0.5 0.6 0.7 Normalized time (orbits)
0.8
0.9
1
Normalized time (orbits)
• momentum wheel inertia (kg m2 ): J = 0.01; ¯ = −200. • nominal wheel speed (rad/s): Ω
10
We recall that σ , which is actually a parameter of the optimization problem, has been introduced to weight the control effort in the cost function (21) and thus to limit the actuator energy. As a matter of fact, however, the selection of σ is not straightforward. Therefore, the approach discussed so far has been used (with N = 540, rmax = 0.995 and µ = 10−3 ) to generate three controllers Kσ , corresponding to σ = 60, 120 and 240. In order to compare the obtained controllers, the following quantities have been computed at the termination of the optimization algorithm:
5
0 0
0.1
0.2
0.3
Fig. 2. Norm of the orientation error δ qR and of the magnetic dipole m for σ = 60 (blue), σ = 120 (green), σ = 240 (red).
8482
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
quaternion) are strongly coupled and their transients correspond to the response of the stabilised precession and nutation modes of the open loop dynamics, while the pitch axis (q4 component of the attitude quaternion) is stabilised much more efficiently thanks to the above mentioned favourable controllability properties. 1
1 q
that, for increasing σ , the rejection of the disturbance torque is weaker, see Figure 3, while a better disturbance rejection does not require an additional effort in the control signal. This is somehow confirmed by the data in Table 1. In fact, as one can see, as the weight on the control effort σ increases, the H∞ norm of the frequency response from the disturbance torque to the attitude error, i.e. c1 , rapidly increases, as well. On the other hand, the sensitivity of the H2 norm of the frequency response from the disturbance torque to the magnetic dipole, namely c2 , is less remarkable, justifying the behaviour in Figure 3.
0.9998 0.9996
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
0.02 q2
The simulations presented in the following have been carried out using an object-oriented environment for satellite dynamics (see Lovera [2006] for details) developed using the Modelica language (Fritzson and Bunus [2002]). More precisely, a full nonlinear simulation of the coupled rigid body orbital and attitude dynamics has been performed and the following models of the space environment have been implemented: the JGM3 spherical expansion for the geopotential as a gravitational model, the Harris-Priester model for the atmospheric density distribution (see Montenbruck and Gill [2000] for details) and the International Geomagnetic Reference Field (IGRF, see Wertz [1978]) for the Earth’s magnetic field (up to order 10). Disturbance torques due to gravity gradient (including J2 effects), magnetic residual dipole (assuming a residual dipole of 1 A m2 along each spacecraft body axis) and solar radiation pressure (computed using the solar coordinates formulas given in Montenbruck and Gill [2000]) have been taken into account in the simulation. Finally, all the controllers have been implemented in digital form, using a conventional sample and hold scheme, using a sampling period Ts = 1 s.
0 −0.02
q3
0.05 0 −0.05
0 −3
4
5
x 10
q
0 −5
0 −3
x 10
2 ω1 [rad/s]
1 0 −1 −2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
−3
x 10
4 ω2 [rad/s]
2 0 −2 −4
0 −5
x 10
5 ω [rad/s]
The simulation results obtained for the controllers corresponding to σ = 60 and σ = 240 are presented, respectively, in Figures 4 and 5. As can be seen from the Figures, both controllers damp out the effect of the initial angular rate perturbation (31) in less than one orbit (an adequate performance for a magnetic attitude control scheme), with the controller for σ = 60 providing a settling time of about a quarter of the orbit period. The system is then brought to its steady state response under the effect of the cyclic external disturbance torques. As expected, the yaw and roll axes (q2 and q3 components of the attitude
3
0 −5 −10
0
2
mx [A m ]
2 0 −2 −4
2
my [A m ]
2
−3
x 10 8
1 0 −1 −2
6 2
m [A m ]
10 5 0
z
4 2
−5 −10
0 1
1.5
2
2.5
3
3.5
4
4.5
5
Fig. 4. Quaternion, angular rates and control dipole moments for σ = 60.
Normalized time (orbits)
1.5
1
Finally, a sensitivity analysis of the closed-loop performance to the choice of the sampling period in the implementation of the control laws has been carried out. The detailed results have been omitted for brevity, but the sampling period can be increased up to Ts = 20 s without any significant performance degradation.
0.5
0 1
1.5
2
2.5
3
3.5
4
4.5
5
Normalized time (orbits)
6. CONCLUSIONS Fig. 3. Norm of the orientation error δ qR and of the magnetic dipole m for σ = 60 (blue), σ = 120 (green), σ = 240 (red).
The problem of designing attitude controllers for magnetically actuated spacecraft has been considered. An approach to the
8483
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
q
1
1 0.9998 0.9996
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
0.02 q2
0 −0.02 −0.04 0.04 q3
0.02 0 −0.02
q4
0.01 0 −0.01
0 −3
ω1 [rad/s]
2
x 10
1 0 −1 −2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
−3
ω2 [rad/s]
4
x 10
2 0 −2 −4
0 −5
ω3 [rad/s]
5
x 10
0 −5 −10
0
m [A m2]
1
x
0 −1 −2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
0.5 m [A m2]
0
y
−0.5 −1 −1.5
z
m [A m2]
2 0 −2 −4
Fig. 5. Quaternion, angular rates and control dipole moments for σ = 240. tuning of projection based controllers has been proposed, relying on robust periodic state feedback control techniques. The performance of the proposed control algorithm has been illustrated in a simulation study. REFERENCES Bittanti, S. and Colaneri, P. (2009). Periodic systems - filtering and control. Springer. Boyd, S. and Barratt, C. (1991). Linear controller design: limits of performance. Prentice Hall. Corno, M. and Lovera, M. (2009). Spacecraft attitude dynamics and control in the presence of large magnetic residuals. Control Engineering Practice, 17(4), 456–468. Fritzson, P. and Bunus (2002). Modelica - a general objectoriented language for continuous and discrete-event system modelling and simulation. In Proceedings of the 35th IEEE Annual Simulation Symposium.
Hablani, H. (1995). Comparative stability analysis and performance of magnetic controllers for bias momentum satellites. Journal of Guidance, Control and Dynamics, 18(6), 1313– 1320. Lewis, A. (2007). Nonsmooth optimization and robust control. Annual Reviews in Control, 31, 167–177. Lovera, M. (2006). Control-oriented modelling and simulation of spacecraft attitude and orbit dynamics. Journal of Mathematical and Computer Modelling of Dynamical Systems, Special issue on Modular Physical Modelling, 12(1), 73–88. Martel, F., Pal, P.K., and Psiaki, M. (1988a). Active magnetic control system for gravity gradient stabilised spacecraft. In 2nd Annual AIAA/USU Conference on Small Satellites, Logan (Utah), USA. Martel, F., Pal, P., and Psiaki, M. (1988b). Active magnetic control system for gravity gradient stabilised spacecraft. In 2nd Annual AIAA/USU Conference on Small Satellites, Logan (Utah), USA. Montenbruck, O. and Gill, E. (2000). Satellite orbits: models, methods, applications. Springer. Psiaki, M. (2001). Magnetic torquer attitude control via asymptotic periodic linear quadratic regulation. Journal of Guidance, Control and Dynamics, 24(2), 386–394. Pulecchi, T., Lovera, M., and Varga, A. (2010). Optimal discrete-time design of three-axis magnetic attitude control laws. IEEE Transactions on Control Systems Technology, 18(3), 714–722. Silani, E. and Lovera, M. (2005). Magnetic spacecraft attitude control: A survey and some new results. Control Engineering Practice, 13(3), 357–371. Stickler, A.C. and Alfriend, K.T. (1976). An elementary magnetic attitude control system. Journal of Spacecraft and Rockets, 13(5), 282–287. Tempo, R., Calafiore, G., and Dabbene, F. (2005). Randomized algorithms for analysis and control of uncertain systems. Springer. Vigano, L., Bergamasco, M., Lovera, M., and Varga, A. (2010). Optimal periodic output feedback control: a continuous-time approach and a case study. International Journal of Control, 83(5), 897–914. Wereley, N.M. and Hall, S.R. (1990). Frequency response of linear time periodic systems. In IEEE International Conference on Decision and Control, CDC. Wertz, J. (1978). Spacecraft attitude determination and control. D. Reidel Publishing Company. Wood, M., Chen, W.H., and Fertin, D. (2006). Model predictive control of low earth orbiting spacecraft with magnetotorquers. In IEEE International Conference on Control Applications, Munich, Germany. Yan, H., Ross, M., and Alfriend, K. (2005). Three-axis magnetic attitude control using pseudospectral control law. In Proc. of the 2005 AAS/AIAA Astrodynamics Specialists Conference, Lake Tahoe, USA. Paper AAS 05-417. Zhigljavsky, A. and Zilinskas, A. (2008). Stochastic global optimization. Springer. Zhou, J. and Hagiwara, T. (2002a). Existence conditions and properties of the frequency response operators of continuoustime periodic systems. SIAM Journal on Control and Optimization, 40(6), 1867–1887. Zhou, J. and Hagiwara, T. (2002b). H2 and H∞ norm computations of linear continuous-time periodic systems via the skew analysis of frequency response operators. Automatica, 38(8), 1381–1387.
8484
H∞ attitude control of magnetically actuated satellites Andrea Maria Zanchettin ∗ Marco Lovera ∗ ∗
Dipartimento di Elettronica e Informazione, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133, Milano, Italy (e-mail: {zanchettin, lovera}@elet.polimi.it)
Abstract: The problem of attitude regulation for magnetically actuated spacecraft is considered and a randomized approach for the solution of the robust periodic state-feedback control problem is proposed. The algorithm is applied to the tuning of "projection-based" attitude controllers. The performance of the designed controllers is illustrated in a detailed simulation study. Keywords: Guidance, navigation and control of vehicles; Autonomous systems. 1. INTRODUCTION Magnetic torquers are an effective and reliable technology for the attitude control of small satellites in low Earth orbit. The principle of operation of such actuators is the interaction between the magnetic dipole generated by the torquers and the magnetic field of the Earth. The main difficulty in the design of attitude control laws based on magnetic torquers is that the torques they generate are instantaneously constrained to lie in the plane orthogonal to the local direction of the geomagnetic field vector. This implies that the attitude control problem is formulated over a time-varying model. Furthermore, the open-loop attitude dynamics is usually marginally stable (i.e., undamped) or unstable and the time-variability of the dynamics introduced by the magnetic actuators is characterised by a very long period. In recent years this control problem has been studied extensively and two main lines of work can be identified in the literature (at least as far as linear regulation near a nominal attitude is concerned). Methods based on averaged models have been proposed in, e.g., Stickler and Alfriend [1976], Hablani [1995]: the time-varying model of the magnetically actuated spacecraft is replaced with an approximate time-invariant one obtained using averaging techniques, thus leading to a much simplified design problem. The main drawbacks of this approach are that the designer has to verify a posteriori the achieved stability degree and performance level and that averaging implies limitations in closed-loop performance. Methods based on periodic models, on the other hand, exploit the quasi-periodic variability of the geomagnetic field. Most of the recent work on this problem has focused on the use of optimal and robust periodic control theory for the design of state and output feedback regulators (see, e.g., Psiaki [2001], Silani and Lovera [2005], Yan et al. [2005], Wood et al. [2006]). While periodic control design methods have the advantage of guaranteeing closed-loop stability a priori, they lead to the synthesis of time-periodic regulators, which are difficult to implement and operate in practice for a number of reasons (see, e.g., Silani and Lovera [2005] for a detailed discussion). In view of this, we will focus on the problem of tuning so-called "projection based" controllers for magnetic attitude regulation (first proposed in Martel et al. [1988b]), a ⋆ This work has been supported by the MIUR project "‘New Methods for Identification and Adaptive Control for Industrial Systems"’.
978-3-902661-93-7/11/$20.00 © 2011 IFAC
controller structure widely used in the practice of magnetic attitude control. Previous work on this problem has relied on optimisation-based tuning of fixed structure LQ-like controllers (see Pulecchi et al. [2010], Vigano et al. [2010]. Optimizationbased design is a very active research area aiming at the development of computationally efficient methods for the solution of challenging control issues (e.g. static output feedback) in terms of constrained optimization problems. Among others, randomized methods, see, e.g., Tempo et al. [2005], Zhigljavsky and Zilinskas [2008] and non-smooth optimization techniques, see Boyd and Barratt [1991] and Lewis [2007], appear today as the most promising paths to the development of viable design methods applicable to practical engineering problems. In the light of the above discussion, the aim of this paper is to present the results obtained in the design of control laws parameterised via constant gains with guaranteed nominal stability and H∞ performance, using an optimisation-based method. 2. BACKGROUND ON LINEAR PERIODIC SYSTEMS In this Section some basic definitions associated with LTP systems are provided (see, e.g., Bittanti and Colaneri [2009] for an in-depth treatment of this topic) and a framework extending frequency domain analysis to continuous-time LTP systems is reviewed. Consider the LTP system x˙ = A (t) x (1) where matrix A (·) is T -periodic, namely A (t + T ) = A (t). Then, given the initial condition x (t = τ ) = x (τ ) at time τ the system (1) has the solution (2) x (t) = φA (t, τ ) x (τ ) where φA (·, ·) is the so-called transition matrix (compare with φA (·, ·) = eA(t−τ ) for LTI systems). The transition matrix can be computed solving the ODE system φ˙A (t, τ ) = A (t) φA (t, τ ) , (3) φA (τ , τ ) = I. Consider ψA (τ ) = φA (τ + T, τ ), the so called monodromy matrix, then the following similarity holds:
ψA (τ1 ) = φA (τ1 , τ2 ) ψA (τ2 ) φA (τ1 , τ2 )−1 (4) from which it follows that the eigenvalues of ψA (·) are timeinvariant quantities. Such eigenvalues are called characteristic
8479
10.3182/20110828-6-IT-1002.02751
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
• Earth Centered Inertial reference axes (ECI). The Earth’s center is the origin of these axes. The positive X-axis points in the vernal equinox direction. The Z-axis points in the direction of the North Pole. The Y-axis completes the right-handed orthogonal triad. • Orbital Axes (X0 , Y0 , Z0 ). The origin of these axes is at the satellite center of mass. The X-axis points to the Earth’s center; the Y-axis points in the direction of the orbital velocity vector. The Z-axis is normal to the satellite orbit plane. • Satellite Body Axes. The origin of these axes is at the satellite center of mass; in nominal Earth-pointing conditions the Xb (yaw), Yb (roll) and Zb (pitch) axes are aligned with the corresponding orbital axes.
multipliers. Again, unlike the LTI case, the time-varying eigenvalues of matrix A (·) have nothing to do with the stability of system (1). However, exponential stability can be checked using the characteristic multipliers. The following Theorem holds. Theorem 2.1. System (1) is uniformly exponentially stable if and only if its characteristic multipliers - i.e., the eigenvalues of ψA (τ ) - lie in the open unit disc. As for stable LTI systems, a frequency response operator can be defined for stable LTP systems. The frequency domain representation for LTI systems is related with the fact that sinusoidal inputs are mapped into sinusoidal outputs at the same frequency with possibly different amplitude and phase. This notion can be extended to LTP systems by making reference to a more general class of signals: the set of exponentiallymodulated periodic signals, i.e., h (t) = ∑ hk esk t (5) k∈Z
where sk = s + jkΩ, Ω = 2π /T , which can be arranged as follows T H = . . . hT−2 hT−1 hT0 hT1 hT2 . . . . (6) Consider the following T -periodic stable system x(t) ˙ = A (t) x(t) + B (t) u(t) (7) y(t) = C (t) x(t) + D (t) u(t). The key idea, see Wereley and Hall [1990], further developed and formalized in Zhou and Hagiwara [2002a], is to expand the complex Fourier series of the dynamics matrices A (t) =
∑ Am e jmΩt
(8)
The angular kinematics and dynamics of the spacecraft are T modeled using as state variables the quaternion q = q1 qTR = T [q1 q2 q3 q4 ] describing the attitude of the body axes with respect to the orbital axes, and the inertial angular velocity vector T ω = [ωx ωy ωz ] , with respect to the body axes. With respect to the selected variables, the nominal, Earth-pointing, equilibrium T corresponds to the attitude quaternion q¯ = [1 0 0 0] and to the T angular rate ω¯ = [0 0 −Ω0 ] , where Ω0 is the orbital angular rate. Define the state vector x = [δ qTR δ ω T ]T formed with small displacements of the vector part qR of the attitude quaternion q T from the nominal values q¯R = [0 0 0] and small deviations of the body rates from the nominal values ω¯ x = ω¯ y = 0, ω¯ z = −Ω0 . Then the attitude dynamics can be linearized and a local linear model for the attitude can be defined as in Pulecchi et al. [2010]
m∈Z
and similarly for B (t) ,C (t) and D (t) and perform a harmonic balance approach. It was shown that defining .. . . .. .. . . . . . . . A0 A−1 A−2 . . . A = . . . A1 A0 A−1 . . . (9) . . . A A A . . . 2 1 0 .. .. .. . . . . . .
and similarly B, C and D, the harmonic transfer function G (s) that describes the input-output relation is 1
G (s) = C [sI − (A − N )]−1 B + D (10) where N = blkdiag { jnΩI} , n ∈ Z. Note that the operator defined in (10) is a doubly-infinite dimensional one. On the other hand, it is no longer time-varying and this fact can be exploited to obtain input-output information related with system (7) without numerical integration.
x(t) ˙ = Ax(t) + BT [Tmag (t) + Td (t)]
where Tmag is the magnetic control torque vector and Td is any disturbance torque vector. Taking into account that Tmag can be written as Tmag (t) = m(t) × b(t) = S(b(t))m(t) = −b(t) × m(t) T
equation (11) can be equivalently written as
x(t) ˙ = Ax(t) + Bcm (t)m(t) + BT Td (t),
In the following denoted with G (s) = A(t),B(t),C(t),D(t) .
(14)
where
3. MODELING AND CONTROL OF MAGNETICALLY ACTUATED SATELLITES
1
(12)
where b = [bx by bz ] is the geomagnetic field vector (in body frame), m is the dipole vector of the magnetic torquers and " # 0 bz −by S(b) = −bz 0 bx , (13) by −bx 0
For practical purposes, however, truncation must be, of course, considered. Given a stable system, performance measures such as the H∞ or the H2 norm can be thus computed using a suitably truncated frequency response operator, see, e.g., Zhou and Hagiwara [2002a] or Zhou and Hagiwara [2002b] for details.
In order to represent the attitude motion of an Earth-pointing spacecraft on a circular orbit the following reference systems are adopted:
(11)
0 Ω0 0 A= 0 0 0 0 0 0 BT = I −1 xx 0 0
−Ω0 0 0.5 0 0 0 0 0 0.5 0 0 0 0.5 0 0 0 0 0 Wx 0 Wy 0 0 −6ky Ω20 0 0 +6kz Ω20 0 0 0 0 0 0 0 0 0 , Bcm (t) = BT S(b(t)) 0 0 −1 Iyy 0 −1 0 Izz
(15)
(16)
Ixx −Iyy Izz −Ixx Iyy , kz = Izz , Wx = −kx Ω0 + kwx Ω, J Wy = −ky Ω0 −kwy Ω, kwx = Ixx , kwy = IJyy . Here, Ω is the nominal
and kx =
8480
Iyy −Izz Ixx , ky
=
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
wheel speed 2 . A very common structure for magnetic attitude control laws, which goes back to classical papers such as Martel et al. [1988a], consists of control laws of the kind 1 m(t) = ST (b(t))Tid (t) (17) kb (t) k2 where b(t) is the measurement, at time t, of the geomagnetic field b and Tid (t) is an "ideal" control torque to be determined on the basis of a suitable (usually static) feedback of state or output variables, according to the specific attitude control architecture of the considered spacecraft. The above control law can be readily given a simple geometric interpretation. Recalling (12) we can express the torque generated by the magnetic coils as 1 S(b(t))ST (b(t))Tid (t) (18) Tmag (t) = kb(t)k2 which can be easily interpreted as the projection of vector Tid onto the plane orthogonal to the direction of the magnetic field vector b (hence the name of "projection-based" controllers). 4. H∞ CONTROLLER DESIGN 4.1 Formulation of the design problem Given the satellite attitude dynamics discussed so far and assuming the value of the geomagnetic field b as well as the state vector x are measurable, the attitude dynamics can be written as a LTP system as follows: x˙ (t) = Ax (t) + BT S (b (t)) m (t) + BT Td (t) 1 = Ax (t) + BT S (b (t)) ST (b (t)) u (t) + BT Td (t) kb (t)k2 | {z } B(t)
= Ax (t) + B (t) u (t) + BT Td (t) .
(19) Under the above assumptions, we can design a static state feedback controller Tid = u = Kx. The aim of the controller is to stabilize the system and possibly optimize some given performance criterion. Since the model is obtained through linearization and many effects, such as a magnetic residual dipole, see Corno and Lovera [2009], are not taken into account, the controller should be able to attenuate the effect of any torque disturbance acting on the attitude variables. To this end, we can define the following variable to be used as a performance measure: Z (20) x (t) = Cσ x (t) , σ > 0 z (t) = σK
which implies ρ ψA+B(t)K (0) < 1, where by ρ (X) we denote the spectral radius of the matrix X. In other terms, the tuning of the control law u = Kx can be formalized as follows: σ minkGzT k subject to ρ ψA+B(t)K (0) < 1. (22) d ∞ K
It should be noticed that, since the cost function in (21) is in general non-smooth, gradient-based descent algorithms such as the ones adopted in Pulecchi et al. [2010], Vigano et al. [2010] might fail in converging to its minimum. In order to circumvent this difficulty, a dedicated approach to the solution of this problem is proposed, which is outlined in the following subsection. 4.2 Controller tuning algorithm
A significant amount of work has been carried out to cope with non-smooth optimization, yielding either deterministic or stochastic methods. In this work, we adopted a stochastic method to solve this optimization problem. Many randomized method exists, however, but the key issue is that a control law is optimal if in its neighborhood a better control law cannot be found (or, equivalently, can be found with null probability). From this simple observation, it is possible to design a procedure to optimize a certain criterion. Assume an initial stabilising controller K (i) is available, e.g., tuned using common techniques for LTI systems applied to the linearized time-averaged model, see e.g. Stickler and Alfriend [1976], or as the output of the previous iteration step. To compute a new controller K (i+1) , the idea is to pick-up randomly sampled controllers in the neighborhood of K (i) and select among them the best one in terms of the cost function in (22). When it is no longer possible to compute better controllers, the algorithm ends with the optimal one. The actual algorithm is further detailed in the following. For the sake of brevity, we
define the cost function to be min σ imized as J (K) = GzTd . As for the optimization algorithm, ∞ the controller K is regarded as a column vector.
Tuning algorithm Set the parameter 0 ≪ rmax < 1 to check the convergence, the number N ≫ 1 of M ONTE C ARLO samples to be generated and the initial step size µ .
where Z = [I3 0] and σ has been considered in order to limit the control effort 3 . A natural criterion to achieve an acceptable disturbance rejection is the H∞ norm of the input-output operator from the disturbance Td to the performance measure z, namely:
σ
GzT ≡ sup kzk2 , (21) d ∞ Td ∈L2 kTd k2
σ = [A + B (t) K, B ,C , 0]. where GzT T σ d The controller gain K can be tuned in order to minimize the cost function in (21) subject, of course, to the stability constraint, 2
The assumption of a diagonal inertia matrix is made only for ease of presentation. 3 Notice that zT z = xT Z T Z + σ 2 K T K x = xT Z T Zx + σ 2 uT u.
8481
(1) given a stabilizing controller K (i) generate N controllers using the following M ONTE C ARLO generation, namely for j = 1, . . . , N !
∆K (i)
(i) j j (i)
K ← K + µ K η + (23)
∆K (i) where η ∼ N (018×1, I18 ) and ∆K (i) = K (i) − K (i−1) ; (2) compute the corresponding norms, namely c j ← J K j ; ˆ (3) select the candidate controller K j (24) jˆ ← argmin J K j j
(4) compute an approximate of the rejection ratio r(i+1) o n ˆ (25) r(i+1) ← prob J K j > J K (i) ˆ (5) if ρ ψA+B(t)Kˆ (0) < 1 then set K (i+1) ← K j and i ← i+ 1 else decrease the step µ and goto (1); (6) if r(i+1) > r(i) decrease the step µ , else if r(i+1) < r(i) increase it;
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
• the H∞ norm of the frequency response operator from the disturbance torque Td to the orientation error δ qR = Zx
kδ qR k2 c1 = Gδ qR Td ∞ ≡ sup (28) Td ∈L2 kTd k2
(7) if r(i+1) ≥ rmax then stop and return Kσ ← K (i+1) , else goto (1); Notice that ∆K (i) /∆K (i) can be regarded as an approximation of the steepest descent direction, see Figure 1. Thus, the candidate controllers K j ’s generated within the step (1) are randomly selected by possibly taking into account the knowledge of the steepest descend direction estimated in the previous step. For
• the H2 norm of the frequency response operator from the disturbance torque Td to the magnetic dipole vector m
kmk∞ c2 = GmTd 2 ≡ sup (29) kT d k∞ Td ∈L∞ • the spectral radius of the characteristic multipliers of A + B (t) Kσ ρ = max λi ψA+B(t)Kσ (0) (30) i
which are reported in Tab. 1 as functions of the selected values of σ . The design model in (19) has been then simulated in Input-output norms c1 (σ ) c2 (σ )
σ 60 120 240
5. SIMULATION STUDY The considered spacecraft is of the type described in Section 3 and operates in a near polar orbit (87 deg inclination) with an altitude of 450 km and a corresponding orbital period of T = 5614.8 s = 2π /Ω0. For the given orbit, the following nominal (periodic) model for the components of the geomagnetic field has been considered: # " 7 cos(Ω0t) + 48 sin(Ω0t) −6 23 cos(Ω0t) − 2 sin(Ω0t) T. (26) b(t) = 10 5 The numerical values for the parameters of the linearized model are as follows: • satellite inertia tensor (kg m2 ): I = diag (Ixx , Iyy , Izz ) = diag (35, 17, 25)
ρ (σ )
12072.45 11039.75 9348.98
0.00283 0.00345 0.00598
Table 1. Computed norms for the closed-loop system.
Fig. 1. One iteration of the tuning algorithm for two decision variables K1 and K2 . the sake of completeness, note that the computed controller is actually a random matrix. To check whether the algorithm converges to the best controller, the procedure has been implemented and run several times, for each given value of σ , ending with a standard deviation in the cost function of about ±0.7%. In other words, the randomized operation of the optimization algorithm is such that the minimality of the cost function is guaranteed within a very small tolerance.
239.19 283.96 399.32
closed-loop using the control law u = Kσ x. The following initial perturbation has been applied to the angular rate vector: T δ ω (0) = 1 × 10−3 2 × 10−3 0 rad/s (31)
while a residual dipole of 1 Am2 along each axis has been considered as an unmeasurable disturbance, i.e., the following disturbance in (19) has been taken into account T Td (t) = S (b (t)) [1 1 1] .
(32) q Time histories of the norm of the orientation error δ qTR δ qR = √ √ xT Z T Zx as well as norm of the control dipole mT m are shown during the first orbit, i.e. t ∈ [0, T ], and during the following four orbits, t ∈ [T, 5T ], in Fig. 2, and Fig. 3, respectively. Figure 2 shows that for decreasing σ , the transient related with a non-zero initial condition is faster and more damped. On the other hand, increasing σ the required control effort increases, as well. As for the effects of disturbances, it can be noticed 0.03
0.02
0.01
0 0
(27)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.4 0.5 0.6 0.7 Normalized time (orbits)
0.8
0.9
1
Normalized time (orbits)
• momentum wheel inertia (kg m2 ): J = 0.01; ¯ = −200. • nominal wheel speed (rad/s): Ω
10
We recall that σ , which is actually a parameter of the optimization problem, has been introduced to weight the control effort in the cost function (21) and thus to limit the actuator energy. As a matter of fact, however, the selection of σ is not straightforward. Therefore, the approach discussed so far has been used (with N = 540, rmax = 0.995 and µ = 10−3 ) to generate three controllers Kσ , corresponding to σ = 60, 120 and 240. In order to compare the obtained controllers, the following quantities have been computed at the termination of the optimization algorithm:
5
0 0
0.1
0.2
0.3
Fig. 2. Norm of the orientation error δ qR and of the magnetic dipole m for σ = 60 (blue), σ = 120 (green), σ = 240 (red).
8482
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
quaternion) are strongly coupled and their transients correspond to the response of the stabilised precession and nutation modes of the open loop dynamics, while the pitch axis (q4 component of the attitude quaternion) is stabilised much more efficiently thanks to the above mentioned favourable controllability properties. 1
1 q
that, for increasing σ , the rejection of the disturbance torque is weaker, see Figure 3, while a better disturbance rejection does not require an additional effort in the control signal. This is somehow confirmed by the data in Table 1. In fact, as one can see, as the weight on the control effort σ increases, the H∞ norm of the frequency response from the disturbance torque to the attitude error, i.e. c1 , rapidly increases, as well. On the other hand, the sensitivity of the H2 norm of the frequency response from the disturbance torque to the magnetic dipole, namely c2 , is less remarkable, justifying the behaviour in Figure 3.
0.9998 0.9996
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
0.02 q2
The simulations presented in the following have been carried out using an object-oriented environment for satellite dynamics (see Lovera [2006] for details) developed using the Modelica language (Fritzson and Bunus [2002]). More precisely, a full nonlinear simulation of the coupled rigid body orbital and attitude dynamics has been performed and the following models of the space environment have been implemented: the JGM3 spherical expansion for the geopotential as a gravitational model, the Harris-Priester model for the atmospheric density distribution (see Montenbruck and Gill [2000] for details) and the International Geomagnetic Reference Field (IGRF, see Wertz [1978]) for the Earth’s magnetic field (up to order 10). Disturbance torques due to gravity gradient (including J2 effects), magnetic residual dipole (assuming a residual dipole of 1 A m2 along each spacecraft body axis) and solar radiation pressure (computed using the solar coordinates formulas given in Montenbruck and Gill [2000]) have been taken into account in the simulation. Finally, all the controllers have been implemented in digital form, using a conventional sample and hold scheme, using a sampling period Ts = 1 s.
0 −0.02
q3
0.05 0 −0.05
0 −3
4
5
x 10
q
0 −5
0 −3
x 10
2 ω1 [rad/s]
1 0 −1 −2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
−3
x 10
4 ω2 [rad/s]
2 0 −2 −4
0 −5
x 10
5 ω [rad/s]
The simulation results obtained for the controllers corresponding to σ = 60 and σ = 240 are presented, respectively, in Figures 4 and 5. As can be seen from the Figures, both controllers damp out the effect of the initial angular rate perturbation (31) in less than one orbit (an adequate performance for a magnetic attitude control scheme), with the controller for σ = 60 providing a settling time of about a quarter of the orbit period. The system is then brought to its steady state response under the effect of the cyclic external disturbance torques. As expected, the yaw and roll axes (q2 and q3 components of the attitude
3
0 −5 −10
0
2
mx [A m ]
2 0 −2 −4
2
my [A m ]
2
−3
x 10 8
1 0 −1 −2
6 2
m [A m ]
10 5 0
z
4 2
−5 −10
0 1
1.5
2
2.5
3
3.5
4
4.5
5
Fig. 4. Quaternion, angular rates and control dipole moments for σ = 60.
Normalized time (orbits)
1.5
1
Finally, a sensitivity analysis of the closed-loop performance to the choice of the sampling period in the implementation of the control laws has been carried out. The detailed results have been omitted for brevity, but the sampling period can be increased up to Ts = 20 s without any significant performance degradation.
0.5
0 1
1.5
2
2.5
3
3.5
4
4.5
5
Normalized time (orbits)
6. CONCLUSIONS Fig. 3. Norm of the orientation error δ qR and of the magnetic dipole m for σ = 60 (blue), σ = 120 (green), σ = 240 (red).
The problem of designing attitude controllers for magnetically actuated spacecraft has been considered. An approach to the
8483
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
q
1
1 0.9998 0.9996
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
0.02 q2
0 −0.02 −0.04 0.04 q3
0.02 0 −0.02
q4
0.01 0 −0.01
0 −3
ω1 [rad/s]
2
x 10
1 0 −1 −2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
−3
ω2 [rad/s]
4
x 10
2 0 −2 −4
0 −5
ω3 [rad/s]
5
x 10
0 −5 −10
0
m [A m2]
1
x
0 −1 −2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5 Orbits
3
3.5
4
4.5
5
0.5 m [A m2]
0
y
−0.5 −1 −1.5
z
m [A m2]
2 0 −2 −4
Fig. 5. Quaternion, angular rates and control dipole moments for σ = 240. tuning of projection based controllers has been proposed, relying on robust periodic state feedback control techniques. The performance of the proposed control algorithm has been illustrated in a simulation study. REFERENCES Bittanti, S. and Colaneri, P. (2009). Periodic systems - filtering and control. Springer. Boyd, S. and Barratt, C. (1991). Linear controller design: limits of performance. Prentice Hall. Corno, M. and Lovera, M. (2009). Spacecraft attitude dynamics and control in the presence of large magnetic residuals. Control Engineering Practice, 17(4), 456–468. Fritzson, P. and Bunus (2002). Modelica - a general objectoriented language for continuous and discrete-event system modelling and simulation. In Proceedings of the 35th IEEE Annual Simulation Symposium.
Hablani, H. (1995). Comparative stability analysis and performance of magnetic controllers for bias momentum satellites. Journal of Guidance, Control and Dynamics, 18(6), 1313– 1320. Lewis, A. (2007). Nonsmooth optimization and robust control. Annual Reviews in Control, 31, 167–177. Lovera, M. (2006). Control-oriented modelling and simulation of spacecraft attitude and orbit dynamics. Journal of Mathematical and Computer Modelling of Dynamical Systems, Special issue on Modular Physical Modelling, 12(1), 73–88. Martel, F., Pal, P.K., and Psiaki, M. (1988a). Active magnetic control system for gravity gradient stabilised spacecraft. In 2nd Annual AIAA/USU Conference on Small Satellites, Logan (Utah), USA. Martel, F., Pal, P., and Psiaki, M. (1988b). Active magnetic control system for gravity gradient stabilised spacecraft. In 2nd Annual AIAA/USU Conference on Small Satellites, Logan (Utah), USA. Montenbruck, O. and Gill, E. (2000). Satellite orbits: models, methods, applications. Springer. Psiaki, M. (2001). Magnetic torquer attitude control via asymptotic periodic linear quadratic regulation. Journal of Guidance, Control and Dynamics, 24(2), 386–394. Pulecchi, T., Lovera, M., and Varga, A. (2010). Optimal discrete-time design of three-axis magnetic attitude control laws. IEEE Transactions on Control Systems Technology, 18(3), 714–722. Silani, E. and Lovera, M. (2005). Magnetic spacecraft attitude control: A survey and some new results. Control Engineering Practice, 13(3), 357–371. Stickler, A.C. and Alfriend, K.T. (1976). An elementary magnetic attitude control system. Journal of Spacecraft and Rockets, 13(5), 282–287. Tempo, R., Calafiore, G., and Dabbene, F. (2005). Randomized algorithms for analysis and control of uncertain systems. Springer. Vigano, L., Bergamasco, M., Lovera, M., and Varga, A. (2010). Optimal periodic output feedback control: a continuous-time approach and a case study. International Journal of Control, 83(5), 897–914. Wereley, N.M. and Hall, S.R. (1990). Frequency response of linear time periodic systems. In IEEE International Conference on Decision and Control, CDC. Wertz, J. (1978). Spacecraft attitude determination and control. D. Reidel Publishing Company. Wood, M., Chen, W.H., and Fertin, D. (2006). Model predictive control of low earth orbiting spacecraft with magnetotorquers. In IEEE International Conference on Control Applications, Munich, Germany. Yan, H., Ross, M., and Alfriend, K. (2005). Three-axis magnetic attitude control using pseudospectral control law. In Proc. of the 2005 AAS/AIAA Astrodynamics Specialists Conference, Lake Tahoe, USA. Paper AAS 05-417. Zhigljavsky, A. and Zilinskas, A. (2008). Stochastic global optimization. Springer. Zhou, J. and Hagiwara, T. (2002a). Existence conditions and properties of the frequency response operators of continuoustime periodic systems. SIAM Journal on Control and Optimization, 40(6), 1867–1887. Zhou, J. and Hagiwara, T. (2002b). H2 and H∞ norm computations of linear continuous-time periodic systems via the skew analysis of frequency response operators. Automatica, 38(8), 1381–1387.
8484