Open access peer-reviewed chapter

# Modeling and Attitude Control of Satellites in Elliptical Orbits

Written By

Espen Oland

Submitted: March 1st, 2018 Reviewed: July 20th, 2018 Published: November 5th, 2018

DOI: 10.5772/intechopen.80422

From the Edited Volume

## Applied Modern Control

Edited by Le Anh Tuan

Chapter metrics overview

View Full Metrics

## Abstract

The attitude determination and control system (ADCS) for spacecraft is responsible for determining its orientation using sensor measurements and then applying actuation forces to change the orientation. This chapter details the different components required for a complete attitude determination and control system for satellites moving in elliptical orbits. Specifically, the chapter details the orbital mechanics; perturbations; controller design; actuation methods such as thrusters, reaction wheels, and magnetic torquers; actuation modulation methods such as bang-bang, pulse-width modulation, and pulse-width pulse-frequency; as well as attitude determination using vector measurements combined with mathematical models. In sum, the work describes in a tutorial manner how to put everything together to enable the design of a complete satellite simulator.

### Keywords

• attitude control
• attitude estimation
• thrusters
• reaction wheels
• magnetic torquers
• elliptical orbits
• PWM
• PWPF
• bang-bang
• Sun vector model
• magnetic field model
• sliding surface control
• quaternions
• angular velocity

## 1. Introduction

The problem of developing attitude determination and control systems (ADCS) has received much attention in the last century with general books such as [1, 2, 3, 4], as well as description of individual ADCS designs for different satellites in works such as [5, 6, 7, 8]. While Refs. [1, 2] can be considered excellent foundation books within the topic of spacecraft modeling and control, there is a need for a more concise presentation of the attitude control problem and how this can be modeled in a simple manner, both not only as a tutorial for new researchers but also to give insight into the different components required for ADCS design drawing on ideas and results from previous works as mentioned above.

This chapter is an extension of [9] and builds on much of the previous work, as well as the research done through the HiNCube project as presented in [10, 11, 12]. This work considers the problem of designing a complete ADCS system comprising all the required components. Figure 1 shows the control structure and the required signal paths, giving an overview of the contents in this chapter, as each block is described in detail to put the reader in position to design their own ADCS system. The required sensors for this system are a magnetometer to measure the Earth’s magnetic field, a gyro to measure the angular velocity of the satellite, as well as Sun sensors for measuring the direction toward the Sun. Further, the mathematical models used together with sensor measurements to determine the attitude of the satellite requires a real-time clock, as the time and date are required to know the direction toward the Sun as well as what the magnetic field looks like at a given day and time. Comparing the sensor measurements with the mathematical models allows for the determination of the attitude of the satellite, something that is done using the Madgwick filter as presented in [13]. With an estimated attitude obtained using the Madgwick filter, the attitude can be controlled to point a sensor onboard the satellite in a desired direction, something that is solved in this chapter using a PD+ attitude controller, calculating the desired torques required in order to make the attitude and angular velocity errors go to zero. In order to create the desired moments, this chapter presents how this can be achieved using a number of different actuators, namely, magnetic torquers, reaction wheels, and thrusters. The orbital mechanics block describes how the satellite moves in its elliptical orbit, while the perturbing forces and moments block describe how the different perturbations affect the satellite. Simulations show the performance of the different methods and should put the reader in a position to simulate and design new attitude determination and control systems for satellites in elliptical orbits.

## 2. Mathematical modeling

### 2.1. Notation

This subsection is similar to the author’s previous works, e.g., [9, 14]. Let ẋ=dx/dtdenote the time derivative of a vector, while the Euclidean length is defined as x=xΤx. Superscript denotes frame of reference for a given vector. The rotation matrix is denoted RacSO3=RR3×3:RΤR=IdetR=1, which rotates a vector from frame ato frame cand where Idenotes the identity matrix. The inverse rotation is found by taking its transpose, such that Rca=RacΤ. The angular velocity of frame crelative to frame areferenced in frame eis denoted ωa,ce, and angular velocities can be added together as ωa,fe=ωa,ce+ωc,fe(cf. [15]). The derivative of the rotation matrix is defined as Ṙac=RacSωc,aawhere Sdenotes the cross product operator, which is defined such that for two vectors v1,v2R3, Sv1v2=v1×v2, Sv1v2=Sv2v1, Sv1v1=0, and v1ΤSv2v1=0.

Quaternions can be used to parameterize the rotation matrix, where qc,aS3=qR4:qΤq=1denotes the quaternion representing a rotation from frame ato frame cthrough the angle of rotation ϑc,aaround the axis of rotation kc,a. The inverse quaternion is defined as qa,c=ηc,aεc,aΤΤ, also sometimes denoted as q. A quaternion comprises a scalar and a vector part, where ηc,adenotes the scalar part, while εc,aR3denotes the vector part. This allows the rotation matrix to be constructed using quaternions as Rac=I+2ηc,aSεc,a+2S2εc,a. Composite quaternions can be found through quaternion multiplication as qc,e=qc,aqa,e=Tqc,aqa,ewith

Tqc,a=ηc,aεc,aΤεc,aηc,aI+Sεc,a.E1

The use of the quaternion product ensures that the resulting quaternion maintains the unit length property. The quaternion kinematics is defined as

q̇c,a=12qc,a0ωc,aa=12Tqc,a0ωc,aa.E2

For attitude control, several different frames are needed:

Inertial: The Earth-centered inertial (ECI) has its origin in the center of the Earth, where the xiaxis points toward the vernal equinox and the zipoints through the North Pole, while yicompletes the right-handed orthonormal frame. The inertial frame is denoted by Fi.

Orbit: The orbit frame has its origin in the center of mass of the satellite (cf. [16], p. 479). The eraxis coincides with the radius vector riR3, which goes from the center of the Earth to the center of mass in the satellite. The ehaxis is parallel to the orbital angular momentum vector, pointing in the normal direction of the orbit. The eθcompletes the right-handed orthonormal frame where the vectors can be described as er=riri, eθ=eh×er, and eh=hhwhere h=ri×ṙi. The orbit frame is denoted by Fo.

Body: The body frame has its origin in the center of mass of the satellite, where its axes coincide with the principal axes of inertia. The body frame is denoted by Fb.

Desired: The desired frame can be defined arbitrarily to achieve any given objective (cf. [17]). The desired frame is denoted by Fd.

### 2.2. Orbital mechanics

This section describes how the orbit frame can be related to the inertial frame through the six classical orbital parameters, and for more details, the reader is referred to [1]. Specifically, the objective with this subsection is to find the radius, velocity, and acceleration vector of the orbit, as well as its angular velocity and acceleration. From well-known orbital mechanics, the six classical parameters can be defined as the semimajor axis a, the eccentricity e, the inclination i, the right ascension of the ascending node Ω, the argument of the perigee ω, and the mean anomaly M. The distance to the apogee and perigee from the center of the Earth can be defined, respectively, as raand rp, allowing the semimajor axis to be found as a=ra+rp2and the eccentricity of the orbit as e=rarpra+rp, while the mean motion can be found from n=μa3. Here, μ=GMEarth, where Gis the gravitational constant, while MEarthis the mass of the Earth. With knowledge of the mean motion, the mean anomaly can be found as M=ntt0=ψesinψwhere ψis the eccentric anomaly, tis the time, and t0is the time when passing the perigee. To properly describe where in the orbit the satellite is located, the true anomaly can be found as θ=cos1cosψe1ecosψ, while its derivative can be found as θ̇=n1+ecosθ21e232([18], p. 42). When running a simulation, it is desirable to have a continuously increasing true anomaly, while the direct method will map the angle between 0and 180°. Instead, by integrating the derivative overtime, a smooth true anomaly can be found that increases continuously. The eccentric anomaly, however, cannot be found analytically, but can be found through an iterative algorithm as described in ([1], p. 26) ψkt=Mt+esinψk1t, where kis the iteration number. This algorithm is valid as long as 0<e<1, which holds for elliptical orbits. From these calculations, the rotation matrix from inertial frame to orbit frame can now be constructed as

Rio=cosω+θcosΩcosisinω+θsinΩcosω+θsinΩ+sinω+θcosicosΩsinω+θsinisinω+θcosΩcosisinΩcosω+θsinω+θsinΩ+cosω+θcosicosΩcosω+θsinisinisinΩsinicosΩcosi.E3

The radius, velocity, and acceleration vector can be defined in the orbit frame, respectively, as ([1], pp. 26–27)

ro=acosψaeasinψ1e20Τ,E4
vo=a2nrsinψa2nr1e2cosψ0Τ,E5
ao=a3n2r2cosψa3n2r21e2sinψ0Τ,E6

where r=riis the length of the radius vector. Each of these vectors can be rotated to the inertial frame using Eq. (3), such that ri=Rioro, vi=Roivo, and ai=Roiao. The angular velocity of the orbit frame relative to the inertial frame can be found as ωi,oi=ri×viriΤri, while the angular acceleration can be found through differentiation as

ω̇i,oi=ri×airiΤri2ri×viviΤririΤri2.E7

In order to implement the orbital mechanics in, e.g., a Simulink framework, the required input to the subsystem would be the time (t). Further, the orbital parameters must be defined as given in Table 1 and can be changed depending on the orbit. These constants allow for the calculations of the eccentricity (e), the semimajor axis (a), and mean motion (n). With the mean motion, the mean anomaly (M) can be found and used to approximate the eccentric anomaly (ψ) using the iterative algorithm presented above. The rate of change of the true anomaly (θ̇) can also be found and by integration enables the calculation of true anomaly (θ). All these values allow for calculating the rotation matrix from the orbit frame to the inertial frame (Roi), the radius vector (ro), the velocity vector (vo), the acceleration vector (ao), angular velocity (ωi,oi), and angular acceleration vector (ω̇i,oi). Hence, all the outputs from this subsystem can easily be found following this procedure and will be used in several other subsystems.

ParameterValueUnit
G6.674081011m3kg1s2
MEarth5.97421024kg
raRe+1200km
rpRe+800km
Re6378km
i75°
Ω0°
ω0°

### Table 1.

Parameters required for calculation of the orbital dynamics.

## 3. Attitude dynamics and control

### 3.1. Attitude dynamics

The attitude dynamics can be derived with the basis in Euler’s moment Equation ([1], p. 95). The angular momentum of a rigid body in the body frame is given as

hb=Jωi,bb,E8

where JR3x3is the inertia matrix, while ωi,bbR3is the angular velocity of the body frame relative to the inertial frame. The angular momentum can be found in the inertial frame as

hi=Rbihb.E9

The rate of change of angular momentum is equal to the total torque, such that τi=ḣi. Hence, by differentiating Eq. (9), it is obtained that

τi=ḣi=RbiSωi,bbhb+Rbiḣb,E10

which can be written in the body frame by using Eq. (8) as τb=Sωi,bbJωi,bb+Jω̇i,bb, where the inertia matrix is assumed to be constant. Decomposing the total torque into an actuation component and a perturbing component, τb=τab+τpb, allows the rotational dynamics to be written as

Jω̇i,bb=Sωi,bbJωi,bb+τab+τpb,E11

where τabR3denotes the actuation torques (e.g., output from reaction wheels), while τpbR3denotes the perturbing torques (e.g., gravity torque). Further, by using quaternion representation, the update law for the quaternion representing the attitude of the body frame relative to the inertial frame can be written as

q̇i,b=12Tqi,b0ωi,bb.E12

Hence, Eqs. (11) and (12) serve as governing equations describing the attitude and angular velocity of the satellite. The inputs that affect these values are the perturbation and actuation torques, where the latter will be found in the following sections.

### 3.2. Error dynamics

From Euler’s moment equation, the angular acceleration is defined relative to the inertial frame. For attitude control, it is often more interesting controlling the attitude and angular velocity relative to the orbit frame. For that reason, the angular velocity of the body frame relative to the orbit frame can be found as ωo,bb=ωi,bbRibωi,oi, which can be differentiated as

Jω̇o,bb=Sωi,bbJωi,bb+τab+τpb+JSωi,bbRibωo,iiJRibω̇i,oiE13

giving a description of the attitude dynamics relative to the orbit frame. It is also possible to find the error dynamics, to enable tracking of a desired attitude and angular velocity. Let qo,d,ωo,dd,ω̇o,ddLdenote a desired quaternion, angular velocity, and acceleration; then, the quaternion and angular velocity error can be found as

qd,b=qd,oqo,b,E14
ωd,bb=ωo,bbRdbωo,dd,E15

with the kinematics as

η̇d,b=12εd,bΤωd,bb,E16
ε̇d,b=ηd,bI+Sεd,bωd,bb.E17

The angular acceleration error can be found by differentiating Eq. (15) as

Jω̇d,bb=Sωi,bbJωi,bb+τab+τpb+JSωi,bbRibωo,iiJRibω̇i,oi+JSωo,bbRdbωo,ddJRdbω̇o,dd.E18

Hence, the control objective can be defined as that of making qd,bωd,bb00, which will make the satellite point in a desired direction and move with a desired angular velocity.

### 3.3. PD+ attitude controller

Takegaki and Arimoto [19] proposed in 1981a simple method for position control of robots, something that was extended by [20] to enable trajectory tracking. The so-called PD+ controller has been applied for spacecraft by [21, 22] showing good results. The author has applied this method in previous works such as [23, 24].

In order to design a PD+ attitude controller, let a Lyapunov function candidate be chosen as V=12ωd,bbΤJωd,bb+kp1ηd,b2+kpεd,bΤεd,bwhere kpis a positive scalar gain. The derivative is found by using Eqs. (16)(18) as

V̇=kpεd,bΤωd,bb+ωd,bbΤSωi,bbJωi,bb+τab+τpb+JSωi,bbRibωo,iiJRibω̇i,oi+JSωo,bbRdbωo,ddJRdbω̇o,dd.E19

A PD+ attitude control law can now be chosen as

τdb=JRdbω̇o,ddJSωo,bbRdbωo,dd+JRibω̇i,oiJSωi,bbRibωo,iiτpb+Sωi,bbJωi,bbkpεd,bkdωd,bb,E20

where kdis another positive scalar gain and τdbdenotes the desired torque required to make the attitude and angular velocity errors go to zero. Assuming no actuator dynamics, i.e., τab=τdb, and then by inserting Eq. (20) into Eq. (19), it is obtained that V̇kdωd,bb2, which is negative semidefinite. By applying the Matrosov theorem (cf. [24]), it can be shown that the origin (εd,b,ωd,bb) = (0,0) is uniformly asymptotically stable.

The inputs to the control law (Eq. 20) are the desired states qo,d, ωo,dd, and ω̇o,dd, which are to be defined by the reader, e.g., as part as a guidance block depending on the mission objective. The inertia matrix (J) is assumed to be known, while the angular velocity vector between the body frame and orbit frame (ωo,bb) can be found as described above. The other angular velocities are found from the orbital mechanics, while the rotation matrices are found as composite rotations, e.g., Rib=RobRio, or by using the relationship between the quaternions and rotation matrices directly (cf. Section 2A). The error quaternion and angular velocity are found from Eqs. (14) and (15), while the perturbing torques will be described in the following section.

## 4. Perturbing torques

There are different kinds of perturbing torques, such as gravity torque, aerodynamic torque, magnetic field due to the electronics inside the satellite, as well as solar radiation torque. This section only considers the gravity torque. In [16], p. 147, the gravity torque is defined as

τgb=3GMEarthr5ri×Jri,E21

where the terms have previously been defined. As can be seen from this equation, non-diagonal inertia matrices will induce gravitational torques to align the satellite with the gravity field. This is also sometimes used for passive control, using e.g., a gravity boom to ensure that one side of the satellite is always facing the Earth. For this chapter, the perturbing toque is set equal to the gravity torque such that τpb=τgb.

## 5. Actuators

The control signal must be mapped to an actuator that must generate the desired torque. With limitations in actuation, the saturation must be modeled in order to obtain realistic results when simulating attitude control. This section considers three types of actuators commonly used for spacecraft attitude control: magnetic torquers, reaction wheels, and thrusters.

### 5.1. Magnetic torquers

Magnetic torquers operate by creating a local magnetic field that interacts with the Earth’s magnetic field. In simple terms, magnetic torques can be explained as that of a compass needle. By applying current through a coil, a local magnetic field is created, which will try to align itself with the Earth’s magnetic field. This allows the attitude of a spacecraft to be changed and is a very popular approach for small satellites, e.g., cubesats. One of the drawbacks or challenges with magnetic actuation lies in the fact that the Earth’s magnetic field goes from the North Pole to the South Pole as shown in Figure 2.1 As can be seen, when the satellite crosses the North Pole, there will be mainly a downward magnetic field component, reducing the possibility of actuation to only two axes, and similarly along the equator. This subsection is based on [12] and will describe how to model magnetic torquers and how it can be applied for attitude control. A magnetic torquer produces a magnetic torque by applying a current through a coil, which can be expressed as [2].

τmb=Smbbb,E22

where mb=mxmymzΤ=NAixiyizΤis the induced magnetic field, Nis the number of turns of the coils, iis the current around a given axis, and Ais the area of the coils. The Earth’s magnetic field is represented through the vector bb, meaning that to use magnetic actuation for attitude control, a good model of the Earth’s magnetic field is needed.

From a control point of view, the physical parameters Aand Nare defined when the spacecraft is designed, such that the controller needs to dictate which currents that must be sent to the torquers in order to get a desired torque. This means that Eq. (22) must be inverted with regard to mb, which is not straightforward due to the cross product, meaning that you obtain rank 2when inverting the right-hand side, losing information about one of the axes. To that end, an approximation to inverting this equation is given in [25].

mb=Sbbτdbbb2,E23

enabling the currents to be found as

ixiyizΤ=1NAmxmymzΤ.E24

It is here assumed that all three torquers are identical, but depending on satellite configuration, there might be differences in the number of turns and areas. Hence, the desired torque τdbcan be used to find the magnetic moment mbin Eq. (23) and solved for the currents and applied resulting in the actuation torque in Eq. (22). Hence, the control law (Eq. 20) can be mapped to a desired magnetic moment (Eq. 23), which then can be used to find the desired current to each of the three coils. Then, the limits in current will dictate the maximum magnetic moment that can be generated.

Consider the HiNCube satellite as shown in Figure 3. The cubesat comprises three orthonormal magnetic torquers with an area A=0.00757m2and with a maximum current of 47.27 mA and N=100turns. This gives a maximum magnetic moment of mmax=0.03578mA2. Hence, an implementation of using magnetic torquers for attitude control would encompass mapping the output from the control law to a desired magnetic moment using Eq. (23) and then imposing the maximum magnetic moment on each axis, before finding the resulting actuation torque using Eq. (22). Note that to ensure sign correctness due to the projection, the actuation torque can be found as

τx,aτy,aτz,aΤ=signτx,dτx,msignτy,dτy,msignτz,dτz,mΤ,E25

where τab=τx,aτy,aτz,aΤ, τdb=τx,dτy,dτz,dΤ, and τmb=τx,mτy,mτz,mΤ.

To show the performance of magnetic torquers, consider again the HiNCube satellite, which had an inertia matrix of J=1.6710−3Ikg m2. Consider the problem of making rotating 90° from an initial quaternion qo,b=0001Τto qo,d=1000Τ. The gains for the PD+ controller are set kp=1105and kd=5103, and the satellite is assumed to have an orbit of rp=500km and ra=600km, with inclination of 75°. Figure 4 shows the attitude, angular velocity, and actuation torque. It is evident that magnetic torquers produce very low torque, such that it takes a very long time to change the attitude of the spacecraft (about 1h). To some extents, this can be improved by being in a lower orbit where the magnetic field is stronger or by using larger coils with higher currents. Also, note that the actuation signal varies in strength as a function of time, depending on the orbital position of the satellite.

### 5.2. Reaction wheels

Another way of changing the attitude of a satellite is through reaction wheels. Reaction wheels are based on the principle of Newton’s third law: When one body exerts a force on a second body, the second body simultaneously exerts a force equal in magnitude and opposite in direction on the first body. This means that by spinning a reaction wheel in one direction, the satellite will rotate in the other direction. Mounting three reaction wheels in an orthogonal configuration enables three-axis attitude control of spacecraft. From Newton’s third law, the momentum generated by the reaction wheels will have opposite sign of the momentum of the satellite, such that ḣi=ḣwiwhere ḣwiis the momentum production by the reaction wheels, ḣiis the momentum acting on the satellite, and τwb. By employing Euler’s moment equation similarly as in Section 3, the torque generated by the reaction wheels can be found by differentiating hwi=Rbihwbwith hwb=Jwwwbwhere ωwbis the angular velocity of the reaction wheels and Jwdenotes their inertia. This gives

τab=ḣwbSωi,bbhwb,E26

where τwb=ḣwbis the torque generated by the reaction wheels. Now, consider a set of three orthonormal reaction wheels, where one produces torques around the x-axis, one around the y-axis, and one around the z-axis of the body frame. Then, the PD+ control law dictates a desired torque, τdb, which shall be achieved by the reaction wheels. To that end, the torque by the reaction wheel can be rewritten as τwb=τdbSωi,bbhwb, where τwbmust be bounded by the motor torque limit, while the angular momentum will be bounded as a function of maximum rotational speed. After imposing the torque and speed constraints, the angular momentum of the reaction wheels is found by integrating ḣwballowing the actuation torque to be calculated using Eq. (26).

Consider the HiNCube satellite again, where it is possible to use three small reaction wheels as described in [11] where the main idea is to place most of the mass away from the center as shown in Figure 5. The inertia of an individual reaction wheel was found to be Jw=1.46105kg m2, and by assuming a maximum rotation speed of 13,700 rpm with maximum torque of τmax=0.0047Nm, the maximum momentum generated by the reaction wheels is found as hmax=Jwωw=1.46105137002π60=1.52389106. Now, consider the same simulation as when using the magnetic torquers, where the gains for the PD+ controller is changed to kp=kd=2and the reaction wheels has the limits as defined above. Figure 6 shows the simulation results, where it is obvious that by using reaction wheels, the satellite is able to change its orientation after about 80s. To some extent, this can be credited to the higher gains, but it lies mainly with the better actuation system that is able to produce higher torque than the reaction wheels. From the figure, the reaction wheels quickly go into saturation of 13,700RPM, where the angular velocity also goes into saturation. As the quaternion error goes toward zero, the reaction wheel despin, reducing the angular velocity and the control objective, is completed.

### 5.3. Thrusters

The third kind of actuator that will be studied is using reaction control thrusters. This section presents how to map the control signal (Eq. 20) to four thrusters used for attitude control. Let the location of each thruster be denoted by rib=rxryrzΤ, and let them have an azimuth and an elevation angle described by χand γ. Then, the torque produced by a given thruster can be found as [1], p. 262

τib=rib×fib=rysinχcosγrzsinγrzcosγcosχrxcosγsinχrxsinγrycosγcosχfi,E27

where fidenotes the total thrust from the ith thruster. Given the thruster configuration defined in Table 2, let the vector of thruster signals be denoted u=f1f2f3f4Τ, and then the torque can be found as τab=Buwith

B=252525252424242424242424.E28
ThrusterElevation (γ)Azimuth (χ)rxryrz
f14590−0.5−0.45−0.05
f213590−0.5−0.450.05
f3−4590−0.50.45−0.05
f4−13590−0.50.450.05

### Table 2.

Thruster configuration.

Given a desired torque from the PD+ control law, it must be mapped to the desired thruster firings, such that the combination of thrusters produces the desired torque. To that end, there are several different modulation methods that can be applied, ranging from a simple bang-bang modulation to more sophisticated pulse-width pulse-frequency modulation. This section will give an introduction to the different methods and detail how they can be implemented. In general the desired torque can be mapped to the desired thruster firings as ud=Bτdbwhere denotes the Moore-Penrose pseudoinverse and ud=u1u2u3u4Τdenotes the magnitude of each of the thrusters.

1. Bang-bang modulation:The easiest approach to thruster firings is bang-bang modulation, where the thruster is fully actuated as long as the ith signal of udis above zero, such that

fi=fmaxifui>00ifui0,E29

where fmaxdenotes the maximum available thrust from the ith thruster. After applying bang-bang modulation, the vector ucan be constructed allowing the actuator torque to be found as τab=Bu.

2. Bang-bang modulation with dead zone:One of the major drawbacks by using simple bang-bang modulation is when the tracking error has converged to zero, where the thruster firings will continue to maintain the desired attitude. Sensor noise is another source that leads to continuous firings, quickly spending all the propellant. To that end, bang-bang modulation can be augmented with a dead zone, giving

fi=fmaxifui>D0ifuiD,E30

where D>0denotes the dead zone. By properly selecting a suitable dead zone enables the thrusters to avoid firing when close to the equilibrium point.

3. Pulse-width modulation:Another approach that is often used for thruster firings is by using pulse-width modulation (PWM), where an analogue signal (desired torque) can be mapped to discrete signals using PWM. Instead of changing the thrust level, the duration of the pulses can be changed, leading to a pulse that is proportional to the torque command from the PD+ controller. A simple way of achieving this is by using the intersective technique, which uses a sawtooth signal that is compared to the control signal. When the sawtooth is less than the control signal, the PWM signal is in a high state and otherwise in a low state. This makes it possible to go from continuous control signal to a discrete representation which can be used for thruster firings. Figure 7 shows how to achieve the PWM signal, enabled through a simple comparison of the two signals.

4. Pulse-width pulse-frequency modulation:In addition to controlling the width of the pulse as in PWM, it is also possible to control the frequency of the pulse—something that is done through pulse-width pulse-frequency (PWPF) modulation ([1], p. 265) (Figure 8). The modulation approach comprises a lag filter and a Schmitt trigger as shown in Figure 9. As long as the input to the Schmitt trigger is below Uon, the output is kept at zero and must be larger than UonKto produce an output, where Kis a DC gain, τis the time constant, Uonand Uoffare the on and off limits for the Schmitt trigger, while Umis the maximum output. Much research has been performed on improving PWPF modulation, and in [26], the authors propose the following settings (cf. Figure 9): 2<K<6, 0.1<τ<0.5, Uon>0.3, Uoff<0.8Uon, and Um=1.

5. Simulations of thruster modulations:Consider a satellite with an inertia matrix as

J=0.50.20.10.20.50.20.10.20.5,E31

where the objective is to perform a yaw maneuver of 90°using four thrusters with 0.1N force, with a specific impulse of 200s. Figure 10 shows the attitude and angular velocity vectors when using bang-bang modulation, where the satellite is able to make the errors go to zero. However, due to the modulation, the thrusters will continue firing as shown in Figure 11. To that end, consider the bang-bang modulation with dead zone. Let the dead zone be chosen as D=0.05, and then the satellite obtains an accuracy as shown in Figure 12 where there is a small offset from the origin which is proportional to the dead zone. On the other hand, the thruster firings are much less prone to do continuous firings as shown in Figure 13.

Now, consider pulse-width modulation. Let the sawtooth signal have an amplitude of 1and a frequency of 1Hz. Then, the attitude and angular velocity error is obtained as shown in Figure 14, while the thruster firings are shown in Figure 15. It is possible to tune on sawtooth frequency to improve the performance.

The final scenario is using PWPF modulation, where the parameters are chosen as K=3, τ=0.2, Uon=0.35, and Uoff=0.28. Figure 15 shows the attitude and angular velocity, which go close to zero, while the thruster firings are shown in Figure 16, which is able to constrain the amount of thruster firings, and therefore propellant.

For satellite control using thrusters, propellant is a critical resource that must not be wasted. To that end it is desirable to limit the amount of propellant while at the same time obtain good pointing accuracy (Figure 17). With the basis in Tsiolkovsky’s rocket equation, the propellant consumption during thruster firings can be found as mpropellant=0tfIspgdtwhere fis the force from one of the thrusters and Ispis the specific impulse, while g=9.81m/s2is the acceleration due to gravity. Figure 18 shows a comparison between the different modulation methods, where it is evident that the PWPF method allows for the least amount of propellant while obtaining close to acceptable performance. The bang-bang modulation will continue spending propellant until running out of fuel but on the other hand obtains the best tracking performance.

## 6. Attitude determination

As a preliminary step before trying to estimate the attitude of the satellite, some knowledge of measurement vectors must be known, i.e., what is the direction toward the Sun and how does the magnetic field vector look like at a given position. There are several other quantities that can be measured to obtain the attitude, where star trackers are known to be the most accurate. For the reader to obtain insight into using multiple measurements and combine it to find the attitude, this work presents a Sun vector model and a simplified magnetic field model that can be used for simulation purposes.

### 6.1. Sun vector model

To find the direction toward the Sun, there are several models that can be applied. The simplest would be to divide a circle into 365days and have a vector always point toward the Sun. Then, by knowing which day it is, it is straightforward to find the direction toward the Sun. This approach would be coarse, such that more accurate models exist. For example, the Sun vector model in [3], pp. 281–282, has an accuracy of 0.01and is valid until 2050. First, the time and date must be converted into the Julian date as [3], p. 189.

JD=367yrINT7yr+INTmo+9124+INT275mo9+d+1,721,013.5+s60+min60+h24,E32

where a real truncation is denoted by INTand the year, month, day, hour, minute, and second are denoted by yr,mo,d,h,min,s. If the day contains a leap second, 61s should be used instead of 60. This gives the Sun vector model as

TUT1=JD2,451,545.036,525,E33
λM=280.460+36,000.771TUT1,E34
M=357.5277233+35,999.05034TUT1,E35
λecliptic=λM+1.914666471sinM+0.019994643sin2M,E36
ε=23.4392910.0130042TUT1,E37
so=Riocosλeclipticcosεsinλeclipticsinεsinλecliptic.E38

Here, the number of Julian centuries is denoted by TUT1, the mean longitude of the Sun is denoted by λM, the mean anomaly for the Sun is denoted by M, the ecliptic longitude is denoted by λecliptic, and the obliquity of the ecliptic is denoted by ε, while the Sun vector in orbit frame is denoted by so.

### 6.2. Magnetic field model

Several different geomagnetic models can be applied for attitude determination in conjunction with a magnetometer. The most basic are simple dipole models [27], while more advanced are, e.g., the chaos model or the 12th generation IGRF model [28], which is the most commonly used model for attitude determination. This section presents the simple dipole model by [27], which can be described by the magnetic field vector in orbit frame as

mo=μfa3cosω0tsinicosi2sinω0tsiniΤ,E39

where the time measured from passing the ascending node of the orbit relative to the geomagnetic equator is denoted by tand the dipole strength is denoted μf=7.91015Wb-m, while ω0=ωi,oidenotes the angular speed of the orbit. For a real application, the reader is recommended to apply the IGRF model, which is available in Simulink inside the Aerospace Toolbox, as C++ implementation2 or using Python.3

### 6.3. Attitude determination using the Madgwick filter

The objective of attitude determination is to find what direction the satellite is pointing. In its core, it mainly requires two vector measurements and two mathematical models that can be compared and used to find the attitude. There are several different kinds of filters applied for attitude estimation, such as the Triad method [29], the Kalman filter [30], or the Mahony filter [31]. The Madgwick filter by [13] has shown good results in attitude estimation based on IMU measurements and is commonly applied in drone applications. The main idea by the filter is to use gradient descent in combination with the complementary filter to fuse sensor data together to produce an estimated quaternion. This section presents an application of the Madgwick filter by using measurements of the Sun vector and the magnetic field vector as well as the acceleration vector (gravity) and shows how to fuse that data together to estimate the attitude of a satellite in an elliptical orbit.

Let the quaternion estimate be denoted by q̂=q1q2q3q4Τ. The measured acceleration, Sun vector, and magnetic field vectors can be defined, respectively, as ab, sb, and mband can be combined with the mathematical models of the acceleration, Sun vector, and magnetic field vector given in Eqs. (6), (38), and (39) to estimate the attitude. Here, the current estimate is denoted by subscript k, while the previous estimate is denoted by k1. Let the objective function be

F=fq̂k1,aoabfq̂k1,sosbfq̂k1,momb,E40

where the objective is found in an estimated quaternion that minimizes this function, something that can be achieved by using gradient descent. The Jacobian matrix can be found as

J=Jqq̂k1aoJqq̂k1soJqq̂k1moE41

and allows the gradient to be found as f=JΤF. Now, let a vector in the orbit frame obtained from a mathematical model be denoted by zo=oxoyozΤand a vector in the body frame obtained through measurement be denoted by zb=bxbybzΤ. Then, the functions fq̂k1zozband Jqq̂k1zoare given by

fq̂k1,zozb=2ox0.5q32q42+2oyq1q4+q2q3+2ozq2q4q1q3bx2oxq2q3q1q4+2oy0.5q22q42+2ozq1q2+q3q4by2oxq1q3+q2q4+2oyq3q4q1q2+2oz0.5q22q32bz,E42
Jqq̂k1zo=2oyq42ozq32oyq3+2ozq44oxq3+2oyq22ozq14oxq4+2oyq1+2ozq22oxq4+2ozq22oxq34oyq2+2ozq12oxq2+2ozq42oxq14oyq4+2ozq32oxq32oyq22oxq42oyq14ozq22oxq1+2oyq44ozq32oxq2+2oyq3.E43

Given the gyro measurement ωgyrob(relative to inertial frame), the angular velocity relative to orbit frame can be found as

ωo,gyrob=0ωgyrobRobq̂k1Rioωi,oiR4,E44

where the rotation matrix from orbit to body frame is constructed using the estimated quaternion and denoted by Robq̂k1. The Madgwick filter can now be presented as

ω̂kb=2Tq̂k1ff,E45
ωbias,kb=ωbias,k1b+ζω̂kbΔT,E46
ωo,bb=Hωo,gyrobωbias,kb,E47
q̂̇k=12Tq̂k10ωo,bbβff,E48
q̂k=q̂k1+q̂̇kΔT,E49
q̂k=q̂kq̂k,E50

where βand ζare gains, the time step is denoted by ΔT, the estimated angular velocity based on vector measurements is denoted ω̂kbR4, and the estimated gyro bias is denoted by ωbias,kbR4, while the angular velocity of the body frame relative to the orbit frame is denoted ωo,bbR3(expected output) and the estimated quaternion is denoted q̂kdescribing the body frame relative to the orbit frame. Note that the quaternion must be normalized to ensure unit length and that the first elements of ω̂kb,ωbias,kbR4are enforced to zero. To map the angular velocity from four to three elements, the projection matrix is defined as H=0IR3×4, which has a column vector of zeros followed by the identity matrix such that ωo,bbR3.

### 6.4. Simulation

Let a satellite have the inertia matrix as given in Eq. (43), which contains non-diagonal terms which therefore will create perturbing moments due to the gravity. Furthermore, let the satellite have the following initial conditions: qo,b0=0.50.50.50.5Τand ωo,bb=0.10.20.3Τwith q̂o,b0=1000Τand ω̂o,bb=000Τ. The desired quaternion can be chosen as qo,d=1000Τ, while ωo,dd=ω̇o,dd=0. To model noise in the sensor measurements, the quaternion is converted into Euler angles, where noise is added to the different sensors. Then, creating the rotation matrix from the noisy Euler angles allows the Sun vector, acceleration, and magnetometer models in the orbit frame to be rotated to the body frame, where the measurements now contain noise. The step size of the simulation is 0.01, while the sensors are sampled every 0.1s.

The quaternion and angular velocity error of the satellite are shown in Figure 18. After about 50s, the objective of making the attitude error and angular velocity error go to zero is achieved. Since the attitude is not measured directly, the Madgwick filter is used for attitude estimation as shown in Figure 19. Both the quaternion error (estimated truth) and angular velocity error converge close to zero.

From the PD+ controller, the desired torque is mapped to the desired thrust firings using bang-bang modulation (Figure 20). Figure 21 shows the thruster firings required to maintain the attitude error close to zero.

## 7. Conclusion

This chapter has presented all the components required to create an attitude determination and control system for satellites in elliptical orbits. With this as basis, it is the hope by the author that the work can help in developing new results within attitude determination and control systems, ranging from nonlinear controllers to new understanding in orbital mechanics, attitude determination, new sensors, and new actuation methods and strategies.

## References

1. 1. Sidi MJ. Spacecraft Dynamics & Control. New York: Cambridge University Press; 1997
2. 2. Wertz J, editor. Spacecraft Attitude Determination and Control. D. Dordrecht, Holland: Reidel Publishing Company; 1978
3. 3. Vallado DA. Fundamentals of Astrodynamics and Applications. 3rd ed. El Segundo, California: Microcosm Press and Dordrecht, Boston, London: Kluwer Academic Publishers; 2007
4. 4. Markley FL, Crassidis JL. Fundamentals of Spacecraft Attitude Determination and Control. El Segundo, California: Microcosm Press and New York, Heidelberg, Dordrecht London: Springer; 2014
5. 5. Alarcon J, Örger N, Kim S, Soon L, Cho M. Aoba VELOX-IV attitude and orbit control system design for a LEO mission applicable to a future lunar mission. In: Proceedings of the 67th International Astronautical Congress; Guadalajara City, Mexico; 2016
6. 6. Dechao R, Tao S, Lu C, Xiaoqian C, Yong Z. Attitude control system design and on-orbit performance analysis of nano-satellite - tian Tuo 1. Chinese Journal of Aeronautics. 2014;27(3):593-601
7. 7. Nakajima Y, Murakami N, Ohtani T, Nakamura Y, Hirako K, Inoue K. SDS-4 attitude control system: In-flight results of three axis attitude control for small satellites. IFAC Proceedings Volumes. 2013;16(19):283-288
8. 8. Fritz M, Shoer J, Singh L, Henderson T, McGee J, Rose R, Ruf C. Attitude determination and control system design for the CYGNSS microsatellite. In: Proceedings of the IEEE Aerospace Conference; Big Sky, Montana; 2015
9. 9. Oland E. Attitude determination and control system for satellites in elliptical orbits—A complete solution. In: Proceedings of the 8th International Conference on Recent Advances in Space Technology (RAST); Istanbul, Turkey; 2017
10. 10. Oland E, Houge T, Vedal F. Norwegian student satellite program—HiNCube. In: Proceedings of the 22nd Annual AIAA/USU Conference on Small Satellites; Utah, USA; 2008
11. 11. Oland E, Schlanbusch R. Reaction wheel design for cubesats. In: Proceedings of the 4th International Conference on Recent Advances in Space Technologies (RAST); Istanbul, Turkey; 2009
12. 12. Schlanbusch R, Oland E, Nicklasson PJ. Modeling and simulation of a cubesat using nonlinear control in an elliptic orbit. In: Proceedings of the 4th International Conference on Recent Advances in Space Technologies (RAST); Istanbul, Turkey; 2009
13. 13. Madgwick S, Harrison A, Vaidyanathan R. Estimation of IMU and MARG orientation using a gradient descent algorithm. In: Proceedings of the 2011 IEEE International Conference on Rehabilitation Robotics; Zurich, Switzerland; 2011
14. 14. Oland E. Autonomous inspection of the International Space Station. In: Proceedings of the 8th International Conference on Mechanical and Aerospace Engineering; Prague, Czech; 2017
15. 15. Egeland O, Gravdahl JT. Modeling and Simulation for Automatic Control. Trondheim, Norway: Marine Cybernetics; 2002. ISBN: 82-92356-01-0
16. 16. Schaub H, Junkins JL. Analytical Mechanics of Space Systems. AIAA American Institute of Aeronautics & Ast; 2003. pp. 1-19
17. 17. Oland E, Andersen TS, Kristiansen R. Subsumption architecture applied to flight control using composite rotations. Automatica. 2016;69:195-200
18. 18. Kristiansen R. Dynamic synchronization of spacecraft [PhD thesis]. Trondheim, Norway: Norwegian University of Science and Technology; 2008
19. 19. Takegaki M, Arimoto S. A new feedback method for dynamic control of manipulators. ASME Journal of Dynamic Systems, Measurement, and Control. 1981;103(2):119-125
20. 20. Paden B, Panja R. Globally asymptotically stable ‘PD+’ controller for robot manipulators. International Journal of Control. 1988;47(6):1697-1712
21. 21. Kristiansen R, Nicklasson PJ, Gravdahl JT. Spacecraft coordination control in 6DOF: Integrator backstepping vs. passivity-based control. Automatica. 2008;44(11):2896-2901
22. 22. Schlanbusch R, Loría A, Kristiansen R, Nicklasson PJ. PD+ based output feedback attitude control of rigid bodies. IEEE Transactions on Automatic Control. 2012;57(8):2146-2152
23. 23. Oland E. A command-filtered backstepping approach to autonomous inspections using a quadrotor. In: Proceedings of the 24th Mediterranean Conference on Control and Automation; Athens, Greece; 2016
24. 24. Hahn W. Stability of Motion. Berlin/Heidelberg/New York: Springer-Verlag; 1967. ISBN: 978-3-642-50087-9
25. 25. Wiśniewski R. Satellite attitude control using only electromagnetic actuation [PhD thesis]. Aalborg University; 1996
26. 26. Hu Q, Ma G. Flexible spacecraft vibration suppression using PWPF modulated input component command and sliding mode control. Asian Journal of Control. 2007;9(1):20-29
27. 27. Psiaki M. Magnetic torquer attitude control via asymptotic periodic linear quadratic regulation. Journal of Guidance, Control, and Dynamics. 2001;24(2):386-394
28. 28. Thébault E et al. International geomagnetic reference field: The 12th generation. Earth, Planets and Space. 2015;67(79):1-19
29. 29. Shuster MD, Oh SD. Three-axis attitude determination from vector observations. Journal of Guidance and Control. 1981;4(1):70-77
30. 30. Farrell JL. Attitude determination by kalman filtering. Automatica. 1970;6(3):419-430
31. 31. Mahony R, Hamel T, Pflimlin J. Nonlinear complementary filters on the special orthogonal group. IEEE Transactions on Automatic Control. 2008;54(5):1203-1217

## Notes

• Figure created using the MATLAB script “international geomagnetic reference field (IGRF) model” by Drew Compston.
• https://github.com/JDeeth/MagDec
• https://github.com/scivision/pyigrf12

Written By

Espen Oland

Submitted: March 1st, 2018 Reviewed: July 20th, 2018 Published: November 5th, 2018