# Multiple Spacecraft Formation Flying Control around Artificial Equilibrium Point Using Propellantless Approach

Authors: Lei Zhao; Changqing Yuan; Xiaoming Li; Jingjiu He
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

This paper demonstrates a detailed analysis of the feasibility for compact formation system around an L2-type artificial equilibrium point by means of continuous low-thrust propulsion in the hybrid form of solar sail and Coulomb force propulsion. Firstly, in view of non-ideal solar sail, the position of L2-type artificial equilibrium point and numerical periodic orbits around L2 utilized as leader’s nominal trajectory are given. Secondly, considering the external disturbances in the deep space environment, the nonlinear dynamic model of the spacecraft formation system based on the circular restricted three-body problem (CRTBP) is derived, under the assumption that the leader covers the nominal trajectory and each follower adjusts its propulsive acceleration vector (that is, both its sail attitude and electrostatic charge) in order to track a desired relative trajectory. Thirdly, based on a new double power combination function reaching law, a fast integral terminal sliding mode control methodology (MFITSM) is ameliorated to achieve orbital tracking rapidly, which has better robustness against external disturbances and the buffeting effect during spacecraft propulsion simultaneously. To properly allocate control inputs, a novel optimal allocation scheme is designed to calculate the charge product of the spacecrafts and sail attitude angles, which can make the magnitude of the acceleration required from the Coulomb propulsion system minimum and avoid formation geometry instabilities by balancing electrostatic interaction between adjacent spacecraft. Finally, several numerical examples are conducted to validate the superiority of the proposed control algorithm.


---

## Body

## 1. Introduction

Formation flight of spacecraft has always been a hot topic in space research. Spacecraft formation is generally composed of a group of small space-distributed spacecrafts orbiting each other, and each satellite in the formation forms a virtual spacecraft through information interaction and cooperative work, in order to achieve the functions of large satellites or the missions that cannot be achieved by traditional single large satellite. Some of its main applications include outer space exploration, space scientific experiments, monitoring the earth and its surrounding atmosphere, and deep space imaging or exploration. Compared with single spacecraft, the formation flight of spacecraft can get longer flight duration, stronger capacity of disturbance rejection, more flexibility for execution of the mission, and lower cost. Controling the relative position of spacecraft accurately is crucial for spacecraft formation system. Conventional formation control is performed using chemical thruster, which consume energy and can generate caustic plumes to contaminate onboard equipments. Therefore, the formation control by means of propellantless propulsion is of great significance. As an emerging technology, propellantless control methods were proposed utilizing the solar sail [1], Coulomb forces [2], electromagnetic forces [3], tethered systems [4], momentum exchange [5], etc. In addition, when some exploration missions are carried out in low orbit (LEO), it is possible to use the aerodynamic forces [6] and Lorentz force [7] as the control resource. In this paper, we propose a novel propellantless propulsion technology by using the combination of Coulomb force and solar radiation pressure which would be widely used deep space exploration missions. The merits of this propulsion technology over the conventional propulsion method include no fuel consumption, and the contamination from plumes is avoided, higher propulsive efficiency and longer flight duration.Formation flight of spacecraft around the libration points to observe heavenly bodies has been extensively studied. For example, the cosmic microwave background can be effectively measured by deploying a formation of spacecraft in a periodic orbit around L2 point. Since sunlight hitting the orbits around L1 point is not blocked by the Earth, spacecraft formation in such orbits can give warning of interstellar storms [8]. To set aside enough reaction time for safety protection equipment, the spacecrafts are considered to be deployed in the orbit around the artificial equilibrium point [9], which is closer to the Sun than the natural equilibrium point. Since the orbits around the artificial equilibrium point is unstable, a continuous control force is vital to realize the long-term bounded relative motion of spacecraft.Although the solar radiation pressure provided by the solar sail around the artificial equilibrium point can offset the gravitational imbalance between the Sun and the Earth, orbit motion is still unstable, so active orbit control is needed. At present, the research results on orbit control of solar sail around the artificial equilibrium point are abundant. Lawrence and Piggott [10] and Zhang et al. [11] proposed LQR controller with attitude angle amplitude constraint. Bookless and McInnes [12] realized the transfer of a solar sail spacecraft from Earth orbit to an artificial equilibrium point orbit using an LQR controller. All the above methods take the two attitude angles of the sail surface as the control input to control the position of the solar sail in the three-dimensional space, which leads to the underactuated case of the control system. To that end, some solutions have been proposed for the underactuated problem. Gong and Li [13] proposed a variable reflectivity solar sail model (RCD) based on the reflectivity control principle used by the IKAROS solar sail, and designed an artificial equilibrium point orbit keeping controller based on Floquet theory. But the design of variable reflectivity solar sail model is too ideal, such as continuous treatment of discrete RCD, ignoring the power supply film covered by sail, and the assumption of high surface to mass ratio, which are difficult to be applied in engineering practice. Soldini et al. [14] proposed a variable structure solar sail model, taking the sail area as the third control input, and designed a Hamiltonian structure holding controller. However, once solar sail is applied in practice, it is usually difficult to achieve flexible structural changes, so the scheme’s practicality is questionable. Zhang [15] adopted the ideal reflective solar power sail structure, ignored the influence of solar cells, and designed an active disturbance rejection track keeping controller. This scheme is also too ideal, and its electric propulsion fuel consumption is large, which cannot meet the long-term mission requirements. At present, the studies of formation around AEP mainly focus on large-scale formation, with formation spacing ranging from several thousand meters to tens of kilometers [16, 17], and the research on extremely close formation is few (below 100 meters). On the one hand, the extremely close formation requires higher control precision; a slight mistake may lead to the loss of configuration and collision damage. On the other hand, ion propulsion and chemical propulsion, which are commonly used in large-scale formation and cannot be used in extremely close formation due to plume pollution. Recognizing these open issues of the previous work, we propose a propellantless propulsion method using the combination of Coulomb force and solar radiation pressure for deep space exploration missions.King et al. [18] firstly proposed the idea of using electrostatic force for spacecraft formation. Recently, considering that the constraint conditions mainly include maximum current saturation limit and maximum spacecraft charging value, Felicetti and Palmerini [19] investigated three spacecrafts formation configuration control problem in GEO by means of electrostatic force. And an optimal charge allocation strategy [20] is proposed and discussed, that is, the charge is distributed to all spacecrafts starting from the calculated charge product. Coulomb force, also known as electrostatic force, is a propellantless propulsion method controlled by electrostatic attraction or repulsion between charged spacecrafts. Solar sail can generate continuous small thrust for long-lived formation missions due to the momentum exchange caused by solar photons hitting the sail membrane, which has been widely applied to spacecraft formation flying around the equilibrium point, halo orbit, displaced orbit, etc. Similarly, the continuous thrust required to maintain a non-Keplerian orbit could be provided by a solar sail [21]. However, the requirements of large size and lightweight materials limit the practical application of solar sail. In addition, the component of the propulsive acceleration produced by the solar sail be pointed towards the Sun [22]. Compared to solar sail, Coulomb force is the inter-spacecraft force, which can be used to adjust the formation geometry and size [18–20]. And Coulomb propulsion is capable of generating propulsion components towards the Sun. As a clean energy source, Coulomb propulsion has been used for some formation missions due to its ability to generate low thrust and high specific impulse. Since Coulomb force is the internal force of formation, Coulomb propulsion cannot be used to change the position of the center of mass and formation orientation. Therefore, the external force of formation is required for the formation maneuver missions, such as electric thruster, chemical thruster, solar sail, E-sail, or magnetic sail.At present, the research on hybrid propulsion for spacecraft formation has seen a lot of interest, which can be roughly divided into two categories: propellant and propellantless propulsion.Hybrid sail propulsion combining solar sail with solar electric propulsion (SEP) has been proposed for non-Kepler-orbit maintenance, and some of the major contributions in this area are described below. Simo and McInnes [23] investigated the hybrid sail spacecraft formation flying around the libration point orbit using feedback linearization methodology. Based on an ADRC technique, Zhang et [24] realized the hybrid sail spacecraft formation control around the heliocentric displaced orbit. Qin et al.[25] established a hybrid propulsion strategy based on SEP and solar sail, conducted simulation experiments for geosynchronous displaced orbit under different lightness numbers, and optimized the control process of solar sail attitude angle and electric propulsion acceleration. Chen et al.[26] proposed a novel sliding mode controller to realize the spacecraft formation control around the heliocentric displaced orbit, which is robust to internal unmodeled dynamics and external unknown disturbances. In addition, the control algorithm has a high control accuracy. Notably, the above hybrid sail propulsion mode belongs to a propellant propulsion due to the fuel requirements of the solar electric propulsion.For the research on other hybrid propulsion methods, Huang et al. [27, 28] introduced a new hybrid propulsion strategy by means of Lorentz force and control force generated by the thruster, which is the part of propellant propulsion system. The formation keeping of spacecraft can be effectively controlled using a propellantless propulsion method [29], which combines Lorentz force and aerodynamic force. An adaptive output feedback control law is introduced to verify the feasibility of the proposed propulsion method against external disturbances and velocity uncertainties. Recently, the designed hybrid propulsion system for the linearized spacecraft formation model around planetary displaced orbits has been presented, which guarantees a given formation geometry by using the combination of ideal solar sail and Coulomb force as a propellantless control method [30].In recent years, various advanced control algorithms have emerged and been applied to spacecraft formation control around artificial equilibrium points. The linear quadratic controller is the first to be used because of its simple structure and easy implementation. Wang et al. [31], Marchand and Howell [32] and Roberts [33] used LQR method to solve the formation control problem around the equilibrium point. Feedback linearization control is another linearization technique to solve nonlinear system control problems, and it also can be applied to formation flight around the equilibrium point. Folta et al.[34] proposed a basic LQR control law to linearize orbits around the libration point. Based on the LTI model, Hamilton et al.[35] designed a LQR control method for formation-keeping and orbital maneuver of SI mission around L2-type libration point. The nonlinear system is approximated to linear system in the above study, and the control effect often cannot meet the requirements of high-precision formation missions, and the system robustness is poor. For the purpose of improving the performance of the control algorithms, some scholars proposed different nonlinear control strategies based on the nonlinear dynamics system of formation flight. Wong and Kapila et al. [36] provide an adaptive controller for the spacecraft formation system around the L2-type libration point. Peng et al. [37] designed an optimal period control law to simultaneously solve the problem of orbit maintenance and formation control in the leader-follower formation configuration. Based on polynomial eigenstructure assignment, Wang et al.[38] designed a nonlinear controller to control formation flying around the L2-type equilibrium point. In addition, nonlinear sliding mode control has been applied many times in equilibrium point formation flight problems. Compared with LQR and other linear control strategies, these nonlinear control strategies achieve higher control accuracy and better system robustness. The above control methods have been verified that they have good control performance for formation keeping. However, some deficiencies still exist and cannot be ignored: the dependence on system model especially poor capability of disturbance rejection. These external disturbances are unavoidable in the execution of formation missions, so they need to be taken into account in the establishment of dynamic models. To ensure system stability and robustness, enhance control accuracy, and get faster convergence rate, the modified fast integral terminal sliding mode control law with a new double power combination function reaching law, known as MFITSMC, is presented in this note.The aim of this note is to investigate the problem of formation flying control around the artificial equilibrium point by means of the combination of Coulomb force and non-ideal solar sail, which can achieve extremely compact formations at distances of 100 meters or less. The hybrid propulsion is a novel formation control method, which has the characteristics of no fuel consumption, high efficiency, the ability to provide continuous small thrust, and so on. Moreover, it overcomes the deficiency that simple solar sail propulsion cannot provide the propulsion component pointing towards the sun, so it is suitable to carry out some specific deep space exploration missions. Compared with the previous hybrid propulsion research, the nonlinear research method and the non-ideal solar sail model are adopted in this paper, which are more suitable for dealing with the strong coupling and nonlinear problems existing in actual spacecraft orbit control. Additionally, the modified fast integral terminal sliding mode controller is designed with a new double power combination function reaching law for spacecraft formation system. Then, in order to allocate control commands more reasonably, the novel optimal allocation scheme is designed by means of particle swarm optimization algorithm and the MATLAB’s fmincon function, which ensures the magnitude of the acceleration required from the Coulomb propulsion system is minimum and avoids formation geometry instabilities by balancing electrostatic interaction between adjacent spacecraft. Finally, numerical simulation results illustrate that the proposed control algorithm has better control performance than the conventional control laws against the perturbations in space environment and system uncertainties.The rest of this paper is organized as follows. The mathematical model of spacecraft formation flying around the L2-AEP is put forward in Section2. The control algorithm is proposed for the nonlinear dynamics model in Section 3. Numerical simulations are analyzed in Section 4. Finally, the paper is summarized in Section 5.


## 2. Mathematical Model

This paper considers a mission scenario in which multiple spacecrafts flying around an L2-type artificial equilibrium point to explore the surface of the Sun and deep space. In the process of formation establishment, the leader moves along the nominal trajectory and each follower adjusts its thrust vector (that is, the combination of solar radiation pressure and Coulomb force) in order to track a desired relative trajectory. In this section, the halo orbit around L2 point is designated as the nominal trajectory and the nonlinear relative dynamic model of the spacecraft formation system is derived. The research in this paper is assumed to be carried out in the Circle Restricted Three-Body Problem model. The following are the contents about artificial equilibrium points, nominal trajectory, spacecraft relative dynamic model, and the external disturbance model.


### 2.1. Generation of the Artificial Equilibrium Point and Nominal Trajectory of the Leader

The base case of the CRTBP used in this paper is described in Figure1, which can be assumed that two celestial bodies (Sun and Earth) rotate around their center of gravity in circular orbits with a constant angular velocity. The mass of the spacecraft can be ignored with respect to the mass of the two celestial bodies. Therefore, the effects of the third body spacecraft on the two central bodies are negligible. Three collinear libration points L1, L2, and L3 are also shown in Figure 1. For convenience, the reference frames used in the paper are given below.Figure 1
Relevant coordinate systems.RI - Sun-Earth barycenter synodic reference frame: the origin of this reference is located at the Sun-Earth’s center-of-mass, with plane XY coinciding to the ecliptic. The X unit vector is pointed from the Sun towards the Earth, the Z-axis is normal to the plane of the celestial bodies’ orbit, and the Y-axis completes the right hand orthogonal frame system and is thus defined as normal to the X-axis in the plane of the celestial bodies’ orbit and along the prograde rotational direction.RB - body-fixed reference frame located at the spacecraft’s center of mass: with origin at the spacecraft center of mass o, which also coincides with the solar sail center of mass. The unit vectors i, j, and k are assumed to lie along the principal axes of inertia of the (ideal) unwrinkled sail.In order to facilitate the subsequent calculation, the system parameters are normalized. The total system’s mass is defined as:(1)m=m1+m2,where m1 and m2 represent the masses of the celestial bodies (Sun and Earth), while μ is the mass ratio of the planetary system.
(2)μ=m2m1+m2.To make it easier to normalize variables, the following assumptions need to be made:(1)
The total mass of the celestial bodies is normalized to one(2)
The distance between the celestial bodies is normalized to one(3)
The angular velocity of the central bodies rotating around their center of mass is normalized to oneAccordingly, the period of one Earth–Sun revolution is2π, the Sun-[Earth+Moon] reference distance is rSE=1AU, the dimensionless mass of the [Earth+Moon] system is μ=m2/m1+m2=3.05425×10−6, the O-Sun distance is rS=μrSE=μ, the O-[Earth+Moon] distance is rE=1−μrSE=1−μ, and the dimensionless mass of the Sun can be given as m1=1−μm1+m2=1−μ.The propulsive acceleration vectoras of non-ideal solar sail can be presented as:
(3)as=β1−μrSs2rSs⋅nb1rSs+b2rSs⋅n+b3n,where rSs is the Sun-solar sail vector, β is referred to as the light number, which defines the (dimensionless) ratio of the maximum propulsive acceleration magnitude to the local Sun’s gravitational acceleration, which is taking values between [0,1]. n is the unit normal vector of solar sail. The terms {b1, b2, b3} are the sail force coefficients [39] defined as
(4)b1=121+sρ,b2=12ρ1−sBf−ρ−1εfBf−εbBbεb+εf,b3=121−sρ.The values of the optical parameters of a non-ideal solar sail membrane, involved in the calculation of {b1, b2, b3}, can be obtained from the experimental data collected during the NASA scout mission [40].
(5)ρ=0.91,s=0.89,εf=0.0025,εb=0.27,Bf=0.79,Bb=0.67.The sail normal unit vectorn can be expressed as a function of two attitude angles θ,φ, written as
(6)n=cosθcosφ,cosθsinφ,sinθT,where θ∈−π/2,π/2 rad is the angle between the direction of n and the i,j plane, while φ∈−π/2,π/2 rad is the angle between the direction of i and the projection of n onto i,j, as shown in Figure 1.Consider wrinkled and imperfect solar sails, which move within the Sun-[Earth+Moon] system. Combining equations (3) and (6), the dimensionless propulsive acceleration vector as of the non-ideal solar sail becomes
(7)asx=β1−μrSs2cosθcosφb1+b2cos2θcos2φ+b3cosθcosφasy=β1−μrSs2cosθcosφcosθsinφb2cosθcosφ+b3asz=β1−μrSs2cosθcosφ−sinθb2cosθcosφ+b3.The equations of orbital motion of solar sail are formulated based on Newton’s law of motion and law of universal gravitation:(8)r¨=−Gm1rSs3rSs−Gm2rEs3rEs+as.Therefore, the non-dimensional equations of motion is(9)x¨=2y.+x−1−μx+μrSs3−μx−1+μrEs3+asxy¨=−2x.+y−1−μyrSs3−μyrEs3+asyz¨=−1−μzrSs3−μzrEs3+asz,where rEs is the Earth-solar sail vector and r=xyzT is the position vector of the solar sail. Due to the existence of solar sail, the balance of forces on spacecraft in space has changed. When the solar radiation pressure balances both the celestial body gravitational pull and the centrifugal force along the Sun-[Earth+Moon] line, a new family of collinear AEPs is generated. Compared with classical equilibrium points, artificial equilibrium points are closer to the sun, which enables spacecraft to perform better in missions, such as observing the Sun and relaying deep space communications.Let us assume that the sail membrane is perpendicular to the direction of solar radiation, according to the definition of equilibrium point(10)x¨=y¨=z¨=0x.=y.=z.=0y=z=0θ=φ=0x=xAEPX^.Substituting Eq. (10) into Eq. (9), the result is a single scalar equation
(11)xAEP−1−μ1−βb1+b2+b3xAEP+μ2−μxAEP+μ−12=0.The solution of Eq. (11) gives the location of the L2-type AEP along the X-axis. Increasing the solar sail lightness number has the effect of shifting the position of L2 towards the Sun. Hence, the required non-ideal sail lightness number β necessary to maintain an L2-type AEP is
(12)β=xAEP+μ21−μ/xAEP+μ2+μ/xAEP+μ−12−xAEP1−μb1+b2+b3.Figure2 shows the relationship curve between the position of the artificial equilibrium point and sail lightness number.Figure 2
Non-ideal solar sail lightness numberβ as a function of dimensionless Sun-AEP distance xAEP.As shown in Figure2, the larger β is, the collinear artificial equilibrium point of solar sail is closer to the central celestial body (Sun). Based on this situation, the solar sail can achieve more accurate observation of the Sun and deep space communication missions by being placed on a non-Keplerian orbit around the AEP, which can hardly be realized by conventional spacecraft.For the purpose of generating quasi-halo orbits around the artificial equilibrium points, the differential corrections methodology is adopted in this paper, which is generally used to calculate unstable long-periodic orbits [41], such as displaced non-Keplerian orbits and J2 invariant relative orbits. For detailed research on differential corrections algorithm, please refer to the works of Howell [41] and Ross [42]. The differential correction algorithm is not self-starting and requires an initial approximate solution. Richardson expansion equation obtained according to Lindstedt-Poincare method [43, 44] can be used as an effective initial value:
(13)x=a21Ax2+a22Az2−Axcosωt+a23Ax2−a24Az2cos2ωt+a31Az3−a32AxAz2cos3ωt,y=Ay+b33Ax3+b34−b35AxAz2sinωt+b21Ax2−b22Az2sin2ωt+b31Ax3−b32AxAz2sin3ωt,z=−3d21AxAz+Azcosωt+d21AxAzcos2ωt+d32AzAx2−d31Az3cos3ωt,where the values of the parametersa, b, d, and ω are available from Zhu [16]. After the initial value is determined, the equation is iterated and corrected numerically. Relevant simulation parameter was selected as β=0.035. As shown in Figure 3, the differential correction results of the halo orbit around the L2-type AEP can be presented, which is considered the nominal trajectory of the leader spacecraft.Figure 3
The quasi-halo orbit near L2-type AEP created by the differential corrector.


### 2.2. Spacecraft Relative Motion Dynamic Model

Having analyzed the generation of the AEP and the nominal trajectory of the leader spacecraft, we are now in a position to investigate the relative motion of spacecraft formation system. We purpose a relative motion scheme where charged spacecrafts (followers) move about a leader that produces an electric field. In addition, both spacecrafts are equipped with several solar sails, which can be controlled to produce solar radiation pressure. Moreover, the leader and followers all have the ability to modulate its electrostatic charge. Therefore, the spacecrafts are subject to the gravitational force of the two primary bodies (Sun and Earth), solar radiation pressure and the Coulomb force induced by each other. A relative dynamic model is established below, which considers the influences of the Coulomb force and solar radiation pressure.The leader-follower spacecraft formation is shown in Figure4. The position vectors of the leader with respect to the coordinate origin O and the two central bodies are rL=xL,yL,zLT, r1L=x1L,y1L,z1LT, and r2L=x2L,y2L,z2LT, respectively. The position vectors of the follower i from the coordinate origin O and the two primaries are rFi=xFi,yFi,zFiT, r1Fi=x1Fi,y1Fi,z1FiT, and r2Fi=x2Fi,y2Fi,z2FiT, respectively. To describe the relative motion of spacecraft in formation, the motion of leader and followers is firstly required to be expressed in the synodic reference frame. The dynamic equations for leader and follower in the hybrid propulsion formation are presented by
(14)r¨L=−Gm1r1L3r1L−Gm2r2L3r2L+asL+acL+wL,(15)r¨Fi=−Gm1r1Fi3r1Fi−Gm2r2Fi3r2Fi+asFi+acFi+wFi.Figure 4
Leader-follower spacecraft formation.HereacL and acFi are the acceleration of the leader and follower due to Coulomb force, respectively; see Eq. (16). asL and asFi are the acceleration of the leader and follower generated by solar sail, respectively, which can be calculated by (7), while wL and wFi represent external disturbances;
(16)acL=∑i=1,2,3−kc⋅QLFimLrLFi31+rLFiλdexp−rLFiλd⋅rLFiacFi=∑i≠j;j=13kc⋅QFiFjmFirFiFj31+rFiFjλdexp−rFiFjλd⋅rFiFj+kc⋅QLFjmFirLFi31+rLFiλdexp−rLFiλd⋅rLFi,where mL and mFi are the mass of leader and follower, QLFi and QFiFj represent charge product between the spacecrafts, kc=8.988×109N⋅m2/C2 represents the electrostatic constant, rLFi is the distance between leader and follower, and rFiFj is the distance between follower i and follower j. The Debye length is λd=24m in this paper.The relative position vector of theith spacecraft can be denoted as
(17)ρi=rFi−rL=xi,yi,ziT.The dimensionless relative equation of motion of the followeri can be written as
(18)x¨i=2y.i+xi+1−μxL+μr1L3−xL+xi+μr1L+rLFi3+μxL−1+μr2L3−xL+xi−1+μr2L+rLFi3+asxi+acxi+wxiy¨i=−2x.i+yi+1−μyLr1L3−yL+yir1L+rLFi3+μyLr2L3−yL+yir2L+rLFi3+asyi+acyi+wyiz¨i=1−μzLr1L3−zL+zir1L+rLFi3+μzLr2L3−zL+zir2L+rLFi3+aszi+aczi+wzi,where asi=asFi−asL=asxi,asyi,asziT represents the difference of solar radiation pressure between the follower i and leader, aci=acFi−acL=acxi,acyi,acziT denote the relative acceleration is produced by the Coulomb force, while wi=wFi−wL=wxi,wyi,wziT represent external disturbances; and the external disturbance models are introduced below.


### 2.3. External Disturbances

When the formation performs its mission in deep space environment, it will suffer from various perturbations, such as the effects arise from the eccentric nature of the Earth’s orbit, the perturbations of the moon and other planets. The most significant disturbances are the effects due to the eccentricity of Earth’s orbit and Lunar gravitation, denoted aswe and wm, respectively. Note that the effects that arise from the moon have been included in the CR3BP.Let us assume that the effects that arise from the eccentric nature of Earth’s orbit are mainly in thex direction. The non-dimensional expression of disturbances can be presented by
(19)we=μ1rc2−1re2,0,0T,where rc is the distance between the center of the earth’s circular orbit and the L2-Type equilibrium point. re is distance between the equilibrium point and Earth at the pericenter of its elliptical orbit.Meanwhile, the unknown uncertain interference items on the spacecraft system during the exploration mission mainly include the unknown effects of plasma in solar wind particles on hybrid propulsion system and the vibration of spacecrafts during acceleration, which vary with the state of the formation system. And the magnitude of these disturbances is small, which can be known as system uncertainties.All these disturbances are collectively referred to as external disturbances. It must be emphasized that external disturbances are only used to demonstrate the robustness of the controller, which is unnecessary to model all the disturbances precisely in this note. To better control formation flight while overcoming external disturbances, the effective control algorithms are designed below.


## 2.1. Generation of the Artificial Equilibrium Point and Nominal Trajectory of the Leader

The base case of the CRTBP used in this paper is described in Figure1, which can be assumed that two celestial bodies (Sun and Earth) rotate around their center of gravity in circular orbits with a constant angular velocity. The mass of the spacecraft can be ignored with respect to the mass of the two celestial bodies. Therefore, the effects of the third body spacecraft on the two central bodies are negligible. Three collinear libration points L1, L2, and L3 are also shown in Figure 1. For convenience, the reference frames used in the paper are given below.Figure 1
Relevant coordinate systems.RI - Sun-Earth barycenter synodic reference frame: the origin of this reference is located at the Sun-Earth’s center-of-mass, with plane XY coinciding to the ecliptic. The X unit vector is pointed from the Sun towards the Earth, the Z-axis is normal to the plane of the celestial bodies’ orbit, and the Y-axis completes the right hand orthogonal frame system and is thus defined as normal to the X-axis in the plane of the celestial bodies’ orbit and along the prograde rotational direction.RB - body-fixed reference frame located at the spacecraft’s center of mass: with origin at the spacecraft center of mass o, which also coincides with the solar sail center of mass. The unit vectors i, j, and k are assumed to lie along the principal axes of inertia of the (ideal) unwrinkled sail.In order to facilitate the subsequent calculation, the system parameters are normalized. The total system’s mass is defined as:(1)m=m1+m2,where m1 and m2 represent the masses of the celestial bodies (Sun and Earth), while μ is the mass ratio of the planetary system.
(2)μ=m2m1+m2.To make it easier to normalize variables, the following assumptions need to be made:(1)
The total mass of the celestial bodies is normalized to one(2)
The distance between the celestial bodies is normalized to one(3)
The angular velocity of the central bodies rotating around their center of mass is normalized to oneAccordingly, the period of one Earth–Sun revolution is2π, the Sun-[Earth+Moon] reference distance is rSE=1AU, the dimensionless mass of the [Earth+Moon] system is μ=m2/m1+m2=3.05425×10−6, the O-Sun distance is rS=μrSE=μ, the O-[Earth+Moon] distance is rE=1−μrSE=1−μ, and the dimensionless mass of the Sun can be given as m1=1−μm1+m2=1−μ.The propulsive acceleration vectoras of non-ideal solar sail can be presented as:
(3)as=β1−μrSs2rSs⋅nb1rSs+b2rSs⋅n+b3n,where rSs is the Sun-solar sail vector, β is referred to as the light number, which defines the (dimensionless) ratio of the maximum propulsive acceleration magnitude to the local Sun’s gravitational acceleration, which is taking values between [0,1]. n is the unit normal vector of solar sail. The terms {b1, b2, b3} are the sail force coefficients [39] defined as
(4)b1=121+sρ,b2=12ρ1−sBf−ρ−1εfBf−εbBbεb+εf,b3=121−sρ.The values of the optical parameters of a non-ideal solar sail membrane, involved in the calculation of {b1, b2, b3}, can be obtained from the experimental data collected during the NASA scout mission [40].
(5)ρ=0.91,s=0.89,εf=0.0025,εb=0.27,Bf=0.79,Bb=0.67.The sail normal unit vectorn can be expressed as a function of two attitude angles θ,φ, written as
(6)n=cosθcosφ,cosθsinφ,sinθT,where θ∈−π/2,π/2 rad is the angle between the direction of n and the i,j plane, while φ∈−π/2,π/2 rad is the angle between the direction of i and the projection of n onto i,j, as shown in Figure 1.Consider wrinkled and imperfect solar sails, which move within the Sun-[Earth+Moon] system. Combining equations (3) and (6), the dimensionless propulsive acceleration vector as of the non-ideal solar sail becomes
(7)asx=β1−μrSs2cosθcosφb1+b2cos2θcos2φ+b3cosθcosφasy=β1−μrSs2cosθcosφcosθsinφb2cosθcosφ+b3asz=β1−μrSs2cosθcosφ−sinθb2cosθcosφ+b3.The equations of orbital motion of solar sail are formulated based on Newton’s law of motion and law of universal gravitation:(8)r¨=−Gm1rSs3rSs−Gm2rEs3rEs+as.Therefore, the non-dimensional equations of motion is(9)x¨=2y.+x−1−μx+μrSs3−μx−1+μrEs3+asxy¨=−2x.+y−1−μyrSs3−μyrEs3+asyz¨=−1−μzrSs3−μzrEs3+asz,where rEs is the Earth-solar sail vector and r=xyzT is the position vector of the solar sail. Due to the existence of solar sail, the balance of forces on spacecraft in space has changed. When the solar radiation pressure balances both the celestial body gravitational pull and the centrifugal force along the Sun-[Earth+Moon] line, a new family of collinear AEPs is generated. Compared with classical equilibrium points, artificial equilibrium points are closer to the sun, which enables spacecraft to perform better in missions, such as observing the Sun and relaying deep space communications.Let us assume that the sail membrane is perpendicular to the direction of solar radiation, according to the definition of equilibrium point(10)x¨=y¨=z¨=0x.=y.=z.=0y=z=0θ=φ=0x=xAEPX^.Substituting Eq. (10) into Eq. (9), the result is a single scalar equation
(11)xAEP−1−μ1−βb1+b2+b3xAEP+μ2−μxAEP+μ−12=0.The solution of Eq. (11) gives the location of the L2-type AEP along the X-axis. Increasing the solar sail lightness number has the effect of shifting the position of L2 towards the Sun. Hence, the required non-ideal sail lightness number β necessary to maintain an L2-type AEP is
(12)β=xAEP+μ21−μ/xAEP+μ2+μ/xAEP+μ−12−xAEP1−μb1+b2+b3.Figure2 shows the relationship curve between the position of the artificial equilibrium point and sail lightness number.Figure 2
Non-ideal solar sail lightness numberβ as a function of dimensionless Sun-AEP distance xAEP.As shown in Figure2, the larger β is, the collinear artificial equilibrium point of solar sail is closer to the central celestial body (Sun). Based on this situation, the solar sail can achieve more accurate observation of the Sun and deep space communication missions by being placed on a non-Keplerian orbit around the AEP, which can hardly be realized by conventional spacecraft.For the purpose of generating quasi-halo orbits around the artificial equilibrium points, the differential corrections methodology is adopted in this paper, which is generally used to calculate unstable long-periodic orbits [41], such as displaced non-Keplerian orbits and J2 invariant relative orbits. For detailed research on differential corrections algorithm, please refer to the works of Howell [41] and Ross [42]. The differential correction algorithm is not self-starting and requires an initial approximate solution. Richardson expansion equation obtained according to Lindstedt-Poincare method [43, 44] can be used as an effective initial value:
(13)x=a21Ax2+a22Az2−Axcosωt+a23Ax2−a24Az2cos2ωt+a31Az3−a32AxAz2cos3ωt,y=Ay+b33Ax3+b34−b35AxAz2sinωt+b21Ax2−b22Az2sin2ωt+b31Ax3−b32AxAz2sin3ωt,z=−3d21AxAz+Azcosωt+d21AxAzcos2ωt+d32AzAx2−d31Az3cos3ωt,where the values of the parametersa, b, d, and ω are available from Zhu [16]. After the initial value is determined, the equation is iterated and corrected numerically. Relevant simulation parameter was selected as β=0.035. As shown in Figure 3, the differential correction results of the halo orbit around the L2-type AEP can be presented, which is considered the nominal trajectory of the leader spacecraft.Figure 3
The quasi-halo orbit near L2-type AEP created by the differential corrector.


## 2.2. Spacecraft Relative Motion Dynamic Model

Having analyzed the generation of the AEP and the nominal trajectory of the leader spacecraft, we are now in a position to investigate the relative motion of spacecraft formation system. We purpose a relative motion scheme where charged spacecrafts (followers) move about a leader that produces an electric field. In addition, both spacecrafts are equipped with several solar sails, which can be controlled to produce solar radiation pressure. Moreover, the leader and followers all have the ability to modulate its electrostatic charge. Therefore, the spacecrafts are subject to the gravitational force of the two primary bodies (Sun and Earth), solar radiation pressure and the Coulomb force induced by each other. A relative dynamic model is established below, which considers the influences of the Coulomb force and solar radiation pressure.The leader-follower spacecraft formation is shown in Figure4. The position vectors of the leader with respect to the coordinate origin O and the two central bodies are rL=xL,yL,zLT, r1L=x1L,y1L,z1LT, and r2L=x2L,y2L,z2LT, respectively. The position vectors of the follower i from the coordinate origin O and the two primaries are rFi=xFi,yFi,zFiT, r1Fi=x1Fi,y1Fi,z1FiT, and r2Fi=x2Fi,y2Fi,z2FiT, respectively. To describe the relative motion of spacecraft in formation, the motion of leader and followers is firstly required to be expressed in the synodic reference frame. The dynamic equations for leader and follower in the hybrid propulsion formation are presented by
(14)r¨L=−Gm1r1L3r1L−Gm2r2L3r2L+asL+acL+wL,(15)r¨Fi=−Gm1r1Fi3r1Fi−Gm2r2Fi3r2Fi+asFi+acFi+wFi.Figure 4
Leader-follower spacecraft formation.HereacL and acFi are the acceleration of the leader and follower due to Coulomb force, respectively; see Eq. (16). asL and asFi are the acceleration of the leader and follower generated by solar sail, respectively, which can be calculated by (7), while wL and wFi represent external disturbances;
(16)acL=∑i=1,2,3−kc⋅QLFimLrLFi31+rLFiλdexp−rLFiλd⋅rLFiacFi=∑i≠j;j=13kc⋅QFiFjmFirFiFj31+rFiFjλdexp−rFiFjλd⋅rFiFj+kc⋅QLFjmFirLFi31+rLFiλdexp−rLFiλd⋅rLFi,where mL and mFi are the mass of leader and follower, QLFi and QFiFj represent charge product between the spacecrafts, kc=8.988×109N⋅m2/C2 represents the electrostatic constant, rLFi is the distance between leader and follower, and rFiFj is the distance between follower i and follower j. The Debye length is λd=24m in this paper.The relative position vector of theith spacecraft can be denoted as
(17)ρi=rFi−rL=xi,yi,ziT.The dimensionless relative equation of motion of the followeri can be written as
(18)x¨i=2y.i+xi+1−μxL+μr1L3−xL+xi+μr1L+rLFi3+μxL−1+μr2L3−xL+xi−1+μr2L+rLFi3+asxi+acxi+wxiy¨i=−2x.i+yi+1−μyLr1L3−yL+yir1L+rLFi3+μyLr2L3−yL+yir2L+rLFi3+asyi+acyi+wyiz¨i=1−μzLr1L3−zL+zir1L+rLFi3+μzLr2L3−zL+zir2L+rLFi3+aszi+aczi+wzi,where asi=asFi−asL=asxi,asyi,asziT represents the difference of solar radiation pressure between the follower i and leader, aci=acFi−acL=acxi,acyi,acziT denote the relative acceleration is produced by the Coulomb force, while wi=wFi−wL=wxi,wyi,wziT represent external disturbances; and the external disturbance models are introduced below.


## 2.3. External Disturbances

When the formation performs its mission in deep space environment, it will suffer from various perturbations, such as the effects arise from the eccentric nature of the Earth’s orbit, the perturbations of the moon and other planets. The most significant disturbances are the effects due to the eccentricity of Earth’s orbit and Lunar gravitation, denoted aswe and wm, respectively. Note that the effects that arise from the moon have been included in the CR3BP.Let us assume that the effects that arise from the eccentric nature of Earth’s orbit are mainly in thex direction. The non-dimensional expression of disturbances can be presented by
(19)we=μ1rc2−1re2,0,0T,where rc is the distance between the center of the earth’s circular orbit and the L2-Type equilibrium point. re is distance between the equilibrium point and Earth at the pericenter of its elliptical orbit.Meanwhile, the unknown uncertain interference items on the spacecraft system during the exploration mission mainly include the unknown effects of plasma in solar wind particles on hybrid propulsion system and the vibration of spacecrafts during acceleration, which vary with the state of the formation system. And the magnitude of these disturbances is small, which can be known as system uncertainties.All these disturbances are collectively referred to as external disturbances. It must be emphasized that external disturbances are only used to demonstrate the robustness of the controller, which is unnecessary to model all the disturbances precisely in this note. To better control formation flight while overcoming external disturbances, the effective control algorithms are designed below.


## 3. Orbital Maneuver Control Strategy

In this section, first, we design a modified fast integral terminal sliding mode control law with a new double power combination function reaching law for formation tracking of the desired relative trajectories. This control law combines the advantages of both the reduced-order disturbance observer and sliding mode control: the reduced-order disturbance observer is able to estimate the unknown disturbances, and the novel FITSMC with a new double power combination function reaching law can guarantee the tracking errors of the spacecraft formation system convergence rapidly and improving the robustness and disturbance rejection capability. Then, the basic PID control law and fast terminal sliding mode control law are proposed to compare with the MFITSMC method. Finally, based on the control input calculated, a novel optimal allocation strategy to determine the control commands is investigated.


### 3.1. The MFITSMC with Reduced-Order Disturbance Observer

Definex1=ρi, x2=ρ.i as the state vector. xr=ρr is being the desired trajectory of relative position, then the state-space equation can be obtained as follows:
(20)x.1=x2x.2=fx1,x2+ui+wi,where(21)ui=aci+asi,(22)fx1,x2,wi=2y.i+xi+1−μxL+μr1L3−xL+xi+μr1L+rLFi3+μxL−1+μr2L3−xL+xi−1+μr2L+rLFi3−2x.i+yi+1−μyLr1L3−yL+yir1L+rLFi3+μyLr2L3−yL+yir2L+rLFi31−μzLr1L3−zL+zir1L+rLFi3+μzLr2L3−zL+zir2L+rLFi3.The aim was to drive the outputx1 to track its desired trajectory xr, and the tracking error can be written as e=x1−xr. From Eq. (20) the following can be obtained:
(23)e¨=fx1,x2+wi+ui−x¨r.The fast integral terminal sliding mode surface is written as:(24)s=e.+αs1+ec1/c2s.1=e.+ηe.m/ns=sx,sy,szT.Whereα>0, η>0 are positive constants. c1, c2, m, and n are odd integers and meet 1<c1/c2<2 and c1/c2>m/n.A novel double power combination function reaching law with the capability of disturbance rejection and rapid convergence is given by(25)s.=−k1fals,a,δ−k2sbsgns.Herea=1+γ, b=1−γ, δ=1, 0<γ<1, and the nonlinear exponential composite function fal⋅ can be given by
(26)fals,a,δ=sasgns,sδ1−a,s>δs≤δ.Consider the spacecraft formation flying system given in Eq. (20) and the proposed fast integral terminal sliding mode surface, the control law is designed as
(27)ui=−k1fals,a,δ−k2sbsgns−αe.−αηe.m/n−c1c2ec1/c2−1⋅e.−fx1,x2−w^i+x¨r.In view of the external perturbations in spacecraft formation system, the reduced-order disturbance observer can be designed below,(28)w^i=v+ξ1∫0te¨dτ−ξ2∫0tsignΘdτξ1>0,ξ2>ς,(29)v.=−ξ1w^i+ξ1−f−a+x¨r,where(30)Θ=w^i−e¨−f−a+x¨r.The disturbance vectorwi can be rewritten wi=e¨−f−a+x¨r. Taking the time derivative of equation (28) yields
(31)w^i=−ξ1w^+ξ1−fx1,x2−a+x¨r+e¨−ξ2signΘ.The disturbance estimation error can be denoted asw~i=w^i−wi, and we can obtain
(32)w^.i=−ξ1w~i−ξ2signΘ.In fact,Θ=w~i, combining with Eq. (31), the derivative of the disturbance estimation error w~i can be obtained as
(33)w~.i=w^.i−w.i=−ξ1w~i−ξ2signw~i−w.i.Select a Lyapunov function for stability analysis, as shown below(34)Vw=12w~iTw~i.Taking the time derivative of equation (34) yields:
(35)V.w=w~Tw~.=−ξ1w~2−ξ2w−w~Tw.≤−ξ1w~2−ξ2−ςw~,where ξ=ξ2−ς, select appropriate variable ξ>0, then Eq.(35) can be concluded as follows,
(36)V.w+2ξ1Vw+2ξVw1/2≤0.In conclusion, the disturbance estimation errorw~ of the formation system can tend towards zero in finite time tw.
(37)tw≤1ξ1ln2ξ1Vw1/20+2ξ2ξ,where Vw0 is the initial value of Vw.Next part is to demonstrate the relative position errors can tend towards 0 in finite time, consider another Lyapunov function,(38)V=12sTs.Taking the derivative of Equation (38), and then it can be given by:
(39)V.=sTs.=sT−k1fals,a,δ−k2sbsgns+w^i−wi≤−sk1fals,1+γ,δ−k2s2−γ+sw^i−wi≤s≤1−k2s2−γ−sk1s−D−k1s2−sk2s1−γ−Ds>1−k2s2−γ−sk1s1+γ−D−k1s2+γ−sk2s1−γ−D.When1≥s≥D/k1, according to Eq. (39), V.≤−k2s2−γ=−21−γ/2k2V1−γ/2≤0.1≥s≥D/k21/1−γ, according to Eq. (39), V.≤−k1s2=−2k1V≤0.s>D/k11/1+γ>1, according to Eq. (39), V.≤−k1s2−γ=−21−γ/2k2V1−γ/2≤0.s>D/k21/1‐γ>1, according to Eq. (39), V.≤−k1s2+γ=−21+γ/2k1V1+γ/2≤0.According to the Lyapunov stability theory, the closed-loop system has global asymptotic stability and can tend towards zero in finite time.Figure5 depicts the block diagram of the MFFTSMC controller with the disturbance observer. In order to test the performance of the proposed MFITSMC, two classical control methods are designed for the purpose of comparison, which are PID controller and fast terminal sliding mode controller (FTSMC). The control law of the system is defined as follows, respectively.Figure 5
Block diagram of the MFITSMC controller for the follower spacecraft.FTSMC:A classical FTSM manifold is chosen as:(40)s=e.+λe+βeq/p,where λ and β are positive constants, and p and q are positive odd integers satisfying p>q, while eq/p=exq/peyq/pezq/p.For the sake of suppressing the chattering phenomenon of sliding mode control, the reach law is designed as(41)s.=−k3s−k4tanhs.The fast terminal sliding mode controller can be designed as:(42)uFTSM=−k4tanhs−λk3⋅e−k3βeq/p−βγeq/p−1⋅e.−fx1,x2−wi+x¨r.The control law for the PID controller is given by(43)uPID=−kpe−ki∫0teτdτ−kd⋅e.−fx1,x2−wi+x¨r.


### 3.2. Control Commands Allocation

Here, the control inputui consisted of the acceleration due to Coulomb force and the acceleration due to solar sail. However, the control commands for the proposed propulsion system are the charge products of spacecrafts and the attitude angles of the solar sails. Therefore, it is needful to distribute the control input to obtain actual control commands. Since the mapping between control commands and control inputs is highly nonlinear and the number of control command variables is more than that of control inputs, the determination of control commands is a nonlinear programming problem.According to Eq. (16), ac is determined by leader-follower charge products QLFii=1,2,3 and follower-follower charge products QFiFji≠j. For mission lifetime extension, the sail attitude angles θ and φ should be optimized so that the magnitude of the acceleration required due to the Coulomb force is minimum. To avoid possible system instabilities and violent chattering in system, it is rational to limit the attitude angles within certain range, such as θ∈30,60 deg and φ∈30,60 deg. Since the charge and discharge of the Coulomb spacecraft needs to consume a fraction of energy, it is firstly distributed according to the principle of optimal energy. A particle swarm optimization (PSO) method is used here to solve the constrained optimization problem. Consequently, the constrained optimization problem can be described below:
(44)θ∗,φ∗=argminθ,φaci=argminθ,φui−asi,where the total required acceleration ui is given by Eq. (27). θ∗ and φ∗ represent the optimal solutions of solar sail attitude angles in the actual control input, respectively.As shown in equation (16), the magnitude of the coulomb acceleration input aci is determined by three different charge products. Take follower 1 for example, the three different charge products are QLF1, QF1F2, and QF1F3. Another allocation mechanism was required to be designed to determine the appropriate charge product of the spacecraft as the control command for the formation system. This control allocation can also be regarded as an optimization problem with nonlinear constraints. The objective is to minimize the difference between the expected acceleration input a^ci due to Coulomb force and the available acceleration aciQLF1,QF1F2,QF1F3. Hence, the optimization problem can be written below:
(45)minimizeJ=a^ci−aciQLF1,QF1F2,QF1F3.Subject to:(46)Qmin≤Qj≤Qmax,Qj−Q^j≤Qlim,j=1,2,3,(47)Q0=Q^,where Q^ are the current charge products. Q0 are the initial values of the charge product. Qmin and Qmax are the minimal and maximal available charge products, respectively. Qlim limit the variations of charge products between control intervals.


## 3.1. The MFITSMC with Reduced-Order Disturbance Observer

Definex1=ρi, x2=ρ.i as the state vector. xr=ρr is being the desired trajectory of relative position, then the state-space equation can be obtained as follows:
(20)x.1=x2x.2=fx1,x2+ui+wi,where(21)ui=aci+asi,(22)fx1,x2,wi=2y.i+xi+1−μxL+μr1L3−xL+xi+μr1L+rLFi3+μxL−1+μr2L3−xL+xi−1+μr2L+rLFi3−2x.i+yi+1−μyLr1L3−yL+yir1L+rLFi3+μyLr2L3−yL+yir2L+rLFi31−μzLr1L3−zL+zir1L+rLFi3+μzLr2L3−zL+zir2L+rLFi3.The aim was to drive the outputx1 to track its desired trajectory xr, and the tracking error can be written as e=x1−xr. From Eq. (20) the following can be obtained:
(23)e¨=fx1,x2+wi+ui−x¨r.The fast integral terminal sliding mode surface is written as:(24)s=e.+αs1+ec1/c2s.1=e.+ηe.m/ns=sx,sy,szT.Whereα>0, η>0 are positive constants. c1, c2, m, and n are odd integers and meet 1<c1/c2<2 and c1/c2>m/n.A novel double power combination function reaching law with the capability of disturbance rejection and rapid convergence is given by(25)s.=−k1fals,a,δ−k2sbsgns.Herea=1+γ, b=1−γ, δ=1, 0<γ<1, and the nonlinear exponential composite function fal⋅ can be given by
(26)fals,a,δ=sasgns,sδ1−a,s>δs≤δ.Consider the spacecraft formation flying system given in Eq. (20) and the proposed fast integral terminal sliding mode surface, the control law is designed as
(27)ui=−k1fals,a,δ−k2sbsgns−αe.−αηe.m/n−c1c2ec1/c2−1⋅e.−fx1,x2−w^i+x¨r.In view of the external perturbations in spacecraft formation system, the reduced-order disturbance observer can be designed below,(28)w^i=v+ξ1∫0te¨dτ−ξ2∫0tsignΘdτξ1>0,ξ2>ς,(29)v.=−ξ1w^i+ξ1−f−a+x¨r,where(30)Θ=w^i−e¨−f−a+x¨r.The disturbance vectorwi can be rewritten wi=e¨−f−a+x¨r. Taking the time derivative of equation (28) yields
(31)w^i=−ξ1w^+ξ1−fx1,x2−a+x¨r+e¨−ξ2signΘ.The disturbance estimation error can be denoted asw~i=w^i−wi, and we can obtain
(32)w^.i=−ξ1w~i−ξ2signΘ.In fact,Θ=w~i, combining with Eq. (31), the derivative of the disturbance estimation error w~i can be obtained as
(33)w~.i=w^.i−w.i=−ξ1w~i−ξ2signw~i−w.i.Select a Lyapunov function for stability analysis, as shown below(34)Vw=12w~iTw~i.Taking the time derivative of equation (34) yields:
(35)V.w=w~Tw~.=−ξ1w~2−ξ2w−w~Tw.≤−ξ1w~2−ξ2−ςw~,where ξ=ξ2−ς, select appropriate variable ξ>0, then Eq.(35) can be concluded as follows,
(36)V.w+2ξ1Vw+2ξVw1/2≤0.In conclusion, the disturbance estimation errorw~ of the formation system can tend towards zero in finite time tw.
(37)tw≤1ξ1ln2ξ1Vw1/20+2ξ2ξ,where Vw0 is the initial value of Vw.Next part is to demonstrate the relative position errors can tend towards 0 in finite time, consider another Lyapunov function,(38)V=12sTs.Taking the derivative of Equation (38), and then it can be given by:
(39)V.=sTs.=sT−k1fals,a,δ−k2sbsgns+w^i−wi≤−sk1fals,1+γ,δ−k2s2−γ+sw^i−wi≤s≤1−k2s2−γ−sk1s−D−k1s2−sk2s1−γ−Ds>1−k2s2−γ−sk1s1+γ−D−k1s2+γ−sk2s1−γ−D.When1≥s≥D/k1, according to Eq. (39), V.≤−k2s2−γ=−21−γ/2k2V1−γ/2≤0.1≥s≥D/k21/1−γ, according to Eq. (39), V.≤−k1s2=−2k1V≤0.s>D/k11/1+γ>1, according to Eq. (39), V.≤−k1s2−γ=−21−γ/2k2V1−γ/2≤0.s>D/k21/1‐γ>1, according to Eq. (39), V.≤−k1s2+γ=−21+γ/2k1V1+γ/2≤0.According to the Lyapunov stability theory, the closed-loop system has global asymptotic stability and can tend towards zero in finite time.Figure5 depicts the block diagram of the MFFTSMC controller with the disturbance observer. In order to test the performance of the proposed MFITSMC, two classical control methods are designed for the purpose of comparison, which are PID controller and fast terminal sliding mode controller (FTSMC). The control law of the system is defined as follows, respectively.Figure 5
Block diagram of the MFITSMC controller for the follower spacecraft.FTSMC:A classical FTSM manifold is chosen as:(40)s=e.+λe+βeq/p,where λ and β are positive constants, and p and q are positive odd integers satisfying p>q, while eq/p=exq/peyq/pezq/p.For the sake of suppressing the chattering phenomenon of sliding mode control, the reach law is designed as(41)s.=−k3s−k4tanhs.The fast terminal sliding mode controller can be designed as:(42)uFTSM=−k4tanhs−λk3⋅e−k3βeq/p−βγeq/p−1⋅e.−fx1,x2−wi+x¨r.The control law for the PID controller is given by(43)uPID=−kpe−ki∫0teτdτ−kd⋅e.−fx1,x2−wi+x¨r.


## 3.2. Control Commands Allocation

Here, the control inputui consisted of the acceleration due to Coulomb force and the acceleration due to solar sail. However, the control commands for the proposed propulsion system are the charge products of spacecrafts and the attitude angles of the solar sails. Therefore, it is needful to distribute the control input to obtain actual control commands. Since the mapping between control commands and control inputs is highly nonlinear and the number of control command variables is more than that of control inputs, the determination of control commands is a nonlinear programming problem.According to Eq. (16), ac is determined by leader-follower charge products QLFii=1,2,3 and follower-follower charge products QFiFji≠j. For mission lifetime extension, the sail attitude angles θ and φ should be optimized so that the magnitude of the acceleration required due to the Coulomb force is minimum. To avoid possible system instabilities and violent chattering in system, it is rational to limit the attitude angles within certain range, such as θ∈30,60 deg and φ∈30,60 deg. Since the charge and discharge of the Coulomb spacecraft needs to consume a fraction of energy, it is firstly distributed according to the principle of optimal energy. A particle swarm optimization (PSO) method is used here to solve the constrained optimization problem. Consequently, the constrained optimization problem can be described below:
(44)θ∗,φ∗=argminθ,φaci=argminθ,φui−asi,where the total required acceleration ui is given by Eq. (27). θ∗ and φ∗ represent the optimal solutions of solar sail attitude angles in the actual control input, respectively.As shown in equation (16), the magnitude of the coulomb acceleration input aci is determined by three different charge products. Take follower 1 for example, the three different charge products are QLF1, QF1F2, and QF1F3. Another allocation mechanism was required to be designed to determine the appropriate charge product of the spacecraft as the control command for the formation system. This control allocation can also be regarded as an optimization problem with nonlinear constraints. The objective is to minimize the difference between the expected acceleration input a^ci due to Coulomb force and the available acceleration aciQLF1,QF1F2,QF1F3. Hence, the optimization problem can be written below:
(45)minimizeJ=a^ci−aciQLF1,QF1F2,QF1F3.Subject to:(46)Qmin≤Qj≤Qmax,Qj−Q^j≤Qlim,j=1,2,3,(47)Q0=Q^,where Q^ are the current charge products. Q0 are the initial values of the charge product. Qmin and Qmax are the minimal and maximal available charge products, respectively. Qlim limit the variations of charge products between control intervals.


## 4. Numerical Simulations

To testify the effectiveness of the proposed control algorithm, the nonlinear relative motion equations defined by Eq. (18) are numerically simulated along with the perturbative force described in Eq. (19) and the control inputs defined by Eq. (27). In this note, we consider a hybrid propulsion formation system that consists of a leader and three followers. To better observe the leader-follower spacecraft formation, consider the leader and followers maintaining the required distance in the formation to form a specific shape. In particular, the spacecrafts in the formation are assumed to have a mass mL=mF1=mF2=mF3=100kg. A periodic Halo orbit, with a vertical displacement of 300,000 km around the L2-type AEP, is considered nominal trajectory. The leader revolves on the nominal orbit, while the follower spacecrafts trajectories follow a projected circular orbit with a radius of 40 m relative to the leader, as shown in Figure 1. Parameters in the proposed orbital controllers are given in Table 1, and the simulation time is considered 1000 days.Table 1
Controller parameters for numerical simulation.
ControllerParameterMFITSMCk1=diag4,4,4, k2=diag1,1,1, γ=0.5, δ=1, α=diag0.1,0.1,0.1, η=diag0.01,0.01,0.01, c1/c2=25/13, m/n=19/13, ξ1=diag25,25,25, ξ2=diag0.5,0.5,0.5FTSMCλ=diag1,1,1×102,k3=diag40,40,40, k4=diag0.01,0.01,0.01,β=diag0.04,0.04,0.04, q=3, p=5PIDkp=diag0.00612,0.00612,0.00612, ki=diag0.15,0.15,0.15×10−3, kd=diag3,3,3It is assumpted that the three followers are in the same position as the leader at the initial time. The initial orbit parameters are given in Table2. Note that the lightness number of the sail is chosen as 0.035. In the chief rotating reference frame Rc, the parametric solution of the i-th follower trajectory (relative to the leader) is assumed to have the following algebraic form.
(48)ρr=0,40cosnrt+ϕi,40sinnrt+ϕiTm,where nr is the angular velocity of the relative motion, and ϕi=2i−1π/3 is the phase angle.Table 2
Initial position and velocity of the followers in the rotation coordinate system.
x0my0mz0mx.0m/sy.0m/sz.0m/sF10000.5×10−5m/s−3.5×10−5m/s1.5×10−5m/sF2000−0.3×10−5m/s2.5×10−5m/s−0.8×10−5m/sF30000.2×10−5m/s−3×10−5m/s−1.5×10−5m/sAt present, the available maximal specific charge is about 0.03 C/kg [27], and the minimal and maximal charge products are chosen as Qmin=−9C2 and Qmax=9C2. Qlim are set as 1×10‐5C2. To prove the controllers’ disturbance rejection capacity, an unmodeled disturbance w¯i is set up during the simulation:
(49)w¯i=w¯icos2ωsub−L2t,sinωsub−L2t,cosωsub−L2tsinωsub−L2tT,where w¯i=1×10−4AU/TU2, absolutely larger than we, which is introduced in Section 2.3. w¯i is large enough to cover all these uncertain disturbances. ωsub−L2 is the frequency of halo orbit. In the normalized system of units, ωsub−L2=1.2023.Some simulation parameters are as follows:The numerical simulation results are as follows:The phase portrait of the trajectoryρi of the followers relative to the leader is presented in Figures 6–8. It demonstrates the desired orbit tracking can be achieved with the controllers using the hybrid propulsion system. The solar sails are primarily used to change the orientation of formation system, which makes the whole formation system move around a halo orbit. Since Coulomb force belongs to the internal force of formation, it is mainly responsible for changing formation geometry and size. To analyze the simulation results more conveniently, only the results of first 1000 days are presented below. From Figures 9(a) and 9(b), it can be seen from the simulation that it takes nearly 50 days and 35 days for the PID and FTSMC controllers to drive the relative position errors gradually approach the prescribed state as time progresses. In contrast, it takes less than 20 days for the MFITSMC controller to reach steady state, as shown in Figure 9(c). The order of magnitude of position errors at steady state using the MFITSMC controller is 0.01 m, which is superior to the other two controllers. The reason is that the unmodeled disturbance is estimated and compensated, and the steady-state errors remain at a low level. In the velocity error responses of the MFITSMC with a new double power combination function reaching law, the nonlinear PID control law and FTSMC are presented in Figures 10(a)–10(c), respectively. It can also be seen that the relative velocity errors of formation system tend towards zero and fluctuate within little range, while the convergence rate of the proposed control algorithm is faster than the others, whose steady-state errors of the relative velocity are smaller than that of other control algorithms. As the followers move near the leader, the control input has a small order of magnitude, as shown in Figures 11(a)–11(c). We can see that the thrusts of each controller at reconfiguration and orbital maintenance stage are at the mN level, which is achievable by the hybrid propulsion system. Figures 12 and 13 show the time histories of the control commands for the follower, including the charge product and attitude angle. In the early stage of simulation, the value of charge products are positive, and the mutually exclusive Coulomb force and solar radiation pressure contribute to the separation of spacecrafts together. After the reconfiguration task is completed, the charge products keep fluctuating near the value of 0, which is responsible for maintaining the stability of formation configuration. As shown in Figure 13, solar sails participated in the spacecraft reconfiguration process. After that, the angle control variables change slowly within little range, and the variations are closely related to the inherent periods of the halo orbit.Figure 6
The trajectories of the followers relative to leader (PID).Figure 7
The trajectories of the followers relative to leader (FTSMC).Figure 8
The trajectories of the followers relative to leader (MFITSMC).
(a)(b)(c)(d)Figure 9
Components of the position error (Follower 1) ((a) PID; (b) FTSMC; and (c) MFITSMC).
(a)(b)(c)Figure 10
Components of the velocity error (Follower 1) ((a) PID; (b) FTSMC; and (c) MFITSMC).
(a)(b)(c)Figure 11
Time histories of the control input (Follower 1) ((a) PID; (b) FTSMC; and (c) MFITSMC).
(a)(b)(c)Figure 12
The time histories of charge products (Follower 1) ((a) Leader-Follower 1; (b) Follower 1-Follower 3; and (c) Follower 1-Follower 2).
(a)(b)(c)Figure 13
The time histories of the attitude angles (Follower 1).It can be seen that all the controllers can make all errors converge towards 0 as time progresses, which proves the effectiveness of the propellantless propulsion strategy. The control performances of the PID controller, FTSMC, and the MFITSMC are shown in Table3. The results show that there are some differences in the average variation of control commands among the three controllers. Since the observer continuously estimates the first derivative of the relative position, the MFITSMC requires much more control cost in terms of angle variations of the solar sails and average charge products. Meanwhile, the MFITSMC possesses much smaller state errors and the faster convergence rate in the presence of external disturbances, which illustrates the superiority of the MFITSMC controller.Table 3
The control performance of control algorithms (take follower 1 as an example).
ControllerTerminal relative position errors max (x,y,z) [m]Terminal relative velocity errors max (x,y,z)[m/s]Average variation of angle [deg]Average variation of charge product [C2]MFITSMC(0.012,0.018,0.022)0.21,1.46,1.24×10−7(0.0285,0.0101)5.754885480953×10‐12FTSMC(0.15,0.118,0.038)4.32,1.81,2.44×10−6(0.0244,0.0078)5.1569581477773×10‐12PID(0.95,0.78,0.42)1.42,1.15,0.39×10‐5(0.0084,0.0042)1.366201555087×10‐12Aforementioned simulation results verify the superiority and efficiency of the proposed MFITSMC algorithm with a new double power combination function reaching law, which demonstrates enhanced performances on the convergence rate, steady-state errors, percent overshoot, and the capacity of disturbance rejection.


## 5. Conclusion

The present work proposed a methodology to perform spacecraft formation orbit control around an L2-type artificial equilibrium point using the hybrid propulsion strategy due to Coulomb force and non-ideal solar sail. Coulomb force is produced by the modulation of the leader and follower spacecraft’s electric charge, and solar radiation pressure accelerations are generated by the rotations of the solar sail’s plates attached to the satellites. To handle the control problem of the formation reconfiguration against external disturbances, a sliding mode observer is introduced to estimate and compensate for the perturbations in space environment and system uncertainties, and an observer-based modified fast integral terminal sliding mode controller with a new double power combination function reaching law is designed. The proposed optimal control allocation scheme can determine control commands properly, which include the charge products of spacecrafts and attitude angles of solar sails. The asymptotic stability of closed-loop system is proved by using Lyapunov stability theorem. Numerical simulation results illustrated the validity and superiority of the proposed controller in the presence of external disturbances.It is obvious that the non-ideal solar sail model is still hypothetical in this note because the flexible characteristics of the solar sail, the structural characteristics brought by the large size structure, and the sail degradation mechanism are not considered. In the meantime, the shielding effect of other celestial bodies or objects on the sail surface is ignored. As such, more work should be done to obtain more useful data to perfect the solar sail model.Future work will focus on a detailed discussion on the process of using small solar arrays to charge the spacecraft through energy conversion and using the electron gun on board to release the unnecessary charge. For different spacecraft, building up and modulating the electric charge are complex processes, which needs further research.



---

*Source: 8719645-2022-05-02.xml*