Adaptive Terminal Sliding Mode Control of Electromagnetic Spacecraft Formation Flying in Near-Earth Orbits

Abstract

An adaptive terminal sliding mode control for six-degree-of-freedom electromagnetic spacecraft formation flying (EMFF) in near-Earth orbits is presented. By using terminal sliding mode (TSM) technique, the output tracking error can converge to zero in finite time, and strong robustness with respect to disturbance forces can be guaranteed. Based on a rotated frame F ^r and the adaptive TSM controller, the special magnetic moment of the steerable magnetic dipole is computed. The angular momentum management strategy (AMM) is implemented in a periodically switching fashion, by which the angular momentum buildup was limited. Illustrative simulations of EMFF are conducted to verify the effectiveness of the proposed controller.

1. Introduction

Spacecraft formation flying (SFF) represents the concept of distributing the functionality of large spacecraft among smaller, less-expensive, cooperative spacecraft [1, 2]. Specifically, NASA and the U.S. Air Force have identified spacecraft formation flying as an enabling technology for future missions. The practical implementation of the SFF concept relies on the accurate control of the relative positions and orientations between the participating spacecraft for formation configuration. The conventional thruster-based schemes may require continuous expenditure of fuel to maintain formation geometry that can contaminate the sensitive sensors on board and mission lifetime also becomes dependent on the fuel available [3, 4].

To alleviate these concerns, several propellant-free formation flying methods have been proposed in the literature. The propulsive conducting tethers and spin-stabilized tether systems have been proposed in place of on-board propulsion systems to form and maintain satellite formations [5, 6]. King et al. [7] have presented Coulomb force approaches to maintain a formation. The flux pinning technology has been applied to achieve passively stable configurations by HTS electromagnetics [8]. LaPointe [9] has presented the microwave scattering formation flight method. Radiation forces on the order of 10-9 N/W may be generated using electromagnetic gradient forces or scattering forces; microwave beam powers of 10 kW can thus produce restoring forces of approximately 10 μN, which are sufficient to correct a number of orbital perturbations. Miller et al. [10, 11] address the novel concept of electromagnetic formation flying (EMFF) in which high temperature superconducting (HTS) wire technology is used to create magnetic dipoles on each spacecraft that can be used to maintain and reconfigure the spacecraft formation. However, Since magnetic force on each spacecraft in the formation can be applied in any arbitrary direction which can be easily created by steerable magnetic dipoles, EMFF has advantages in terms of controllability [4].

A critical component of EMFF is an effective formation flying control. The dynamics and control problem associated with EMFF become highly challenging, due to the nonlinear nature of the magnetic forces. Ahsun and Miller [3] have presented a hybrid adaptive control scheme in which translation control is implemented in a centralized fashion with a decentralized attitude control. Elias et al. [12] designed a linear optimal controller based on the linearized dynamics. Kong et al. [13] addressed the use of electromagnetic dipoles for relative position and orientation maintenance as needed for the terrestrial planet finder (TPF). Reference [14] derived the dynamics of an N-spacecraft EMFF (in 2D) for deep space missions and discussed a nonlinear control law using potential functions.

The orientation of a dipole obviously depends on the orientation of the body axes in the inertial space, and changing the dipole on one spacecraft affects actuation on all other spacecraft in EMFF. Therefore, general asymptotical stability may not deliver fast enough convergence to meet EMFF control for high-precision situation. The recently developed terminal sliding mode (TSM) control enables convergence to the desired state in finite time [15–18]. This technique has been used successfully in some control designs, such as robotic manipulators [17] and mobile target tracking [18]. The physical interpretation of finite time convergence lies in the fact that the convergence rate of TSM grows exponentially when the state is near equilibrium.

In this Paper, we confine our attention to adaptive TSMC design for EMFF in near-Earth orbits. The remainder of this paper is organized as follows. Section 2 presents a detailed 6-DOF dynamic modeling for EMFF in near-Earth orbits. In Section 3, a position/attitude tracking controller based on TSM is developed for EMFF, which ensures error convergence in finite time and strong robustness of the bounded disturbances. Based on a rotated frame F ^r and the adaptive TSM controller, the special magnetic moments of the Steerable Magnetic Dipole (SMD) are computed. The AMM strategy is implemented in a periodically switching fashion. Simulation results are presented in Section 4. Finally, Section 5 concludes this paper.

2. System Model

Spacecraft formation constitutes an N-body mechanics problem. Figure 1 shows a typical formation system including N spacecrafts orbiting a central body (Earth). The spacecraft 0's designated as the “leader” spacecraft. We make the following considerations. (1) The ECI frame F ^I is defined with its origin at the center of Earth, its X-axis points toward vernal equinox, Z-axis points toward celestial north pole, and Y-axis completes a right-handed axis system. (2) The Hill frame F ^H is used to visualize the relative motion of each follower spacecraft with respect to the leader spacecraft. The x-axis aligned in the radial (zenith) direction, the z-axis is perpendicular to the orbital plane and points in the direction of the angular momentum vector, and the y-axis completes the right-hand system. (3) The orbital frame F ^Oi is used in the attitude control. The reference frame is rotating about the y-axis with respect to ECI frame at orbital rate. The roll axis x is in the fight direction, the pitch axis y is perpendicular to the orbital plane, and the yaw axis z points toward the Earth. (4) A body frame F ^Bi attached to the body of ith spacecraft with the center of mass of the spacecraft is also defined to describe the orientation of each spacecraft in the inertial space. (5) R _i ∈ ℝ³ denotes the position vector from the origin of ECI frame to ith spacecraft. r _ij ∈ ℝ³ denotes the relative position vector from ith spacecraft to the follower jth spacecraft. ${\mathit{\rho }}_i\triangleq {\left[\begin{array}{ccc}{x}_i& y_i& z_i\end{array}\right]}^T\in {ℝ}^{\mathrm{3}}$ is the expression of the relative position vector r_0i from the origin of the leader spacecraft coordinates system to the follower spacecraft i in the Hill frame F ^H . (6) The rotated frame F ^ri is used in a simplified algebraic form of the magnetic force equation. The x-axis is aligned with vector rij [3].

Figure 1:

Schematic representation of the EMFF system.

2.1. Translational Dynamics

Designating spacecraft 0 as the leader of the formation, the relative dynamics of ith spacecraft in the Hill reference frame F ^H can be represented as

${\ddot{\mathit{\rho }}}_i+{\mathbf{C}}_{iT}\left({\omega }_O,{\dot{\mathit{\rho }}}_i\right)+{\mathbf{N}}_{iT}\left({\mathit{\rho }}_i,{\omega }_O,{\mathbf{R}}_{\mathrm{0}}\right)=\widehat{{\mathbf{F}}_{di}}+\widehat{{\mathbf{F}}_{ci}},$

(1)

where ω _O is the orbital angular velocity of the leader. C _iT (·)∈ℝ³ is a nonlinear term defined as

${\mathbf{C}}_{iT}\left({\omega }_O,{\dot{\mathit{\rho }}}_i\right)\triangleq \mathrm{2}{\omega }_O{\left[\begin{array}{ccc}-{\dot{y}}_i& {\dot{x}}_i& \mathrm{0}\end{array}\right]}^T.$

(2)

N _iT (·)∈ℝ³ is nonlinear term defined as

${\mathbf{N}}_{iT}({\mathit{\rho }}_i,{\omega }_O,{\mathbf{R}}_{\mathrm{0}})\triangleq \left[\begin{array}{c}\mu \left(\frac{x_i+\parallel {\mathbf{R}}_{\mathrm{0}}\parallel }{{\parallel {\mathbf{R}}_{\mathrm{0}}+{\mathit{\rho }}_i\parallel }^{\mathrm{3}}}-\frac{\mathrm{1}}{{\parallel {\mathbf{R}}_{\mathrm{0}}\parallel }^{\mathrm{2}}}\right)-{\omega }_O^{\mathrm{2}}{x}_i-{\dot{\omega }}_O{y}_i\\ \mu \frac{y_i}{{\parallel {\mathbf{R}}_{\mathrm{0}}+{\mathit{\rho }}_i\parallel }^{\mathrm{3}}}-{\omega }_O^{\mathrm{2}}{y}_i+{\dot{\omega }}_O{x}_i\\ \mu \frac{z_i}{{\parallel {\mathbf{R}}_{\mathrm{0}}+{\mathit{\rho }}_i\parallel }^{\mathrm{3}}}\end{array}\right],$

(3)

where μ is gravitational constant of earth.

Note that $\widehat{{\mathbf{F}}_{di}}$ denotes the parametric uncertainty and external disturbances, such as differential J ₂ and higher-order terms, differential solar pressure, and differential drag. The specific relative magnetic force $\widehat{{\mathbf{F}}_{ci}}$ provided by the coils needs to be strong enough to cancel these disturbances and at the same time provide additional force for trajectory following. Consider the following:

$\widehat{{\mathbf{F}}_{ci}}=\frac{\mathrm{1}}{m_i}{\mathbf{F}}_{cmi}-\frac{\mathrm{1}}{m_{\mathrm{0}}}{\mathbf{F}}_{cm\mathrm{0}},$

(4)

where

${\mathbf{F}}_{cmi}=\sum _{j=\mathrm{0},\hspace{0.17em}j\ne i}^{N-\mathrm{1}}{\mathbf{F}}_{ij}^m.$

(5)

By approximating the coils on each spacecraft to SMD, the magnetic field due to the jth spacecraft can be written as [3, 10]

${\mathbf{B}}_j\left({\mathbf{r}}_{ij}\right)=\frac{{\mu }_{\mathrm{0}}}{\mathrm{4}\pi }\left(\frac{\mathrm{3}{\mathit{\mu }}_j·{\mathbf{r}}_{ij}}{r_{ij}^{\mathrm{5}}}{\mathbf{r}}_{ij}-\frac{{\mathit{\mu }}_j}{r_{ij}^{\mathrm{3}}}\right)$

(6)

and the force between spacecrafts can be written as

$\begin{array}{l}{\mathbf{F}}_{ij}^m\triangleq f\left({\mathit{\mu }}_i,{\mathit{\mu }}_j,{\mathbf{r}}_{ij}\right)\\ =\frac{\mathrm{3}{\mu }_{\mathrm{0}}}{\mathrm{4}\pi }(\phantom{\frac{\left({\mathit{\mu }}_i·{\mathbf{r}}_{ij}\right)\left({\mathit{\mu }}_j·{\mathbf{r}}_{ij}\right)}{{\mathbf{r}}_{ij}^{\mathrm{7}}}}-\frac{{\mathit{\mu }}_i·{\mathit{\mu }}_j}{r_{ij}^{\mathrm{5}}}{\mathbf{r}}_{ij}-\frac{{\mathit{\mu }}_i·{\mathbf{r}}_{ij}}{r_{ij}^{\mathrm{5}}}{\mathit{\mu }}_j-\frac{{\mathit{\mu }}_i·{\mathbf{r}}_{ij}}{r_{ij}^{\mathrm{5}}}{\mathit{\mu }}_i+\mathrm{5}\frac{\left({\mathit{\mu }}_i·{\mathbf{r}}_{ij}\right)\left({\mathit{\mu }}_j·{\mathbf{r}}_{ij}\right)}{r_{ij}^{\mathrm{7}}}{\mathbf{r}}_{ij}),\end{array}$

(7)

where μ₀ = 4π × 10⁻⁷ T·m/A is the permeability constant. μ _i and μ _j are the dipole strength of the ith satellite and the jth satellite, respectively. Note that (7) gives the force on dipole i (located on ith spacecraft) due to dipole j (located on jth spacecraft), which depends on the distance between the two dipoles and the orientation of both dipoles in the inertial space. It is the dependence on the orientation of the dipoles that gives rise to the complexity of the expression for the force, since the orientation of a dipole obviously depends on the orientation of the body axes in the inertial space [3].

2.2. Attitude Dynamics

To avoid singular points, the Euler parameter is chosen to describe the attitude of the spacecraft. Let q i represent the Euler parameters corresponding to the attitude of the body frame relative to the orbital frame F ^Oi [19],

$\begin{array}{c}{\mathbf{q}}_i={\left[\begin{array}{cc}{\mathbf{\epsilon }}_i^T& {\eta }_i\end{array}\right]}^T\in {ℝ}^{\mathrm{4}},\\ {\mathbf{\epsilon }}_i\equiv {\left[\begin{array}{ccc}{\epsilon }_{i\mathrm{1}}& {\epsilon }_{i\mathrm{2}}& {\epsilon }_{i\mathrm{3}}\end{array}\right]}^T\in {ℝ}^{\mathrm{3}\times \mathrm{1}},\hspace{1em}{\eta }_i\in {ℝ}^{\mathrm{1}\times \mathrm{1}}.\end{array}$

(8)

The Euler parameters, which are equivalent to the coefficients of unit quaternion, have unit norm by definition; hence

${\mathit{\epsilon }}_i^T{\mathit{\epsilon }}_i+{\eta }_i^{\mathrm{2}}=\mathrm{1}.$

(9)

Attitude kinematics and dynamics of the spacecraft are governed by [19]

${\mathbf{\Omega }}_i=\mathrm{2}{{\mathbf{T}}_i}^{-\mathrm{1}}{\dot{\mathbf{\epsilon }}}_i+{\mathbf{\omega }}_{Oi},$

(10)

${\mathbf{J}}_i{\dot{\mathbf{\Omega }}}_i+{\mathbf{\Omega }}_i^{\times }\left({\mathbf{J}}_i{\mathbf{\Omega }}_i+\mathbf{h}\right)={\mathit{\tau }}_{mi}+{\mathit{\tau }}_{cwi}+{\mathit{\tau }}_{di},$

(11)

where ${\mathbf{T}}_i\equiv {\eta }_i\mathbf{I}+{\mathbf{\epsilon }}_i^{\times }\in {ℝ}^{\mathrm{3}\times \mathrm{3}}$ , the superscript “×” denotes skew-symmetric matrix of a vector, and Ω _i = ω _Oi + ω _bi , ωbi is angular velocity of the spacecraft relative to the orbital frame F ^Oi . ωOi is the orbital angular velocity of ith spacecraft. Ji is the moment of inertia of the spacecraft. h is the angular momentum of the flywheel. The magnetic torque term in (11) can be written as

${\mathit{\tau }}_{mi}=\sum _{j=\mathrm{0},\hspace{0.17em}i\ne j}^{N-\mathrm{1}}{\mathit{\tau }}_{ij}^m,$

(12)

where ${\mathit{\tau }}_{ij}^m$ is the torque exerted on the SMD of ith spacecraft due to that of jth spacecraft [3, 10],

${\mathit{\tau }}_{ij}^m=\frac{{\mu }_{\mathrm{0}}{\mathit{\mu }}_i}{\mathrm{4}\pi }\times \left(\frac{\mathrm{3}{\mathbf{r}}_{ij}({\mathit{\mu }}_j·{\mathbf{r}}_{ij})}{r_{ij}^{\mathrm{5}}}-\frac{{\mathit{\mu }}_j}{r_{ij}^{\mathrm{3}}}\right)$

(13)

τ _cwi and τ _di are the control and disturbance torques, respectively. τ _di includes the gravity gradient torques τ _gi and the Earth magnetic torques τ _ei . Consider the following:

${\mathit{\tau }}_{di}={\mathit{\tau }}_{gi}+{\mathit{\tau }}_{ei}.$

(14)

A very important factor presented in the disturbance torque term is the torque that acts on the spacecraft due to the Earth's magnetic field. This torque can be written as [4]

${\mathit{\tau }}_{ei}={\mathit{\mu }}_i\times {\mathbf{B}}_e\left({\mathbf{R}}_i\right),$

(15)

where Be is Earth's magnetic field strength at the location of the spacecraft.

Substituting (10) into (11), the attitude dynamics can be written compactly as

${\mathbf{M}}_{iA}\left({\mathbf{J}}_i,{\mathbf{q}}_i\right){\ddot{\mathbf{\epsilon }}}_i+{\mathbf{C}}_{iA}\left({\mathbf{J}}_i,{\mathbf{\Omega }}_i,\mathbf{h}\right)+{\mathbf{N}}_{iA}\left({\mathbf{J}}_i,{\mathbf{q}}_i,{\mathbf{\omega }}_{bi},{\mathbf{\omega }}_{Oi},{\dot{\mathbf{\omega }}}_{Oi}\right)={\mathit{\tau }}_{cwi}+{\mathit{\tau }}_{mi}+{\mathit{\tau }}_{di},$

(16)

where

$\begin{array}{c}{\mathbf{M}}_{iA}\left({\mathbf{J}}_i,{\mathbf{q}}_i\right)=\mathrm{2}{\mathbf{J}}_i{{\mathbf{T}}_i}^{-\mathrm{1}},\\ {\mathbf{C}}_{iA}\left({\mathbf{J}}_i,{\mathbf{\Omega }}_i,\mathbf{h}\right)={\mathbf{\Omega }}_i^{\times }\left({\mathit{J}}_i{\mathbf{\Omega }}_i+\mathbf{h}\right),\hspace{1em}\hspace{1em}{\mathbf{X}}_i={\left[{\mathit{\rho }}_i^T,{\mathbf{\epsilon }}_i^T\right]}^T,\\ {\mathbf{N}}_{iA}\left({\mathbf{J}}_i,{\mathbf{q}}_i,{\mathbf{\omega }}_{bi},{\mathbf{\omega }}_{Oi},{\dot{\mathbf{\omega }}}_{Oi}\right)=\mathrm{2}{\mathbf{J}}_i{{\mathbf{T}}_i}^{-\mathrm{1}}{f}_i,\hspace{1em}{f}_i=-{\mathbf{\omega }}_{bi}^{\times }{\mathbf{\omega }}_{Oi}+{\dot{\mathbf{\omega }}}_{Oi}.\end{array}$

(17)

2.3. Combined Attitude and Translational Dynamics

Combining the attitude dynamics in (1) and the translational dynamics in (16), the following 6-DOF dynamics equation for formation flying is obtained:

${\mathbf{M}}_i{\ddot{\mathbf{X}}}_i+{\mathbf{C}}_i+{\mathbf{N}}_i={\mathbf{u}}_i+{\mathbf{D}}_i,$

(18)

where

$\begin{array}{c}{\mathbf{M}}_i=\left[\begin{array}{cc}{\mathbf{I}}_{\mathrm{3}\times \mathrm{3}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{M}}_{iA}\end{array}\right]\in {ℝ}^{\mathrm{6}\times \mathrm{6}},\hspace{1em}\hspace{1em}{\mathbf{C}}_i=\left[\begin{array}{c}{\mathbf{C}}_{iT}\\ {\mathbf{C}}_{iA}\end{array}\right]\in {ℝ}^{\mathrm{6}\times \mathrm{1}},\\ {\mathbf{N}}_i=\left[\begin{array}{c}{\mathbf{N}}_{iT}\\ {\mathbf{N}}_{iA}\end{array}\right]\in {ℝ}^{\mathrm{6}\times \mathrm{1}},\hspace{1em}{\mathbf{u}}_i=\left[\begin{array}{c}\widehat{{\mathbf{F}}_{ci}}\\ {\mathit{\tau }}_{cwi}\end{array}\right],\\ {\mathbf{D}}_i\triangleq {\left[\begin{array}{cccccc}{d}_{\mathrm{1}i}& d_{\mathrm{2}i}& d_{\mathrm{3}i}& d_{\mathrm{4}i}& d_{\mathrm{5}i}& d_{\mathrm{6}i}\end{array}\right]}^T=\left[\begin{array}{c}\widehat{{\mathbf{F}}_{di}}\\ {\mathit{\tau }}_{mi}+{\mathit{\tau }}_{di}\end{array}\right].\end{array}$

(19)

3. Adaptive Terminal Sliding Mode Control Design

In this section, an adaptive terminal sliding mode controller is designed for the follower spacecraft based on the dynamic model in (18) and the TSM technique. With this controller, the follower spacecraft can track the desired attitude and relative position trajectories simultaneously.

3.1. Error Dynamics Equation

The trajectory tracking errors of the follower spacecraft are defined as

${\mathit{\rho }}_i^e\triangleq {\mathit{\rho }}_i-{\mathit{\rho }}_i^d,\hspace{1em}\hspace{1em}{\dot{\mathit{\rho }}}_i^e\triangleq {\dot{\mathit{\rho }}}_i-{\dot{\mathit{\rho }}}_i^d,$

(20)

where ${\mathit{\rho }}_i^d,{\dot{\mathit{\rho }}}_i^d\in {ℝ}^{\mathrm{3}\times \mathrm{1}}$ are the relative position/velocity of the desired trajectory with respect to the leader. Differentiating (20) and substituting it into (1), one can obtain that

${\ddot{\mathit{\rho }}}_i^e+{\mathbf{C}}_{iT}\left({\omega }_O,{\dot{\mathit{\rho }}}_i\right)+{\mathbf{N}}_{iT}^e\left({\mathit{\rho }}_i,{\omega }_O,{\mathbf{R}}_{\mathrm{0}},{\ddot{\mathit{\rho }}}_i^d\right)=\widehat{{\mathbf{F}}_{di}}+\widehat{{\mathbf{F}}_{ci}},$

(21)

where ${\mathbf{N}}_{iT}^e({\mathit{\rho }}_i,{\omega }_O,{\mathbf{R}}_{\mathrm{0}},{\ddot{\mathit{\rho }}}_i^d)={\mathbf{N}}_{iT}+{\ddot{\mathit{\rho }}}_i^d$ .

The error Euler parameter is chosen to describe the attitude error of the spacecraft with respect to the reference attitude. The error Euler parameter is determined as follow [19]:

${\mathbf{\epsilon }}_i^e={\mathbf{U}}_i^{d^T}{\mathbf{q}}_i,$

(22)

${\eta }_i^e={\mathbf{q}}_i^{d^T\mathrm{\hspace{0.17em}\hspace{0.17em}}}{\mathbf{q}}_i,$

(23)

where ${\mathbf{U}}_i^d\equiv {\left[\begin{array}{cc}{\mathbf{T}}_i^{d^T}& -{\mathbf{\epsilon }}_i\end{array}\right]}^T\in {ℝ}^{\mathrm{4}\times \mathrm{3}}$ , ${\mathbf{T}}_i^d\equiv {\eta }_i^d\mathbf{I}+{\mathbf{\epsilon }}_i^{\times d}\in {ℝ}^{\mathrm{3}\times \mathrm{3}}$ , and ${\mathbf{q}}_i^d\equiv {\left[\begin{array}{cc}{\mathbf{\epsilon }}_i^{d^T}& {\eta }_i^d\end{array}\right]}^T\in {ℝ}^{\mathrm{4}\times \mathrm{1}}$ is the desired attitude.

Differentiating (22) with respect to time and utilizing (10), the attitude kinematics equations with error Euler parameter are given by

${\mathbf{\Omega }}_{bi}^{}=\mathrm{2}{{\mathbf{T}}_i}^{e-\mathrm{1}}{\dot{\mathbf{\epsilon }}}_i^e+{\mathbf{A}}_i^e{\mathbf{\omega }}_{bi}^d+{\mathbf{A}}_{bIi}{\mathbf{\omega }}_{Oi}^{},$

(24)

where ${\mathbf{A}}_i^e$ is the transformation matrix from the reference attitude frame to the body frame.

Differentiating (24), one can obtain that

${\dot{\mathbf{\Omega }}}_{bi}^{}=\mathrm{2}{\mathbf{T}}_i^{e^{-\mathrm{1}}}{{\ddot{\mathbf{\epsilon }}}_i}^e+{\mathbf{f}}_i^e,$

(25)

where

${\mathbf{f}}_i^e={\mathbf{f}}_i+\frac{{\mathbf{\epsilon }}_i^e}{\mathrm{2}{\eta }_i^{e\mathrm{2}}}\left({\mathbf{\omega }}_{bi}^{}-{\mathbf{A}}_i^e{\mathbf{\omega }}_{bi}^d\right)\left({\mathbf{\omega }}_{bi}-{\mathbf{A}}_i^e{\mathbf{\omega }}_{bi}^d\right)+{\mathbf{A}}_i^e{\dot{\mathbf{\omega }}}_{bi}^d-{\mathbf{\omega }}_{bi}^{\times }{\mathbf{A}}_i^e{\mathbf{\omega }}_{bi}^d$

(26)

Substituting (24) and (25) into (11), the error attitude dynamics are derived as follows [19]:

${\mathbf{M}}_{iA}^e\left({\mathbf{J}}_i,{\mathbf{q}}_i\right){\ddot{\mathbf{\epsilon }}}_i+{\mathbf{C}}_{iA}\left({\mathbf{J}}_i,{\mathbf{\Omega }}_i,\mathbf{h}\right)+{\mathbf{N}}_{iA}^e\left({\mathbf{J}}_i,{\mathbf{q}}_i^e,{\mathbf{\omega }}_{bi},{\mathbf{\omega }}_{bi}^d,{\dot{\mathbf{\omega }}}_{bi}^d\right)={\mathit{\tau }}_{cwi}+{\mathit{\tau }}_{mi}+{\mathit{\tau }}_{di},$

(27)

where ${\mathbf{N}}_{iA}^e=\mathrm{2}{\mathbf{J}}_i{\mathbf{T}}_i^{e^{-\mathrm{1}}}{\mathbf{f}}_i^e$ , ${\mathbf{M}}_{iA}^e=\mathrm{2}{\mathbf{J}}_i{\mathbf{T}}_i^{e^{-\mathrm{1}}}$ .

Let ${\mathbf{e}}_i={\left[\begin{array}{cc}{\mathit{\rho }}_i^{e^T}& {\mathbf{\epsilon }}_i^{e^T}\end{array}\right]}^T$ , and utilizing (18), (21), and (27), the error dynamics of the whole system can be obtained

${\mathbf{M}}_i^e{\ddot{\mathbf{e}}}_i+{\mathbf{C}}_i+{\mathbf{N}}_i^e={\mathbf{u}}_i+{\mathbf{D}}_i,$

(28)

where

${\mathbf{M}}_i^e=\left[\begin{array}{cc}{\mathbf{I}}_{\mathrm{3}\times \mathrm{3}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{M}}_{iA}^e\end{array}\right]\in {ℝ}^{\mathrm{6}\times \mathrm{6}},\hspace{1em}\hspace{1em}{\mathbf{N}}_i^e=\left[\begin{array}{c}{\mathbf{N}}_{iT}^e\\ {\mathbf{N}}_{iA}^e\end{array}\right]\in {ℝ}^{\mathrm{6}\times \mathrm{1}}.$

(29)

3.2. Controller Design

The terminal switching planes can be describe as [16–18]

${\mathbf{S}}_i={\dot{\mathbf{e}}}_i+{\mathit{\alpha }}_i{\mathbf{e}}_i+{\mathit{\beta }}_i{\mathbf{e}}_i^{{\gamma }_i},$

(30)

where ${\mathbf{S}}_i\triangleq {\left[\begin{array}{cccccc}{s}_{\mathrm{1}i}& s_{\mathrm{2}i}& s_{\mathrm{3}i}& s_{\mathrm{4}i}& s_{\mathrm{5}i}& s_{\mathrm{6}i}\end{array}\right]}^T$ , ${\mathit{\alpha }}_i=\mathrm{diag}(a_{1i},a_{2i},\dots ,a_{6i})$ , β _i = diag⁡(b_1i,b_2i,…, b_6i), a _n >0, b _n >0 (n = 1,2,…,6), and 0<γ _i <1. The unknown disturbances Di in (18) are bounded but unknown. σ _ni ^* is the upper bound of d _ni . Here, it can be estimated by

${\dot{\widehat{\sigma }}}_{ni}=\frac{\mathrm{1}}{{\vartheta }_{ni}}\left|s_{ni}\right|,\hspace{1em}\left(n=\mathrm{1,2},\dots ,\mathrm{6}\right),$

(31)

where ϑ _ni is the adaptive gain.

Utilizing (28), (30), and (31), an adaptive terminal sliding mode control law is designed as follows:

${\mathbf{u}}_i\triangleq \left[\begin{array}{c}{\mathbf{u}}_{iT}\\ {\mathbf{u}}_{iA}\end{array}\right]={\mathbf{C}}_i+{\mathbf{N}}_i^e+{\mathbf{M}}_i^e\left(-{\mathit{\alpha }}_i{\dot{\mathbf{e}}}_i-{\gamma }_i{\mathit{\beta }}_i\mathrm{diag}\left({{\mathbf{e}}_i}^{{\gamma }_i-\mathrm{1}}\right){\dot{\mathbf{e}}}_i-\mathrm{diag}\left({\widehat{\mathit{\sigma }}}_i\right)\mathrm{sign}\left({\mathbf{S}}_i\right)\right),$

(32)

where ${\widehat{\mathit{\sigma }}}_i\triangleq {\left[\begin{array}{cccccc}{\widehat{\sigma }}_{1i}& {\widehat{\sigma }}_{2i}& {\widehat{\sigma }}_{3i}& {\widehat{\sigma }}_{4i}& {\widehat{\sigma }}_{5i}& {\widehat{\sigma }}_{6i}\end{array}\right]}^T$ and the vector signum function sign⁡(S _i ) is a column of signum functions

$\mathrm{sign}\left({\mathbf{S}}_i\right)\triangleq {\left[\begin{array}{cccccc}\mathrm{sign}\left(s_{\mathrm{1}i}\right)& \mathrm{sign}\left(s_{\mathrm{2}i}\right)& \mathrm{sign}\left(s_{\mathrm{3}i}\right)& \mathrm{sign}\left(s_{\mathrm{4}i}\right)& \mathrm{sign}\left(s_{\mathrm{5}i}\right)& \mathrm{sign}\left(s_{\mathrm{6}i}\right)\end{array}\right]}^T.$

(33)

Substituting (32) into (18) produces the closed-loop dynamics

${\ddot{\mathbf{e}}}_i=-{\mathit{\alpha }}_i{\dot{\mathbf{e}}}_i-{\gamma }_i{\mathit{\beta }}_i\mathrm{diag}\left({{\mathbf{e}}_i}^{\gamma -\mathrm{1}}\right){\dot{\mathbf{e}}}_i-\mathrm{diag}\left({\widehat{\mathit{\sigma }}}_i\right)\mathrm{sign}\left({\mathbf{S}}_i\right)+{\mathbf{D}}_i.$

(34)

Consider a Lyapunov function as follows:

$V=\frac{\mathrm{1}}{\mathrm{2}}{\mathbf{S}}_i^T{\mathbf{S}}_i+\frac{\mathrm{1}}{\mathrm{2}}{\tilde{\mathit{\sigma }}}_i^T\mathrm{diag}\left({\mathit{\vartheta }}_i\right){\tilde{\mathit{\sigma }}}_i,$

(35)

where ${\mathit{\vartheta }}_i\triangleq {\left[\begin{array}{cccccc}{\vartheta }_{\mathrm{1}i}& {\vartheta }_{\mathrm{2}i}& {\vartheta }_{\mathrm{3}i}& {\vartheta }_{\mathrm{4}i}& {\vartheta }_{\mathrm{5}i}& {\vartheta }_{\mathrm{6}i}\end{array}\right]}^T$ . The estimation error is defined as ${\tilde{\mathit{\sigma }}}_i={\widehat{\mathit{\sigma }}}_i-{\mathit{\sigma }}_i^{*}$ , ${\mathit{\sigma }}_i^{\mathit{*}}\triangleq {\left[\begin{array}{cccccc}{\sigma }_{1i}^{\mathit{*}},& {\sigma }_{2i}^{\mathit{*}},& {\sigma }_{3i}^{\mathit{*}},& {\sigma }_{4i}^{\mathit{*}},& {\sigma }_{5i}^{\mathit{*}},& {\sigma }_{6i}^{\mathit{*}}\end{array}\right]}^T$ .

The first derivative of (30) can be expressed as

${\dot{\mathbf{S}}}_i={\ddot{\mathbf{e}}}_i+{\mathit{\alpha }}_i{\dot{\mathbf{e}}}_i+{\gamma }_i{\mathit{\beta }}_i\mathrm{diag}\left({{\mathbf{e}}_i}^{\gamma -\mathrm{1}}\right){\dot{\mathbf{e}}}_i.$

(36)

Differentiating (35) and utilizing (36), one can obtain that

$\begin{array}{l}\dot{V}={\mathbf{S}}_i^T{\dot{\mathbf{S}}}_i+{\tilde{\mathit{\sigma }}}_i^T\mathrm{diag}\left({\mathit{\vartheta }}_i\right){\dot{\tilde{\mathit{\sigma }}}}_i\\ ={\mathbf{S}}_i^T\left({\ddot{\mathbf{e}}}_i+{\mathit{\alpha }}_i{\dot{\mathbf{e}}}_i+{\gamma }_i{\mathit{\beta }}_i\mathrm{diag}\left({{\mathbf{e}}_i}^{\gamma -\mathrm{1}}\right){\dot{\mathbf{e}}}_i\right)+{\tilde{\mathit{\sigma }}}_i^T\parallel {\mathbf{S}}_i\parallel ,\end{array}$

(37)

where $\parallel {\mathbf{S}}_i\parallel \triangleq {\left[\begin{array}{cccccc}\left|s_{\mathrm{1}i}\right|& \left|s_{\mathrm{2}i}\right|& \left|s_{\mathrm{3}i}\right|& \left|s_{\mathrm{4}i}\right|& \left|s_{\mathrm{5}i}\right|& \left|s_{\mathrm{6}i}\right|\end{array}\right]}^T$ . Substituting (31) and (34) into (37) produces that

$\dot{V}=-{\mathbf{S}}_i^T\left(\mathrm{diag}\left({\widehat{\mathit{\sigma }}}_i\right)\mathrm{sign}\left({\mathbf{S}}_i\right)-{\mathbf{D}}_i\right)+\left({\widehat{\mathit{\sigma }}}_i^T-{\mathit{\sigma }}_i^{*T}\right)\parallel {\mathbf{S}}_i\parallel \le -\sum _{n=\mathrm{1}}^{\mathrm{6}}\left|s_{ni}\right|\left({\sigma }_{ni}^{*}-\left|d_{ni}\right|\right)<\mathrm{0}\hspace{1em}\text{for}\hspace{0.17em}{s}_{ni}\ne \mathrm{0}.$

(38)

It means that the switching planes s _ni (n = 1, 2, …, 6) converge to zero [20]. On the other hand, in the TSM s _ni = 0 (n = 1, 2, …, 6), the system state will reach zero in finite time [17].

Remark 1. For the purpose of eliminating chattering, a common practice is to replace the signum function of (32) with a continuous saturation function

$\mathrm{sat}\left(s_{ni},{\xi }_{ni}\right)=\{\begin{array}{ccc}\frac{s_{ni}}{{\xi }_{ni}}& \text{if}\hspace{0.17em}\left|s_{ni}\right|<{\xi }_{ni}& \\ & & n=\mathrm{1,2},\dots ,\mathrm{6},\\ \mathrm{sign}\left(s_{ni}\right)& \text{if}\hspace{0.17em}\left|s_{ni}\right|\ge {\xi }_{ni},& \end{array}$

(39)

where ξ _ni >0 is the width of the boundary layer. Equation (32) can be written as

${\mathbf{u}}_i={\mathbf{C}}_i+{\mathbf{N}}_i^e+{\mathbf{M}}_i^e(-{\mathit{\alpha }}_i{\dot{\mathbf{e}}}_i-{\gamma }_i{\mathit{\beta }}_i\mathrm{diag}\left({\mathbf{e}}_i^{{\gamma }_i-\mathrm{1}}\right){\dot{\mathbf{e}}}_i-\mathrm{diag}\left({\widehat{\mathit{\sigma }}}_i\right)\mathrm{sat}\left({\mathbf{S}}_i,{\mathbf{\xi }}_i\right)\phantom{\left({\mathbf{e}}_i^{{\gamma }_i-\mathrm{1}}\right)}),$

(40)

where

$\mathrm{sat}\left({\mathbf{S}}_i,{\mathbf{\xi }}_i\right)={\left[\begin{array}{cccccc}\mathrm{sat}(s_{\mathrm{1}i},{\xi }_{\mathrm{1}i})& \mathrm{sat}(s_{\mathrm{2}i},{\xi }_{\mathrm{2}i})& \mathrm{sat}(s_{\mathrm{3}i},{\xi }_{\mathrm{3}i})& \mathrm{sat}(s_{\mathrm{4}i},{\xi }_{\mathrm{4}i})& \mathrm{sat}(s_{\mathrm{5}i},{\xi }_{\mathrm{5}i})& \mathrm{sat}\left(s_{\mathrm{6}i},{\xi }_{\mathrm{6}i}\right)\end{array}\right]}^T.$

(41)

Remark 2. There exists a possible singularity in sliding mode controller as e _i →0. Since e _i = 0 only approaches along a sliding mode, we observe that for a general choice of γ _i

${\dot{\mathbf{e}}}_i=-{\mathit{\alpha }}_i{\mathbf{e}}_i-{\mathit{\beta }}_i\mathrm{diag}\left({\mathbf{e}}_i^{\gamma }\right)$

(42)

while sliding and that the component in (32)

${\gamma }_i{\mathit{\beta }}_i\mathrm{diag}\left({\mathbf{e}}_i^{\gamma -\mathrm{1}}\right){\dot{\mathbf{e}}}_i⟶-{\gamma }_i{\mathit{\beta }}_i\left({\mathit{\alpha }}_i\mathrm{diag}\left({\mathbf{e}}_i^{{\gamma }_i}\right)+{\mathit{\beta }}_i\mathrm{diag}\left({\mathbf{e}}_i^{\mathrm{2}{\gamma }_i-\mathrm{1}}\right)\right).$

(43)

Consequently there will be a singularity in (32) unless γ _i is chosen so that 2γ _i >1. To satisfy this requirement we set γ _i = 3/5 for the examples to follow.

3.3. Compute the Special Magnetic Moments of SMD

The adaptive sliding mode controller uiT gives the desired special forces that can be used as input for a thruster-based system. However, for EMFF, the desired special forces are produced by SMD. The SMD of each spacecraft is a complicated function of current position and orientation [3, 10]. For EMFF, the control variables consist of the dipoles of individual spacecraft. Utilizing (32) and (4), one can obtain that

${\mathbf{F}}_{cm\mathrm{1}}=\frac{m_{\mathrm{0}}{m}_{\mathrm{1}}}{m_{\mathrm{0}}+m_{\mathrm{1}}}{\mathbf{u}}_{iT},\hspace{1em}\hspace{1em}{\mathbf{F}}_{cm\mathrm{0}}=-\frac{m_{\mathrm{0}}{m}_{\mathrm{1}}}{m_{\mathrm{0}}+m_{\mathrm{1}}}{\mathbf{u}}_{iT}.$

(44)

Substituting (5) and (7) into (44) produces

$\sum _{j=\mathrm{0},\hspace{0.17em}j\ne i}^{N-\mathrm{1}}{f}_{ij}({\mathit{\mu }}_i,{\mathit{\mu }}_j,{\mathbf{r}}_{ij})=\frac{m_i{m}_j}{m_i+m_j}{\mathbf{u}}_{iT}.$

(45)

Note that (7) gives the force on the ith SMD due to the jth SMD and it depends on the distance between the two SMDs and the orientation of both SMDs in the inertial space. It rises to the complexity of the expression for the magnetic. Here, we consider a two-spacecraft electromagnetic formation flying (EMFF) array. By defining a rotated frame F ^r (see Section 2), a simplified algebraic form of (7) is obtained

$f_{ij}^r\left({\mathit{\mu }}_i,{\mathit{\mu }}_j,{\mathbf{r}}_{ij}\right)=\frac{\mathrm{3}{\mu }_{\mathrm{0}}}{\mathrm{4}\pi r_{ij}^{\mathrm{4}}}\left[\begin{array}{c}\mathrm{2}{\mu }_{ix}{\mu }_{jx}-{\mu }_{iy}{\mu }_{jy}-{\mu }_{iz}{\mu }_{jz}\\ -{\mu }_{ix}{\mu }_{jy}-{\mu }_{iy}{\mu }_{jx}\\ -{\mu }_{ix}{\mu }_{jz}-{\mu }_{iz}{\mu }_{jx}\end{array}\right].$

(46)

Utilizing (45) and (46) produces that

$\left[\begin{array}{c}\mathrm{2}{\mu }_{ix}{\mu }_{jx}-{\mu }_{iy}{\mu }_{jy}-{\mu }_{iz}{\mu }_{jz}\\ -{\mu }_{ix}{\mu }_{jy}-{\mu }_{iy}{\mu }_{jx}\\ -{\mu }_{ix}{\mu }_{jz}-{\mu }_{iz}{\mu }_{jx}\end{array}\right]=\frac{\mathrm{4}\pi r_{ij}^{\mathrm{4}}}{\mathrm{3}{\mu }_{\mathrm{0}}}\frac{m_i{m}_j}{m_i+m_j}{\mathbf{A}}_{rIi}{\mathbf{A}}_{Ibi}{\mathbf{u}}_{iT}^{},$

(47)

where ArIi is the transformation matrix from the inertial ECI reference frame F ^I to the rotated frame F ^r . AIbi is the transformation matrix from the body frame F ^Bi of ith spacecraft to ECI frame.

Let ${\breve{\mathbf{u}}}_{iT}^r\triangleq {[{\breve{u}}_{iTx}^r,{\breve{u}}_{iTy}^r,{\breve{u}}_{iTz}^r]}^T=(4\pi r_{12}^4/3{\mu }_0)(m_{\mathrm{1}}{m}_{\mathrm{2}}/(m_{\mathrm{1}}+m_{\mathrm{2}})){\mathbf{A}}_{rIi}{\mathbf{A}}_{Ibi}{\mathbf{u}}_{iT}^{}$ , μ _i = −μ _j , and utilizing (47), one can obtain that

$\begin{array}{l}{\mu }_{ix}=\frac{1}{2}\sqrt{\sqrt{{\breve{u}}_{iTx}^{r2}+2{\breve{u}}_{iTy}^{r2}+2{\breve{u}}_{iTz}^{r2}}-{\breve{u}}_{iTx}^{r2}},\\ {\mu }_{iy}=\frac{{\breve{u}}_{iTy}^{r2}}{\sqrt{\sqrt{{\breve{u}}_{iTx}^{r2}+2{\breve{u}}_{iTy}^{r2}+2{\breve{u}}_{iTz}^{r2}}-{\breve{u}}_{iTx}^{r2}}},\\ {\mu }_{iz}=\frac{{\breve{u}}_{iTz}^{r2}}{\sqrt{\sqrt{{\breve{u}}_{iTx}^{r2}+2{\breve{u}}_{iTy}^{r2}+2{\breve{u}}_{iTz}^{r2}}-{\breve{u}}_{iTx}^{r2}}}.\end{array}$

(48)

The special magnetic moments of SMD can be computed using (48), by which the follower spacecraft can track the desired relative position trajectories.

3.4. Angular Momentum Management

For EMFF in near-Earth orbits, a constant disturbance torque may act on each spacecraft due to the Earth's magnetic field that causes the reaction wheels on each spacecraft to quickly become saturated. In order to avoid this situation, the angular momentum management strategy herein can be utilized, which is based on a unique nonlinearity of the magnetic torques [3, 4]. The force acting between a pair of dipoles depends on the product of the individual dipole values. By switching the polarity of all dipoles in the EMFF, the torque acting on the spacecraft due to Earth's magnetic field changes sign, but the torques and forces due to the other spacecraft in the system do not. As can be seen from (15), it results in a net cancellation of the effect of the Earth's magnetic field on the average sense.

4. Numerical Simulation

The adaptive terminal sliding mode controller equation (40) was simulated for a two-spacecraft formation flying control. Considering the nonlinear dynamics with disturbance D₁, the effectiveness of the proposed controller was verified. The orbital parameters of leader and follower spacecraft in the simulations are listed in Table 1. The parameters of spacecraft are summarized in Table 2. The dynamics includes the second harmonic of the gravitational field as well as the Earth's magnetic field, of which the vector is calculated using the World Magnetic (WMM2005) Model block in the paper.

Table 1:

Orbital parameters.

Spacecraft	Semimajor axis	Inclination	Right ascension of ascending node	Argument of perigee	Eccentricity	Mean anomaly
LS	6878 km	1.106538745764 rad	0.523598775598 rad	0 rad	0.003	0 rad
FS	6878 km	1.106538745764 rad	0.523598071512 rad	0.000000315261 rad	0.002999636522	0 rad

Table 2:

Parameter values of model.

Spacecraft	Mass	Moment of inertia
Ls	100 kg	Diag (19, 19, 32) kg·m2
FS	100 kg	Diag (19, 19, 32) kg·m2

The follower was commanded to move around the leader in an elliptic orbit. The desired trajectory was generated by solving nonlinear equation (1) numerically (set $\widehat{{\mathbf{F}}_{d\mathrm{1}}}$ and $\widehat{{\mathbf{F}}_{c\mathrm{1}}}$ equal to zero) with the following initial condition:

$\begin{array}{c}{\mathit{\rho }}_{d\mathrm{1}}\left(\mathrm{0}\right)=\left[\begin{array}{ccc}\mathrm{2.499998640217},& \mathrm{0.000000000217},& \mathrm{4.317138211336}\end{array}\right]\left[\text{m}\right],\\ {\dot{\mathit{\rho }}}_{d\mathrm{1}}\left(\mathrm{0}\right)=\left[\begin{array}{ccc}\mathrm{0.000000000000},& -\mathrm{0.000005559083863},& -\mathrm{0.000000000758}\end{array}\right]\left[\text{m}/\text{s}\right].\end{array}$

(49)

The relative position/velocity of the follower was initialized to

$\begin{array}{c}{\mathit{\rho }}_{\mathrm{1}}\left(\mathrm{0}\right)=\left[\begin{array}{ccc}\mathrm{2.599999999855},& \mathrm{0.722841614669},& \mathrm{4.503332099822}\end{array}\right]\left[\text{m}\right],\\ {\dot{\mathit{\rho }}}_{\mathrm{1}}\left(\mathrm{0}\right)=\left[\begin{array}{ccc}-\mathrm{0.000000000000},& -\mathrm{0.005557083863},& -\mathrm{0}.\mathrm{000000000758}\end{array}\right]\left[\text{m}/\text{s}\right].\end{array}$

(50)

The parameters of the adaptive sliding mode controller are then

$\begin{array}{c}{\alpha }_1=\mathrm{diag}\left(\mathrm{2,2},\mathrm{3,250,80,80}\right)\times {\mathrm{10}}^{-4},\\ {\beta }_1=\mathrm{diag}\left(\mathrm{10,15,15,1200,800,800}\right)\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}},\\ {\gamma }_{\mathrm{1}}=\frac{\mathrm{3}}{\mathrm{5}},\hspace{1em}\hspace{1em}\xi =\left[\mathrm{5,5},\mathrm{5,1},\mathrm{1,1}\right]\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}},\\ {\vartheta }_{\mathrm{1}}=\left[\mathrm{15,22,22,80,125,125}\right].\end{array}$

(51)

In this simulation, a switching time of 100 second was used for AMM. Figures 2–6 show the simulation results, which are obtained by simulating the adaptive terminal sliding mode controller. The phase portrait of the trajectory ρ₁(t) of the follower spacecraft relative to the leader spacecraft is illustrated in Figure 3, where * represents the leader spacecraft at the origin. Figure 3 depicts the position tracking error of the follower spacecraft. Figure 4 gives attitude error of the follower spacecraft. Figure 5 shows the magnetic dipole strength of the follower spacecraft. Figure 6 gives the reaction wheel angular momentum buildup of the follower spacecraft, and it is clear from this figure that the angular momentum buildup was limited by AMM. In a stable tracking situation, simulation results show that $\parallel {\mathit{\rho }}_{\mathrm{1}}^e\parallel$ is smaller than 6 × 10⁻³ m and $\parallel {\mathbf{\epsilon }}_{\mathrm{1}}^e\parallel$ is smaller than 0.075° which meets the required tracking accuracy.

Figure 2:

Actual trajectory ρ₁(t) of the follower spacecraft (* the leader spacecraft.)

Figure 3:

Position tracking error.

Figure 4:

Attitude tracking error.

Figure 5:

Magnetic dipole strength in F^B1 frame.

Figure 6:

Reaction wheel angular momentum.

5. Conclusions

In this paper, we consider the control problem for six-degree-of-freedom electromagnetic spacecraft formation flying (EMFF) in near-Earth orbits. Firstly, the combined attitude and translational dynamics are presented. Secondly, we confine our attention mainly to the position/attitude tracking control using the terminal sliding mode and adaptive control technique. Based on a rotated frame F ^r and the adaptive TSM controller, the special magnetic moments of SMD are computed. The AMM strategy is implemented in a periodically switching fashion. Through numerical simulation, it is clear that the output tracking error can converge to zero in finite time by using the proposed controller and that the angular momentum buildup was limited by AMM. At the same time, strong robustness with respect to bounded disturbances can be guaranteed (in the simulation, the J2 disturbance, the gravity gradient torques τ_g1, the Earth magnetic torques τ_e1, and the magnetic torque term τ_m1 have been considered).