Lyapunov-Based Control via Atmospheric Drag for Tetrahedral Satellite Formation

Abstract

The problem of small tetrahedral satellite formation maintenance in a Low Earth Orbit is being considered. The main purpose is to develop a simple algorithm for tetrahedron control via atmospheric drag. To design a controller, the direct Lyapunov method is used. The control obtained is suitable for tetrahedral formation maintenance, with an average distance of about 1 km. During the controlled motion, the geometric characteristics of the tetrahedron are preserved.

1. Introduction

The concept of four satellites orbiting around the Earth while forming a tetrahedron is useful for studying the dynamic characteristics of the Earth’s electromagnetic field. It is evident that the minimum number of satellites required to obtain continuous measurements that make it possible to calculate the spatial dynamic characteristics of the electromagnetic field is four [1]. As a result, the development and launch of a mission to study, for example, the geomagnetic field, the density of ionized particles in the Earth’s atmosphere, etc., inevitably face the need to deploy, build, maintain, and manage a tetrahedral formation of satellites. The relevant research on tetrahedral satellite formations, including the tasks of constructing, maintaining, and controlling their motion, can be roughly divided into three groups: publications related to the Cluster/Cluster II mission, publications related to the Magnetospheric MultiScale (MMS) mission, and other publications.

The Cluster mission (ESA project), was proposed at the end of 1982 and launched (as Cluster II) in 2000 [2,3]. In this mission, four satellites orbit the Earth in highly elliptical polar orbits (perigee to apogee: 20–120 thousand km). The main goal of the mission is a detailed study of the plasma of the Earth’s magnetosphere, in particular, study of the solar wind and the formation of the bow shock wave; study of the magnetopause; study of polar cusps, auroras, etc. [4,5,6]. The main analytical and numerical results obtained during the preparation of the mission, along with the developed methods, are presented in report [7].

The concept of the MMS mission (NASA project) was initially outlined in [8]; a review of the main scientific problems to be solved within the framework of this mission is presented in [9]. The mission was launched in 2015, and the MMS mission working group released a series of works on the construction of a tetrahedral formation suitable for the implementation of this mission. The general information about the construction of a tetrahedron is presented in [10,11]; more sophisticated information about the modeling of the tetrahedral configuration and the requirements for relative orbits is obtainable in [12,13].

Other works on tetrahedral configurations include the construction of tetrahedrons that span the Earth [14], and the study of signal delay properties for global positioning satellites [15]. The maintenance of a tetrahedral formation using a solar sail is considered in papers [16,17]. The problem of launching and maintaining a tetrahedral formation of a small spacecraft in a highly elliptical orbit was solved using semi-analytical and numerical methods in [18]. A publication related to The Auroral Lites mission deserves special mention [19].

Quite a lot of missions (perhaps using formations, but not tetrahedral ones) have been developed and implemented to study the ionosphere of the Earth. Among them are AEROS [20], CHAMP [21], Ørsted [22], Swarm [23], Chibis-M [24], Demeter [25], ICON [26]. The DICE mission showed the capability of using CubeSat-type nanosatellites to study the ionosphere [27]. Moreover, the study aiming at the usage of nanosatellites for studying the ionosphere can be found in [28]. The complexity of atmospheric research lies in the fact that the atmosphere is a non-stationary environment; the processes occurring in it depend on the state of the Earth’s magnetosphere, climate, atmospheric movement, time of year, solar activity, and even anthropogenic factors [29]. The study of such a heterogeneous environment allows us to better understand the physics of phenomena occurring in the atmosphere, as well as evaluate and predict the influence of various natural and man-made factors on each other. In addition, the heterogeneity of the ionosphere reduces the accuracy of satellite navigation using GPS constellations and causes a delay in signals when providing communications [30]. This implies the need to take into account the influence of the ionosphere and make amendments to theoretical and empirical models of radio signal propagation in the ionosphere [31].

This research is a continuation of the authors’ work presented in [32]. Here, we briefly present the relevant results of the previous work for the convenience of the reader. For additional information, we kindly ask the reader to refer to previous articles [32,33]. In these papers, we found special relative motions that ensure the stationary quality and volume of the tetrahedron. Note that this result was obtained for a linear model of motion. Hence, when we take into account non-linear terms and additional disturbances (e.g., J2 harmonics), the quality of the tetrahedron rapidly degrades. Therefore, we need some kind of active relative motion control. Since we consider the mission based on small satellites, the utilization of thrusters is hardly possible, and fuel-less approaches are desirable. In Low Earth Orbits (LEOs), we can use atmospheric drag for relative motion control. Control law is obtained using the well-known Lyapunov-based approach, which allows us to ensure the asymptotic stability of the required motion [34,35].

The remaining part of the paper is organized as follows. In Section 2, we briefly highlight the main results of our previous papers and provide a description of the reference tetrahedral formation. In Section 3, the relative motion equations and the effect of atmospheric drag are discussed. Section 4 is dedicated to the control law derivation, and Section 5 provides several numerical examples that demonstrate the feasibility of the suggested relative motion control approach.

2. Preliminaries

The object of the research is a tetrahedral satellite formation, such as the one in the MMS mission, but in a low, near-circular orbit. We consider a small formation of four spacecrafts that orbit around the Earth. The goal is to obtain a reference orbit for a formation in which it conserves its volume and shape. The concrete definitions of volume and shape are presented in aforementioned paper [32]; here, we only present the results.

The problem statement is as follows:

• Four satellites move passively in near-circular orbits, forming a tetrahedron.
• The quality of the tetrahedron at any particular moment in time is defined as

$$\mathbb{Q}=12\frac{{\left(3\mathbb{V}\right)}^{2/3}}{\mathbb{L}}$$

(1)

where $\mathbb{V}$ is the volume and $\mathbb{L}$ is the sum of squares of all six tetrahedron edge lengths.

In general, the quality $\mathbb{Q}$ depicts how close the tetrahedron is to the regular one: the quality equals one for any regular tetrahedron and equals zero for any degenerate tetrahedron.

In previous works, we found rather interesting reference motions that preserved the volume and quality of the tetrahedron. This result was obtained using linearized equations of motion derived by Clohessy and Wiltshire (CW equations) [36].

$$\begin{array}{l}\ddot{x}-2n\dot{y}-3n^{2}x=0,\\ \ddot{y}+2n\dot{x}=0,\\ \ddot{z}+n^{2}z=0.\end{array}$$

(2)

Here, $n=\sqrt{\mu /r_c^3}$ is the mean motion, $\mu$ is the Earth’s gravitational parameter, $r_c$ is the radius of the reference circular orbit. The reference frame here is the Local Vertical Local Horizon (LVLH) system of one of the satellites.

The bounded solutions are then described using equations

$$\begin{array}{l}{x}_i(t)=a_i\sin \nu +b_i\cos \nu ,\\ y_i(t)=2a_i\cos \nu -2b_i\sin \nu +c_i,\\ z_i(t)=d_i\sin \nu +e_i\cos \nu ,\end{array}$$

(3)

where $\nu =nt$. The constants $a_i,b_i,c_i,d_i,e_i$ depend on the initial (at value $t=0$) parameters of the motion; index $i$ attains values 1, 2, and 3. The exact dependence is described using

$$a=\frac{{\dot{x}}_0}{n}, b=-\frac{2{\dot{y}}_0}{n}-3x_0, c=y_0-\frac{2{\dot{x}}_0}{n}, d=\frac{{\dot{z}}_0}{n}, e=z_0, \frac{{\dot{y}}_0}{n}+2x_0=0.$$

(4)

The last equation is the standard CW trajectory boundedness condition. The motion is presented w.r.t. satellite number 0, which stays at the origin.

The one particular case studied in [32] is of interest in the scope of this paper. The (unique) family of initial values for three satellites (the fourth spacecraft is at the origin) simultaneously achieve the following:

• in the CW motion model, the tetrahedron preserves quality and volume;
• the constant quality is maximal;
• one of the three satellites moves along the same orbit (in inertial space) as the reference satellite, i.e., there is a phase shift between them;
• two remaining satellites move in LVLH somehow;

The above is given by the following equations:

$$\begin{array}{l}\left(\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right)=K\left(\begin{array}{c}0\\ \frac{\sqrt{6}}{3}\cos \phi +\frac{\sqrt{3}}{3}\sin \phi \\ \frac{\sqrt{6}}{3}\cos \phi -\frac{\sqrt{3}}{3}\sin \phi \end{array}\right),\left(\begin{array}{c}{e}_1\\ e_2\\ e_3\end{array}\right)=-\sqrt{5}\left(\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right),\\ \left(\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right)=K\left(\begin{array}{c}0\\ -\frac{\sqrt{3}}{3}\cos \phi +\frac{\sqrt{6}}{3}\sin \phi \\ \frac{\sqrt{3}}{3}\cos \phi +\frac{\sqrt{6}}{3}\sin \phi \end{array}\right),\left(\begin{array}{c}{d}_1\\ d_2\\ d_3\end{array}\right)=\sqrt{5}\left(\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right),\\ \left(\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right)=K\sqrt{\frac{5}{3}}\left(\begin{array}{c}2\\ 1\\ 1\end{array}\right),\end{array}$$

(5)

where $\phi$ is an arbitrary initial phase, and $K$ is an arbitrary size factor. Figure 1 depicts such a formation.

It is evident that the motion with initial conditions (4) and (5) in a more realistic non-linear motion model does not preserve the quality and/or volume of the tetrahedron. Even worse, the motion appears to be unstable, so even numerical optimization of the reference orbit is of little use. In addition, the motion simulation, which takes into account the error in the initial conditions, shows poor results in terms of maintaining the quality of the tetrahedron. The exceptional sensitivity of passive motion to the initial conditions implies that in order to maintain the constellation, a control algorithm is required.

3. Controlled Motion Model

3.1. Curvilinear Coordinates

The construction of a simple control algorithm for a tetrahedral formation is carried out using a linear motion model. However, the utilization of standard Descartes coordinates leads to large errors, even for rather small formations. We can improve the accuracy of the model of motion by introducing curvilinear coordinates [37]. The controlled CW equations in curvilinear coordinates then have the following form:

$$\begin{array}{l}\frac{d^2}{d{t}^2}\rho -2n\frac{d}{dt}\left(r_c\theta \right)-3n^2\rho =u_x,\\ \frac{d^2}{d{t}^2}\left(r_c\theta \right)+2n\frac{d}{dt}\rho =u_y,\\ \frac{d^2}{d{t}^2}\left(r_c\phi \right)+n^2\left(r_c\phi \right)=u_z.\end{array}$$

(6)

The aforementioned equations have the same form as CW equations, so the uncontrolled solutions are the same:

$$\begin{array}{l}\rho =2C_{\mathrm{drift}}+C_{\mathrm{inplane}}\sin \left(\nu -{\sigma }_{\mathrm{inplane}}\right),\\ r_c\theta =2C_{\mathrm{inplane}}\cos \left(\nu -{\sigma }_{\mathrm{inplane}}\right)+C_{\mathrm{shift}}-3C_{\mathrm{drift}}\nu ,\\ r_c\phi =C_{\mathrm{outplane}}\sin \left(\nu -{\sigma }_{\mathrm{outplane}}\right).\end{array}$$

(7)

The geometrical meaning of $C$ with various subscripts is similar to amplitudes in (3), i.e., $\sqrt{a^2+b^2}$ is similar to $C_{\mathrm{inplane}}$, though two values are not exactly equal due to the difference between rectangular and curvilinear coordinates. Hereinafter, we neglect this difference when considering the tetrahedron shape and quality.

To be effective in maintaining the tetrahedral formation, the control law we seek must preserve the following quantities:

• the oscillation amplitudes in and out of the orbital plane;
• the shift of each satellite along the orbit;
• drift (which must be equal to zero);
• the relative motion plane orientation;
• phase difference of oscillations between the pair of satellites that move outside the orbital plane.

During passive motion, the above parameters change, but over short time intervals, one can expect them to change slowly. Let us make a transition to new geometrical variables [38]:

$$\begin{array}{l}\rho =A\sin \eta +2C,\frac{d}{dt}\rho =An\cos \eta ,\\ r_c\theta =2A\cos \eta +D,\frac{d}{dt}\left(r_c\theta \right)=-2An\sin \eta -3Cn,\\ r_c\phi =B\sin \lambda ,\frac{d}{dt}\left(r_c\phi \right)=Bn\cos \lambda .\end{array}$$

(8)

The new variables describing relative motion are $A,B,C,D,\eta ,\lambda$. The first two correspond to in-plane and out-plane amplitudes, $C$ to the drift rate, $D$ to the shift, and angle variables to the initial phases.

With (6), we obtain

$$\begin{array}{ll}\dot{C}=\frac{1}{n}{u}_y,\hfill & \dot{D}=-\frac{2}{n}{u}_x-3Cn,\hfill \\ \dot{A}=\frac{1}{n}\left(\cos \eta \cdot u_x-2\sin \eta \cdot u_y\right),\hfill & \dot{B}=\frac{1}{n}\cos \lambda \cdot u_z,\hfill \\ \dot{\eta }=n-\frac{1}{nA}\left(\sin \eta \cdot u_x+2\cos \eta \cdot u_y\right),\hfill & \dot{\lambda }=n-\frac{1}{nB}\sin \lambda \cdot u_z.\hfill \end{array}$$

(9)

To preserve quality, due to (5), the following must be true:

• for the satellite shifted along the reference orbit,

  $$A_1^{ref}=B_1^{ref}=C_1^{ref}=0,D_1^{ref}=2\sqrt{\frac{5}{3}}K$$

  (10)
• for the two remaining satellites,

  $$A_{2,3}^{ref}=K,B_{2,3}^{ref}=\sqrt{5}K,C_{2,3}^{ref}=0,D_{2,3}^{ref}=\sqrt{\frac{5}{3}}.$$

  (11)

Moreover, phases must comply with restrictions

$$\begin{array}{l}{\eta }_3-{\eta }_2={\lambda }_3-{\lambda }_2=\arccos \frac{1}{3},\\ {\eta }_3-{\lambda }_3={\eta }_2-{\lambda }_2=\frac{\pi }{2}\end{array}$$

(12)

3.2. Actuator Model Choice

We need some kind of orbital control instrument. The first obvious choice is to use reactive thrusters: they can provide the necessary maneuvers and maintain the orbital height and formation with high accuracy. However, if we consider small satellites (e.g., 3U-6U CubeSats), it is almost impossible to install such actuators. When constructing a control algorithm for a satellite formation at a height of 400 km, it is quite natural to make a choice in favor of atmospheric drag control. We use the specular–diffuse model [39], where some of the air molecules are reflected from the surface absolutely elastically, while the rest are reflected diffusely with a Maxwellian distribution. The acceleration of a rectangular-shaped satellite in this case is

$$a_{atm}=-{\rho }_{atm}\frac{S}{m}\left(\delta v,n\right)\cdot \left((1-ϵ)\cdot \delta v+2\cdot ϵ\left(\delta v,n\right)n+(1-ϵ)\alpha \cdot \delta v\cdot n\right).$$

(13)

Here, ${\rho }_{atm}$ is the atmospheric density, $S$ is the surface area, $m$ is the satellite mass, $\delta v$ is the satellite velocity with respect to the atmosphere, $\delta v=‖\delta v‖$ is its absolute value, $n$ is the unit normal to the surface, $ϵ$ is the specular reflection coefficient, $\alpha$ is the diffuse reflection coefficient, the round brackets depict inner product. Note that this model allows us to introduce the lift force; however, since $\alpha$ and $ϵ$ are close to zero, its effect is rather small. Since atmospheric drag depends on cross-section area, we can actually control its value by rotating the satellite, e.g., for the 3U CubeSat, the ratio between the minimum and maximum value of atmospheric drag reaches three, which makes it feasible for utilization in relative motion control.

4. Lyapunov Direct Method

We synthesize the control algorithm using the direct Lyapunov method. Our goal is to ensure that relative amplitudes $A_k,B_k$, relative drifts and shifts $C_k,D_k$, and phases ${\lambda }_k,{\eta }_k$ correspond to the reference values defined by Equations (10)–(12). Let us again note that in case of aerodynamic drag forces, out of the three components $(u_x,u_y,u_z)$ of the control acceleration, the value of the force component $u_y$ is at least an order of magnitude greater than the values $u_x$ and $u_z$. Due to that fact, we neglect the values $u_x$ and $u_z$ while synthesizing the control in the $O_y$ direction. That allows us to split the whole control into three different parts.

4.1. Shift and Drift Control

From (9), due to $u_x\approx 0$, it follows that

$$\dot{D}=-3Cn,\hspace{1em}\dot{C}=\frac{1}{n}{u}_y.$$

(14)

The reference values for shift and drift are as follows:

$$\begin{array}{l}{D}_{ref,1}=0,D_{ref,2}=D_{ref,3}=\sqrt{\frac{5}{3}}K,D_{ref,4}=2\sqrt{\frac{5}{3}}K,\\ C_{ref,1}=C_{ref,2}=C_{ref,3}=C_{ref,4}=0.\end{array}$$

(15)

We use the candidate Lyapunov function

$$V_{CD,i}=C_i{}^2+k_D{(D_i-D_{ref,i})}^2,\hspace{1em}{k}_D>0.$$

(16)

Further, the index $i$ is omitted to improve readability; the candidate Lyapunov function and the corresponding control are constructed separately for each satellite.

The function $V_{CD}$ is non-negative and equal to zero if and only if $C=C_{ref}=0, D=D_{ref}$. The derivative of this function due to (14) is

$${\dot{V}}_{CD}=2C\left(\frac{u_y}{n}-3n{k}_D\left(D-D_{ref}\right)\right).$$

(17)

In order for the candidate Lyapunov function to satisfy the Barbashin–Krasovsky–LaSalle theorem [40], it is necessary that its derivative, due to the equations of motion, be non-positive. We require that

$$\left(\frac{u_y}{n}-3n{k}_D\left(D-D_{ref}\right)\right)=-k_{C}C,\hspace{1em}\hspace{1em}{k}_C>0,$$

(18)

which implies

$${\dot{V}}_{CD}=-2k_C{C}^2\le 0.$$

(19)

That in turn implies that the control given by

$$u_y=3n^2{k}_D\left(D-D_{ref}\right)-n{k}_{C}C$$

(20)

ensures the asymptotic stability of the reference motion.

4.2. In-Plane Motion Control

The main troublesome feature of tetrahedron quality control is that the relative positions must be maintained with high precision. According to Equation (5), the phase difference is responsible for the relative motion plane orientation, and a constant phase difference corresponds to a constant orientation of the plane. Consequently, maintaining $\eta -\lambda =\frac{\pi }{2}$ is crucial for relative orbit maintenance. Consider the candidate Lyapunov function

$$V_{A,\varphi }={\left(A_2-A_{ref,2}\right)}^2+{\left(A_3-A_{ref,3}\right)}^2+k_{\varphi }{\left({\eta }_2-{\lambda }_2-\frac{\pi }{2}\right)}^2+k_{\varphi }{\left({\eta }_3-{\lambda }_3-\frac{\pi }{2}\right)}^2,$$

(21)

here, as before, $k_{\varphi }$ is a positive constant. As in the previous step, the derivative is equal to

$$\begin{array}{r}{\dot{V}}_{A,\varphi }=u_{y,2}\frac{4}{n}\left(-k_{\varphi }\left({\eta }_2-{\lambda }_2-\frac{\pi }{2}\right)\frac{\cos {\eta }_2}{A_2}-\left(A_2-A_{ref,2}\right)\sin {\eta }_2\right)\\ +u_{y,3}\frac{4}{n}\left(-k_{\varphi }\left({\eta }_3-{\lambda }_3-\frac{\pi }{2}\right)\frac{\cos {\eta }_3}{A_3}-\left(A_3-A_{ref,3}\right)\sin {\eta }_3\right),\end{array}$$

(22)

and therefore the control

$$\begin{array}{l}{u}_{y,2}=-k_A\left(-k_{\varphi }\left({\eta }_2-{\lambda }_2-\frac{\pi }{2}\right)\frac{\cos {\eta }_2}{A_2}-\left(A_2-A_{ref,2}\right)\sin {\eta }_2\right),\\ u_{y,3}=-k_A\left(-k_{\varphi }\left({\eta }_3-{\lambda }_3-\frac{\pi }{2}\right)\frac{\cos {\eta }_3}{A_3}-\left(A_3-A_{ref,3}\right)\sin {\eta }_3\right),\end{array}$$

(23)

where $k_A$ is positive, maintains the necessary in-plane configuration.

4.3. Out-of-Plane Motion Control

Outside the orbital plane, the motion equations have the form

$$\begin{array}{l}\dot{B}=\frac{1}{n}\cos \lambda \cdot u_z,\\ \dot{\lambda }=n-\frac{1}{nB}\sin \lambda \cdot u_z.\end{array}$$

(24)

The reference values for the second and third satellites are

$$\begin{array}{l}{B}_{ref,2}=B_{ref,3}=\sqrt{5}K,\\ {\lambda }_3-{\lambda }_2=\Delta {\lambda }_{ref}=\arccos \frac{1}{3}.\end{array}$$

(25)

The same procedure with function

$$V_B={\left(B_2-B_{ref,2}\right)}^2+{\left(B_3-B_{ref,3}\right)}^2+k_{\lambda }{\left({\lambda }_3-{\lambda }_2-\Delta {\lambda }_{ref}\right)}^2$$

(26)

produces the control

$$\begin{array}{l}{u}_{z,2}=-k_B\left(k_{\lambda }\left({\lambda }_3-{\lambda }_2-\Delta {\lambda }_{ref}\right)\frac{\sin {\lambda }_2}{B_2}+2\left(B_2-B_{ref,2}\right)\cos {\lambda }_2\right),\\ u_{z,3}=-k_B\left(-k_{\lambda }\left({\lambda }_3-{\lambda }_2-\Delta {\lambda }_{ref}\right)\frac{\sin {\lambda }_3}{B_3}+2\left(B_3-B_{ref,3}\right)\cos {\lambda }_3\right).\end{array}$$

(27)

4.4. The Choice of Control Law

Since relative motion parameters can be arbitrary, in general, control laws (20) and (23) cannot be applied simultaneously; their goal is to preserve different geometric parameters, so the necessary control values would be different. Moreover, these two types of control are designed to preserve different geometric characteristics of the tetrahedron, and therefore, while maintaining the tetrahedron, the priority of these laws is different:

• if the current relative drift is nonzero, ensuring the amplitudes and/or phases of satellite oscillations is meaningless because the relative motion is unbounded, so the tetrahedron degrades in a couple of revolutions;
• on the other hand, if the shift of the satellites in the group is close to the required one and the drift is close to zero, control law (23) can be used to maintain the amplitudes and phases of the oscillations.

The control law is constructed as follows. Some cutoffs of shift and drift are selected:

$$\delta C_{\mathrm{lower}}<\delta C_{\mathrm{upper}},\hspace{1em}\hspace{1em}\delta D_{\mathrm{lower}}<\delta D_{\mathrm{upper}}.$$

(28)

If both current values are less than lower cutoffs,

$$\left|C\right|\le \delta C_{\mathrm{lower}}\hspace{1em}\mathrm{and}\hspace{1em}|D-D_{\mathrm{ref}}|\le \delta D_{\mathrm{lower}},$$

(29)

control law (23) is used. If, on the other hand,

$$\left|C\right|\ge \delta C_{\mathrm{upper}}\hspace{1em}\mathrm{or}\hspace{1em}|D-D_{\mathrm{ref}}|\ge \delta D_{\mathrm{upper}},$$

(30)

the control switches to shift and drift correction mode (20). Two different cutoffs prevent the algorithm from rapidly switching modes.

4.5. Atmospheric Drag Implementation

Since we utilize aerodynamic drag forces for orbital motion control, the suggested control law must be modified. Firstly, aerodynamic drag is significantly limited in magnitude: let $u_{\max }$ be the maximum possible value for the control applied. Secondly, the drag is always pointing in the direction opposite to satellite velocity. Since the entire reference motion model is relative to one of the satellites of the group, we assume that this satellite is under the influence of a constant aerodynamic drag force, $u_{y,4}=\frac{u_{max}}{2}.$ If so, the control $u_{y,i}\in [0,u_{max}/2)$ leads to a negative value of the relative control force between the satellite number $i$ and the reference satellite, and $u_{y,i}\in (u_{max}/2,u_{max}]$ leads to a positive value. Thirdly, as was already mentioned, the magnitude of the force along the $z$ axis is much less than along the $y$ axis. In addition, the formation requires constant adjustment in a direction perpendicular to the orbital plane. In this regard, we always maximize control output in the direction governed by law (27).

Given the necessary control output $u_y,u_z$, we determine the normal to the surface of the satellite:

$$\left(\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right)=\left(\begin{array}{c}\sin \vartheta \sin \zeta \\ \cos \zeta \\ \cos \vartheta \sin \zeta \end{array}\right),$$

(31)

here, $\zeta$ is the angle between the normal and the satellite velocity, and $\vartheta$ is the angle that determines the rotation angle of the satellite around the velocity direction:

$$\begin{array}{l}\zeta =\arccos \left(\frac{u_y-\frac{u_{max}}{2}}{u_{max}}\right),\\ \vartheta =\left(1-\mathrm{sign}(u_z)\right)\cdot \frac{\pi }{2}.\end{array}$$

(32)

Figure 2 depicts the control scheme we have developed.

5. Numerical Simulation

We use the following simulation parameters:

• Initial values are as in (4), (5), where $\phi =0$, and the value $K$ is specified below.
• The altitude of the circular reference orbit is 400 km. The inclination is equal to $56°$.
• The motion of each satellite is integrated independently via the RK4 method in IRF.
• Initial values contain normally distributed noise.
• The Earth’s gravitational field is modelled up to terms of order (10, 10).
• The atmosphere is considered to rotate with the Earth; its density is calculated according to the Russian GOST model (it takes into account solar activity, day/night, and summer/winter variations).
• The satellite is a 5 kg square plate with an area of 0.1 m². Its interaction with the atmosphere is in accordance with (13), where $\alpha =ϵ=0.1$.
• The position and velocity in inertial space are known precisely.
• It is assumed that the necessary attitude is always achievable. Attitude change occurs instantly.

5.1. Small Formation

Figure 3 depicts simulation results of a small formation with $K=100 \mathrm{m}$. It is evident that the tetrahedron quality remains almost the same for a long period of time.

5.2. Possible Insufficiency of Control

For bigger formation of size $K=1000 \mathrm{m}$, the simulation results vary. Figure 4 depicts a typical result of a simulation.

As can be seen, the quality of the tetrahedron is within 0.4 over approximately 200 revolutions, and within 0.2 over 600 revolutions. The tetrahedron remains non-degenerate for about 1000 revolutions (about 60 days). It might be explained by increased values of disturbances, both from non-linear terms in relative motion equations and additional harmonics in the Earth gravity model.

Table 1. The probability of control sufficiency for K = 1000 m. Different variances for normally distributed noise are considered. One hundred simulations with independent initial condition errors were carried out.

$\mathit{\sigma }$	1 m	5 m	10 m	20 m
0.5 cm/s	97%	92%	82%	41%
1 cm/s	91%	86%	73%	32%
2 cm/s	76%	62%	48%	23%

summarizes the dependence of the formation degradation rate on the initial conditions’ error. In each specific case, the control is either sufficient to counter the initial disturbance or insufficient. In the former case, the tetrahedron retains nonzero quality for at least 600 revolutions. In the latter case, the tetrahedron becomes uncontrollable and degrades (i.e., the quality becomes zero) within the first month.

As can be seen, the initial velocity error significantly reduces the chances of successful control. That is due to the low possible thrust outcome of the aerodynamic-based control and, furthermore, to the built-in unpredictability of the exact value of aerodynamic force. The position error affects the algorithm less. According to the results, it may be necessary to use additional control in the first phase of the mission to achieve high accuracy in the initial motion parameters.

The numerical simulations presented above show that the proposed control can maintain the quality of the tetrahedron at an acceptable level for several weeks.

6. Discussion and Conclusions

The presented paper provides a control law construction algorithm that is able to significantly increase the lifetime of a near-Earth tetrahedral satellite formation.

The proposed simple atmospheric drag-based approach for controller construction works rather well for a small formation. The MMS mission, though using a slightly different approach to quality definition, had a minimal requirement of a quality of 0.7. This was achieved using high-quality optimization techniques and fuel consumption. Our approach works not quite as well but surely can be used as an auxiliary algorithm to reduce the delta-V necessary for formation maintenance. Furthermore, for missions in LEO, which would study Earth’s ionosphere, time and space intervals between measurements are naturally smaller compared to the magnetospheric scale. This allows us to hope that the proposed algorithm alone will cope with maintaining the tetrahedral formation, possibly with some loss of measurement accuracy.

For larger formations, the proposed algorithm performs significantly worse. The resulting loss in quality entails the need for additional control. However, the proposed technique can save a significant amount of fuel, especially when providing suitable initial formation data. A quality of 0.2–0.4 is unlikely to allow for high-quality measurements, which, however, can be offset to some extent by the extended duration of the experiment. Specific requirements for the quality of the formation depend on the mission in question; we note here again that any nonzero quality corresponds to a non-degenerate tetrahedron and therefore allows us, in principle, to make 3D measurements.