Spacecraft Magnetic Control Using Dichotomous Coordinate Descent Algorithm with Box Constraints

In this paper we present magnetic control of a spacecraft using the Dichotomous Coordinate Descent (DCD) algorithm with box constraints. What is common for most work on magnetic spacecraft control is the technique for solving for the control variables of the magnetic torquers where a cross product is included which is well known to be singular. The DCD algorithm provides a new scheme which makes it possible to use a general control law and then adapt it to work for magnetic torquers including restrictions in available magnetic moment, instead of designing a specialized controller for the magnetic control problem. A non-linear passivity-based sliding surface controller is derived for a fully actuated spacecraft and is then implemented for magnetic control by utilizing the previous mentioned algorithm. Results from two simulations are provided, the first comparing the results from the DCD algorithm with older results, and the second showing how easily the derived sliding surface controller may be implemented, improving our results.


Introduction
As mass, space and power consumption are restricted, magnetic torquers are often used as control actuators on board low Earth orbiting spacecraft.Magnetic spacecraft control has a number of implications which makes the control problem different than using other types of actuators.For magnetic torquers it is not possible to provide three independent control torques at each time instant, and the behavior is time dependent because of the variations in the magnetic field along the spacecraft orbit.Therefore other actuators are used to compensate for this such as reaction wheels or a gravity boom.This increases either mass, space and/or power consumption in addition to heightening the complexity of the overall system.Therefore it is desirable to design a control system using only magnetic torquers giving reduced controllability, but capable of fulfilling the mission requirements regarding pointing accuracy of the spacecraft, while saving mass.
A considerable amount of work has been dedicated over the years to solve the nature of magnetic control in the nonlinear case.In (Wiśniewski, 1996) a continuous sliding mode controller was derived and proven asymptotically stable, while backstepping design was utilized in (Wang et al., 1998).Control law for gravity-gradient stabilized spacecraft was derived in (Wiśniewski and Blanke, 1999), and a nonlinear lowgain PD-like control law was proposed in (Lovera and Astolfi, 2001) using only magnetic coils.Proving global uniform asymptotic stability without assuming periodicity of the geomagnetic field was performed by (Gravdahl et al., 2003), and almost global stabilization results were achieved in (Lovera and Astolfi, 2004) for feedback control without rate measurements.This is further pursued in (Reyhanoglu and Drakunov, 2008) doi:10.4173/mic.2010.4.1 c 2010 Norwegian Society of Automatic Control using passivity-based techniques.Attitude tracking control is naturally inspired by literature on tracking control of robot manipulatorscf.(Kelly et al., 2005).The passivity-based approach to robot control have gained much attention, which, contrary to computed torque control, coupes with the robot control problem by exploiting the robots' physical structure (Berghuis and Nijmeijer, 1993).A classic in robot control literature is the PD+ controller of Paden and Panja -cf.(Paden and Panja, 1988) which, together with the Slotine and Li controller - (Slotine and Li, 1987), was the first algorithm for which global asymptotic stability was demonstrated, which has been widely utilized for different control applications such as in (Sira-Ramírez and Siguerdidjane, 1996) for space vehicle stabilization.
In this paper we utilize an iterative technique called the Dichotomous Coordinate Descent (DCD) algorithm (cf.(Zakharov et al., 2008)) which is used to solve for the magnetic moment in the magnetic control equation, which is well known to be singular.We also propose a passivity based sliding surface controller, reminiscent of the classical robot control law (Slotine and Li, 1987) adapted to the topology of SO(3).Simulation results are presented, showing how this method is behaving compared to earlier results.First we utilize an equal controller as derived in (Lovera and Astolfi, 2004), using similar simulation parameters, and compare results.Then the sliding surface controller is implemented using the aforementioned algorithm, giving satisfactory results which are presented through performance functionals to give meaningful comparison.We stress the point that the closed-loop system for the magnetically controlled spacecraft is not uniformly asymptotically stable as for a fully actuated spacecraft.However, in this paper we show how easily an advanced control law can be implemented for magnetic spacecraft control using the presented scheme, leading to increased performance.

Modeling
In the following, we denote by ẋ the time derivative of a vector x, i.e. ẋ = dx/dt, and moreover, ẍ = d 2 x/dt 2 .The cross product operator a × b is denoted S(a)b, ω c b,a is the angular velocity of frame a relative frame b, expressed in frame c, R b a is the rotation matrix from frame a to frame b, and • denotes the 2 -norm.Coordinate reference frames are denoted by F (•) , where the superscript denotes the frame.When the context is sufficiently explicit, we may omit the arguments of a function, vector or matrix.

Cartesian Coordinate Frames
The coordinate reference frames used throughout the paper are defined as follows: Earth-Centered Inertial Frame: The Earthcentered inertial (ECI) frame is denoted F i , and has its origin in the center of the Earth.The axes are denoted x i , y i , and z i , where the z i axis is directed along the axis of rotation of the Earth toward the celestial North Pole, the x i axis is pointing in the direction of the vernal equinox, Υ, which is the vector pointing from the center of the sun toward the center of the Earth during the vernal equinox, and finally the y i axis complete the right handed orthonormale frame.
Earth-Centered, Earth-Fixed Frame: The Earth-Centered, Earth-Fixed (ECEF) frame is denoted F e , and has its origin in the center of the Earth.The axes are denoted x e , y e , and z e , where the z e axis is directed along the axis of rotation of the Earth toward the celestial North Pole, x e intersects the sphere of the Earth at 0 • latitude and longitude and the y e axis complete the right handed orthonormale frame.The ECEF frame is therefore fixed to the earth with an angular rate of rotation of ω e = 7.2921 • 10 −5 rad/s relative to the ECI frame.
North-East-Down Frame: The North-East-Down frame (NED), denoted F n , is defined relative to the Earths' reference ellipsoid and is defined as the tangent plane of the surface of the Earth.The x n axis is pointing toward true north, y n towards true east, and z n points downwards normal to the surface of the Earth.
Spacecraft Orbit Reference Frame: The orbit frame, denoted F o , has its origin located in the center of mass of the spacecraft.The e r axis in the frame coincide with the vector r ∈ R 3 from the center of the Earth to the spacecraft center of mass, and the e h axis is parallel to the orbital angular momentum vector, pointing in the orbit normal direction.The e θ axis completes the right-handed orthonormale frame.The basis vectors of the frame can be defined as Spacecraft Body Reference Frame: The body frame of the spacecraft is denoted F b , and is located at the center of mass of the spacecraft, and its basis vectors are aligned with the principle axis of inertia.

Quaternions
The attitude of a rigid body is represented by a rotation matrix R ∈ SO(3) fulfilling which is the special orthogonal group of order three.Quaternions, often referred to as Euler parameters, are used to parameterize members of SO(3), where the unit quaternion is defined as q = [η, ] ∈ S 3 = {x ∈ R 4 : x x = 1}, where η ∈ R is the scalar part and ∈ R 3 is the vector part.The rotation matrix may be described by (Egeland and Gravdahl, 2002) (1) The inverse rotation can be performed by using the inverse conjugated of q as q = [η, − ] .The set S 3 forms a group with quaternion multiplication, which is distributive and associative, but not commutative, and the quaternion product of two arbitrary quaternions q 1 and q 2 is defined as (Egeland and Gravdahl, 2002) The rotation from F n to F e is expressed as (Fossen, 2002) where l is the longitude while µ is the latitude.F e rotates with respect to F i at a rate ω e about the z i axis.This rotation may be written as where t is the time since F e and F i were aligned.

Kinematics
The time derivative of eq. ( 1) can be written as (Egeland and Gravdahl, 2002) where ω a a,b ∈ R 3 is the angular velocity of frame F b relative to frame F a .The kinematic differential equations may be expressed as (Egeland and Gravdahl, 2002) where

Dynamics
The dynamical model of the satellite can be described by a differential equation for angular velocity, and is deduced from Euler's moment equation.This equation describes the relationship between applied torque and angular momentum in a rigid body as (Sidi, 1997) where τ ∈ R 3 is the total torque working on the body frame, and ḣi indicates the derivative in the inertial frame, while ḣb indicates the derivative in the rotating body frame, and J ∈ R 3×3 is the moment of inertia.
The torque working on the body is derived from two parameters, where τ b d is the disturbance torque, and actuator (control) torque τ b a , such as The dynamical model is derived from eqs. ( 2) and (3) expressed as where ω b i,b is the angular velocity of the body frame relative to the inertial reference system, and the angular velocity of the body frame relative to the orbit frame, ω b o,b is expressed as where

Magnetic Field
The Earth is surrounded by a magnetic field which may be used for spacecraft control purposes.Although this approach has its limitations and challenges, it is a simple and cheap way of achieving adequate stability performance where high pointing accuracy is not required.
The magnetic field of Earth can for simplicity be seen as a perfect dipole (Psiaki, 2001) or a more advanced model such as (Wertz, 1978) [g n,m cos (mφ) + h n,m sin (mφ)] P n,m (θ) where g n,m and h n,m are Gaussian coefficients, P n,m (θ) is the Gauss function of colatitude only, θ and φ are geographic colatitude and longitude, respectively, a is Earth radius, and r is the Earth orbit radius.B r is the outward positive radial component of the field, B θ is the south positive coelevation component, and B φ is the east positive azimuthal component.The field may be expressed in the NED coordinates such as where = µ − δ < 0.2 • , µ d is the geodetic latitude, and δ = 90 • − θ is the declination.To express the magnetic field in the body frame we apply (Fossen, 2002)

Magnetic Torquers
Magnetic torquers have been used for attitude control of spacecraft since the mid-sixties.A magnetic coil produces a magnetic dipole when current is applied, and by the influence of magnetic field of the Earth, a torque is produced which may be expressed as (Sidi, 1997) where m b = [m x , m y , m z ] is the magnetic dipole moment produced by each coil, and where N is the number of coil windings, i c is the current flowing through the coil, and A is the coil area.

Dichotomous Coordinate Descent (DCD) Algorithm
The DCD Algorithm (Zakharov and Tozer, 2004) is one among many iterative techniques for solving the linear Least-Square (LS) problem.According to (Golub and Van Loan, 1996) such techniques typically requires O(N 2 ) to O(N 3 ) operations per sample, including multiplication and division.As processing power on board a spacecraft is highly restricted we want to keep number of processes to a minimum.Another interesting property of this algorithm is called the box constraint which helps limiting the solution range within given bounds.This is particulary interesting for spacecraft control because of the limitations in available torque.
The objective is to solve the equation where Z ∈ R N ×N , h, d ∈ R N .We start by multiplying (7) by Z T , thus obtaining where A = Z Z ∈ R N ×N , and β = Z d ∈ R N .We define a cost function for the LS problem as where h i is the elements of vector h, and H > 0 is the bounded constraint of our solution vector h.The algorithm is presented in Table 1 (Zakharov et al., 2008), where N is the number of unknowns, M b is the number of bits within the amplitude range [−H, H], and R (p) is the p-th column of the appropriate matrix.The algorithm starts an iterative approximation of the solution vector h from the most significant bit.Once the most significant bit has been found for all vector elements, the algorithm starts updating the next less significant bit, and so on.If a bit update happens, we call it a successful update and β is also updated.To limit the complexity we use a predefined number for successful iterations N u .If this number is large enough, the accuracy of the solution is 2 −M b H.As the algorithm only requires P ≤ N (2N u + M b − 1) + N u number of additions it is well suited for solving the magnetic control problem in real-time.

Controller Design
For control of the spacecraft we incorporate a passivity based sliding surface controller, similar to the one if q > N u , the algorithm stops 12 if flag=1, then repeat step 2 derived in (Slotine and Li, 1987), adapted for quaternion feedback attitude control.For the control law it is assumed that the spacecraft has available information of its orbit position r, orbit velocity v, attitude q, and angular velocity ω b i,b .Then the implementation of the control law is shown by use of the DCD algorithm.

Problem Formulation
The control problem is to design a controller that makes the state q(t) converge towards the generated reference specified as q d satisfying the kinematic equation qd = T(q d )ω d , and acts as a solution to the dynamic model presented in (4).The error quaternion q = [η, ˜ ] is found by using the quaternion product q = q ⊗ qd = 1 2 and the error dynamic can according to (Fjellstad, 1994) be expressed as Due to the redundancy in the quaternion representation, q and −q represent the same physical attitude but mathematically it differs by a 2π rotation about an arbitrary axis.Therefore we are not able to achieve global representation since the term global refers to the whole state space R n according to (Hahn, 1967).Since both equilibrium points represent the same physical representation we choose the equilibrium point which require the shortest rotation, minimizing the path length.Hence q+ = [1, 0 ] is chosen if η(t 0 ) ≥ 0, and q− = [−1, 0 ] is chosen if η(t 0 ) < 0. We then apply a coordinate transformation such that the stable equilibrium point is located in the origin.For the positive equilibrium point, the attitude error vector is chosen as (-cf.(Kristiansen, 2008)) ẽq+ = [1 − η, ˜ ] , while for the negative equilibrium point, the error vector is chosen as ẽq− = [1 + η, ˜ ] .The angular velocity error vector is chosen as e ω = ω − ω d .

Control Law
In the following e q = e q+ , and we define a reference trajectory as where the desired trajectory for pointing is defined as q d (t) and ω d (t) with the relationship using equation ( 9), Γ = Γ > 0 is a feedback gain matrix, and using the rotation error expressed in (10).We then define a sliding surface by applying eq. ( 11) which leads to where ω = ω b i,b .By using a control law expressed as where K p = K p > 0 and K d = K d > 0 are feedback gain matrixes, we obtain the closed-loop system by inserting ( 12) into (4) A radial unbounded, positive definite Lyapunov function candidate is defined as and by differentiating eq. ( 14) results in V = s Jṡ + e q K p ėq , and by inserting (13) we end up with where the first term in eq. ( 15) is zero because S (Jω) is a skew-symmetric matrix, which leads to By employing Lyapunov arguments (cf.(Khalil, 2002)), we find that the closed-loop system given by eq. ( 13) is uniformly asymptotically stable (UAS) in the equilibrium point (e q , e ω ) = (0, 0).The proof for the negative equilibrium point e q− is performed in the same way leading to a similar result.Hence, it follows that the dual equilibrium points (e q± , e ω ) = (0, 0) are UAS, thus we do not obtain global results.

Implementation
The torques produced by the magnetic torquers do not act about their mounted axes but vary according to the orientation of the magnetic field (5).Since the control variable for the actuators is the current flowing through the coils (6), we need to solve (5) for m b .Since the matrix S(•) is skew-symmetric, the equation cannot be solved analytically.In previous work this has typically been solved by introducing τ a = Γ(t)τ c where τ a is the actuator torque used in the dynamical model, τ c is the control law, and Γ(t) = S(B)S (B)/ B 2 .Instead we now rewrite (5) as and by applying the DCD algorithm on ( 16) we observe that A(t) = S (B)S(B) is a symmetric time varying matrix with eigenvalues λ 1 = 0, and λ 2 = λ 3 = B 2 x + B 2 y +B 2 z , which according to (Horn and Johnson, 1985) leads to positive semi-definiteness.Since this is true for all t the matrix is uniformly positive semi-definite.It is no surprise that λ 1 = 0 because rank(A) < 3. On the right hand side of (16) β = −S (B)τ , and we end up on the form (8) where m = h.
Again we stress that even if the closed-loop system is proven UAS for the fully actuated spacecraft, it does not guarantee UAS when the control law is implemented by using the DCD algorithm presented in this section.

Simulation Results
In the following, simulation results are presented to illustrate the performance of the magnetic control scheme using the DCD algorithm.In all simulations we have been using similar initial values and parameters as in (Lovera and Astolfi, 2004) for comparative reasons.The simulations were performed in Simulink using a variable sample-time Runge-Kutta ODE45 solver, with tolerance of 1•10 −6 .The moments of inertia were given as J = diag{27, 17, 25} kgm 2 , and the spacecraft was chosen to operate in a near polar orbit with inclination at 87 • in a circular earth orbit with an altitude of 450 km, and the argument of perigee and the right ascension of the ascending node at 0 • .The maximum magnetic moment from the magnetic torquers was chosen as 8 Am 2 , and gravity gradient perturbation was introduced according to (Schaub and Junkins, 2003), but not accounted for in the controllers.
To evaluate and give meaningful comparison of the performance of the controllers, we use the performance functionals where t 0 and t f define the start and end of the simulation window, respectively.The functionals J q and J eq describe the integral functional error of the attitude between body and desired frame, and body and estimated frame, respectively, while J p describes the integral of the applied control torque.
Figure 1 shows the results of a simulation performed using a PD controller similar to the one in (Lovera and Astolfi, 2004) except that the saturation is accounted for in the DCD algorithm instead of in the control law.The controller gains were chosen as K p = diag{1 × 10 −3 I} and K d = diag{5 × 10 −2 I}.As the results show, the spacecraft is faster stabilized, in less than two orbit periods, compared to about three periods in (Lovera and Astolfi, 2004), but at a higher cost of control dipole moments.Now we utilize a more advanced controller such as the one derived in Section 4.2, using equal gains and magnetic dipole moment constraint, and Γ = I, we get the results as pictured in Figure 2. It should be noted that ωd = ω d = 0 for sake of comparison since a tracking controller can be expected to perform better than a regulator.The spacecraft is stabilized in less than one orbit period, which is a significant improvement compared the PD controller, which is evident by looking at the performance functionals in Table 2.Note that even as the sliding surface controller is faster, it also consumes less energy compared to the PD regulator.

Conclusions
In this paper we have applied the Dichotomous Coordinate Descent (DCD) algorithm with box constraints for magnetic control of a spacecraft.Using this scheme we do not have to design a controller especially to work for magnetic torquers but can use any stable controller  and then adapt it to work for magnetic control.Simulation results were presented comparing our results with earlier results on magnetic control, and it was shown how different controllers easily can be implemented, providing satisfactory results of convergence.As the algorithm is used for finding the solution of ( 16) the solution doesn't necessarily provide the correct answer of (5).This is because going from ( 16) to ( 5) has resemblance to zero division, but even so, good results are acquired, especially for controller torques not close to zero.
= e h × e r and e h := h h , where h = r× ṙ is the angular momentum vector of the orbit, h = |h| and r = |r|.This frame is also known as the Local Vertical/Local Horizontal (LVLH) frame.

Figure 1 :
Figure1: Attitude quaternion, angular velocity and control dipole moments for attitude acquisition using a PD controller with saturation and disturbance torque.

Figure 2 :
Figure 2: Attitude quaternion, angular velocity and control dipole moments for attitude acquisition using a passivity-based sliding surface controller with saturation and disturbance torque.129

Table 2 :
Values of performance functionals for attitude