Hybrid propulsion spacecraft formation control around the planetary displaced orbit

Abstract

This paper aims to investigate the feasibility of using the combination of solar radiation pressure and Coulomb force as a propellantless control method for spacecraft formation around the planetary displaced orbit. Firstly, the dynamical equation of spacecraft formation is derived and linearized. Based on the linearized dynamic model, an integral sliding mode controller (ISMC) is designed. Aimed to stabilize the spacecraft formation, the control method is proposed to adjust the product of the charge and the attitude angles of two spacecrafts. Finally, numerical simulations are conducted and the results show that the controller can make the formation achieve the desired configuration with favorable control performances.

Introduction

Spacecraft formation flying has been considered as an attractive technology because of the improvement it offers over large monolithic spacecraft in terms of system reliability and launch flexibility. Controlling the relative motion between spacecrafts is traditionally conducted via propulsion system. As an emerging technology, propellantless control has attracted the attention of many researchers. Propellantless control methods including using the combination of Coulomb force and solar sail would be widely used space exploration missions, which can save propellant and avoid thruster plume impingement.

The Coulomb formation was first proposed by King et al. ¹ in 2002. It uses the electrostatic attraction or repulsion between spacecrafts as a control force without fuel consumption. Solar sail is a propellantless propulsion system that can generate continuous thrust over a long period by reflecting solar photons. The concept of solar sail propulsion was first proposed by Tsiolkovsky in 1921, and it has been widely studied in deep space exploration missions around libration point orbit, halo orbit, displaced orbit, etc. Much research has focused on using solar sailing as the primary propulsion system on a spacecraft to maintain NKOs. ² However, there are still some challenges to overcome, such as the difficulties of designing and building large and lightweight membrane structures. In addition, the solar sail is unable to generate a thrust component toward the Sun. ³ Compared to solar sailing, Coulomb propulsion of free-flying vehicles is to control the spacecraft formation shape and size using the inter-spacecraft forces created by electrostatically charging the spacecraft to different potentials. ¹ It can produce a relatively low thrust with high specific impulse. It can be successfully used in several space missions. However, as the electrostatic forces are internal to the formation, Coulomb force cannot be used to reorient a full formation to a new orientation. Coulomb force cannot be used to control the center of mass of the formation, therefore, to reorient a spacecraft formation required for a special mission, external forces such as thrusters or solar radiation pressure must be used.

There are several prior references with regard to hybrid propulsion system in the literature, which can be divided into two categories, propellant control method and propellantless control method.

Considering the complementarity of solar sailing and SEP, the concept (hybrid sail) of hybridizing a relatively small, near-term solar sail and SEP has been previously proposed.^4,5 In this field, research is flourishing, aiming to maximize the potential of the hybrid sail propulsion system. Proposed concepts include optimal transfers from Earth to Venus and Mars.^6,7 Based on the feedback linearization method, Simo and McInnes ⁸ realized the maintenance and control of the hybrid sail around the libration point orbit. Zhang et al. ⁹ designed a self-anti-interference control method that can realized the formation control of the hybrid propulsion spacecraft around the displaced solar orbit. Gong^4,10–12 has carried out further research on the solar sail spacecraft, and classified the displaced orbits according to different orbit periods. Based on the hybrid sail formation around the geosynchronous displaced orbit, Qin et al. ¹³ devised an optimal propulsion strategy around the geostationary displaced orbit. Chen et al. ¹⁴ used a high-precision sliding mode controller to control the formation of the hybrid sail spacecraft around the displaced solar orbit. Lou ¹⁵ built the improved model of solar radiation pressure and the thrust produced by the solar sail with variable reflectivity, and solved the underactuation of solar radiation pressure propulsion system. Obviously, the existence of SEP would lead to the requirement of the propulsion system, which indicates that the aforementioned hybrid formulations are propellant methods rather than propellantless ones.

In other field, Saaj et al. ¹⁶ proposed a propellant control method using Coulomb forces and electric/ion propulsion for spacecraft formation flying in GEO and other high Earth orbits. A tetrahedron formation scenario is used to demonstrate the effectiveness of the hybrid propulsion technique. Hybrid formulations combining Lorentz force with impulsive thrusts based on both orbital elements-based model and Cartesian formulation were proposed.^17,18 Huang et al.^19,20 proposed a hybrid formulation that consists of the specific charge and the thruster-generated control acceleration of Lorentz spacecraft. The resulting strategy is a thruster-assisted one but not a propellantless one. Sun et al. ²¹ suggested a propellantless control method for spacecraft formation-keeping problem by using the combination of the Lorentz force and aerodynamic force and designed adaptive output feedback control algorithm can guarantee the effectiveness of the proposed propellantless control approach in the presence of external perturbations and unavailability of velocity measurements. Zuo et al. ²² proposed a new propellantless control method using the combination of Coulomb force and solar sail, and designed the LQR controller to control the formation around the displaced solar orbit.

At present, the control algorithms applied to formation control are developing rapidly. Zhou et al. ²³ addressed the formation control of unmanned surface vehicles and designed adaptive fuzzy backstepping-based control algorithm to solve unknown nonlinear items in the USV. The effectiveness of the proposed approach is demonstrated by the experimental results. Huang et al. ¹⁹ designed a adaptive fast nonsingular terminal sliding mode control law, based on which an adaptive controller is designed for Lorentz-augmented spacecraft relative motion to deal with the uncertain parameters and perturbations and proposed a fuel-optimal distribution law of Lorentz force and traditional chemical propulsion. Numerical simulations substantiate the feasibility and validity of the proposed controller for Lorentz-augmented relative orbital control. In this paper, an integral sliding mode control algorithm is used to solve complex spacecraft dynamics problems, and simulation results show that the control algorithm has good performance.

The main contribution of the paper is to propose a propellantless strategy for full-dimensional spacecraft formation control using the combination of Coulomb force and solar radiation pressure force. The propellantless control spacecraft is a new type of continuous small thrust spacecraft, on the one hand, it has the dual characteristics of solar sail propulsion without energy consumption and Coulomb force propulsion with high efficiency. On the other hand, it overcomes the defect that solar sail cannot provide the propulsion component pointing to the Sun, and is suitable for carrying out complex orbital missions. Additionally, the proposed integral sliding mode controller can guarantee the robustness of the proposed propellantless control approach in the presence of initial disturbance.

Dynamic equations of spacecraft formation flying around planetary displaced orbits

This section establishes the mathematical model of the equation of motion. It is assumed that both spacecraft, namely, the chief and deputy, are equipped with several solar sails. The satellites have the ability to modulate their electrostatic charge. A relative dynamic model is established, which considers the influences of the Coulomb force and solar radiation pressure. Two reference frames are defined. $\mathrm{OXYZ}$ is an planet-centered inertial frame and O is the center of the planet. The X axis is along the angular velocity of planet’s orbit around the Sun, the Z axis is along the Sun-planet line, and the Y axis forms a right-handed coordinate system with X and Z axes. The displaced orbit is a circular orbit in the anti-sun direction, as shown in Figure 1. $o\xi \varsigma \eta$ is a rotating frame associated to the displaced orbit. o is the center of chief spacecraft, the $\xi$ axis points from the origin to the center of the displaced orbit, the $\eta$ axis is parallel to the X axis, and the $\varsigma$ axis forms a right-handed coordinate system with $\xi$ and $\eta$ axes (as shown in Figure 1). The solar radiation pressure is assumed to be uniform because the variation near the planet is very small. The angular velocity of the planet is far less than the angular velocity of the spacecraft around the planet. Therefore, the angular velocity of the planet around the Sun can be ignored. Meanwhile, solar sails are assumed to be ideal full-reflection sails. The chief spacecraft evolves on a displaced orbit and the deputy spacecraft is on a nearby orbit.

Figure 1.

Planetary displaced orbit.

As shown in Figure 1, ${\mathit{r}}_{\mathit{m}}$ and ${\mathit{r}}_{\mathit{s}}$ are the position vectors of the chief spacecraft and the deputy spacecraft relative to the planet, respectively. $h$ is the height of the displaced orbit, and $\tilde{\rho }$ is the radius of the displaced orbit. In the frame $o\xi \varsigma \eta$ , the position vector of the deputy spacecraft can be written as $\mathit{\rho }={\mathit{r}}_{\mathit{s}}-{\mathit{r}}_{\mathit{m}}={[x,y,z]}^T$ .

In the frame $\mathrm{OXYZ}$ , the equation of motion the spacecrafts can be given by

$$\frac{d^2{\mathit{r}}_{\mathit{m}}}{d{t}^2}=-{\mathit{f}}_{\mathit{g}}({\mathit{r}}_{\mathit{m}})+{\mathit{f}}_{\mathit{s}}({\mathit{r}}_{\mathit{m}})+{\mathit{f}}_{\mathit{c}}({\mathit{r}}_{\mathit{m}})$$

(1)

$$\frac{d^2{\mathit{r}}_{\mathit{s}}}{d{t}^2}=-{\mathit{f}}_{\mathit{g}}({\mathit{r}}_{\mathit{s}})+{\mathit{f}}_s({\mathit{r}}_{\mathit{s}})+{\mathit{f}}_{\mathit{c}}({\mathit{r}}_{\mathit{s}})$$

(2)

where ${\mathit{f}}_{\mathit{g}}(\mathit{r})=\frac{{\mu }_{\mathrm{e}}}{r^3}·\mathit{r}$ is the planet’s gravitational acceleration, ${\mathit{f}}_{\mathit{s}}(\mathit{r})=\beta \frac{\mu }{{({\mathit{r}}_a+\mathit{r})}^2}·{\cos }^2\alpha ·\mathit{n}$ is solar radiation pressure acceleration, ${\mathit{f}}_{\mathit{c}}\left({\mathit{r}}_{\mathit{m}}\right)=-k_c·\frac{Q}{m_m{L}^3}(1+\frac{L}{{\lambda }_d})\exp (-\frac{L}{{\lambda }_d})\mathit{\rho }$ is Coulomb force acceleration acting on the chief spacecraft, ${\mathit{f}}_{\mathit{c}}\left({\mathit{r}}_{\mathit{s}}\right)=k_c·\frac{Q}{m_s{L}^3}(1+\frac{L}{{\lambda }_d})\exp (-\frac{L}{{\lambda }_d})\mathit{\rho }$ is Coulomb force acceleration acting on the deputy spacecraft, $\beta$ is the lightness number of the spacecraft. $a$ is defined as the angle between the unit normal vector of the solar sail $\left(\mathit{n}\right)$ and sunlight, ${\mathit{r}}_{\mathit{a}}$ is planet’s position vector with respect to the sun. Because the solar radiation pressure is assumed to be uniform, the sunlight direction can be given by $\mathit{l}={\left[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right]}^T$ in both reference frames above. The solar radiation pressure can be written as:

$${\mathit{f}}_{\mathit{s}}(\mathit{r})=\beta \frac{\mu }{{{\mathit{r}}_{\mathit{a}}}^2}·{\cos }^2\alpha ·\mathit{n}=\kappa ·{\cos }^2\alpha ·\mathit{n}=\kappa ·(\mathit{n}·\mathit{l}{)}^2·\mathit{n}$$

$\kappa$ is the characteristic acceleration, $\mu$ is the solar gravitational constant, and $k_c$ is electrostatic constant. ${\lambda }_d$ is the length of Debye. $Q=q_1·q_2$ is the product of the charge of two spacecrafts. Combining equations (1) and (2), the equation of relative motion in OXYZ frame can be written as

$$\begin{array}{c}\frac{d^2\mathit{\rho }}{d{t}^2}=\frac{d^2{\mathit{r}}_{\mathit{s}}}{d{t}^2}-\frac{d^2{\mathit{r}}_{\mathit{m}}}{d{t}^2}={\mathit{f}}_{\mathit{g}}({\mathit{r}}_{\mathit{m}})-{\mathit{f}}_{\mathit{g}}({\mathit{r}}_{\mathit{s}})+{\mathit{f}}_{\mathit{s}}({\mathit{r}}_{\mathit{s}})-{\mathit{f}}_{\mathit{s}}({\mathit{r}}_{\mathit{m}})\hfill \\ +k_c\frac{Q}{L^3}(\frac{1}{m_m}+\frac{1}{m_s})(1+\frac{L}{{\lambda }_d})\exp (-\frac{L}{{\lambda }_d})\mathit{\rho }\hfill \end{array}$$

(3)

The relationship of derivatives in two reference frames can be given by:

$$\begin{array}{c}\multicolumn{1}{c}{\frac{d^2\mathit{\rho }}{d{t}^2}=\frac{d^2{\mathit{r}}_{\mathit{s}}}{d{t}^2}-\frac{d^2{\mathit{r}}_{\mathit{m}}}{d{t}^2}=\stackrel{··}{\mathit{\rho }}+2\mathit{\omega }\times \mathit{\rho }+\mathit{\omega }}\\ \multicolumn{1}{c}{\times (\mathit{\omega }\times \mathit{\rho })+\stackrel{·}{\mathit{\omega }}\times \mathit{\rho }}\end{array}$$

(4)

Where $\frac{d}{\mathrm{dt}}$ denotes derivative with respect to the OXYZ frame, and $\stackrel{·}{\mathit{\rho }}$ denotes derivative with respect to the $o\xi \varsigma \eta$ frame. The rotating frame moves around the center of the displaced orbit at a constant speed, the angular velocity is a constant vector, $\stackrel{·}{\mathbf{\omega }}=0$ . Therefore, the equation of relative motion in the rotating frame can be written as:

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{··}{\mathit{\rho }}+2\mathit{\omega }\times \mathit{\rho }+\mathit{\omega }\times (\mathit{\omega }\times \mathit{\rho })=\mathrm{\Delta }{\mathit{f}}_{\mathit{g}}+\mathrm{\Delta }{\mathit{f}}_{\mathit{s}}}\\ \multicolumn{1}{c}{+k_c\frac{Q}{L^3}(\frac{1}{m_m}+\frac{1}{m_s})(1+\frac{L}{{\lambda }_d})\exp (-\frac{L}{{\lambda }_d})\mathit{\rho }}\end{array}$$

(5)

Assume that the distance between two spacecraft be very small compared with $\left|{\mathit{r}}_{\mathit{m}}\right|$ , $\left|{\mathit{r}}_{\mathit{s}}\right|$ . Then the Coulomb force term can be linearized. The acceleration term due to the Coulomb force can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{f}}_{\mathit{c}}=k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}}\\ \multicolumn{1}{c}{·[-\exp (-\frac{L_0}{{\lambda }_d})r·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]\mathit{\rho }}\end{array}$$

(6)

To maintain the deputy spacecraft in the vicinity of the chief spacecraft, the accelerations of the solar radiation pressure need to be similar. Therefore, ${\mathit{f}}_{\mathit{s}}({\mathit{r}}_{\mathit{s}})$ can be expanded around $({\mathit{n}}_{\mathit{m}},{\kappa }_m)$ :

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{f}}_{\mathit{u}}={\mathit{f}}_{\mathit{s}}({\mathit{r}}_{\mathit{s}})-{\mathit{f}}_{\mathit{s}}({\mathit{r}}_{\mathit{m}})=\frac{\partial {\mathit{f}}_{\mathit{s}}}{\partial \mathit{n}}({\mathit{n}}_{\mathit{s}}-{\mathit{n}}_{\mathit{m}})+\frac{\partial {\mathit{f}}_{\mathit{s}}}{\partial \kappa }({\kappa }_s-{\kappa }_m)}\\ \multicolumn{1}{c}{={\cos }^{2}a·{\mathit{n}}_{\mathit{m}}·\delta \kappa +\kappa {\cos }^{2}a·E_3·\delta \mathit{n}}\end{array}$$

(7)

The relative motion equation can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{··}{\mathit{\rho }}+2\mathit{\omega }\times \mathit{\rho }+\mathit{\omega }\times (\mathit{\omega }\times \mathit{\rho })=\frac{{\mu }_e}{{r_m}^3}[{\mathit{r}}_{\mathit{m}}-{\left(\frac{r_m}{r_s}\right)}^3{\mathit{r}}_{\mathit{s}}]}\\ \multicolumn{1}{c}{+{\mathit{f}}_{\mathit{c}}+{\mathit{f}}_{\mathit{u}}}\end{array}$$

(8)

$\mathit{\rho }$ and $\mathit{\omega }$ can be projected in the rotating coordinate system as $\mathit{\rho }={\left[xyz\right]}^{\mathbf{T}}$ , $\mathit{\omega }={\left[00\mathrm{n}\right]}^{\mathbf{T}}$ , where $n=\sqrt{\frac{{\mu }_e}{{\rho }^3}}$ is. The position vectors of the two spacecrafts relative to the planet are given by ${\mathit{r}}_{\mathit{m}}={\left[r_{m}00\right]}^{\mathbf{T}}$ and ${\mathit{r}}_{\mathit{s}}={[r_m+xyz]}^{\mathbf{T}}$ , respectively. The simplification of the gravitational acceleration is similar to the derivation of the C-W equation.

Thus, the term of $r_m$ / $r_s$ in the gravitational acceleration can be linearized as

$$\begin{array}{c}\multicolumn{1}{c}{r_s=\sqrt{{(r_m+x)}^2+y^2+z^2}=\sqrt{{\rho }^2+{r_m}^2+2x{r}_m}}\end{array}$$

(9)

$$\begin{array}{c}\multicolumn{1}{c}{{\left(\frac{r_m}{r_s}\right)}^3={\left(\frac{{\rho }^2+{r_m}^2+2x{r}_m}{r_m}\right)}^{-\frac{3}{2}}={[1+{\left(\frac{\rho }{r_m}\right)}^2+\frac{2x}{r_m}]}^{-\frac{3}{2}}}\\ \multicolumn{1}{c}{\approx {(1+\frac{2x}{r_m})}^{-\frac{3}{2}}}\end{array}$$

(10)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{r}}_{\mathit{m}}-{\left(\frac{r_m}{r_s}\right)}^3{\mathit{r}}_{\mathit{s}}\approx {\mathit{r}}_{\mathit{m}}-(1-3\frac{x}{r_m})({\mathit{r}}_{\mathit{m}}+\mathit{\rho })\approx 3\frac{x}{r_m}{\mathit{r}}_{\mathit{m}}-\mathit{\rho }}\\ \multicolumn{1}{c}{={[2x-y-z]}^{\mathbf{T}}}\end{array}$$

(11)

Using equations (8)–(11), equation (12) can be expressed as:

$$\begin{array}{c}\stackrel{··}{x}-2n\stackrel{·}{y}-n^{2}x-2x·{\mu }_e/{r_m}^3={\mathit{f}}_{\mathit{cx}}+{\mathit{f}}_{\mathit{ux}}\hfill \\ \stackrel{··}{y}+2n\stackrel{·}{x}-n^{2}y+y·{\mu }_e/{r_m}^3={\mathit{f}}_{\mathit{cy}}+{\mathit{f}}_{\mathit{uy}}\hfill \\ \stackrel{··}{z}+z·{\mu }_e/{r_m}^3={\mathit{f}}_{\mathit{cz}}+{\mathit{f}}_{\mathit{uz}}\hfill \end{array}$$

(12)

Where six right-hand items are the components of the Coulomb acceleration ${\mathit{f}}_{\mathit{c}}$ and the solar radiation pressure acceleration ${\mathit{f}}_{\mathit{u}}$ , respectively. Because the sunlight direction around the planetary displaced orbit is assumed to be parallel to the $\eta$ axis in the rotating frame.

In the rotating frame, the unit normal vector of the solar sail can be given as (as shown in Figure 2):

$$\mathit{n}=\sin \varphi \sin \theta ·\mathit{i}+\sin \varphi \cos \theta ·\mathit{j}+\cos \varphi ·\mathit{k}$$

(13)

Figure 2.

Description of solar sail normal vector in rotating frame.

In the rotating frame, the normal vector of the deputy spacecraft can be expressed as

$${\mathit{n}}_{\mathit{s}}=\left[\begin{array}{ccc}\sin \varphi \frac{x}{\sqrt{x^2+y^2}}\hfill & \sin \varphi \frac{y}{\sqrt{x^2+y^2}}\hfill & \cos \varphi \hfill \end{array}\right]$$

(14)

$\delta \mathit{n}$ can be approximated in the position feedback form by linearizing ${\mathit{n}}_{\mathit{s}}$ around ${\mathit{n}}_{\mathit{m}}$ , given by

$$\delta \mathit{n}={\mathit{n}}_{\mathit{s}}-{\mathit{n}}_{\mathit{m}}={\frac{\partial \mathit{n}}{\partial \mathit{r}}|}_{r_m}\mathit{\rho }=\left[\begin{array}{ccc}0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \sin \varphi (1/\tilde{\rho })\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}x\hfill \\ \multicolumn{1}{c}{y}\\ \multicolumn{1}{c}{z}\end{array}\right]$$

(15)

Substitution of equation (15) into equation (12) gives

$$\begin{array}{c}\stackrel{··}{x}-2n\stackrel{·}{y}-n^{2}x-2x·{\mu }_e/{r_m}^3={\cos }^2\varphi ·\sin \varphi ·\sin \theta ·\delta \kappa +\hfill \\ k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}\hfill \\ \multicolumn{1}{c}{·[-\exp (-\frac{L}{{\lambda }_d})·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]·x}\\ \stackrel{··}{y}+2n\stackrel{·}{x}-n^{2}y+y·{\mu }_e/{r_m}^3=\kappa ·{\cos }^2\varphi \hfill \\ \multicolumn{1}{c}{·\sin \varphi (1/\tilde{\rho })+{\cos }^2\varphi ·\sin \varphi ·\cos \theta ·\delta \kappa +}\\ k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}\hfill \\ \multicolumn{1}{c}{·[-\exp (-\frac{L}{{\lambda }_d})·(r\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]·y}\\ \stackrel{··}{z}+z·{\mu }_e/{r_m}^3={\cos }^2\varphi ·\cos \varphi ·\delta \kappa +k_c·Q\hfill \\ \multicolumn{1}{c}{·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}·}\\ [-\exp (-\frac{L}{{\lambda }_d})·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})r]·z\hfill \end{array}$$

(16)

Semi-natural formation is mainly studied in this paper, the sails of two spacecraft are designed with the same capacity of lightness number, namely $\delta \beta =0$ , (formation flying near the planet, two spacecrafts are approximately the same distance from the Sun).

The equations of relative motion can be written as

$$\begin{array}{c}\stackrel{··}{x}-2n\stackrel{·}{y}-n^{2}x-2x·{\mu }_e/{r_m}^3=k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}·\hfill \\ [-\exp (-\frac{L}{{\lambda }_d})·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]·x\hfill \\ \stackrel{··}{y}+2n\stackrel{·}{x}-n^{2}y+y·{\mu }_e/{r_m}^3=\kappa ·{\cos }^2\varphi ·\sin \varphi (1/\tilde{\rho })\hfill \\ +k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}·\hfill \\ [-\exp (-\frac{L}{{\lambda }_d})·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]·y\hfill \\ \stackrel{··}{z}+z·{\mu }_e/{r_m}^3=k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}·\hfill \\ [-\exp (-\frac{L}{{\lambda }_d})·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]·z\hfill \end{array}$$

(17)

The components of the hybrid propulsion acceleration are given by

$$\begin{array}{c}{\mathit{a}}_{\mathit{x}}=k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}·\hfill \\ [-\exp (-\frac{L}{{\lambda }_d})·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]·x\hfill \\ {\mathit{a}}_{\mathit{y}}=\kappa ·{\cos }^2\varphi ·\sin \varphi (1/\tilde{\rho })+\hfill \\ k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}·\hfill \\ [-\exp (-\frac{L}{{\lambda }_d})·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]·y\hfill \\ {\mathit{a}}_{\mathit{z}}=k_c·Q·(\frac{1}{m_m}+\frac{1}{m_s})·\frac{1}{{\lambda }_d}·\hfill \\ [-\exp (-\frac{L}{{\lambda }_d})·(\frac{3{\lambda }_d}{L^4}+\frac{3}{L^3}+\frac{1}{{\lambda }_d·L^2})]·z\hfill \end{array}$$

(18)

When the chief spacecraft is on the planetary displaced orbit, the reference attitude angles can be obtained from Li et al. ²⁴ In the second attitude control scheme of semi-natural formation, the angle $\varphi$ of the deputy spacecraft is the same as that of the chief spacecraft. Therefore, the controller mainly adjusts the attitude angle $\varphi$ of the chief spacecraft and the charge product $Q$ of the spacecrafts to realize the control of formation configuration.

Controller design

The integral sliding mode controller is an improved form of traditional sliding mode controller. Traditional sliding mode control includes approaching and sliding processes. In the approach process, the state of the system will reach the sliding mode surface, then the state will move along the sliding mode surface until it reaches the target state. The integral sliding mode control eliminates the approach process, and ensures that the state of the system can reach the sliding mode surface quickly. Therefore, the controller has a strong robust stability during the whole operational process.

The state variables are defined as $x_1=(x,y,z{)}^T$ , $x_2=(\stackrel{·}{x},\stackrel{·}{y},\stackrel{·}{z}{)}^T$ . The equations of motion is written as

$$\{\begin{array}{c}{\stackrel{·}{x}}_1=x_2\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{x}}_2=f(x_1,x_2)+u+d}\end{array}$$

(19)

Where

$$f(x_1,x_2)=\{\begin{array}{c}2n\stackrel{·}{y}+n^{2}x+2x·{\mu }_e/{r_m}^3\hfill \\ \multicolumn{1}{c}{-2n\stackrel{·}{x}+n^{2}y-y·{\mu }_e/{r_m}^3}\\ \multicolumn{1}{c}{-z·{\mu }_e/{r_m}^3}\end{array}$$

(20)

$d$ in equation (19) are the external disturbance acceleration. The controller can make the actual relative distance between the chief and deputy spacecraft approach the reference distance gradually. The relative distance error is defined as $e(t)=x_1(t)-x_{1d}(t)$ , and the control goal is $e(t)\to 0$ .

The sliding mode surface of the system is defined as

$$s(t)=\stackrel{·}{e}(t)-\stackrel{·}{e}(0)+k_p[e(t)-e(0)]+k_i{\int }_0^{t}e(\tau )d\tau$$

(21)

Where $k_p=diag[240,240,240]$ and $k_i=[0.6,0.6,0.6]$ are positive definite matrices and coefficient vector, respectively. To suppress the chatter phenomenon of sliding mode control, the exponential approach law is written as

$$\stackrel{·}{s}(t)=-\mathrm{Ps}(t)-Q\tanh [s(t)]$$

(22)

Where $p=diag[0.01,0.01,0.01]$ , $Q=diag[0.08,0.08,0.08]$ . The sliding mode control law is finally described as

$$\begin{array}{c}\multicolumn{1}{c}{u(t)=-[f(x_1,x_2)-{\stackrel{··}{x}}_{1d}(t)]-k_p\stackrel{·}{e}(t)-k_{i}e(t)}\\ \multicolumn{1}{c}{-\mathrm{Ps}(t)-Q\tanh [s(t)]}\end{array}$$

(23)

Consider a Lyapunov function as follows:

$$V=\frac{1}{2}{\mathit{S}}^T\mathit{S}$$

(24)

Differentiating equation (24), and utilizing equation (21), one can obtain that

$$\begin{array}{c}\stackrel{·}{V}={\mathit{S}}^{\mathit{T}}\stackrel{·}{\mathit{S}}\hfill \\ ={\mathit{S}}^{\mathit{T}}(\stackrel{··}{\mathit{e}}(t)+k_p\stackrel{·}{\mathit{e}}(t)+k_i\mathit{e}(t))\hfill \end{array}$$

(25)

Where $\stackrel{··}{e}(t)={\stackrel{··}{x}}_1-{\stackrel{··}{x}}_{1d}$ . The error dynamics equation can be written as:

$$\stackrel{··}{e}(t)=f(x_1,x_2)+u+d-{\stackrel{··}{x}}_{1d}$$

(26)

Substituting equations (26) (23) into equations (25) produces that

$$\stackrel{·}{V}={\mathit{S}}^{\mathit{T}}(d-Q\tanh [s(t)])$$

(27)

For a satellite around the planetary displaced orbit, planetary perturbation is the dominant perturbation. The perturbation acceleration $d$ due to space perturbation can be assumed to be bounded. $F$ is the upper bound of $d$ .

$$\stackrel{·}{V}=\sum _{i=1}^3(\left|S_i\right|F_i-Q_i\left|S_i\right|)\le -\sum _{i=1}^3\left|S_i\right|(Q_i-F_i)\le 0$$

(28)

Where $\Vert \mathit{S}\Vert ={\left[\begin{array}{ccc}\left|S_1\right|\hfill & \left|S_2\right|\hfill & \left|S_3\right|\hfill \end{array}\right]}^T$ . $\stackrel{·}{V}\le 0$ can be made by controlling the gain $Q_i\ge F_i$ . According to Lyapunov’s theorem, the system state is stable and can converge to zero in finite time.

In order to demonstrate the superiority of the designed controller, contrasting it with traditional sliding mode control (SMC). The sliding mode surface of the system is defined as: $s(t)=\lambda e(t)+\stackrel{·}{e}(t)$ .

The exponential approach law is written as: $\stackrel{·}{s}(t)=-\mathrm{ks}(t)-\epsilon \tanh [s(t)]$ , the sliding mode control law is finally described as:

$$\begin{array}{c}\multicolumn{1}{c}{u(t)=-[f(x_1,x_2)-{\stackrel{··}{x}}_{1d}(t)]-(\lambda +k)·\stackrel{·}{e}(t)}\\ \multicolumn{1}{c}{-k\lambda ·e(t)-\epsilon ·\tanh [s(t)]}\end{array}$$

(29)

$x_{\mathrm{ref}}$ , $y_{\mathrm{ref}}$ , $z_{\mathrm{ref}}$ are the components of the ideal relative distance between the chief and deputy spacecraft, respectively. $\mathrm{ex}$ , $\mathrm{ey}$ , $\mathrm{ez}$ are the error components of the relative distance.

$\varphi$ , $Q$ are the attitude angle and the charge product of the spacecrafts, respectively. Using equation (18), the input acceleration used to maintain the deputy spacecraft in the reference relative orbit can be obtained. In this paper, the input acceleration is mainly provided by the attitude angle and charge product of the spacecrafts.

Numerical simulations

This section examines the performance of the proposed formation-maintenance strategy through numerical simulations. The simulation scenario is set to perform long-period on-orbit missions in which two spacecrafts form master-slave formation around the planetary displaced orbit. It is assumed that the chief spacecraft is in a planetary displaced orbit of the first type (The period of the displaced orbit is the same as the period of the polar orbit of radius $\tilde{\rho }$ , given by $\omega =\sqrt{{\mu }_e/{\tilde{\rho }}^3}$ ). According to the results in Gong et al., ¹¹ the stability region of the first kind of planetary displaced orbit under the second attitude control scheme is $\tilde{\rho }\ge 2.2635h$ , therefore, $\tilde{\rho }$ is set $\tilde{\rho }=2.5h$ .

In the rotating frame $o\xi \varsigma \eta$ , deputy spacecraft is evolving in a circular relative orbit with a radius of 100 m around the chief spacecraft. The orbital period is chosen as 1/100 of the lunar revolution period. Initial error $\delta x=10\mathrm{m}$ is given at the initial time, the integral sliding mode controller can be used to design the coulomb force and the solar radiation pressure force to make the spacecraft converge to the desired orbit.

Some simulation parameters are shown in Table 1 and Table 2.

Table 1.

Model parameters.

Parameter name	Parameter	Numerical	Unit
The quality of the chief spacecraft	$m_m$	$100$	$\mathrm{Kg}$
The quality of the deputy spacecraft	$m_s$	$100$	$\mathrm{Kg}$
Orbital altitude	$h$	$50$	$\mathrm{Km}$
The ideal value of relative distance	$L_{\mathrm{ref}}$	$100$	$\mathrm{m}$
Light pressure factor	$\beta$	$0.20$
The position vector of the Earth relative to the Sun	$r_a$	$1$	$\mathrm{AU}$
Initial error	$\delta x$	$10$	$\mathrm{m}$
Initial error	$\delta \varphi$	$1$
Static dielectric constants	$k_c$	$8.99\times {10}^9$	$\mathrm{N}{\mathrm{m}}^2/{\mathrm{C}}^2$
The size of a solar sail	$s_1\times s_2$	$20\times 20$	${\mathrm{m}}^2$
Earth radius	$R_e$	$6.371\times {10}^6$	$\mathrm{m}$
Orbit radius	$\tilde{\rho }$	$125$	$\mathrm{Km}$
Geocentric gravitational constant	${\mu }_e$	$3.986\times {10}^{14}$
Debye length	${\lambda }_d$	$25$	$\mathrm{m}$

Table 2.

Controller parameters.

Control method
ISM	$k_p=\mathrm{diag}[240,240,240]$ , $k_i=[0.6,0.6,0.6]$ $p=\mathrm{diag}[0.01,0.01,0.01]$ , $Q=\mathrm{diag}[0.08,0.08,0.08]$
SMC	$\lambda =\mathrm{diag}[0.02,0.02,0.02]$ , $k=\mathrm{diag}[0.08,0.08,0.08]$ , $\epsilon =\mathrm{diag}[0.01,0.01,0.01]$

The numerical simulation results are as follows:

The phase portrait of the trajectory $\rho (t)$ of the deputy spacecraft relative to the chief spacecraft is illustrated in Figure 3, where reference orbit tracking can be rapidly achieved with the designed controller using a hybrid propulsion system. The solar sails are used to control the center of mass in the planetary displaced orbit plane. As the electrostatic forces are internal to the formation, Coulomb force is responsible to control the spacecraft formation shape and size. Figures 4 to 8 show the simulation results, which are obtained by simulating the integral sliding mode controller and sliding mode controller. Figure 4 depicts the time histories of attitude angle for the solar sail. After 5 days, the attitude angle of the duty spacecraft becomes stable around a fixed angle, and small angle fluctuation continuously provides small solar radiation pressure force control (ISMC). On the contrary, it is clear from Figure 4 that the attitude angle converges slowly to the steady state by using traditional sliding mode control. The charge product of the spacecraft is shown in Figure 5. By charging and discharging, the controller adjusts the magnitude of coulomb force. The error caused by the initial disturbance needs to be adjusted about 20d (ISMC). After that, the charge product basically remains stable, and the subsequent change curve will not be shown. Time histories of the relative position errors and relative distance are shown in Figures 6 and 8, respectively, from which it is clear that both controllers eliminate the initial relative state errors and guarantee the system robustness against the external disturbances, but the proposed ISMC presents distinctively faster convergence rate. Take the relative position errors for example, it takes about 6d for ISMC to drive the system to the desired relative states while SMC requires a duration of nearly 10d. Furthermore, terminal relative position accuracies of both ISMC and SMC are, respectively, on the order of 10^-3 and 10^-2 m. Figure 7 depicts the time histories of the control inputs. As can be seen, the control inputs are nonsingular with no significant chattering. Aforementioned simulation results substantiate the validity and efficiency of the proposed ISMC. The ISMC presents enhanced performances on the convergence rate, mean error distance, and percent overshoot.

Figure 3.

Relative trajectory of the deputy spacecraft.

Figure 4.

Attitude angle change.

Figure 5.

Charge product’s change.

Figure 6.

Relative distance error and enlarge the image after stabilization.

Figure 7.

Control input.

Figure 8.

Relative distance.

In summary, it can be seen that the controller can make the relative position errors converge to the neighborhood of zero rapidly, which verifies the robustness and speed ability of the proposed integral sliding mode controller.

Conclusion

This paper mainly studied the relative motion dynamics and control of the hybrid propulsion spacecraft formation around the planetary displaced orbit. A propellantless control method by using the coulomb force and solar radiation pressure is presented. Firstly, the formation’s nonlinear dynamics model is established. Then, through the linearization, the linearized dynamic model with weak coupling is obtained. Finally, numerical simulations are provided to substantiate the validity of the proposed control strategy in the presence of initial perturbation. And then the simulation results show that relative motion in the formation can be effectively controlled by the integral sliding mode control method.

In this paper, when modeling the solar sail spacecraft, the flexible characteristics of the solar sail and the structural characteristics brought by the large size structure are not considered, but the solar sail spacecraft is only considered as a particle. The flexibility and structural characteristics of solar sails are rarely analyzed at home and abroad. These characteristics can be taken into account in subsequent studies to obtain a more complete model profit. When considering the solar sail, it is regarded as an ideal solar sail model, that is, there is no loss of light utilization and the shielding effect of other celestial bodies or objects on the sail surface is not considered. It is difficult to provide such an ideal situation in practical application, so practical factors can be added as limiting conditions in subsequent studies.