Robust attitude control of a spacecraft with flexible dynamics

(1)

STOCKHOLM SWEDEN 2019,

Robust attitude control of a

spacecraft with flexible dynamics

FRANCESCO GIULIANO

KTH ROYAL INSTITUTE OF TECHNOLOGY

(2)

List of Figures

1.1 Artistic representation of a full-electric spacecraft with orientable thrust . 2

1.2 Star tracker. Courtesy of OHB Sweden. . . 3

1.3 Thruster used for PRISMA mission. Courtesy of OHB Sweden. . . 5

1.4 Spacecraft body frame . . . 9

1.5 Max and min singular values of G(s) . . . 10

1.6 RGA number for decentralize controller on the main diagonal . . . 12

1.7 PID model in Simulink . . . 13

1.8 Standard closed-loop form used for design . . . 13

1.9 Current design. Reference response (top) and torque command (bottom) for the x-axis. . . 15

1.10 Current design. Disturbance rejection (top) and torque command (bottom) for the x-axis . . . 15

1.11 Current design. Reference response (top) and torque command (bottom) for the y-axis. . . 16

1.12 Current design. Disturbance rejection (top) and torque command (bottom) for the y-axis . . . 16

1.13 Current design. Reference response (top) and torque command (bottom) for the z-axis. . . 17

1.14 Current design. Disturbance rejection (top) and torque command (bottom) for the z-axis . . . 17

1.15 Stability analysis considering one loop at the time. Inputs are w1, w2. Outputs are z1, z₂. . . 19

2.1 Representation of system A with input w and output z. . . . 23

2.2 Standard closed-loop form. . . 24

2.3 Generalized system description with weights . . . 25

2.4 Characteristic loci for a real system. . . 26

2.5 Example of WP and WT. On top we can see the unfeasible formulation of weights where at the same frequency ¯σ(S) + ¯σ(T ) ≥ 1. At the bottom a correct formulation. . . 28

2.6 Use of different γ on closed-loop transfer functions. On the left γ ≤ 1 means that the weight will always be respected all over the spectrum. On the right γ ≥ 1 shows how the closed-loop transfer function overflows the weight for at least a frequency. . . 29

(5)

2.8 Design trade-off for L(jω) . . . 36

2.9 Closed-loop between the controller and the shaped plant . . . 37

2.10 Left coprime factorization of the perturbed set of plants Gp . . . 39

2.11 Example of closed-loop configuration with prefiltering Fr . . . 41

2.12 Representing uncertainty of system G(jω) in Nyquist plot . . . 43

2.13 Complex region (continuous line) enclosing the real uncertainty (dashed line) 44 2.14 Complex bounded regions contour the nominal plant with discs of radius |W_a(jω)| . . . 45

2.15 Plant frequency response for 2 setpoints. When SADM is 0°(parallel to X). When SADM is 90°(parallel to Z). Notice axis swapat 90°. . . 48

2.16 Multiplicative input uncertainty realization for plant G(s) . . . 49

2.17 Uncertainty weight Wi(jω) covering the relative error. After the first peak, the weights copes for unmodelled uncertainty . . . 49

2.18 Uncertainty weight Wi(jω) magnitude plot. Notice bandwidth limitations for X-Z axis at 0.36 rad/s and Y axis at 2.5 rad/s. . . 50

3.1 Generalized plant form for mixed sensitivity . . . 53

3.2 σ(W_T⁻¹) vs σ(W_i⁻¹). Note that condition kWik_∞<kW_Tk_∞ is respected. . 54

3.3 Singular values of W_u⁻¹ . . . 55

3.4 Singular values of S . . . 56

3.5 Singular values of T . . . 56

3.6 Singular values of KS . . . 57

3.7 Reference step response on 3-axis . . . 57

3.8 Input disturbance step response on 3-axis . . . 58

3.9 Torque command on 3 axis. No saturation. . . 58

3.10 Singular values of KS against W_u⁻¹ . . . 59

3.11 Singular values of T against WT . . . 60

3.12 Singular values of G(s) . . . 61

3.13 Second order filter with low damping ratio - Magnitude plot . . . 63

3.14 Second type Chebyshev filter - Magnitude plot . . . 63

3.15 Shaped plant Gs(jω) . . . 64

3.16 Singular values of S (Glover-McFarlane) . . . 65

3.17 Singular values of T (Glover-McFarlane) . . . 65

3.18 singular values of KS (Glover-McFarlane) . . . 66

3.19 Reference step response on 3-axis . . . 67

3.20 Input disturbance step response on 3-axis . . . 68

3.21 Torque command on 3 axis. Max control input of 0.86 Nm on X-axis. . . 68

3.22 Robust stability condition is respected since kT Wik∞<1 . . . 69

3.23 N∆ structure for robust stability . . . 71

3.24 Scaled version of M∆ structure . . . 72

3.25 µ for Robust controller from Mixed Sensitivity . . . 76

3.26 µ for Robust controller from Glover-McFarlane . . . 76

4.1 Simple closed-loop model in Simulink, with saturation of actuators . . . . 78

(6)

4.2 Reference response on the x-axis - PID . . . 80

4.3 Disturbance rejection on the x-axis - PID . . . 80

4.4 Reference response on the y-axis - PID . . . 81

4.5 Disturbance rejection on the y-axis - PID . . . 81

4.6 Reference response on the z-axis - PID . . . 82

4.7 Disturbance rejection on the z-axis - PID . . . 82

4.8 Reference response on the x-axis - Mixed Sensitivity . . . 83

4.9 Disturbance rejection on the x-axis - Mixed Sensitivity . . . 83

4.10 Reference response on the y-axis - Mixed Sensitivity . . . 84

4.11 Disturbance rejection on the y-axis - Mixed Sensitivity . . . 84

4.12 Reference response on the z-axis - Mixed Sensitivity . . . 85

4.13 Disturbance rejection on the z-axis - Mixed Sensitivity . . . 85

4.14 Reference response on the x-axis - Glover McFarlane . . . 86

4.15 Disturbance rejection on the x-axis - Glover McFarlane. . . 86

4.16 Reference response on the y-axis - Glover McFarlane . . . 87

4.17 Disturbance rejection on the y-axis - Glover McFarlane. . . 87

4.18 Reference response on the z-axis - Glover McFarlane . . . 88

4.19 Disturbance rejection on the z-axis - Glover McFarlane . . . 88

4.20 Example of integration for the robust celestial controller in Simulink. . . . 90

4.21 Control error in Spacecraft Body frame. Notice the computer swap at 10 min, for 45 s. PID controller . . . 92

4.22 Absolute Pointing Error. Focus After Reboot. PID controller . . . 93

4.23 Reaction Wheels torque on Spacecraft Body. PID controller . . . 94

4.24 Control error in Spacecraft Body frame. Notice the computer swap at 10 min, for 45 s. Robust MS controller . . . 95

4.25 Absolute Pointing Error. Focus After Reboot. Robust MS controller . . . 96

4.26 Reaction Wheels torque on Spacecraft Body. Robust MS controller . . . . 97

4.27 Control error in Spacecraft Body frame. Notice the computer swap at 10 min, for 45 s. Robust GM controller . . . 98

4.28 Absolute Pointing Error. Focus After Reboot. Robust GM controller . . . 99

4.29 Reaction Wheels torque on Spacecraft Body. Robust GM controller . . . 100

A.1 Representation of z = Aw . . . 110

A.2 Which system norm depending on the input. . . 111

B.1 Control error in Spacecraft Body frame. Focus after reboot. PID controller113 B.2 Absolute Pointing Error. Complete simulation. PID controller . . . 114

B.3 Torque commands. PID controller . . . 115

B.4 Estimated torque applied on Spacecraft Body. PID controller . . . 116

B.5 Control error in Spacecraft Body frame. Focus after reboot. Robust MS controller . . . 117

B.6 Absolute Pointing Error. Complete simulation. Robust MS controller . . 118

B.7 Torque commands. Robust MS controller . . . 119

B.8 Estimated torque applied on Spacecraft Body. Robust MS controller . . . 120

(7)

controller . . . 121 B.10 Absolute Pointing Error. Complete simulation. Robust GM controller . . 122 B.11 Torque commands. Robust GM controller . . . 123 B.12 Estimated torque applied on Spacecraft Body. Robust GM controller . . . 124

(8)

Acronyms

FES Full-Electric Spacecraft

NASA National Administration for Astronautics & Aeronautics ESA European Space Agency

LSS Large Space Strucure

LTI Linear Time Invariant system DPS Distributed Parameters System LTI Linear Time-Invariant system MoI Moment of Inertia

SADM Solar Array Driver Mechanism LMI Linear Matrix Inequalities SGT Small Gain Theorem EP Electric Propulsion SSV Structured Singular Value LFT Lower Fractional Transformation UFT Upper Fractional Transformation RS Robust Stability

AOCS Attitude and Orbital Control System APE Absolute Pointing Error

SSC Swedish Space Corporation

(9)

When this Thesis started I never imagined It would have been such a gratifying journey.

The approach to Multivariable feedback control requested me to pause, ponder and investigate the subject in such a way that let me to understand also the context that brought to develop these new theories, about 25 years ago. This inspired me consistently when I tried to apply robust control to the attitude problem in space.

Nonetheless, this work has demanded a lot of time, efforts and trust, I was very luck to have had my two supervisors Per and Elling on my side that advised me during this year and had always kept their office door open. Furthermore I want to thank the great AOCS team of OHB Sweden that, from the first moment, made me feel in a warm environment.

They have always been open to discussion and gave me precious suggestions on how to design a controller for space applications. Thank you: Tristan, Sten, Nicolas, Per, Erik and Kristian. I am also grateful to my friends that have always been patients whenever I had last minute work and could not join them. Last but not least I want to thank Alessia for being my constant in an highly nonlinear world.

(10)

Abstract

The employment of Electrical Propulsion, to reduce the cost of satellites, has also increased the use of large solar arrays. Such solar panels are associated with low-frequency resonant modes, and to ensure fine pointing of the satellite, the attitude control should not excite these modes. This is a difficult task in such complex systems, whose exact a priori knowledge is largely unavailable during the early design phase.

The primary scope of this thesis is to synthesize a robust algorithm for attitude control that deals with modeling errors.

The first part is dedicated to understanding the principal causes of uncertainty that will inevitably limit the design. Once a description of the uncertainty has been obtained, two different approaches to robust control algorithms have been investigated. The outcome is the syntheses of two multivariable controllers with different performance and stability margins. Since this project has been carried out during the development of a new full electric spacecraft concept by OHB, it has been possible to validate the two controllers in the on-board software simulator of the spacecraft. The results confirm the robust controllers as possible alternatives to the current design.

(11)

Införandet av Elektrisk Framdrivning för att reducera kostnaden för satelliter har också ökat användningen av stora solpaneler. Sådana stora solpaneler är associerade med lågfrekventa resonanta moder och för att säkerställa satellitens prestanda får attitydkon- trollen inte exitera dessa moder. Detta är en svår uppgift i sådana komplicerade system för vilka exakt a priori kunskap till stora delar inte finns tillgänglig under projektets tidiga designfas.

Det här examensarbetet handlar i första hand om att syntetisera en robust algoritm för attitydkontroll vilken hanterar modelleringsfel. Första delen av arbetet avser att förstå vilka de primära osäkerhetskällor som begränsar systemets prestanda. Efter att en osäkerhetsbeskrivning har tagits fram undersöks två olika tillvägagångssätt för att ta fram robusta styralgoritmer. Resultatet är synteser av två olika multivariabla regulatorer med olika prestanda och stabilitetsmarginaler. Eftersom projektet har genomförts samtidigt som OHB har utvecklat en satellitplattform baserad helt på Elektrisk Framdrivning, har det varit möjligt att validera de två regulatorerna i ombordmjukvarans simulatormiljö.

Resultaten bekräftar att de robusta regulatorerna är möjliga alternativ till nuvarande design.

(12)

Introduction

The advent of space exploration has widened the employment of large satellites.

If on one side, companies like Bigelow Space are redefining the concept of inflatable structures, other players of space panorama, like SpaceX, are planning their way to Mars by means of huge launchers, like BFR (Big Falcon Rocket). Moreover,NASA andESA are considering to invest in a lunar gateway: a series of space stations that will serve as staging posts for mission to Moon and Mars.

At the present, private companies are committing to Electric Propulsion (EP) as a cheaper way to allocate a satellite into orbit. The challenges of using this kind of propulsion derives from two factors: the uncertainty of low thrust, from an orbit determination point of view, and the significant amount of power to be used by the engine, which requires large solar panels [24]. OHB is currently developing a telecommunication satellite which is part of a new concept of Full-Electric Spacecraft (FES), meaning that transfers, momentum management and station keeping will be performed using EP. Such a large satellite, whose solar array can have more than 50 m², serves as a platform for telecommunication via antennas, that might be bigger than 2 m in diameter.

The main characteristic of a large space structure are the low damped resonant frequencies that, in some cases, represent a major limitation to the attitude and orbit control.

The challenge starts when the spacecraft needs to be controlled in the presence of uncer- tain parameters that are directly connected to the flexible dynamics, such as resonant frequencies themselves, damping ratio, etc. Things complicate even more if one thinks that some mass properties like Moments of Inertia are not well-known at early stage of the project, because either the payload has not been decided yet or the design is likely to change.

How is it possible to design the attitude control of the spacecraft in presence of such uncertainty? The answer stands in Multivariable Feedback control; in particular, to a branch of it called Robust control. Where Robust means that the system is stable and of high performance in presence of uncertainty.

This thesis work has the scope to develop a robust attitude controller for aFES.

The methodology used in this thesis derives from the work of Doyle, Skogestad, Stein, and their research groups in the early nineties, where robust control opened the path to use of Large Space structures [1].

(13)

formal, but still tries to catch the meaning of mathematical entities in their own language.

The work is mainly addressed to a student that has a knowledge in control theory, but interested in space technology. Moreover, I wrote the following pages in first person plural because, despite I did most of the work, it would not have been possible to accomplish it without the help ofAOCS department. So it’s a "we".

The thesis opens with a brief description on telecommunication satellite technology andAOCS hardware. It follows a mathematical model for the spacecraft, that along with the description of the state of art, represents the starting point for analysis and benchmark of results.

In the second chapter, the methodology for designing in H∞ space is provided. It concerns two design methods that lead to two robust controllers, synthesized in the third chapter. Even the uncertainty will be addressed in H∞ space and consequently it will be part of robustness assessment in chapter three.

The report ends with validation of both controllers via Monte Carlo simulations. A last chapter is reserved for conclusions and future work. Eventually, the thesis comes with an appendix that may help the reader. In fact, the intent is to keep the person focused and this is why, in the PDF version, the hyperlink is available.

(14)

Chapter 1

A Full-Electric Spacecraft

Theoretically, a Large Space Strucure (LSS) is described as a Distributed Parameters System (DPS), which is an infinite order system. Because of its nature, a LSS, has many resonant modes at low frequencies, which are lightly damped and appear to be grouped in specific parts of the frequency domain. Since testing and prediction of this behavior is quite difficult on Earth one should take this into account for intrinsic uncertainty of the system when it comes to attitude control. OHB’sFESis one example of LSS, because of its large appendages, like the solar array and antennas.

In the first section a brief introduction to attitude control of telecom satellites and why it is demanding in terms of pointing performance. In the second section we will try to carry out a model reduction ofDPS, in a way that we can obtain a finite order system with its relative resonance frequencies and use them to make aLTIrepresentation of the spacecraft with flexible dynamics.

In the third section, the current attitude control of the satellite will be introduced along with its control objectives. The chapter ends pointing out the shortcomings of using this type of controller and leave the floor to H∞ control.

1.1 The attitude control for a telecommunication satellite

When we turn on the TV and zap to our favorite football match, we do not mind about from how far the signal comes. We do not worry about the precision with which the satellite is pointing to Earth, that is in the order of mdeg. We just care about our football team.

Satellite telecom technology is based on the simple idea that from above, everything appears smaller, and therefore it is easy to irradiate a bigger area than it would be if we were to use an antenna on the ground. Having such a "mirror" in space is as much expensive as challenging.

We can see the problem in two ways, first to put the satellite in geostationary orbit, where its revolution around Earth is synchronous with the planet, and second to ensure a precise pointing for all the mission.

(15)

Usually, a telecom satellite is massive, compared to the ones in Low Earth Orbit¹, it weighs few tons. This is due to the fact that it needs to carry a big payload (mostly reflectors), but, most of all, the propellant used to bring the satellite to geostationary orbit occupies more than 60 % of the total mass. In fact, chemical propulsion produces thrust in the order of 10⁴ N and allows for reaching the orbit in some weeks, but it is not efficient². ([11])

If there are no time constraints on the orbit transfer, another type of propulsion can be considered. It is the case to use electrons instead of hydrazine. A noble gas like Xe is discharged by a strong electromagnetic field, forming plasma, that is then ejected at high velocity. The efficiency increases drastically, but the thrust produced is in the order of millinewtons. Therefore, the mass of the propellant can be also reduced, from thousand of kilos in case of chemical, to hundreds for electrical propulsion. The counter effect is that the requested power for the plasma thruster grows consistently, in some cases up to 5 kW just for the propulsion subsystem. This leads to bigger solar panels that give a major contribution to flexible dynamics and make the attitude control more difficult.

In this thesis, we present a new concept of spacecraft which not only uses electrical propulsion for orbital transfer, but it is also capable to orientate the thrust vector in order to perform orbital manoeuvres in geostationary orbit. To do so, theFEShas 4 Hall Effect thrusters and big solar panels, up to 50m² (bigger than a penalty area in a regular soccer pitch). In fig.1.1and artistic representation of this type of satellite during orbit transfer.

Figure 1.1: Artistic representation of a full-electric spacecraft with orientable thrust

1Low Earth Orbit (LEO) is defined from Earth’s surface in the range 100 - 2000 km. Medium Earth Orbit (MEO) from 2000 km to 35000 km. Geostationary orbit is at 36000 km.

2For electrical propulsion the energy to accelerate the mass is taken from the sun. While for chemical propulsion all the energy is stored in the propellant itself, and therefore carried through the manoeuvre.

(16)

1 – A Full-Electric Spacecraft

The control problem for a spacecraft is not trivial. Without digging into mathematics, a free-falling body is highly non-linear process, because of its governing equations (or Euler equations of motion, see [19]-[16]), control rater and saturation constraint limits. The first missions were actually spin stabilized, like Explorer 1 fromNASAor Odin from Swedish Space Corporation (SSC), while for telecom satellite is mandatory a 3-axis stabilization, where the input is a torque, provided by the actuators, and the measured variable is an angle vector. Usually, a closed-loop control is preferred over an open-loop, since the former can deal with uncertainty of the system. It is common practice to choose a non-linear feedback control for large slews, while having a linear control algorithm for Earth pointing attitude. The latter is the objective of this thesis and therefore will extensively treated further in text. Here we focus on the hardware in the loop, namely sensors and actuators, with which the spacecraft can be controlled.

Sensors

Common sensors for Attitude control are: Star tracker, Gyroscope, Sun Presence sensor, GPS and Magnetometer (only in LEO). The reader is probably familiar with most of them because they are used in many fields. However, the Star tracker and Sun Sensor are a specialty of space sector. The former is basically a camera with and integrated map of the sky. It takes a picture to the stars and compares it with an internal catalogue. Hence, it back calculates the attitude of the spacecraft using distances and angle of constellations. Generally a satellite mounts two Star trackers for redundancy (in Fig.1.2).

The sun sensor is a photocell, that recognizes if the satellite has a panel exposed to sunlight. It consist of one, or more photocells, which give an output current proportional to the intensity of the light falling on it. It is not possible to reconstruct a precise attitude from it, but the sun sensor is very useful during the detumble phase or in case a failure occurs and the satellite should point the sun.

Figure 1.2: Star tracker. Courtesy of OHB Sweden.

(17)

Actuators

In space application, for actuator, we mean a device that can exchange momentum. Com- mon actuators are: Reaction wheels, Control Moment Gyro (CGM), Thrusters, Magnetic Torquesin LEO. Each of these components can provide relatively high or low torque levels, therefore it is not unusual to find more than one in the same satellite. In case ofFES, we have four Reaction Wheels³, Cold Gas Thrusters and 4 Electrical propulsion Thrusters. The reaction wheel is a spinning flywheel that accelerates or decelerates basing on the torque command. It is the principal device used for attitude manoeuvres, and the one that usually has more torque capability. However, there is a limit on the rotational speed, and so to the momentum provided to the spacecraft⁴. Beyond this limit is not possible to control the spacecraft, hence the wheels need to be unloaded. This operation is demanded to the thrusters.

In detumbling phase and in case of emergency, the Cold Gas thrusters are activated.

They eject pressurized gas, like a spray paint, and create a torque around the center of mass (in fig.1.3). During operational life in geostationary orbits, the time of events is much longer respect to the launch phase or transfer, thus it is possible to use a smaller torque but continuously, in order to unload the wheels. This is the purpose of electrical propulsion thrusters. It is more efficient (which translates in decreasing propellant consumption) and also the thrust can be oriented, by means of a boom, to make small in-orbit adjustments.

In the next section we will provide a linear model for attitude control of the spacecraft in geostationary orbit. First we will obtain the resonant frequencies of the body and subsequently use them to derive a simpler transfer function to be used in the control loop.

3One for each axis plus one for redundancy

4It is also possible to think of the reaction wheel on the other way around, as a momentum absorber.

In fact, if the spacecraft is having a momentum around a certain direction, the reaction wheel can start to spin in the opposite direction, making the total momentum of the spacecraft to be 0.

(18)

Figure 1.3: Thruster used for PRISMA mission. Courtesy of OHB Sweden.

(19)

1.2 LTI model

A physical body can be described by means of non-linear partial differential equations (PDE). The system depends on distributed parameters and it is expressed by

m(x)∂²u

∂t²(x, t) + D0

∂u

∂t(x, t) + A0u(x, t) = F (x, t) (1.1) where u(x, t) is the vector of instantaneous displacements off its natural equilibrium points, due to a force F (x, t) that includes both control forces and disturbances. m(x) is the mass density, which is dependent on the structure . A0u is the internal restoring forces which contains the smooth functions that solve the modal eingenvalue problem. Basically it carries information about resonant frequencies and respective modes⁵. D0∂u

∂t is the damping term which contains both the system’s damping ratios and gyroscopic damping if the body is rotating. A state space representation for this problem is given in [1].

It is clear that controlling the infinity number of modes for (1.1) is the best thing can ever happen, but it is unrealistic. An infinite order controller is not realizable. Rather we can try to reduce the model by capturing only the most interesting modes and evaluate residuals such that they will not affect stability.

An effective method to reduce the order of the model is to use Finite Elements (FEM), to derive a linear description of the model by means of spatial discretization of the continuum body. In eq. (1.2) a solution to the DPSproblem is given, whose accuracy⁶ depends on the number and nature of elements N. An extensive description of FEM and DPScan be found in [1], [17] and [14].

u(x, t) =

N

X

k=1

θ_k(x) ˆqk(t) (1.2)

where θk(x) are the shape functions, i.e. hermitian-cubic polynomials connecting consecu- tive elements, while ˆqk(t) is the displacement vector expressed in local coordinates . Notice that a linear combinations of θk(x) approximates the actual system normal modes of vibration. If we substitute (1.2) in (1.1) and let qN(t) = [ ˆq₁(t), ˆq₂(t) . . . ˆq_N(t)]^T, we get

MNq¨N+ DNq˙N + KNqN = B⁰Nf (1.3) Where MN is the mass matrix, KN the stiffness matrix and B⁰Nf the external forces applied to the elements.

Often, it is convenient to visualize (1.3) with a new set of coordinates describing intuitively what is the displacement vector associated with a particular natural frequency.

Namely, we swap to modal coordinates.

Assume that UN is a unitary matrix for all elements N, such that U_N^TMNUN = IN U_N^TKNUN = ΛN

5Mode here means the shape that the body acquires at that specific resonant frequency, for instance bending or torsional mode.

6In terms of relative error

(20)

where Λ is diagonal,. A transform from local coordinates qN to modal ones is uN = U_N^Tq_N, which gives

¨uN + DN˙uN + ΛNuN = BNf

This is an ODE formulation of (1.1), with diagonal ΛN, whose entries are the squared approximate mode eingenfrequencies. The columns of UN, multiplied by the respective shape function, forms the approximate mode shape.

Having a reduced order formulation of a more complex problem, as it is the continuum body (orDPS) is an advantage that we pay with some residual errors for the vibrational modes that we do not consider. However, they will be taken into account during the design of the robust controller as neglected dynamics. Nonetheless, the most important information from the modal analysis is the ΛN, which provides the resonant frequencies of the body. In fact, the next paragraph will use ΛN as input for a simplified version of the transfer function used for attitude control.

Satellite transfer function

For attitude control design, it is convenient to adopt a simplified model of the plant. It is furthermore appropriate to have it as transfer function, so we can easily visualize its properties in terms of resonant modes, gains and uncertainty.

The spacecraft is a system that rotates under an input torque vector τ to give an angle vector ϕ. In Laplace domain this is expressed by

ϕ= G(s)τ (1.4)

where G(s) ∈ R^3×3 and τ ∈ R^3×1 in Nm. Attitude angles can be expressed in different ways. The most common one is to use quaternions, a set of 4 complex parameters which are bounded to the unitary hypersphere⁷, that provide a redundant, nonsingular, attitude description. They are well suited to describe arbitrary, large rotations [19]. However, they cannot be used for computational purposes since they are a set of four, while all the measurements and errors have dimension three. Thus, in order to use quaternions, we need to assume small angle rotations, so that only the vectorial part is taken into account.

This choice limits the computation of control error, especially when the reference is an attitude profile to be given from the ground station.

The satellite is a rigid system, embodied by

G(s) = (J(s)s²)⁻¹ (1.5)

where J ∈ R^3×3 is the Inertia of the spacecraft. It both includes the flexible dynamics and the rigid one; while it is convenient to split their respective contribute in two parts.

Rewrite eq.(1.5) in

G(s) = (Jrs²)⁻¹G_f(s) (1.6)

7A formal definition of quaternion is q(θ) = [1/cos(θ/2) 1/sin(θ/2)ex1/sin(θ/2)ey1/sin(θ/2)ez].

Which corresponds to q = q0+ q1~i+ q2~j+ q3~k. Where e = [exeyez] is the unitary rotation vector.

(21)

where Jr reflects the rigid part, while Gf(s) is the normalize inverse MoI that contains the flexible dynamics. We can show this by

G_f(s) = JrJ(s)⁻¹

In order to model the resonant peaks, given by the modal analysis computed before, we use the second order system:

G_f(s) = (s²+ 2ζfω_fs+ ωf²)[(1 − ρ)s²+ 2ζfω_fs+ ωf²]⁻¹ (1.7) where ωf is the resonant frequency, ζf the damping ratio and ρ it is the element that characterizes the flexible dynamics of appendages Ja respect to the rigidMoI, namely is the ratio

ρ= JaJ_r⁻¹

Eventually, more modes can be captured by means of the following notation

J(s) = Jr 1 −

N

X

k=1

ρ_ks²(s²+ 2ζkω_k+ ωk²)⁻¹

!

(1.8)

System analysis

Consider the spacecraft in its body frame: XSCB (or roll axis) is outward of the East facing panel, YSCB(or pitch axis) is parallel to the Solar Array Driver Mechanism (SADM), ZSCB (or yaw axis) is outward of the Earth facing solar panel and positive towards the top of the satellite. In fig.1.4a representation of this.

(22)

Figure 1.4: Spacecraft body frame

(23)

Depending on the situation, the spacecraft can enter different configurations which are called control modes, for instance if it has been released from the Launch Vehicle or if it is looking at the sun after exiting a failure. We will concentrate our analysis on the normal mode where the satellite will spend most of its time and where it has the most restrictive requirements. It covers when the spacecraft is orbiting the Earth in Geostationary orbit and broadcasts the signal.

Values used in this thesis work are selected to allow publication and correspond to the order of magnitude for numbers that would be used in an early design stage of a real project. Thus, the reader should trust the following text and understand that, despite numbers are similar, they are fabricated ad hoc to not make any reverse engineering on them.

Now, consider the spacecraft in geostationary orbit and at the beginning of life, meaning that the fuel left after orbit raising maneuvers is around one third of initial mass in the tank. An order of magnitude for Jr is

J_r=







27 0.02 −0.024

0.02 4 −0.1

−0.024 −0.1 26





·10³ [Kg · m²]

This matrix is important for two reasons. Firstly, it gives an hint on the contribution of large solar panels on the principal inertia. Then we note how the off-diagonal terms of Jr

are smaller respect the diagonal, meaning that the system is quite symmetrical. In control language we say that the matrix is not ill-conditioned and preliminary we can assume the use of decentralize control. Further in text we provide a mathematical justification for that.

In fig. 1.5, the singular values of G(s) are depicted.

Figure 1.5: Max and min singular values of G(s)

(24)

Around 0.7 rad/s the first resonant peaks occurs. After that more modes with less contribute in terms of gain take place (six in total). In section2.4, we will see how the flexible modes limit the design.

The last consideration to be done concerns the effect of cross coupling on the control- lability of the system both for controller design and uncertainty description. One on the most used technique that provides a measure of coupling interactions is Relative Gain Array (RGA), as described in [2]. Here is a short survey.

Consider uj and yi to be the input-output pair of the system G(s), whatever dimension is considered, and say we want to control yi from uj. Two extreme conditions appear:

• All other loops open: _∂u^∂yⁱ_j

uk=0,∀k /=j = gij

• All other loops closed: _∂u^∂y_jⁱ

yk=0,∀k /=i= ˆg_ij

Here gij = [G]ij represents how much the input ui influences the output yi, whereas ˆg_ij is the inverse of the ji’th element of G⁻¹, namely how much the output yi is influenced by the input uj. The ratio between these two limit cases is the relative gain of the ij’th element, i.e. the input-output relation among channels of G(s).

λij = g_ij ˆ

gij = [G]ij[G⁻¹]ji (1.9)

RGA is an interesting tool to understand which inputs are preferable to be controlled. In fact, a λ close to 1 is desired, while a change of sign in the gain between high frequency and steady state must be avoided. RGA of G shall be computed for all frequencies. Instead of reporting many matrices, here we chose to represent RGA by using the RGA number Γ , which is defined as

Γ(G) =kΛ(G) − Iksum

RGA number is a measure of suitability for the pairing of the controller along the main diagonal and I − Γ indicates a better off-diagonal fit. For the satellite we can see how the principalMoI are disposed along the main diagonal. Thus, it is fair to assume that a decentralized controller could be used. We can see this by computing Γ (G) for the spacecraft. The closer Γ is to 0 the better, see fig.1.6.

Since Γ is relatively small, a decentralized controller on the main diagonal is suggested.

A low RGA number, along with a small condition number⁸, allows for using a multiplicative uncertainty at the input, even though physically would be modeled at the output.

8It is defined as the ratio between the biggest singular value and the smallest one γ = ^¯^σ(G)_σ(G). A big condition number indicates an ill-conditioned system

(25)

10^-2 10^-1 10⁰ 10¹ 10² 10³ 10⁴ 10⁵ Frequency rad/s

0 0.01 0.02 0.03 0.04 0.05 0.06

RGA number

|| - I ||

sum for Diagonal pairing

Figure 1.6: RGA number for decentralize controller on the main diagonal

(26)

Figure 1.7: PID model in Simulink

1.3 The Celestial controller

In geostationary orbit, during its nominal phase, the spacecraft shall point to the Earth, with an accuracy specified by the costumer. To comply these requirements, the spacecraft currently adopts a decentralized controller on the main diagonal, justified by the previous analysis. The controller is a PID⁹, as shown in fig.1.7, whose tunable parameters we omit here.

Not only pointing requirements should be respected, but the satellite shall cope with all disturbances that act on it, both from the environment (like orbit perturbations, solar wind etc.) and internal forces. In fig.1.8, a block diagram used for design is depicted.

− +

+ + ++

W_n G_d

K G

r

d

y

yn

e uk u

n

Figure 1.8: Standard closed-loop form used for design

The biggest disturbance considered during the design phase is originated from the reaction wheel assembly. One can note that there is no output disturbance considered, even though we know that some external forces acting. Nevertheless, for easiness’ sake we neglect those effects because they are very small compared to the input disturbance,

9The transfer function of PID is in the standard form KP ID= Kp+^K_s^I + KDs, where KP, KI and KD are the tuning parameters of proportional, integrative and derivative parts respectively

(27)

represented by a friction phenomenon in the reaction wheel. This is a mechanical problem that is due to the type of lubrication used in space. In fact, liquid lubricant is not indicated at very low pressures, because it creates bubbles. Therefore a solid lubricant should be used. This has a friction characteristic which is function of the angular velocity of the wheel and it is higher when it is approaching zero speed. Since the reaction wheel specifications are usually well known from the manufacturer, it is possible to estimate this friction phenomenon and to compensate for it. The latter is captured by Gd(s) whose is mimicking the estimation transient

G_d=1 − 1 τ s+ 1

I where the time constant τ = 50 s.

The noise is captured by the weight Wn whose transfer function is

Wn= 200s²+ 2 · 1 · 0.2s + 0.2² s²+ 2 · 1 · 40s + 40² One last remark on the controller structure, which is

K = KP ID=







K₁₁ 0 0

0 K₂₂ 0

0 0 K33







The main objective of this thesis work is to design a controller which is both robustly stable and possibly has more performance than the current design.

In a more rigorous way we can summarize the following goals:

• Robust stability

• Bandwidth ωb as high as possible

• Reject in-orbit disturbances

• Minimize actuator energy

The PID has been chosen as benchmark to be compared with the robust controller.

Hereafter is shown the time domain analysis of the closed-loop in Simulink, with a saturation of reaction wheels at 40 mNm (same approach as in fig.4.1). In particular, we have analyzed the system response to a step and an input disturbance of 5 mNm. In figs. (1.9,1.11,1.13) we see the reference step response for the three axis. In figs. (1.10, 1.12,1.14) the disturbance step response. The simulation is performed over 100 cases where small changes in some real parameters of the plant occurs.

(28)

Figure 1.9: Current design. Reference response (top) and torque command (bottom) for the x-axis.

Figure 1.10: Current design. Disturbance rejection (top) and torque command (bottom) for the x-axis

(29)

Figure 1.11: Current design. Reference response (top) and torque command (bottom) for the y-axis.

Figure 1.12: Current design. Disturbance rejection (top) and torque command (bottom) for the y-axis

(30)

Figure 1.13: Current design. Reference response (top) and torque command (bottom) for the z-axis.

Figure 1.14: Current design. Disturbance rejection (top) and torque command (bottom) for the z-axis

(31)

Both X and Z axis have similar performance, both in terms of overshoot in response and disturbance rejection, while the Y axis is faster, since it has less inertia. It is worth mentioning the oscillatory behavior of the torque command that, in a long term, represents a threat to the mechanical functioning of the actuator and therefore should be avoided.

Eventually, the delays of the closed-loop, both considering software and hardware, are around 300 ms and do not represent a major control limitation, because the rate sampling of the controller already occurs at 1 Hz. However, the delay will be taken into account indirectly when a robustness assessment will be derived for the controllers. It will be also included in the closed-loop used for simulations.

Motivation for robustness

We have showed previously that the satellite is a quite symmetrical system, for which the off-diagonal contribute ofMoIis smaller than the principal axis one, at least for low frequencies. Therefore, it is possible to use a decentralized regulator for attitude control, which does not excite the flexible dynamics. It is a common design procedure, for the PID, to choose the tunable parameters as we were designing a SISO per axis, under the assumption that each axis is not interacting with the others.

But, how much is this method sensible to any kind of variation in the system structural parameters? Is stability always guaranteed?

Consider the following example (from [21]):

Example 1. A satellite is stabilized around its principal axis whose transfer function between the torque input and the angular velocity is

G(s) = 1 s²+ a²

"

s − a² a(s + 1)

−a(s + 1) s − a²

#

a= 10 which is translated in a minimal state-space realization:

G=







0 a 1 0

−a 0 0 1

1 a 0 0

−a 1 0 0







from A we see that the plant has two poles at s = ±ja and therefore needs to be stabilized.

Let a decentralized controller be

K= I the closed loop is

A_cl= A − BKC =

"

−1 0 0 −1

#

we see that the eingenvalues of A are two at -1 and the system is nominally stable.

We now want to study robust stability and assume the system to be perturbed in each of its inputs. Those kind of perturbations are always present since they are connected with uncertainty parameters of the actuators. If we break one loop at the input and

(32)

assume that each channel has such small coupling with the others to be considered independent, as in fig.1.15, we have that the transfer function from input w1 to output z1

is t11(s) = _s+1¹ = _1+L^L¹^(s)₁_(s). Because the loop gain is L1(s) = 1/s it allows for an infinite gain margin and a phase margin of 90°. The same results hold for w2 and z2. The system appears to be very robust on each axis, despite the value of a.

G

K w1

z₁

w₂ z₂

+

+ −

−

Figure 1.15: Stability analysis considering one loop at the time. Inputs are w1, w₂. Outputs are z1, z2.

Consider an input uncertainty on the two channels to be

˜

u1= (1 + 1)u1 u˜2 = (1 + 2)u2

Where 1, 2 are real numbers. Because they occur at the input we can rewrite B as B˜ =

"

1 + 1 0 0 1 + 2

#

the corresponding closed-loop ˜A_cl = A − ˜BKC will have the characteristic polynomial given by

det(sI − ˜A_cl) = s²+ 2 + 1+ 2

| {z }

a1

s+ 1 + 1+ 2+ (a²+ 1)1₂

| {z }

a0

which has poles in LHP only if a0, a1 are positive. This condition is verified if the uncertainty gain is positive and just one channel at the time is perturbed.

If we consider simultaneous perturbations in both channels, for example 1 = −2, the system is unstable (a0<0) for

|₁| > √ 1

a²+ 1 ≈0.1

Thus, even for a 10% additive perturbation from the nominal value, the system becomes unstable.

This example has shown how a decentralized controller that is designed upon principal inertia can lead to poor stability in presence of uncertainty at the input. Nonetheless one can argue that the controller can be tuned finer, so that is more conservative and do not enter the instability region for a given set of perturbations. By doing so, we renounce

(33)

to performance, but increase robustness. Eventually, the lack of knowledge in the design phase of the project requires extensive tests, in order to prove robust stability and, in some cases, it forces the engineer to retune the controller. Hence another approach, which consider multivariable feedback control, is needed.

(34)

Chapter 2

Methodology

In this chapter an overview on optimal control design in the H∞ framework will be provided. The first section will open with a brief introduction on historical reasons that led to develop a robust control theory in frequency domain. We will then move to the optimization problem definition and, in particular, we will focus on the choice of weights that make the controller respecting some given performances. In the second section we will investigate another way for designing a robust controller which is closer to the classical rules of loop shaping.

The last section is dedicated to uncertainty in MIMO systems, with focus on two ways for describing it. The chapter ends with an uncertainty characterization for the spacecraft.

2.1 Design in the H

∞

space

As described in [14], dynamics and control of mechanical system were a unique subject during XVII century, but with at the beginning of 1900 they set themselves apart. In fact, what we now call "classical control theory" was born independently during 1930, when pioneers like Nyquist and Bode defined the usage of transfer functions in the frequency domain, or input-output domain. At that time, it was mostly applied to electrical systems and didn’t find much application in control of mechanics & structures.

The prospective changed during the space race that took place in between 1950s - 1970s. Where key figures as Bellman and Kalman developed an optimal control theory in state space, where the problem would have been to minimize signals in time domain.

During this period, some controller like LQR and LQG¹ or Kalman filter made their way.

Yet there were two major flaws in using these theories. Although they were MIMO based, there was not obvious connection of state space with frequency domain, where fundamental limitations and uncertainty were noticeable. On the other hand, it was a common practice to design against robustness by adding more stochastic disturbance to the system, which does not guarantee any stability margins, as concisely stated by Doyle in [4]. To better visualize this concept, here is an easy example.

1Linear Quadratic Regulator and Linear Quadratic Gaussian regulator

(35)

For LQG control we have

y= Gu + w

where w is white noise, as independent variable, that will not appear in the closed loop (I + GK)⁻¹. But in real systems we know that

y= Gu + Eu + Gdd

where E is the model uncertainty and Gdthe disturbance process. The closed loop becomes (I + (G + E)K)⁻¹. From the latter we can see how a certain amount of uncertainty is

able to shift the poles on the RHP, leading to instability.²

It was in the 80’s that the problem to address uncertainty in input-output domain became an hot trend. In 1981 on the IEEE Transactions on Automatic Control, Doyle and Stein arose the issue of how to achieve a good feedback control for MIMO in presence of unstructured uncertainty [5]. Later Glover & McFarlane showed how to apply loop-shaping ideas by means of singular values in frequency domain. [9]. Those theories became very important for years running, especially because they address a problem in frequency domain while solving it in state space, as a signal minimization problem.

This thesis is based on these doctrines.

2.1.1 An optimal control problem formulation

Consider a strictly proper LTI system G(s). The H∞norm is defined as kG(s)k∞= sup

ω ¯σ(G(jω)) (2.1)

It means that the H∞norm is the maximum gain of the system in the worst case direction, i.e. the maximum singular value over frequencies. We can see this concept under different lights. Let us consider two applications of this norm, one in the frequency domain and the other in time domain.

The size of a vector, as discussed in AppendixA, is defined as its 2-norm (or Euclidean norm)

kek₂=^s^X

i

|e_i|²

Let us take two vectors w and z to be respectively the input and the output of matrix A, as in fig.A.1.

The gain is given by the ratio between output and input. However, we are dealing with a matrix, which means that the maximum gain is the one induced by the worst input direction. If we use the Euclidean norm as a measure of w and z, the induced matrix norm will be

kAk₂ = max

w /=0

kAwk₂

kwk₂ = ¯σ(A) (2.2)

2Notice that (I + (G + E)K)⁻¹= (I + L)⁻¹= S

(36)

2 – Methodology

Figure 2.1: Representation of system A with input w and output z.

This is in principle what happens to the gain of system (G(jω)) at any frequency. Then H∞ norm is the peak of ¯σ(G(jω)) all over the spectrum.

Let us now consider an analogy in time domain. In fact, we know that a persistent sinusoidal input w(t) = A sin(ωst+ φi) = w(ωs) ³ will excite the system to give an output z(ω) at t → ∞. The gain at ωs will be kz(ωs)k2/kw(ωs)k2 or ¯σ(G(jωs)), as in (2.2). But in general the following holds

¯σ(G(jω)) = max

w(ω) /=0

kz(ω)k2

kw(ω)k2

Here you can notice the similarity with (2.1), which is actually its representation in frequency domain. However, the most powerful use of the H∞ has to be sought among the properties of norm itself, below defined.

Definition 2.1. A norm of a matrix has to be:

1. Non-negative: kAk > 0 2. Positive: kAk = 0

3. Homogeneous: kα · Ak =|α|·kAk α ∈ C 4. Triangle inequality:

kA+ Bk ≤kAk+kBk 5. Multiplicative property:

kABk ≤kAk·kBk

The last one is actually what makes the H∞ norm to be different from H2 and ρ(A). In fact, because H2 and ρ(A) do not comply with property n.5, they are not considered norms. As a consequence, they cannot be used in the Small Gain Theorem (SGT). This is why we use H∞norm to characterize the size of uncertainty that the system can afford before going into instability. In Control language, we need to ensure

kM(jω)∆(jω)k∞<1 (2.3)

3Written in phasor notation

(37)

Where M is the Closed-loop between exogenous signal inputs w and outputs z, and ∆ is the uncertainty description of the system. We will come back later to this notation.

We now proceed to formulate an optimal control problem which respect to the objectives defined in Ch.1.

For a standard configuration of the plant as in fig.2.2, that here we recall

− +

+ + ++

Wn

G_d

K G

r

d

y

yn

e uk u

n

Figure 2.2: Standard closed-loop form.

In order to respect the objectives, we need that some closed-loop transfer functions will have a particular shape in a certain region of the frequency domain. To do so, we need to include some weights that constrain these functions. There are different approaches to address these kind of problems and all of them share a signal based technique. Here we show the most common practice, called the Mixed sensitivity. In fact, we want to stack our objectives in a minimization matrix such that they are in minimum number, because the more objectives we include in the minimization process the more difficult will be for a numeric solution to converge. Basically, we want a system N, say the closed loop that reflects the objectives, to be

N =

WPS W_TT W_UKS

(2.4)

This system can be obtained from input w = r to output z = [e y u]^T, see fig.2.3.

N incorporates the plant, the objectives and the controller, but in practice we do not have K before solving the problem. So, what we actually compute is a generalized form description of the plant that will be used as baseline in the minimization iteration and whose interconnection with K yields to eq. (2.4).

Consider splitting the model in two parts as in fig.2.3, K ∈ R^u×v is the controller and P ∈ R^z×w is the plant.

(38)

2 – Methodology

Figure 2.3: Generalized system description with weights

Rewrite P by inspection as the transfer functions between [w u]^T and [z v]^T

P =

0 W_UI

0 −W_TG WPI −W_PG

−I −G

(2.5)

the generalized form of (2.5) is equivalent to the partition in state space

z = P11w+ P12u (2.6)

v = P21w+ P22u (2.7)

where the matrix P11 is the pure plant (with no controller) and P22 has dimensions compatible with K, i.e. P22∈ Rⁿ^v^×n^u.

The closed-loop is created by inserting the controller

u= Kv (2.8)

After some algebraic manipulations of equations (2.6-2.8) we obtain

N = P11+ P12K(I − P22K)⁻¹P21= Fl(P, K) (2.9) Here Fl(P, K) stands for Lower Fractional Transformation (LFT) of P , with K as a parameter. By means either Upper or Lower fractional transformation is possible to trace back almost any interconnected system, from the one that has originated it. This is one of the reason why is a common practice to standardize the plant into a generalized form and then close the loop throughLFT.

In a similar fashion we can link the uncertainty to the N system, for both analysis and synthesis. In ch.3 we will be using the Upper Fractional Transformation (UFT) instead, where we will analyze the stability margins of our controller.

(39)

Weights selection

Now that we have an idea on how to translate the design objectives into frequency domain we have one last ingredient missing in order to make seething the H∞ iteration, namely defining the weights.

There is no specific literature that helps us determine WP, W_T and WU, although some material can be found in [10], [18], [15]. In this thesis work we will follow two main rules to define weights. One has to do with the fact that MIMO have natural limitations that cannot be overcome by using any kind of boundaries, the other follows the raison d’être of weight, i.e. a function that reflects the objective addressed by the user in the

design phase.

Starting from the limitations, recall that

e= Sr

Ideally we would like |S(jω)| = 0 ∀ω but this is not possible for real system. In fact the loop gain L(jω) of any stable and proper system will have an excess pole number of at least 2, which translates in crossing the real negative axis of Nyquist plot at least once.

We can better visualize this in fig.2.4.

−1 0

Im

Re

Figure 2.4: Characteristic loci for a real system

Robust attitude control of a spacecraft with flexible dynamics