Attitude Stabilization of a Satellite Having Only Electromagnetic Actuation Using Oscillating Controls

Abstract

We consider the problem of attitude stabilization for a low Earth orbit satellite having only electromagnetic actuation. Such a satellite is not fully actuated, as the control torque is the cross-product of magnetic moment due to magnetorquers and the geomagnetic field. The aim of this work is to study whether oscillating controls can be designed such that a satellite actuated via magnetorquers alone can achieve full three-axis control irrespective of the position of the satellite. To this end, we propose considering oscillating feedback controls which generate the motion of the closed-loop system in the direction of appropriate Lie brackets. Simulation studies show that the proposed control scheme is able to stabilize the considered system.

1. Introduction

In this work, we consider the attitude stabilization problem for a low Earth orbit satellite. The attitude of a satellite is its orientation in the orbital coordinate system and is defined by an attitude matrix which is parameterized by unit quaternions (see [1] for a survey on the application of quaternions for the orientation of rigid bodies). Attitude stabilization can be achieved by applying external torque on the satellite, which is described via Euler’s equations [2]. Besides the torque due to the magnetorquers, the control torque can be provided via gas jets and reaction wheels. The books [3,4] provide an overview of control techniques for the attitude stabilization problem. In this work, we consider a satellite having only electromagnetic actuation. Since such a satellite is actuated via a renewable energy source, it has a longer usable life compared to other actuation techniques. However, it suffers from controllability limitations as the control torque is provided by the cross-product of the geomagnetic field and magnetic moment generated by the satellite’s magnetic coils. Specifically, the control torque is constrained to lie in the plane orthogonal to the geomagnetic field, and therefore, full three-axis control is not possible at all attitudes. However, as the satellite’s position is varying on its orbital plane, correspondingly, the geomagnetic field also varies, and this results in the overall dynamical system being periodic with a frequency depending on the orbital rate of the satellite. Such a nonlinear periodic system is controllable [5,6]. There have been many works on satellites actuated via magnetorquers alone. In [7], full three-axis control is obtained by dividing the dynamics into two loops (inspired by the backstepping method from [8]). A sliding mode controller is used for tracking in the inner loop, and the outer loop is used for stabilization. In [9], the performance of three different control schemes (Detumbling, Linear-Quadratic-Regulator (LQR) and Proportional-Derivative (PD) controls) are compared. The survey paper [10] surveys applications of model-predictive control for such satellites, and the survey paper [11] along with the references therein provide an overview of the present state-of-art in attitude stabilization for satellites having only electromagnetic actuation. To the best of the authors’ knowledge, the current state-of-art in attitude stabilization of satellite relies on the periodic nature of the geomagnetic field to ensure stability given the controllability restriction of a satellite actuated via magnetorquers alone.

In this work, we study whether oscillating controls can be designed such that a satellite actuated via magnetorquers alone can achieve full three-axis control irrespective of the position of the satellite and almost time-invariant geomagnetic field. In order to overcome the aforementioned control theoretic challenges, we consider here oscillating controls in addition to the PD controls discussed in [6]. These oscillating controls are motivated from the literature on the control of underactuated systems [12] and are based on the fact that the considered class of systems admits a time-varying stabilizing feedback law, provided that the local controllability and regularity assumptions as mentioned in [13] are satisfied. The oscillating controls are sinusoidal in nature and generate motion along the direction of Lie brackets of vector fields of the considered dynamical system. The book [14] provides a review on properties of Lie brackets. A generalized design procedure of stabilizing oscillating controls for driftless control-affine systems is presented in [15], and the inclusion of drift is discussed in [16]. The primary aim of this work is to investigate whether oscillating controls can provide sufficient time variation such that full three-axis control can be achieved on a satellite actuated via magnetorquers alone irrespective of the position of the satellite. The control methodology developed subsequently archives full three-axis control in the case of a slowly varying magnetic field. Further work is required to address practical considerations such as magnetorquer saturation, residual magnetic moment and disturbances such as areodynamic drag.

The rest of the paper is organized as follows. In Section 2, we review the mathematical model of the satellite. We rewrite the aforementioned mathematical model in the control-affine form and pose the motivating question in Section 3. Thereafter, we design Lie bracket-based controls in Section 4. Simulation studies are presented in Section 5, and lastly, Section 6 provides the conclusions of this work.

2. Mathematical Model

We begin by stating general coordinate systems (CS) and notations for the satellite.

A pictorial representation of the CS is described in Figure 1 with description of CS given in

Table 1. Coordinate system.

Acronym	Name	Description
BCS	Body (or Control) Coordinate System	CS built on principal axes of inertia of the satellite
OCS	Orbital Coordinate System	Reference CS fixed in the orbit
WCS	World Coordinate System	Inertial right orthogonal CS

. In the sequel, the superscripts ${}^c(·),{}^o(·),{}^w(·)$ represent the vector in parentheses $(·)$ expressed in BCS, OCS and WCS, respectively, and the subscripts ${\mathsf{\Omega }}_{cw},{\mathsf{\Omega }}_{co}$ represent the angular velocity in BCS with respect to WCS and the angular velocity in BCS with respect to OCS.

The equations of motion of a rotating satellite with electromagnetic actuation are governed by kinematics and rigid body dynamics. The kinematics represent the attitude of the rotating satellite and are described via the following equation:

$$\dot{q}=\frac{1}{2}{}^c{\mathsf{\Omega }}_{\mathrm{co}}{q}_4-\frac{1}{2}{}^c{\mathsf{\Omega }}_{\mathrm{co}}\times q,$$

(1)

where q (consisting of three components $q_1$, $q_2$, $q_3$) denotes the vector part of the attitude quaternion, which represents the rotation of the satellite in the BCS with respect to OCS, and $q_4$ is the scalar component of the attitude quaternion. The scalar part $q_4$ is not unique, but is constrained as follows:

$$q_1^2+q_2^2+q_3^2+q_4^2=1.$$

(2)

In this work, we consider only the kinematics of the vector part of the attitude quaternion in (1) and enforce quaternion constraint (2) by considering the space of q as a ball in ${\mathbb{R}}^3$ with $q_4=\pm \sqrt{1-q_1^2-q_2^2-q_3^2}$. The rotation of the satellite in OCS with respect to the BCS is described using the attitude matrix, which is defined as follows:

$$A\left(q\right)=[{}^c{\mathbf{i}}_o,{}^c{\mathbf{j}}_o,{}^c{\mathbf{k}}_o],$$

(3)

where ${}^c{\mathbf{i}}_o,{}^c{\mathbf{j}}_o$ and ${}^c{\mathbf{k}}_o$ are the unit vectors representing the $x,y$ and z axes of the OCS projected on the BCS. The components ${}^c{\mathbf{i}}_o,{}^c{\mathbf{j}}_o$ and ${}^c{\mathbf{k}}_o$ of $A\left(q\right)$ are parameterized by the attitude quaternion q and are defined as follows:

$$\begin{array}{cc}\hfill & {}^c{\mathbf{i}}_o={[q_1^2-q_2^2-q_3^2+q_4^{2}2(q_1{q}_2-q_3{q}_4)2(q_1{q}_3+q_2{q}_4)]}^T,\hfill \\ \hfill & {}^c{\mathbf{j}}_o={[2(q_1{q}_2+q_3{q}_4)q_4^2-q_1^2+q_2^2-q_3^{2}2(q_2{q}_3-q_1{q}_4)]}^T,\hfill \\ \hfill & {}^c{\mathbf{k}}_o={[2(q_1{q}_3-q_2{q}_4)2(q_1{q}_4+q_2{q}_3)q_4^2-q_1^2-q_2^2+q_3^2]}^T.\hfill \end{array}$$

(4)

The dynamics of the satellite represents the time evolution of the angular velocity in OCS. The angular velocity in OCS can be obtained from the angular velocity defined in WCS by

$${}^c{\mathsf{\Omega }}_{\mathrm{co}}={}^c{\mathsf{\Omega }}_{\mathrm{cw}}-{\mathsf{\Omega }}_0{}^c{\mathbf{i}}_o,$$

(5)

where ${\mathsf{\Omega }}_0$ is the orbital rate which will be considered a constant in this work, as the eccentricity of the considered satellite’s orbit is negligible. The Euler equation of the satellite in WCS can be stated as follows:

$$I{}^c{\dot{\mathsf{\Omega }}}_{\mathrm{cw}}=I{}^c{\mathsf{\Omega }}_{\mathrm{cw}}\times {}^c{\mathsf{\Omega }}_{\mathrm{cw}}+N_{ctrl}+N_{gg},$$

(6)

where $N_{ctrl}$ is the control torque, $N_{gg}$ is the torque due to the gravitational gradient, and the remaining terms represent the classical Euler equations as given in [2]. In this work, we have not considered disturbance torque due to the aerodynamic drag and the residual magnetic field as our aim is to study whether oscillating controls can be designed such that a satellite actuated via magnetorquers alone can achieve full three-axis control irrespective of the position of the satellite. Therefore, for the sake of simplicity, we neglect the disturbance torques. Substituting (5) in (6) gives us

$$I({}^c{\dot{\mathsf{\Omega }}}_{\mathrm{co}}+{\mathsf{\Omega }}_0{}^c{\dot{\mathbf{i}}}_o)=I({}^c{\mathsf{\Omega }}_{\mathrm{co}}+{\mathsf{\Omega }}_0{}^c{\mathbf{i}}_o)\times ({}^c{\mathsf{\Omega }}_{\mathrm{co}}+{\mathsf{\Omega }}_0{}^c{\mathbf{i}}_o)+N_{ctrl}+N_{gg}$$

(7)

Simplifying (7) gives us

$$\begin{array}{c}I{}^c{\dot{\mathsf{\Omega }}}_{\mathrm{co}}=(I{}^c{\mathsf{\Omega }}_{\mathrm{co}}\times {}^c{\mathsf{\Omega }}_{\mathrm{co}})+{\mathsf{\Omega }}_0(I{}^c{\mathsf{\Omega }}_{\mathrm{co}}\times {}^c{\mathbf{i}}_o)+{\mathsf{\Omega }}_0(I{}^c{\mathbf{i}}_o\times {}^c{\mathsf{\Omega }}_{\mathrm{co}})+{\mathsf{\Omega }}_0^2(I{}^c{\mathbf{i}}_o\times {}^c{\mathbf{i}}_o)\hfill \\ \hfill +N_{ctrl}+N_{gg}-{\mathsf{\Omega }}_0(I{}^c{\mathsf{\Omega }}_{\mathrm{co}}\times {}^c{\mathbf{i}}_o),\end{array}$$

(8)

The control torque $N_{ctrl}$ in (8) is generated via electromagnetic actuation. Such actuators are known as magnetorquers, and they generate control torque via interaction between magnetic coils stored in the satellites body and the geomagnetic field due to Earth. Mathematically, $N_{ctrl}$ is defined as follows:

$$N_{\mathrm{ctrl}}\left(t\right)=u\left(t\right)\times B\left(t\right),$$

(9)

where $u\left(t\right)$ represents the magnetic moment generated by the satellite, and $B\left(t\right)$ represents the geomagnetic field. The satellite’s magnetic moment is generated via magnetic coils housed inside the satellite body as follows:

$$u\left(t\right)=n_{coil}{i}_{coil}\left(t\right)A_{coil},$$

(10)

where $n_{coil}$ is the number of windings of the coil, $i_{coil}\left(t\right)$ is the current in the coil, and $A_{coil}$ is the area of the coil. The electromagnetic coils are placed perpendicular to the x, y and z axis, and thus, m is a vector with components $m_x,m_y$ and $m_z$. For the rest of the paper, m represents the control signal. For the purpose of simulation studies in this work, the geomagnetic field vector $B\left(t\right)$ is approximated via the dipole model given in [3]. Let $g_1^0$, $g_1^1$, $h_1^1$ represent the empirical Gaussian coefficients of the standard model of Earth’s geomagnetic field. The dipole model in OCS can then be stated as follows:

$$\begin{array}{cc}\hfill & {}^o{B}_x=\frac{a^3{H}_0}{2R^3}\sin {\theta }_m^{\prime }[3\cos (2v^{\prime }-{\alpha }_m)+\cos {\alpha }_m],\hfill \\ \hfill & {}^o{B}_y=\frac{a^3{H}_0}{2R^3}\sin {\theta }_m^{\prime }[3\sin (2v^{\prime }-{\alpha }_m)+\sin {\alpha }_m],\hfill \\ \hfill & {}^o{B}_z=-\frac{a^3{H}_0}{2R^3}\cos {\theta }_m^{\prime },\hfill \end{array}$$

(11)

where a is the radius of Earth, r is the geocentric distance of the satellite, ${\theta }_m^{\prime }$ is the co-elevation of the dipole model, ${\varphi }_m^{\prime }$ is the east longitude of the dipole model, ${\alpha }_m$ is the right ascension of the dipole model, $v^{\prime }$ is the true anomaly measured from the ascending node and $H_0={[g_1^{0^2}+g_1^{1^2}+h_1^{1^2}]}^{1/2}$ is used to calculate the total dipole strength. Furthermore, ${\theta }_m^{\prime }$ and ${\varphi }_m^{\prime }$ are calculated as follows:

$$\begin{array}{cc}\hfill & {\theta }_m^{\prime }=\mathrm{arc}\cos \left(\frac{g_1^0}{H_0}\right),\hfill \\ \hfill & {\varphi }_m^{\prime }=\mathrm{arc}\tan \left(\frac{h_1^1}{g_1^1}\right).\hfill \end{array}$$

(12)

The cosine and sine components of the dipole model have been approximated by a harmonic oscillator. The harmonic oscillator is a 2nd order differential equation, which is described as follows:

$${\ddot{B}}_{\mathrm{mag}}+{\omega }^2{B}_{\mathrm{mag}}=0,$$

(13)

where $\omega$ is frequency of the oscillations, and $B_{mag}$ are the oscillating states representing the cosine and sine components of the dipole model. The solution of (13) approximates the cosine and sine components due to $2v^{\prime }-{\alpha }_m$. Additional harmonic oscillators can be added for more accurate approximations at the cost of adding extra states to the system. The gravity gradient torque is also considered in the model, and it can be stated using the attitude matrix (3) third row k as done in [3], as follows:

$${}^c{N}_{\mathrm{gg}}\left(t\right)=3{\mathsf{\Omega }}_0^2({}^c{\mathbf{k}}_o\times \mathbf{I}{}^c{\mathbf{k}}_o).$$

(14)

In the sequel, $\mathsf{\Omega }$ without subscripts or superscripts represents ${}^c{\mathsf{\Omega }}_{\mathrm{co}}$.

3. Problem Statement

We rewrite system (1)–(8) in the control-affine form as follows:

$$\dot{x}=f_0\left(x\right)+\sum _{i=1}^3{u}_i{f}_i(x,B),x\in X\subset {\mathbb{R}}^6,u\in {\mathbb{R}}^3,$$

(15)

where $x={\left({\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2,{\mathsf{\Omega }}_3,q_1,q_2,q_3\right)}^T$, X is a neighborhood of $x=0$ in ${\mathbb{R}}^6$ such that X does not contain any nontrivial equilibrium of (15) with $u=0$,

$$\begin{gathered} f_0\left(x\right)=\left(\begin{array}{c}\frac{I_2-I_3}{I_1}{\mathsf{\Omega }}_2{\mathsf{\Omega }}_3\\ \frac{I_3-I_1}{I_2}{\mathsf{\Omega }}_1{\mathsf{\Omega }}_3\\ \frac{I_1-I_2}{I_3}{\mathsf{\Omega }}_1{\mathsf{\Omega }}_2\\ \frac{1}{2}({\mathsf{\Omega }}_3{q}_2-{\mathsf{\Omega }}_2{q}_3+{\mathsf{\Omega }}_1{q}_4)\\ \frac{1}{2}({\mathsf{\Omega }}_1{q}_3-{\mathsf{\Omega }}_3{q}_1+{\mathsf{\Omega }}_2{q}_4)\\ \frac{1}{2}({\mathsf{\Omega }}_2{q}_1-{\mathsf{\Omega }}_1{q}_2+{\mathsf{\Omega }}_3{q}_4)\end{array}\right)+\left(\begin{array}{c}3{\mathsf{\Omega }}_0^2{I}^{-1}({}^c{\mathbf{k}}_o\times I{}^c{\mathbf{k}}_o)\\ 0_{3\times 1}\end{array}\right)-\left(\begin{array}{c}{\mathsf{\Omega }}_{0}I{}^c{\mathsf{\Omega }}_{\mathrm{co}}\times {}^c{\mathbf{i}}_o\\ 0_{3\times 1}\end{array}\right)+ \\ \left(\begin{array}{c}{I}^{-1}(I{\mathsf{\Omega }}_0({}^c{\mathsf{\Omega }}_{\mathrm{co}}\times {}^c{\mathbf{i}}_o)+{\mathsf{\Omega }}_0(I{}^c{\mathbf{i}}_o\times {}^c{\mathsf{\Omega }}_{\mathrm{co}})+{\mathsf{\Omega }}_0^2(I{}^c{\mathbf{i}}_o\times {}^c{\mathbf{i}}_o))\\ 0_{3\times 1}\end{array}\right), \end{gathered}$$

$$f_1(x,B)=\left(\begin{array}{c}0\\ \frac{-2(q_1{q}_3+q_2{q}_4)b_1+2(q_2{q}_3-q_1{q}_4)b_2+(q_1^2+q_2^2-q_3^2-q_4^2)b_3}{I_2}\\ \frac{2(q_1{q}_2-q_3{q}_4)b_1+(q_4^2-q_1^2+q_2^2-q_3^2)b_2+2(q_1{q}_4+q_2{q}_3)b_3}{I_3}\\ 0_{3\times 1}\end{array}\right),$$

$$f_2(x,B)=\left(\begin{array}{c}\frac{2(q_1{q}_3+q_2{q}_4)b_1+2(q_2{q}_3-q_1{q}_4)b_2+(q_4^2-q_1^2-q_2^2+q_3^2)b_3}{I_1}\\ 0\\ \frac{(q_2^2+q_3^2-q_1^2-q_4^2)b_1-2(q_1{q}_2+q_3{q}_4)b_2+2(q_2{q}_4-q_1{q}_3)b_3}{I_3}\\ 0_{3\times 1}\end{array}\right),$$

$$f_3(x,B)=\left(\begin{array}{c}\frac{2(q_3{q}_4-q_1{q}_2)b_1+(q_1^2-q_2^2+q_3^2-q_4^2)b_2-2(q_1{q}_4+q_2{q}_3)b_3}{I_1}\\ \frac{(q_1^2-q_2^2-q_3^2+q_4^2)b_1+2(q_1{q}_2+q_3{q}_4)b_2+2(q_1{q}_3-q_2{q}_4)b_3}{I_2}\\ 0\\ 0_{3\times 1}\end{array}\right).$$

For the control design, we consider the following Lyapunov function candidate

$$V\left(x\right)=\frac{1}{2}{\mathsf{\Omega }}^{T}I\mathsf{\Omega }+\frac{3}{2}{\mathsf{\Omega }}_0^2({}^c{\mathbf{k}}_o^{T}I{}^c{\mathbf{k}}_o-I_z)+\frac{1}{2}{\mathsf{\Omega }}_0^2(I_x-{}^c{\mathbf{i}}_o^{T}I{}^c{\mathbf{i}}_o)+2\eta (1-q_4),$$

(16)

where $\eta >0$ is a scalar constant. Then, the directional derivative of $V\left(x\right)$ along the vector field $f_0\left(x\right)$ is

$$L_{f_0}V=\frac{\partial V\left(x\right)}{\partial x}{f}_0\left(x\right)=\eta {\mathsf{\Omega }}^{T}q,$$

(17)

and the time derivative of V along the trajectories of (15) can be written as

$$\dot{V}=\eta {\mathsf{\Omega }}^{T}q+\langle \mathsf{\Omega },u\times \left(A\left(q\right)B\left(t\right)\right)\rangle .$$

(18)

As the scalar triple product $\langle \mathsf{\Omega },u\times \left(A\right(q\left)B\right(t\left)\right)\rangle$ is unchanged under a circular shift of its operands, we have

$$\dot{V}=\alpha \left(x\right)+(u,\beta \left(x\right))$$

(19)

with

$$\alpha \left(x\right)=\eta {\mathsf{\Omega }}^{T}q,$$

$$\beta \left(x\right)=\left(A\right(q\left)B\right(t\left)\right)\times \mathsf{\Omega }.$$

Recall that a positive definite function $V:X\subset {\mathbb{R}}^n\to \mathbb{R}$ of class $C^1$ is said to be a control Lyapunov function (CLF) for a system of the form $\dot{x}=f(x,u)$, $x\in X$, $u\in U\subset {\mathbb{R}}^m$ [17], if, for any $x\in X\backslash \left\{0\right\}$, there exists a vector $u_x\in U$ such that $(\nabla V\left(x\right),f(x,u_x))<0$. Similarly, we call $V\left(x\right)$ a weak CLF if [18], for any $x\in X$, there is a $u_x\in U$ such that $(\nabla V\left(x\right),f(x,u_x))\le 0$. By Artstein’s theorem [19], the existence of a CLF is equivalent to the existence of a stabilizing feedback law of the form $u=v\left(x\right)$, $v\in C(X\backslash \{0\left\}\right)$ for nonlinear control-affine systems. We will see that the design methodology based on Artstein’s theorem (and its extension with a weak CLF [18]) is not directly applicable in our case. Indeed, as one can easily deduce from the expression (19), a necessary and sufficient condition for the function $V\left(x\right)$ to be a weak CLF for (15) is that

$$\mathrm{for}\mathrm{each}x\in X:\mathrm{if}\beta \left(x\right)=0\Rightarrow \alpha \left(x\right)\le 0.$$

(20)

Condition (20) can be equivalently written in the following form: “for each $x\in X$, either $\beta \left(x\right)\ne 0$ or ($\beta \left(x\right)=0$ and $\alpha \left(x\right)\le 0$)”.

We show that condition (20) is violated because the geomagnetic field $B\left(t\right)\ne 0$. Indeed, the condition $\beta \left(x\right)=0$ for $x\ne 0$ is equivalent to

$$\mathsf{\Omega }=\lambda A\left(q\right)B\left(t\right)\mathrm{with}\mathrm{some}\lambda \in \mathbb{R}.$$

(21)

Then, by substituting the above $\mathsf{\Omega }$ into $\alpha \left(x\right)$, we get

$$\alpha \left(x\right)=\lambda \eta \left(A\right(q\left)B\right(t),q),$$

which takes positive values under a suitable choice of $\lambda$ whenever $\left(A\right(q\left)B\right(t),q)\ne 0$. Thus, the only possibility to satisfy the condition (20) for any x and $B\left(t\right)$ at the points $x\ne 0$ is to have $\left(A\right(q\left)B\right(t),q)\equiv 0$. However, we see that the latter property holds only for $B\left(t\right)=0$ due to the structure of $A\left(q\right)$ in (3). Hence, the candidate Lyapunov function $V\left(x\right)$ is not a weak CLF for system (15).

To overcome the above-mentioned constraint for pure Lyapunov-based control design, we will propose a “hybrid” control strategy based on a combination of state feedback controllers with $\epsilon$-periodic oscillating input signals, where $\epsilon >0$ is treated as a small parameter. Such a construction is inspired by the fact that $\dot{V}$ takes positive values at some $x\ne 0$ with any admissible choice of controls u. Thus, the objective of our oscillating control component is to ensure the decreasing of $V\left(x\right(t\left)\right)$ on time intervals of length $\epsilon$ when the time derivative $\dot{V}\left(x\right)$ is not negative.

To be precise, we introduce a parameter $\mu \in (0,1)$ and split the whole state space ${\mathbb{R}}^6$ into the union of two disjoint sets ${\mathcal{A}}_{\mu }\cup {\mathcal{B}}_{\mu }$,

$${\mathcal{A}}_{\mu }=\{x:W_{\mu }\left(x\right)\ge 0\},$$

$${\mathcal{B}}_{\mu }=\{x:W_{\mu }\left(x\right)<0\}.$$

(22)

where

$$W_{\mu }(q,\mathsf{\Omega })=(\mathsf{\Omega },A\left(q\right)B)-\mu \parallel \mathsf{\Omega }\parallel ·\parallel B\parallel .$$

(23)

In the sequel, we will consider the following two cases:

• $x\left(t\right)\in {\mathcal{B}}_{\mu }$;
• $x\left(t\right)\in {\mathcal{A}}_{\mu }$.

Consider first the case $x\left(t\right)\in {\mathcal{B}}_{\mu }$. This implies that $\beta \left(x\right)\ne 0$. We can define the following state feedback law (pre-compensator) aiming to stabilize the trivial solution of system (15)

$$\tilde{u}\left(x\right)=-\gamma \frac{\beta \left(x\right)}{\parallel \beta \left(x\right)\parallel }$$

(24)

where $\gamma >0$ is a design parameter. Note that formula (24) also contains the design parameter $\eta >0$ originating from (16). Then, the time derivative $\dot{V}\left(x\right)$ along the trajectories of the closed-loop system (15), (24) takes the form

$$\dot{V}=\alpha \left(x\right)-\gamma \parallel \beta \left(x\right)\parallel .$$

(25)

The following Lemma proves the existence of a globally stabilizing control law for the case when $x\left(t\right)\in {\mathcal{B}}_{\mu }$. This should be seen as an alternative to the locally stabilizing PD controls discussed in the literature [6,11].

Lemma 1.

For any $\mu \in (0,1)$, there exists a $\gamma >0$ such that the time-derivative of $V\left(x\right)$ (25) along the trajectories of the closed-loop system (15), (24) satisfies

$$\dot{V}\left(x\right)<0\mathrm{for}\mathrm{all}x\in {\mathcal{B}}_{\mu }.$$

Proof. 

We substitute $\tilde{u}\left(x\right)$ from (24) in (18) to obtain

$$\begin{array}{cc}\hfill \dot{V}=\alpha \left(x\right)-\gamma \parallel \beta \left(x\right)\parallel & \le \eta \parallel \mathsf{\Omega }\parallel ·\parallel q\parallel -\gamma \parallel A\left(q\right)\parallel ·\parallel B\parallel ·\parallel \mathsf{\Omega }\parallel ·\mid sin\theta \mid ,\hfill \\ \hfill & \le \eta \parallel \mathsf{\Omega }\parallel -\gamma \parallel B\parallel ·\parallel \mathsf{\Omega }\parallel ·\mid \sqrt{1-{\mu }^2}\mid ,\hfill \\ \hfill & =\parallel \mathsf{\Omega }\parallel ·(\eta -\gamma \parallel B\parallel ·\sqrt{1-{\mu }^2}),\hfill \end{array}$$

(26)

where $\theta$ is the angle between the vectors $A\left(q\right)B$ and $\mathsf{\Omega }$ (recall from the definition of $W_{\mu }$ in (23) that $cos\theta =\frac{(\mathsf{\Omega },A(q\left)B\right)}{\parallel \mathsf{\Omega }\parallel ·\parallel B\parallel }$, and the norm of unit quaternion is 1). By choosing $\gamma$ as

$$\gamma >\frac{\eta }{\parallel B\parallel ·\sqrt{1-{\mu }^2}},$$

(27)

and subsequently $\tilde{u}\left(x\right)$ by substituting (27) into (24), we can ensure negative definiteness of $\dot{V}$ as $x\left(t\right)\in {\mathcal{B}}_{\mu }\Rightarrow \beta \left(x\right)\ne 0$. □

We now consider the case $x\left(t\right)\in {\mathcal{A}}_{\mu }$ and we pose the following question: is it possible to ensure

$$V\left(x\right(\epsilon \left)\right)<V\left(x\right(0\left)\right)$$

for all $x\left(0\right)\in {\mathcal{A}}_{\mu }$ by applying $\epsilon$-periodic oscillating controllers for some $\epsilon >0$? We show in the following section that this is indeed the case if we consider oscillating controls with sine and cosine terms.

4. Lie Bracket-Based Control in ${\mathcal{A}}_{\mu }$

In this section, we construct $\epsilon$-periodic controllers that ensure the decreasing of $V\left(x\right(t\left)\right)$ along the trajectories in ${\mathcal{A}}_{\mu }$. The solutions of system (15) corresponding to the initial data $x\left(0\right)=x^0$ and admissible controls $u_i\left(t\right)$ can be represented by the Chen–Fliess expansion:

$$\begin{array}{cc}\hfill x\left(t\right)=& x^0+\sum _{j_1=0}^3{f}_{j_1}\left(x^0\right){\int }_0^t{u}_{j_1}\left(s\right)ds+\sum _{j_1,j_2=0}^3{L}_{f_{j_2}}{f}_{j_1}\left(x^0\right){\int }_0^t{\int }_0^{s_1}{u}_{j_1}\left(s_1\right)u_{j_2}\left(s_2\right)d{s}_{2}d{s}_1\hfill \\ & +\sum _{j_1,j_2,j_3=0}^3{L}_{f_{j_3}}{L}_{f_{j_2}}{f}_{j_1}\left(x^0\right){\int }_0^t{\int }_0^{s_2}{\int }_0^{s_1}{u}_{j_1}\left(s_1\right)u_{j_2}\left(s_2\right)u_{j_3}\left(s_3\right)d{s}_{3}d{s}_{2}d{s}_1+r\left(t\right),\hfill \end{array}$$

(28)

where $L_{f_i}{f}_j\left(x\right)=\frac{\partial f_j\left(x\right)}{\partial x}{f}_i\left(x\right)$ denotes the directional derivative of $f_j\left(x\right)$ along the vector field $f_i\left(x\right)$, and $r\left(t\right)$ stands for higher-order terms. The Chen–Fliess series are reviewed in [20]. The overall control scheme comprises a time-invariant feedback component with an oscillating component and is similar to the generalized periodic controls presented in [15,16]. Let $\epsilon >0$ be a small parameter which controls the oscillating frequency; then, we define the control functions as follows:

$$\begin{array}{cc}\hfill u_1=& v_1+\frac{2\pi k_{10}{v}_{10}}{\epsilon }sin\left(\frac{2\pi k_{10}t}{\epsilon }\right),\hfill \\ \hfill u_2=& v_2+\frac{2\pi k_{12}\sqrt{2\left|v_{12}\right|}}{\epsilon }sign\left(v_{12}\right)cos\left(\frac{2\pi k_{12}t}{\epsilon }\right),\hfill \\ \hfill u_3=& v_3+\frac{2\pi k_{30}{v}_{30}}{\epsilon }sin\left(\frac{2\pi k_{30}t}{\epsilon }\right),\hfill \end{array}$$

(29)

where $v_1$, $v_2$, $v_3$, $v_{10}$, $v_{30}$, $v_{12}$ are treated as real parameters, and $k_{10},k_{12}$, and $k_{30}$ are nonzero integers such that their magnitudes $|k_{ij}|$ are mutually distinct numbers. We substitute controls (29) into the Chen–Fliess expansion (28) and compute the Taylor series expansion of $x\left(\epsilon \right)$ for small $\epsilon$:

$$\begin{array}{cc}\hfill x\left(\epsilon \right)=& x^0+\epsilon (f_0\left(x^0\right)+v_1{f}_1\left(x^0\right)+v_2{f}_2\left(x^0\right)+v_3{f}_3\left(x^0\right)+[f_1,f_0]\left(x^0\right)v_{10}+[f_3,f_0]\left(x^0\right)v_{30}\hfill \\ & +[f_1,[f_2,f_0]]\left(x^0\right)v_{12})+o\left({\epsilon }^2\right),\hfill \end{array}$$

(30)

where $[f_i,f_j]\left(x\right)=L_{f_i}{f}_j\left(x\right)-L_{f_j}{f}_i\left(x\right)$ represents the Lie brackets of vector fields $f_i\left(x\right)$ and $f_j\left(x\right)$. The idea behind our control scheme is to construct the above $v_1$, $v_2$, $v_3$, $v_{10}$, $v_{30}$, $v_{12}$, depending on the current system state, such that (30) guarantees

$$\left(\nabla V\left(x^0\right),\frac{x\left(\epsilon \right)-x^0}{\epsilon }\right)<0\mathrm{for}\mathrm{sufficiently}\mathrm{small}\epsilon >0\mathrm{and}{x}^0\ne 0,$$

(31)

and thereby ensure $V\left(x\left(\epsilon \right)\right)<V\left(x^0\right)$ along nontrivial trajectories of (15).

If $x^0\in {\mathcal{A}}_{\mu }\setminus \left\{0\right\}$ is chosen in such a way that its $\mathsf{\Omega }$ component is not parallel to $A\left(q\right)B$, then $\beta \left(x^0\right)\ne 0$ because of (21). Thus, by putting $v_{10}=v_{30}=v_{12}=0$ in (29), we see that the condition (31) is satisfied with the control parameters

$$v_i={\tilde{u}}_i\left(x\right),i=1,2,3,$$

(32)

at $x=x^0$, where $\tilde{u}\left(x\right)$ is given by (24) provided that $\gamma >0$ is large enough.

Consider now the case $x^0\in {\mathcal{A}}_{\mu }$ and $\mathsf{\Omega }=\lambda A\left(q\right)B\mathrm{with}\mathrm{some}\lambda \in \mathbb{R}$, where $B\left(t\right)=B$ indicates that the magnetic field in OCS is almost time-invariant (for example, if the satellite is following an equatorial orbit). This implies $\beta \left(x^0\right)=0$, so that $L_{f_1}V\left(x^0\right)=L_{f_2}V\left(x^0\right)=L_{f_3}V\left(x^0\right)=0$, and the property (31) can be ensured by defining

$$\begin{array}{cc}\hfill v_{10}& =-\tilde{\gamma }{L}_{[f_1,f_0]}V\left(x\right),\hfill \\ \hfill v_{30}& =-\tilde{\gamma }{L}_{[f_3,f_0]}V\left(x\right),\hfill \\ \hfill v_{12}& =-\tilde{\gamma }{L}_{[f_1,[f_2,f_0]]}V\left(x\right),\tilde{\gamma }>0,\hfill \end{array}$$

(33)

at $x=x^0$, provided that the gain $\tilde{\gamma }$ is large enough and at least one of the Lie brackets in (33) does not vanish at $x=x^0$. To formalize this scheme, we introduce the following assumption.

Assumption 1.

The vector fields $f_0\left(x\right),f_1(x,B),f_2(x,B),f_3(x,B)$ satisfy the following nonsingularity condition for each $x\in {\mathcal{A}}_{\mu }\setminus \left\{0\right\}$ and $B\in {\mathbb{R}}^3\setminus \left\{0\right\}$:

$$\begin{array}{c}|L_{f_1}V\left(x\right)|+|L_{f_2}V\left(x\right)|+|L_{f_3}V\left(x\right)|+|L_{[f_1,f_0]}V\left(x\right)|+|L_{[f_3,f_0]}V\left(x\right)|\hfill \\ \hfill +|L_{[f_1,[f_2,f_0]]}V\left(x\right)|\ne 0,\end{array}$$

(34)

where we treat B as a parameter in $f_1$, $f_2$, $f_3$, $[f_1,f_0]$, $[f_3,f_0]$, and $[f_1,[f_2,f_0]]$.

From the geometric viewpoint, condition (34) means that the gradient of $V\left(x\right)$ is not perpendicular to the linear span

$$span\{f_1\left(x\right),f_2\left(x\right),f_3\left(x\right),[f_1,f_0]\left(x\right),[f_3,f_0]\left(x\right),[f_1,[f_2,f_0]]\left(x\right)\}$$

for each $x\ne 0$ in ${\mathcal{A}}_{\mu }$ and each admissible value of the geomagnetic field. Then, it is possible to ensure the decreasing of $V\left(x\right(t\left)\right)$ along the trajectories of (15) in ${\mathcal{A}}_{\mu }$ whenever one can implement the motion along the directions $\pm f_1\left(x\right),\pm f_2\left(x\right),\pm f_3\left(x\right),\pm [f_1,f_0]\left(x\right)$,$\pm [f_3,f_0]\left(x\right)$, and $\pm [f_1,[f_2,f_0]]\left(x\right)$.

To sum up, we propose to use the oscillating control strategy (29) with $v_i$ and $v_{ij}$ given by (32) and (33), respectively, when the satellite trajectory is in ${\mathcal{A}}_{\mu }$ and either a globally stabilizing control law (24) or a locally stabilizing PD control law when the satellite trajectory is in ${\mathcal{B}}_{\mu }$.

5. Simulation Results

In this section, we present simulation results on the closed-loop system (15) using controls (29). The geomagnetic field was simulated using (11) with (13). The empirical Gaussian coefficients used in the simulations are obtained from the IGRF model. The inertia (in $kg·m^2$) tensor of the satellite is given as

$$I=\left(\begin{array}{ccc}3.428& 0& 0\\ 0& 2.904& 0\\ 0& 0& 1.280\end{array}\right).$$

In all the simulation results, the frequency parameter $\epsilon$ is $0.001$, and the controller parameters are shown in

Table 2. Controller parameters used in simulations.

	${\mathit{K}}_{{\mathit{P}}_1}$	${\mathit{K}}_{{\mathit{P}}_2}$	${\mathit{K}}_{{\mathit{P}}_3}$	${\mathit{K}}_{{\mathit{D}}_1}$	${\mathit{K}}_{{\mathit{D}}_2}$	${\mathit{K}}_{{\mathit{D}}_3}$	${\mathit{k}}_{10}$	${\mathit{k}}_{12}$	${\mathit{k}}_{30}$
Value	0.0375	0.0375	0.1875	55	60	35	1	2	3

. The form of control law (24) makes it undesirable for practical implementation (as $\parallel \beta \left(x\right)\parallel \to 0$ as $\mathsf{\Omega }\to 0$). We have therefore implemented a locally stable PD controller based on Chapter 4 in [6] when the satellite trajectory is $x\left(t\right)\in {\mathcal{B}}_{\mu }$. The gain parameters for the PD controller were found in the following way. We considered the case when $x\left(t\right)\in {\mathcal{B}}_{\mu }$ (i.e., control scheme (33) is not active) and began by tuning the derivative gains until we stabilized the angular velocities. Thereafter, we tuned the proportional gains such that the kinematics were stabilized.

The altitude of the orbit is 800 km. We ran the simulations with a geomagnetic field that is almost time-invariant as can be seen in Figure 2 and initialized (15) in ${\mathcal{A}}_{\mu }$.

In Figure 3, we can observe the stability of (15) if the satellite was operating in the situation of an almost time-invariant geomagnetic field and (15) were initialized in ${\mathcal{A}}_{\mu }$. During the transient phase (i.e., up to the 4th orbit), we observe high-frequency oscillations up to the 2nd orbit and thereafter, the frequency of oscillations is markedly decreased between the 2nd and the 4th orbit, culminating in stabilization around the 4th orbit. Note that, at the equilibrium point, $q_4$ is equal to $-1$, while $q_1$, $q_2$ and $q_3$ are 0, and this is due to the quaternion constraint (2).

In Figure 4, we observe the angular velocity $\mathsf{\Omega }$ corresponding to Figure 3. The initial angular velocity considered in the simulation can be attained by a damping law (referred to as B-dot law in literature [11]), and attitude maneuvering should be started once the initial velocity is reasonably low. It should be noted that the equilibrium point, $\left\{0\right\}\in {\mathcal{A}}_{\mu }$ as $\beta \left(x\right)\approx 0$ in (19) (since $\mathsf{\Omega }\approx 0$). Consequently, (33) is activated at this point, and therefore, we observe high frequency oscillations at the steady state in Figure 4.

The high-frequency oscillations in the zoomed-in part of Figure 4 are due to the oscillating controls. Figure 5 represents the control torque, and the relatively high torque values in the beginning limit the practical application of this work. In future work, we plan to study these controls with saturation limits included in (29) (imposing saturation limits in an ad hoc manner results in an unstable behavior) to simulate the real-world situation more closely.

We have also simulated the oscillating control scheme with locally stable PD controls for different initial conditions, and the results are shown in Figure 6 and Figure 7.

6. Conclusions and Future Work

In this work, we have proposed an oscillating control scheme for the attitude stabilization of satellite actuated via magnetorquers alone. Simulation results show that the proposed oscillating controls scheme manages to stabilize the attitude and angular velocity in the case of a slowly varying magnetic field. However, the control torque was quite high, and therefore, this technique still needs to be refined. Also, this work does not yet address disturbances due to aerodynamic drag and residual magnetic field. Furthermore, a formal proof is still lacking due to technical challenges in estimating the Chen–Fliess terms of different orders with respect to components of the state vector, particularly in the case of a time-varying geomagnetic field.

Possible future work could be to either find a better candidate Lyapunov function or to use another technique such as combinations of Lyapunov and density functions to prove the stability of the closed-loop system. Another challenging future work is to simulate whether the magnetically actuated satellite with the proposed control scheme can follow an orbit in the geomagnetic equatorial plane. Consequently, $x\left(t\right)\in {\mathcal{A}}_{\mu }$ for all time, and the Lie bracket-based controller (33) will be active always. Unfortunately, due to the limitations of the conventional Ordinary Differential Equations (ODE) solver, we could not simulate this scenario satisfactorily. Specifically, due to the oscillating nature of (29), the ODE solver could not maintain the error tolerances without reducing the step-size to a very low value. Therefore, further investigation is required in this direction, as the small step-size may lead to problems both in the pulse-width modulation implementation and in the residual magnetization of magnetorquers. To summarize, we can conclude that oscillating controls represent a potential research direction in addressing the challenging problem of active three-axis control using magnetorquers only.