Decentralized Differential Aerodynamic Control of Microsatellites Formation with Sunlight Reflectors

Abstract

The paper presents a study of decentralized control for a satellite formation flying mission that uses differential lift and drag to enforce the relative positioning requirements. All spacecraft are equipped with large sunlight reflectors so that, given the appropriate lighting conditions, the formation as a whole can be made visible from the Earth as a configurable pixel image in the sky. The paper analyzes the possibility of achieving a pre-defined lineup of the formation by implementing decentralized aerodynamic-based control through the orientation of sunlight reflectors relative to the incoming airflow. The required relative trajectories are so-called projected circular orbits which ensure the rotation of the image with the orbital period. The choice of the reference trajectory for each satellite is obtained by minimizing the total sum of relative trajectory residuals. The control law is based on the linear-quadratic regulator with the decentralized objective function of reducing the mean deviation of each satellite’s trajectory relative to the other satellites. The accuracy of the required image construction and convergence time depending on the initial conditions and orbit altitude are studied in the paper.

1. Introduction

Small satellite formation missions are widely applied for various important problems and tasks in near-Earth orbits. A number of satellites flying at short relative distances can serve as a distributed instrument, for example, for ionosphere and magnetosphere sampling [1], stereo remote sensing [2], and astrophysical research [3]. Now that Starlink and OneWeb mega constellations are launching satellites in large numbers, satellite flocks have been observed flying over the night sky in different parts of the world, visible to the naked eye [4]. The sunlight reflected by the surface of satellites moving near the terminator line is bright enough to see the flying line of dots. This effect underlies the idea of launching a set of microsatellites equipped with sunlight reflectors to obtain a relative configuration that forms an image visible to an observer on the Earth. Such a formation flying mission can thus function as a space media broadcasting logos or short messages. A recent paper [5] discusses the commercial viability of operating such formations. Surprisingly, the problem of designing a space mission to demonstrate pixel images in the night sky has already been considered in a very practical sense. At least two space advertising projects were discussed in the 20th century. A string of 100 reflectors to form a ring of light, visible throughout the world, could be sent into orbit in 1989 to mark the centennial of the Eiffel Tower [6]. Then, in the 1990s, the city of Atlanta investigated a Space Billboard concept for the 1996 Olympics [7]. As we have already mentioned, the economic viability of space advertising missions was studied in [5] and was found impossible without reactive propulsion systems for frequent image reconfigurations. However, single-event missions such as the Eiffel Tower anniversary or the Olympics opening ceremony which only require a single demonstration might use the control strategy described in this paper, which allows for the omission of the propulsion system from the mission. Another practical aspect of the problem discussed in this paper is the problem statement itself. We believe that the problem of a formation flying mission to form pixel images can be made a benchmark test for any formation control algorithms, because of very clear visibility requirements, which need to be met and can be met on the up-to-date level of technology.

A preliminary feasibility study [5] showed that a formation to meet the visibility requirements under suitable lighting conditions launched at low Earth Sun-synchronous orbits of 400 km can be formed using 12U CubeSats with 2 × 2 m² sails. The minimum distance between any two reflectors should be greater than 600 m to make the two points distinguishable by the human eye. If the attitude of the image is fixed in the orbital reference frame, it requires continuous control using onboard thrusters and leads to excessive fuel consumption. Hence, the preferred approach is to let the image rotate in the orbital reference frame with the orbital period. Ideally, each satellite can be assigned initial conditions that result in a circle motion when projected onto the plane perpendicular to nadir, so that each “pixel” in the image rotates with the same angular velocity. Relative motion control, however, is required to achieve the selected relative trajectories. A hybrid impulsive-continuous control scheme allowing to converge to the required trajectory using continuous thrust was proposed in [8]. This approach relies on a propulsion system and becomes inappropriate if any of the satellites run out of fuel. The recent paper [9] proposes the new transforming equations using topocentric coordinate system allowing to achieve arbitrary advertisement shape, and stable orbits that keep the formation shape without control are found. An approach proposed in this paper is based on using aerodynamic forces, which is acceptable for orbits of 400 km and lower. The sunlight reflectors can function as drag sails. By changing their orientation relative to the incoming airflow, it is possible to employ the differential drag and lift forces for relative motion control. It is assumed that after convergence to the reference relative trajectories, all sunlight reflectors can be synchronously turned and form a graphical image to be observed on Earth at a given point. A preliminary study [10] showed that the linear-quadratic regulator (LQR) approach based on a linear motion model in Cartesian coordinates does not ensure convergence for satellites whose trajectories lie farthest from the image center. The linearization problem is addressed in this paper by the choice of curvilinear coordinates for the relative motion model in the control algorithm. Differential control using aerodynamic forces is widely studied for formation flying applications. The pioneering work of Leonard [11] proposed a proportional-differential controller for cross-section change to control the relative motion of two satellites. Most papers consider two-satellite formation flying with the application of a variety of different control algorithms using a differential drag [12,13,14,15,16]. Differential drag control for multiple satellites is considered in several papers, proposing cyclic and optimal control strategies for cluster flight [17,18], and decentralized control for swarms of satellites with communication restrictions [19]. An application of the differential lift along with the drag for the small satellite rendezvous problem was first proposed by Horsley in [20]. Papers [21,22,23,24] further explored the concept of two-satellite formation flying control using various control approaches. In [25], a decentralized relative motion control is proposed for four 3U CubeSats constructing tetrahedral formation using differential drag and lift forces.

The main contribution of this paper is an improved decentralized control algorithm based on LQR using the mean deviation from the reference relative trajectories expressed in curvilinear relative coordinates. The reference trajectory for each satellite is obtained by minimizing the total sum of relative trajectory residuals. The performance of the proposed control algorithm is investigated by numerical simulations.

The paper has the following structure. In Section 2, the problem statement is formulated and the equations of motion are presented. In Section 3, a detailed description of the proposed decentralized control algorithm is provided, and its implementation using aerodynamic forces is discussed. The numerical simulation results and their analysis are presented in Section 4. The last section summarizes the principal results of the study.

2. Problem Statement and Equations of Motion

A number of small satellites after a cluster launch in LEO are considered. All satellites are equipped with flat sunlight reflectors fixed to the satellite body. An active attitude control system is used to achieve the required angular position relative to the Sun’s direction or to the incoming airflow. Sunlight is reflected to the observer’s position on the Earth’s surface under the required attitude of the reflector relative to the Sun’s direction, and it can be visible to the naked eye in the night sky (see Figure 1). Meanwhile, the sunlight reflector is used as an aerodynamic wing for relative motion control in low-Earth orbits. A number of sunlight reflectors form a specified image, functioning like pixels. It is assumed that all the satellites are provided with information about the relative motion of all the other satellites, i.e., the relative motion is assumed to be known. The main goal of the aerodynamic-force-based relative motion control is to achieve the required configuration to compose the chosen image.

The linear model is used for the reference relative trajectories calculation. It can be described by Hill–Clohessy–Wiltshire equations [26], written in the local-vertical local-horizontal (LVLH) reference frame. The origin O of the LVLH reference frame moves along a circular orbit, the Oz axis is aligned along the vector from the Earth center $O_I$ to O, the Oy axis is directed along the orbital angular momentum vector, and the Ox axis completes the right-handed triad. The HCW equations are as follows:

$$\left\{\begin{array}{l}{\ddot{x}}_{ij}\hspace{0.17em}+\hspace{0.17em}2\omega {\dot{z}}_{ij}\hspace{0.17em}=\hspace{0.17em}{u}_x^{ij},\\ {\ddot{y}}_{ij}\hspace{0.17em}+\hspace{0.17em}{\omega }^2{y}_{ij}\hspace{0.17em}=\hspace{0.17em}{u}_y^{ij},\\ {\ddot{z}}_{ij}\hspace{0.17em}-\hspace{0.17em}2\omega {\dot{x}}_{ij}\hspace{0.17em}-\hspace{0.17em}3{\omega }^2{z}_{ij}\hspace{0.17em}=\hspace{0.17em}{u}_z^{ij}.\end{array}\right.$$

(1)

Here ${[x_{ij},y_{ij},z_{ij}]}^{\mathrm{T}}=r_{ij}=r_j-r_i$ is the difference between the positions of the j-th and i-th satellites in the LVLH reference frame, ω is the orbital angular velocity, and ${[u_x^{ij},u_y^{ij},u_z^{ij}]}^{\mathrm{T}}=u^{ij}=u^j-u^i$ is a vector of control acceleration difference. If the control input in Equation (1) is zero, the system becomes homogeneous and its solution is:

$$\left\{\begin{array}{l}{x}_{ij}(t)=-3C_1^{ij}\omega t+2C_2^{ij}\cos (\omega t)-2C_3^{ij}\sin (\omega t)+C_4^{ij},\\ y_{ij}(t)=C_5^{ij}\sin (\omega t)+C_6^{ij}\cos (\omega t),\\ z_{ij}(t)=2C_1^{ij}+C_2^{ij}\sin (\omega t)+C_3^{ij}\cos (\omega t),\end{array}\right.$$

(2)

where constants $C_1^{ij}-C_6^{ij}$ depend on the initial conditions at time $t_0=0$ as follows:

$$\begin{array}{l}{C}_1^{ij}=\frac{{\dot{x}}_{ij}(0)}{\omega }+2z_{ij}(0),\hspace{1em}{C}_2^{ij}=\frac{{\dot{z}}_{ij}(0)}{\omega },\hspace{1em}{C}_3^{ij}=-2\frac{{\dot{x}}_{ij}(0)}{\omega }-3z_{ij}(0),\\ C_4^{ij}=-2\frac{{\dot{z}}_{ij}(0)}{\omega }+x_{ij}(0),\hspace{1em}{C}_5^{ij}=\frac{{\dot{y}}_{ij}(0)}{\omega },\hspace{1em}{C}_6^{ij}=y_{ij}(0).\end{array}$$

(3)

If all the constants $C_1^{ij}$ are equal to zero, then all the relative trajectories are closed ellipses. For certain sets of the constant values, the relative orbits become projected circular orbits, which means that their projections onto the local horizontal plane Oxy are circular. In this paper, all the reference trajectories for each satellite are considered to be projected circular orbits (PCOs). The constants $C_1^{ij}-C_6^{ij}$ determine the size and the phase parameters of each “pixel” motion relative to the image center.

It is known that the Hill–Clohessy–Wiltshire relative motion Equation (1) expressed in Cartesian coordinates is valid for short relative distances between the satellites and for short-term motion propagation. The comparison of the relative motion model for formation flying was carried out in [27]. It was shown that the greater the distance, the higher the deviation of the free motion trajectory (2) compared to the actual relative trajectory, calculated as the difference of orbital trajectories. These errors are mainly caused by the Cartesian coordinates representation of the relative trajectory. As shown in [8], the satellite formation for image demonstration can be as large as several kilometers in size, and in this case, the control calculation using the Cartesian coordinates can lead to significant errors. Due to these errors, the control algorithm for tracking reference trajectories (2) results in divergence in time [9]. To reduce the control errors, the following relative motion equations in curvilinear coordinates are used in this work:

$$\left\{\begin{array}{l}{\ddot{\phi }}_{ij}\hspace{0.17em}+\hspace{0.17em}2\omega {\dot{\rho }}_{ij}/a_0\hspace{0.17em}=\hspace{0.17em}{w}_{\phi }^{ij},\\ {\ddot{\theta }}_{ij}\hspace{0.17em}+\hspace{0.17em}{\omega }^2{\theta }_{ij}\hspace{0.17em}=\hspace{0.17em}{w}_{\theta }^{ij},\\ {\ddot{\rho }}_{ij}\hspace{0.17em}-\hspace{0.17em}2\omega a_0{\dot{\phi }}_{ij}\hspace{0.17em}-\hspace{0.17em}3{\omega }^2{\rho }_{ij}\hspace{0.17em}=\hspace{0.17em}{w}_{\rho }^{ij}.\end{array}\right.$$

(4)

Curvilinear coordinates ${\phi }_{ij},\hspace{0.17em}\hspace{0.17em}{\theta }_{ij},\hspace{0.17em}\hspace{0.17em}{\rho }_{ij}$ are defined as follows (Figure 2): ${\phi }_{ij}$ is the angle between the radius-vector $R_i$ of the i-th satellite and the projection of the radius-vector $R_j$ of the j-th satellite onto the orbital plane of the i-th satellite, ${\theta }_{ij}$ is the angle between the radius-vector of the j-th satellite $R_j$ and projection of the radius-vector of the j-th satellite onto the orbital plane of the i-th satellite, ${\rho }_{ij}=\hspace{0.17em}\hspace{0.17em}\left|R_j\right|-\left|R_i\right|$, $a_0=\hspace{0.17em}\hspace{0.17em}\left|R_i\right|$ is the distance between the Earth center $O_I$ and the i-th satellite, and $w={\left[w_{\phi }^{ij},\hspace{0.17em}\hspace{0.17em}{w}_{\theta }^{ij},\hspace{0.17em}\hspace{0.17em}{w}_{\rho }^{ij}\right]}^T$ is the control acceleration vector in curvilinear coordinates.

The solution to the free motion Equation (4) without control is similar to Equation (2) and is given by:

$$\left\{\begin{array}{l}{\phi }_{ij}(t)=-3{\tilde{C}}_1^{ij}\omega t+2{\tilde{C}}_2^{ij}\cos (\omega t)-2{\tilde{C}}_3^{ij}\sin (\omega t)+{\tilde{C}}_4^{ij},\\ {\theta }_{ij}(t)={\tilde{C}}_5^{ij}\sin (\omega t)+{\tilde{C}}_6^{ij}\cos (\omega t),\\ {\rho }_{ij}(t)=2{\tilde{C}}_1^{ij}{a}_0+{\tilde{C}}_2^{ij}{a}_0\sin (\omega t)+{\tilde{C}}_3^{ij}{a}_0\cos (\omega t),\end{array}\right.$$

(5)

where constants ${\tilde{C}}_1^{ij}-{\tilde{C}}_6^{ij}$ depend on the initial conditions at time $t_0=0$ similarly to Equation (3) as follows:

$$\begin{array}{l}{\tilde{C}}_1^{ij}=\frac{{\dot{\phi }}_{ij}(0)}{\omega }+2\frac{{\rho }_{ij}(0)}{a_0},\hspace{1em}{\tilde{C}}_2^{ij}=\frac{1}{\omega }\frac{{\dot{\rho }}_{ij}(0)}{a_0},\hspace{1em}{\tilde{C}}_3^{ij}=-2\frac{{\dot{\phi }}_{ij}(0)}{\omega }-3\frac{{\rho }_{ij}(0)}{a_0},\\ {\tilde{C}}_4^{ij}=-\frac{2}{\omega }\frac{{\dot{\rho }}_{ij}(0)}{a_0}+{\phi }_{ij}(0),\hspace{1em}{\tilde{C}}_5^{ij}=\frac{{\dot{\theta }}_{ij}(0)}{\omega },\hspace{1em}{\tilde{C}}_6^{ij}={\theta }_{ij}(0).\end{array}$$

(6)

To evaluate the control algorithm utilizing the curvilinear coordinates representation, a numerical simulation study was conducted. Numerical simulations are based on the nonlinear orbital motion model that takes into account the $J_2$ gravity perturbation and influence of the aerodynamic forces:

$${\ddot{R}}_i=-\frac{{\mu }_g{R}_i}{R_i^3}+f_i^{\mathrm{grav}}+f_i^{\mathrm{aero}},$$

(7)

where ${\mu }_g=3.986\cdot {10}^{14}{\mathrm{m}}^3/{\mathrm{s}}^2$ is the Earth gravitational parameter, $R_i={\left[X_i,Y_i,Z_i\right]}^T$ is the radius-vector of the i-th satellite in the Earth-centered inertial reference frame (ECI), and $f_i^{\mathrm{grav}}$ is the vector of acceleration caused by the second zonal harmonic of Earth gravitational potential:

$$f_i^{\mathrm{grav}}=\left(\begin{array}{c}\frac{\delta }{R_i^5}{X}_i\left(\frac{5Z_i^2}{R_i^2}-1\right)\\ \frac{\delta }{R_i^5}{Y}_i\left(\frac{5Z_i^2}{R_i^2}-1\right)\\ \frac{\delta }{R_i^5}{Z}_i\left(\frac{5Z_i^2}{R_i^2}-3\right)\end{array}\right),\hspace{1em}\delta =\frac{3}{2}{J}_2{\mu }_g{R}_E^2,\hspace{1em}{J}_2=1082.23·{10}^{-6}$$

(8)

$R_E^{}$ is the Earth’s average radius, and $f_i^{\mathrm{aero}}$ is the acceleration vector due to the aerodynamic force, described by the following formula [28]:

$$f_i^{\mathrm{aero}}=-\frac{1}{m_i}\rho (R_i)V_i^2{S}_i[(1-\epsilon )(e_V,n_i)e_V+2\epsilon {(e_V,n_i)}^2{n}_i+(1-\epsilon )\eta (e_V,n_i)n_i],$$

(9)

where $m_i$ is satellite mass, $\rho (R_i)$ is the atmospheric density at the orbital position $R_i$, $V_i$ is the magnitude of satellite’s velocity relative to the incoming airflow, $S_i$ is the area of the reflector, $e_V$ is the unit vector along the velocity, $n_i$ is the unit vector of the normal to the reflector surface, and $\epsilon$ and $\eta$ are parameters of interaction of molecules with the surface, which are estimated by Beletsky in [28]; in LEO, their values are $\epsilon \approx 0.1$, $\eta \approx 0.1$. In aerodynamic force model (9), the first term is the aerodynamic drag, and the second and the third terms describe the lift force. Due to the lift component of the force, which is 10 times lower than the value of drag force, it is possible to control the relative motion in the out-of-orbital plane to achieve the required projected circular relative trajectories.

3. Control Algorithm Description

It is further assumed that an image to be demonstrated by the satellites consists of several English alphabet letters. However, it could be any symbol or drawing. All letters or symbols can be represented as a composition of points with defined relative position. Each point of a letter is a satellite with a sunlight reflector in our case. For image demonstration using sunlight reflectors, the required relative motion of the satellites should be achieved with an acceptable margin. The relative configuration is based on the closed trajectories with zero relative drift between all the satellites. The projected circular orbits are used for image construction to preserve the relative distances in projection onto the horizontal plane, which is suitable for the Earth observer. Instead of Cartesian coordinates, the curvilinear relative coordinates are utilized for the reference relative trajectories construction, which improves the accuracy of the relative motion model. Since the satellites are in low Earth orbits with an altitude of 300–400 km, the difference of aerodynamic forces acting on the satellites is considered as the control source. By changing the attitude of the reflector relative to the incoming airflow, it is possible to obtain the required value and the direction of the aerodynamic drag and lift. The linear-quadratic regulator is applied to the control problem for converging to zero mean deviation of the current trajectory from the reference one for each satellite in the formation. This section provides the details of the control scheme designed to achieve the required relative configuration.

3.1. Reference Trajectories

The free motion trajectory (2) can be rewritten as follows:

$$\left\{\begin{array}{l}{x}_{ij}(t)=-3B_1^{ij}\omega t+2B_2^{ij}\cos (\omega t+{\alpha }_{ij})+B_4^{ij},\\ y_{ij}(t)=B_3\sin (\omega t+{\beta }_{ij}),\\ z_{ij}(t)=2B_1^{ij}+B_2^{ij}\sin (\omega t+{\alpha }_{ij}),\end{array}\right.$$

(10)

where $B_1^{ij}-B_4^{ij}$ are trajectory constants defined using $C_1^{ij}-C_6^{ij}$ constants from Equation (3) as follows:

$$B_1^{ij}=C_1^{ij},\hspace{0.17em}\hspace{0.17em}{B}_2^{ij}=\sqrt{{\left(C_2^{ij}\right)}^2+{\left(C_3^{ij}\right)}^2},\hspace{0.17em}\hspace{0.17em}{B}_3^{ij}=\sqrt{{\left(C_5^{ij}\right)}^2+{\left(C_6^{ij}\right)}^2},\hspace{0.17em}\hspace{0.17em}{B}_4^{ij}=C_4^{ij},$$

(11)

${\alpha }_{ij},{\beta }_{ij}$ are phase angles.

If the relative drift parameter $B_1^{ij}$ and the relative trajectory shift parameter $B_4^{ij}$ equal zero, and the in-plane amplitude $B_2^{ij}$ and the out-of-plane amplitude $B_3^{ij}$ satisfy the relation $B_3^{ij}=2B_2^{ij}$, then the relative trajectories are projected circular orbits. It means that under these initial conditions for any time $t$, the following relation holds true for the trajectory coordinates:

$$x_{ij}^2(t)+y_{ij}^2(t)=x_{ij}^2(0)+y_{ij}^2(0)={\left(B_3^{ij}\right)}^2=const$$

(12)

Projected circular relative trajectories are well-suited for image demonstration by sunlight reflectors because during the motion the relative geometry is maintained for the Earth observer, although during the flight the image will rotate with orbital angular velocity $\omega$ relative to the image center. If all the satellites have the same out-of-plane phase angle ${\beta }_{ij}$, then by choosing the in-plane phase angles ${\alpha }_{ij}$ and the trajectory radius $B_3^{ij}$, it is possible to set the required relative position of each satellite.

Consider, for example, that some satellites are initially placed in two perpendicular lines, and they start to move along the projected circular reference orbits. Figure 3 presents the initial positions of nine satellites on perpendicular red and green lines and the closed relative trajectories in the $Oxy$ plane. Two initially perpendicular lines formed by the satellites remain perpendicular, which means if the letters to be demonstrated are constructed in this way they will not be deformed for the observer at zenith.

The projected circular trajectory in the $Oxy$ plane is a linear approximation of a more complex natural periodic closed trajectory. A model to describe it with greater accuracy uses the solution to Hill–Clohessy–Wiltshire Equation (5) expressed in the curvilinear coordinates. The trajectory (5) can be rewritten into the form similar to Equation (10):

$$\left\{\begin{array}{l}{\phi }_{ij}(t)=-3{\tilde{B}}_1^{ij}\omega t+2{\tilde{B}}_2^{ij}\cos (\omega t+{\tilde{\alpha }}_{ij})+{\tilde{B}}_4^{ij},\\ {\theta }_{ij}(t)={\tilde{B}}_3\sin (\omega t+{\tilde{\beta }}_{ij}),\\ {\rho }_{ij}(t)=2{\tilde{B}}_1^{ij}{a}_0+{\tilde{B}}_2^{ij}{a}_0\sin (\omega t+{\tilde{\alpha }}_{ij}),\end{array}\right.$$

(13)

where ${\tilde{B}}_1^{ij}-{\tilde{B}}_4^{ij},{\tilde{\alpha }}_{ij},{\tilde{\beta }}_{ij}$ are the constants defined using ${\tilde{C}}_1^{ij}-{\tilde{C}}_6^{ij}$ similarly to (11). The condition (12) for the projected circular orbit in the curvilinear coordinates becomes:

$${\left({\phi }_{ij}^{}(t)a_0\right)}^2+{\left({\theta }_{ij}^{}(t)a_0\right)}^2={\left({\tilde{B}}_3^{ij}{a}_0\right)}^2=const.$$

(14)

The trajectory expressed in the curvilinear coordinates with zero relative drift parameter ${\tilde{B}}_1^{ij}$ and zero relative shift parameter ${\tilde{B}}_4^{ij}$ subject to condition (14) has circular projection on the sphere with radius $a_0$, in contrast to circular projection on the $Oxy$ plane for Cartesian coordinates. This kind of closed trajectory is also acceptable for the Earth observer, since at close relative distance the relative geometry is not notably affected. However, such a reference trajectory is a more accurate representation of the natural relative dynamics, and the control algorithm for tracking this trajectory is expected to have a lower error compared to the reference trajectory expressed in the Cartesian coordinates. This aspect of the reference trajectory representation turns out to be essential for a formation flying mission considered in this paper.

3.2. Linear-Quadratic Regulator

In this work, the linear-quadratic regulator (LQR) is used for reference trajectory tracking problems. LQR is a well-known control algorithm, which can applied to the standard linear-quadratic problem in the case of two satellites’ relative motion control. However, in the case of multiple satellites in formation, several different control vectors provided by LQR cannot be implemented simultaneously. Thus, the control implementation strategy to deal with this problem is developed in this paper.

Consider the standard formulation of LQR. The equations of motion (4) can be represented in the following linear form:

$${\dot{x}}_{ij}=A{x}_{ij}+B{w}^{ij}$$

where $x_{ij}$ is the relative state vector consisting of curvilinear coordinates and their derivatives:

$$x_{ij}={\left[\begin{array}{cccccc}{\phi }_{ij}& {\theta }_{ij}& {\rho }_{ij}& {\dot{\phi }}_{ij}& {\dot{\theta }}_{ij}& {\dot{\rho }}_{ij}\end{array}\right]}^{\mathrm{T}},$$

(15)

$A$ is a system matrix, and $B$ is the control matrix of the following form:

$$A=\left(\begin{array}{cc}{O}_{3\times 3}& E_{3\times 3}\\ C& D\end{array}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}B=\left(\begin{array}{c}{O}_{3\times 3}\\ E_{3\times 3}\end{array}\right),$$

where

$$C=\left(\begin{array}{ccc}0& 0& 0\\ 0& -{\omega }^2& 0\\ 0& 0& 3{\omega }^2\end{array}\right),\hspace{1em}D=\left(\begin{array}{ccc}0& 0& -2\omega /a_0\\ 0& 0& 0\\ 2\omega a_0& 0& 0\end{array}\right),$$

$O_{3\times 3}$ is the zero matrix, and $E_{3\times 3}$ is the identity matrix.

Let us designate the trajectory deviation from the reference trajectory $x_{ij}^{\mathrm{des}}$ by $e_{ij}=x_{ij}-x_{ij}^{\mathrm{des}}$. The linear-quadratic regulator is the feedback control [29], minimizing the following function:

$$J={\displaystyle \underset{0}{\overset{+\infty }{\int }}(}{e}_{ij}^{\mathrm{T}}Q{e}_{ij}+w_{ij}^{\mathrm{T}}R{w}_{ij})dt,$$

where $Q$ and $R$ are positive-definite matrices which are assumed to be diagonal in this paper. The control vector is calculated as follows:

$$w_{ij}=-R^{-1}{B}^{T}P{e}_{ij},$$

(16)

where matrix $P$ is the solution to the algebraic Riccati equation:

$$A^{T}P+PA-PB{R}^{-1}{B}^{T}P+Q=0$$

3.3. Calculation of Acceleration Vector in Cartesian Coordinates

According to Equation (16), the control acceleration vector in the curvilinear coordinates $w_{ij}={[{\ddot{\phi }}_{ij}^w\hspace{0.33em}\hspace{0.33em}{\ddot{\theta }}_{ij}^w\hspace{0.33em}\hspace{0.33em}{\ddot{\rho }}_{ij}^w]}^T$ has to be transformed to Cartesian coordinates in order to calculate the required linear acceleration $u_{ij}$ in the LVLH reference frame. To this end, the Lame approach is adopted. The transition between the Cartesian and the curvilinear coordinates is provided by the following relations [30]:

$$\left\{\begin{array}{l}{x}_{ij}=({\rho }_{ij}+a_0)\cos {\theta }_{ij}\sin {\phi }_{ij},\\ y_{ij}=({\rho }_{ij}+a_0)\sin {\theta }_{ij},\\ z_{ij}=({\rho }_{ij}+a_0)\cos {\theta }_{ij}\cos {\phi }_{ij}.\end{array}\right.$$

(17)

For the calculation of the components $u_k$ of control acceleration vector $u$, the following formula is used:

$$u_k=\frac{1}{H_k}\left[\frac{d}{dt}\left(\frac{\partial V^2/2}{\partial {\dot{q}}_k}\right)-\frac{\partial V^2/2}{\partial q_k}\right],$$

(18)

where $q={\left[\begin{array}{ccc}{\phi }_{ij}& {\theta }_{ij}& {\rho }_{ij}\end{array}\right]}^T$ is a vector of curvilinear coordinates, $k=1,\hspace{0.17em}\hspace{0.17em}2,\hspace{0.17em}\hspace{0.17em}3$, and $V$ is the absolute value of the velocity, which is given by:

$$V^2={\displaystyle \sum _{k=1}^3{\dot{q}}_k^2{H}_k^2},$$

(19)

$H_k$ is the Lame coefficient defined by the formula:

$$H_k^2={\left(\frac{\partial x_{ij}}{\partial q_k}\right)}^2+{\left(\frac{\partial y_{ij}}{\partial q_k}\right)}^2+{\left(\frac{\partial z_{ij}}{\partial q_k}\right)}^2.$$

(20)

Using relations (17)–(20), the transition of the control acceleration components calculated in the curvilinear coordinates to control acceleration components expressed in the Cartesian coordinates is as follows:

$$\begin{array}{l}{u}_x^{ij}=2{\dot{\phi }}_{ij}(-({\rho }_{ij}+a_0){\dot{\theta }}_{ij}\sin {\theta }_{ij}+{\dot{\rho }}_{ij}\cos {\theta }_{ij})+{({\dot{\rho }}_{ij}+a_0)}^2{\ddot{\phi }}_{ij}^w\cos {\theta }_{ij},\\ u_y^{ij}=({\rho }_{ij}+a_0)({\ddot{\theta }}_{ij}^w+{\dot{\phi }}_{ij}^2\cos {\theta }_{ij}\sin {\theta }_{ij})+2{\dot{\rho }}_{ij}{\dot{\theta }}_{ij},\\ u_z^{ij}={\ddot{\rho }}_{ij}^w-({\rho }_{ij}+a_0)({\dot{\phi }}_{ij}^2{\cos }^2{\theta }_{ij}+{\dot{\theta }}_{ij}^2).\end{array}$$

(21)

3.4. Relative Control of Multiple Satellites

The control algorithm (16) is aimed at the relative trajectory tracking of two satellites. However, for multiple satellite formation flying when it is required to achieve the reference trajectory relative to all satellites, the control algorithm needs to be modified. In this paper, the decentralized control approach is used, which implies that each satellite implements its own control regardless of the other satellites’ control actions.

In the case of $N$ satellites in the group, each satellite has $N-1$ consensus reference trajectories relative to other satellites. It means that for each $N-1$ pair of satellites, the control acceleration $u_{ij}$ according to (16) and (21) can be calculated. In the decentralized approach, the $j$-th satellite lacks information about the control of the $i$-th, leading the $j$-th satellite to assume that the $i$-th is passive (i.e., $u^i=0$) and it implements the difference in the control acceleration $u^{ij}=u^j-u^i$ by itself. Thus, each $j$-th satellite should implement simultaneously $N-1$ different control accelerations $u^{ij}$, which is impossible. For this reason, in this paper, the mean value of the control acceleration is to be implemented according to the following:

$$u_j=\frac{{\displaystyle \sum _{i=1}^N{u}^{ij}}}{N-1},i\ne j,$$

(22)

assuming that $u^{ii}$ does not exist.

The control acceleration vector $u_j$ is to be implemented by the aerodynamic forces, which has strict limitations on the value according to (9), and it is impossible to implement positive values for the along-track component $u_j^x$. Let us consider the following example. With the averaged relative control (22), it is possible that all satellites except one achieve their respective reference trajectories relative to the majority of other satellites, causing the corresponding control acceleration values $u^{ij}$ to become close to zero. Meanwhile, one satellite still has a large deviation $e_{ij}$ relative to the others. Assume that according to Equation (22), the $j$-th satellite needs to implement a positive control action in $u_j^x$, which is not possible by aerodynamic drag. It means that all other satellites should implement negative values in their control acceleration, since $u^{ij}=u^j-u^i$ by definition. However, other satellites have small values for all but one relative control acceleration, and the resulting sum is divided by $N-1$. Thus, the control actions implemented by the other satellites are much less than the required control $u^{ij}$ relative to the $j$-th satellite that has not converged. In this situation, the convergence time to the reference trajectories for the $j$-th satellite could become considerably large.

To avoid such undesirable cases, the acceptable deviation norm $Err=\left|e_{ij}\right|$ from the reference trajectory is introduced. The control acceleration vector $u^{ij}$ is used for the sum calculation in (22) only when the corresponding deviation norm is less than the acceptable value $\left|e_{ij}\right|>Err$. Thus, the modified control for the $j$-th satellite is as follows:

$$u_j={\displaystyle \sum _{i\in G_j}{u}^{ij}}/N_{G_j},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{G}_j=\{i:\hspace{0.33em}|e_{ij}^{}|\hspace{0.33em}>Err\},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{N}_{G_j}=\hspace{0.33em}|G_j|,$$

(23)

where $G_j$ is a set of satellites with $\left|e_{ij}^{}\right|\hspace{0.33em}>Err$, and $N_{G_j}$ is the number of satellites in $G_j$. With modification (23) the proposed control becomes more sensitive to satellites which are not converged to reference trajectories.

3.5. Reference Trajectories Assignment Problem

Since the number of satellites in the group is $N$, there are $N!$ possible ways to assign to each satellite the reference trajectory according to its position in the pixel image formed by the satellites. Suppose that initial satellite positions in the LVLH reference frame are represented by vectors $r_i\left(t_0\right),\hspace{0.33em}i=1\dots N$, and the corresponding state vectors in curvilinear coordinates are expressed by $x_{ij}(t_0)$. The trajectory assignment is performed to minimize the deviation of the initial state vector from the state vector $x_{ij}^{des}(t_0)$ corresponding to the possible reference trajectory.

First, a satellite closest to the center of the formation at the time $t_0$ is chosen. The selection is performed according to the solution of the following optimization problem:

$$s=\arg \underset{i}{\min }\left|r_i(t_0)-{\displaystyle \sum _{i=1}^N{r}_i(t_0)/N}\right|.$$

(24)

The $s$-th satellite is to be assigned to the $p$-th reference trajectory, which is selected as the center of image rotation. The $p$-th satellite position is the nearest point to the image center according to the following:

$$p=\arg \underset{i}{\min }\left|r_{1i}^{\mathrm{des}}-{\displaystyle \sum _{j=2}^M{r}_{1j}^{\mathrm{des}}/M}\right|,$$

(25)

where $r_{1i}^{\mathrm{des}}$ is the radius-vector of the reference trajectories relative to the first satellite (which is chosen for unambiguity), $M=N-1$.

Next, consider the relative state vectors for all satellites $x_{si}^{}(t_0),\hspace{0.33em}i=1\dots N$ ($x_{ss}^{}(t_0)=0$) and state vectors according to the reference trajectories $x_{pj}^{\mathrm{des}}(t_0),\hspace{0.33em}j=1\dots N$ (taking into account that $x_{pp}^{\mathrm{des}}(t_0)=0$). The following cost function is proposed for each pair of $i$-th satellite and $j$-th reference trajectory:

$${\Delta }_{ij}=\sqrt{{\left(x_{pj}^{\mathrm{des}}(t_0)-x_{si}^{}(t_0)\right)}^{\mathrm{T}}\tilde{Q}\left(x_{pj}^{\mathrm{des}}(t_0)-x_{si}^{}(t_0)\right)},$$

(26)

where $\tilde{Q}$ is a positive definite matrix. The satellite assignment problem then can be formulated as the following optimization problem:

$$j(i)=\arg \underset{j(i)}{\min }{\displaystyle \sum _{i=1}^N{\Delta }_{ij(i)}}.$$

(27)

This problem can be reformulated as a linear programming problem as follows. Consider a matrix $X$ whose elements can be equal to 0 or 1, and there is only one “1” in every line and column. This “1” in the $i$-th line and $j$-th column means that $i$-th satellite is assigned to the $j$-th reference trajectory. Then, the matrix $X$ can be transformed into the vector $\tilde{x}$ according to:

$$\tilde{x}={\left[X_{11},\dots ,X_{1\mathrm{N}},X_{21},\dots ,X_{2\mathrm{N}},\dots \right]}_{N^2}^{\mathrm{T}}$$

The assignment problem (27) is equivalent to the following linear programming problem with two constraints:

$$\tilde{x}=\arg \underset{\tilde{x}}{\min }\hspace{0.17em}\left(c_{}^{\mathrm{T}}\tilde{x}\right),\hspace{0.17em}\hspace{0.17em}\tilde{A}\tilde{x}=b,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{e}_1\le \tilde{x}\le e_2,$$

(28)

where

$$\begin{gathered} c={\left[{\Delta }_{11},\dots ,{\Delta }_{1\mathrm{N}},{\Delta }_{21},\dots ,{\Delta }_{2\mathrm{N}},\dots \right]}_{N^2}^{\mathrm{T}},b={\left[1\dots 1\right]}_{2N}^{\mathrm{T}},\hspace{0.17em}\hspace{0.17em}{e}_1={\left[0\dots 0\right]}_{N^2}^{\mathrm{T}},\hspace{0.17em}\hspace{0.17em}{e}_2={\left[1\dots 1\right]}_{N^2}^{\mathrm{T}},\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{0.33em}\stackrel{N}{\overbrace{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}}}\hspace{0.33em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}\stackrel{N}{\overbrace{\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\hspace{1em}}}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}\stackrel{N}{\overbrace{\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\hspace{1em}}} \\ \tilde{A}=\left(\begin{array}{ccccccccccccc}1& 1& \dots & 1& 0& 0& \dots & 0& \dots & 0& 0& \dots & 0\\ 0& 0& \dots & 0& 1& 1& \dots & 1& \dots & 0& 0& \dots & 0\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0& 0& \dots & 0& 0& 0& \dots & 0& \dots & 1& 1& \dots & 1\\ 1& 0& \dots & 0& 1& 0& \dots & 0& \dots & 1& 0& \dots & 0\\ 0& 1& \dots & 0& 0& 1& \dots & 0& \dots & 0& 1& \dots & 0\\ 0& 0& \ddots & 0& 0& 0& \ddots & 0& \dots & 0& 0& \ddots & 0\\ 0& 0& \dots & 1& 0& 0& \dots & 1& \dots & 0& 0& \dots & 1\end{array}\right)\begin{array}{c}\left.\begin{array}{c}\\ \\ \\ \end{array}\right\}N\\ \left.\begin{array}{c}\\ \\ \\ \end{array}\right\}N\end{array}.\hspace{1em}\hspace{1em}\end{array} \end{gathered}$$

The first constraint $\tilde{A}\tilde{x}=b$ in (28) corresponds to the requirement that only one reference trajectory can be assigned to each satellite, and the second constraint $e_1\le \tilde{x}\le e_2$ is the requirement that the vector $\tilde{x}$ can contain only 0 or 1. The solution vector $\tilde{x}$ can be found numerically by any of the linear programming solver methods [31].

It should be noted that the computational cost of the numerical solution of the assignment problem is high, and it is possible that this problem cannot be solved in orbit with onboard computers. Thus, it is assumed that the problem (28) is to be solved on the ground before the image construction and each time before the reconfiguration. After the assignment of a particular reference trajectory to each satellite, new reassignment is not recalculated during the convergence, even if the sum of the cost functions (26) on the trajectory could eventually be less than the initial value.

3.6. Control Implementation Using Aerodynamic Force

The calculated control (23) is assumed to be implemented by the aerodynamic force (9). Consider the normal vector $n_i$ to the $i$-th satellite flat surface, which is defined in LVLH reference frame by two angles $\vartheta \in [0,\hspace{0.17em}\hspace{0.17em}\pi /2],\hspace{0.17em}\hspace{0.17em}\psi \in [0,\hspace{0.17em}\hspace{0.17em}2\pi ]$, as shown in Figure 4:

$$n_i={[\sin {\vartheta }_i\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\cos {\vartheta }_i\cos {\psi }_i\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\cos {\vartheta }_i\sin {\psi }_i]}^{\mathrm{T}}.$$

Each flat surface has two normal vectors, and the normal vector, which has a nonnegative projection to the along-track component $(n_i,e_x)\ge 0$, is considered. Assuming that the unit vector along the velocity $e_V$ is aligned with the $e_x$ axis, the aerodynamic acceleration vector (9) can be rewritten as follows:

$$f_i=k\left[\begin{array}{c}-2\epsilon {(\sin {\vartheta }_i)}^3+\eta (\epsilon -1){(\sin {\vartheta }_i)}^2+(\epsilon -1)\sin {\vartheta }_i\\ -\cos {\vartheta }_i\sin {\vartheta }_i(\eta -\epsilon \eta +2\epsilon \sin {\vartheta }_i)\cos {\psi }_i\\ -\cos {\vartheta }_i\sin {\vartheta }_i(\eta -\epsilon \eta +2\epsilon \sin {\vartheta }_i)\sin {\psi }_i\end{array}\right],$$

(29)

where $k=\frac{1}{m_i}\rho V_i^2{S}_i$ is a coefficient considered as constant for a given satellite in a given orbit. In this paper, the parameters $\epsilon$ and $\eta$ are set as $\epsilon =0.1$, $\eta =0.1$ according to estimates by Beletsky in [29]. From Equation (29), it can be concluded that the components of the aerodynamic acceleration are bounded. The $f_{}^x$ component cannot be positive. The maximum value of the along-track component $f_{\max }^x$ corresponding to aerodynamic drag at ${\vartheta }_i=0$ is equal to $f_{\max }^x=-1.19k$. At the same time at ${\vartheta }_i=0$, the other two components are zero $f_{}^y=f_{}^z=0$. The maximal value of the projection of the vector $f_i$ onto the $Oyz$ plane is $f_{\max }^{yz}\approx 0.12k$, which is achieved at about ${\vartheta }_i=52\mathrm{deg}$. At ${\vartheta }_i=52\mathrm{deg}$, the value of the $f_{}^x=-0.86k$. The maximum aerodynamic drag value $f_{\max }^x$ is about 10 times higher than the maximum value of the “lift” components $f_{\max }^{yz}$. So, the calculated control vector (23), in the common case, cannot be directly implemented using (29). The following implementation strategy is used in this paper. The implemented control acceleration $f_i$ depends on the calculated control vector $u_i$ as follows [12]:

$$f_i=\left\{\begin{array}{l}{f}_{\max }^x,\hspace{1em}\mathrm{if} u_i^x<f_{\max }^x;\\ f_{\max }^{yz},\hspace{1em}\mathrm{if} \hspace{0.17em}{f}_{\max }^x\le u_i^x<0 \mathrm{and} \sqrt{{(u_i^y)}^2+{(u_i^z)}^2}>f_{\max }^{yz};\\ u_i,\hspace{1em}\mathrm{if} \hspace{0.17em}{f}_{\max }^x\le u_i^x<0 \mathrm{and} \sqrt{{(u_i^y)}^2+{(u_i^z)}^2}\le f_{\max }^{yz};\\ 0,\hspace{1em}\mathrm{if} u_i^x\ge 0;\end{array}\right.$$

(30)

where $f_{\max }^x={[f_{\max }^x\hspace{0.33em}\hspace{0.33em}0\hspace{0.33em}\hspace{0.33em}0]}^T,\hspace{0.33em}{u}_{\max }^{yz}=k{[u_i^x/k\hspace{0.33em}\hspace{0.33em}{u}_i^y/f_{\max }^{yz}\hspace{0.33em}\hspace{0.33em}{u}_i^z/f_{\max }^{yz}]}^T$.

The implementation strategy (30) prioritizes applying the available value of the aerodynamic drag component first. This results in eliminating the relative drift and relative shift of the trajectories. Then, in cases where the required drag component is less than the maximum value, but the required lift component is high, the maximum available value of the lift component is to be implemented as the second priority. If both required lift and drag components are less than the maximum values, then the control acceleration vector $u_i$ can be implemented. In [12], it is shown that even in such cases, it may not be possible to find a pair of angles $\vartheta ,\hspace{0.17em}\hspace{0.17em}\psi$ to exactly implement the vector $u_i$. For this situation, the closest possible vector to the required $u_i$ is implemented. If the component of the calculated control $u_i^x$ is positive, then no control is applied, i.e., the normal vector to the sunlight reflector $n_i$ is perpendicular to the velocity vector.

The required sunlight reflectors’ attitude providing the required values of angles $\vartheta ,\hspace{0.17em}\hspace{0.17em}\psi$ is assumed to be implemented by the onboard attitude control system, for example, by reaction wheels control. In this paper, the attitude-implementation errors are not considered, the required attitude is perfectly realized. However, for attitude calculation, the atmospheric density in the model (29) is assumed to be constant, which results in errors in aerodynamic acceleration implementation, because the GOST upper atmosphere density model [32] is used in the numerical simulations.

4. Numerical Simulation Results

The proposed control of satellite formation flying with sunlight reflectors is based on many assumptions, and there hardly be any analytical proof that the control algorithm would unfailingly result in the achievement of the required relative configuration. The algorithm performance study is conducted using the numerical simulations of the controlled motion of the satellites.

In numerical simulation, the initial conditions for orbital motion are determined by the cluster launch. It is assumed that the separation velocity $V_e$ from the launch vehicle is aligned with orbital velocity, and the satellites are ejected from the launch containers with a fixed time interval. The random launch velocity error $\delta V_e$ is modeled as a normally multivariate distributed vector with a standard deviation of ${\sigma }_{\delta V_e}$ without correlation in component values.

In the equations of motion (7), the $J_2$ disturbance is taken into account and the atmospheric density is given by the GOST model [32]. The integration of Equation (7) is carried out by the 4th-order Runge–Kutta method with a constant integration step. The simulation control loop is presented in Figure 5. The inputs are the initial conditions for the satellite motion, the image, and the control parameters. The image to be demonstrated is discretized and the reference trajectory is assigned to each satellite according to Equation (28). Then, for each pair of satellites, the deviation $e_{ij}$ and the corresponding acceleration vectors $w_{ij}$ by (16) and $u_{ij}$ by (21) are calculated. Then, the modification (23) is applied to calculate $u_i$. Using the implementation strategy (30), the aerodynamic acceleration $f_i^{aero}$ is calculated, which is used for each satellite orbital motion integration.

The following mission scenario is considered. After launching in a near-circular orbit according to the initial conditions, the reference trajectory assignment problem is to be solved. The example image to be constructed is chosen as that consisting of three English alphabet letters. Each letter consists of a certain number of satellites with sunlight reflectors, and the distance between the neighbor satellites is set by 700 m, which is in accordance with the visibility conditions calculated in [5]: The magnitude of the reflectors is about −2.8 and the minimum distance is 380 m at orbit altitude of 350 km. After the convergence to the reference trajectories at a given point with proper lighting conditions, the attitudes of all the satellites are synchronized, and the sunlight is reflected by the satellites to be visible to the Earth observer. After the demonstration of the first image (three letters), the satellite’s relative positions are to be reconfigured to compose another image (three different letters).

The main parameters of the simulated motion scenario are provided in

Table 1. Main parameters of numerical simulations.

Parameter	Value
First letters	“ABC”
Second letters	“DEF”
Number of satellites	60
Distance between satellites	700 m
Mass of satellites	18 kg
Size of sail-reflector	2 × 2 m
Time between separations	20 s
Separation velocity	1.5 m/s
Standard deviation of separation velocity error	0.1 m/s
Initial orbit height	350 km
Orbit inclination	51.7°
Simulation time	30 + 30 h
Integration step	60 s

. The LQR weight matrices are selected as follows:

$$Q=\left(\begin{array}{cccccc}10a_0^2& 0& 0& 0& 0& 0\\ 0& 10a_0^2& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 10a_0^2& 0& 0\\ 0& 0& 0& 0& 10a_0^2& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\hspace{1em}R=\left(\begin{array}{ccc}{10}^{14}{a}_0^2& 0& 0\\ 0& {10}^{15}{a}_0^2& 0\\ 0& 0& {10}^{15}\end{array}\right).$$

Values of elements in these matrices should ensure the absence of saturation in limited aerodynamic control. The drag force is about 10 times greater than the lift force; therefore, the elements of $R$ for the control component ${\ddot{\phi }}^w$ should be less than those for ${\ddot{\theta }}^w$ and ${\ddot{\rho }}^w$. The values of the elements of the matrix $R$ are selected in order to avoid the saturation of the aerodynamic value control. In matrix Q, values corresponding to the deviation of $\phi$ and $\theta$ are greater than those of $\rho$ and its derivatives. Angular coordinates are multiplied by the squared major semiaxis $a_0^{}$. The value of the acceptable deviation $\left|e_{ij}^{}\right|\hspace{0.33em}>Err$ for the modification (23) is set as $Err=100$.

The parameters of the GOST model of atmospheric density are in accordance with the values for the launch date of 1 March 2012, which corresponds to the high solar activity period.

4.1. Example of Numerical Simulation

The simulation example trajectories are visualized in a video available online in [33]. Figure 6 demonstrates the arcs of satellite relative trajectories and the satellite positions forming the letters after the first 30 h of simulation and then 30 h after the reconfiguration. It can be seen that almost all the satellites have achieved the required relative circular trajectories in projection on the $Oxy$ plane in the LVLH reference frame, and the satellites are almost at equal relative distances forming the letters. In Figure 7a, the average trajectory deviations from the relative reference trajectory are presented. It can be concluded that after the first 30 h of simulation of relative motion, almost all the trajectories have converged to the reference ones. Oscillations of about 50 m in amplitude are caused by $J_2$ disturbances, which are not accounted for in the control algorithm, and errors in the aerodynamic force implementation due to uncertainties in the atmospheric density. After 30 h, the new relative reference trajectories have been assigned to the satellites to form the second set of letters “DEF”. This is why the deviations’ values in Figure 7a change abruptly at that moment, although they nearly converge to zero within the following 15 h, which is faster than the initial convergence to “ABC” letters after the launch. Figure 7b demonstrates the particular trajectory deviation of the first satellite relative to the second satellite.

Figure 8 shows the example of the control acceleration vector components calculated by Equation (23) and the aerodynamic force acceleration vector implemented by Equation (30) for the first satellite in the formation. The values are in dimensionless units, divided by the parameter k of the aerodynamic model (29). It can be seen that the along-track component is not saturated during the motion simulation, which is the result of correctly selected elements of the matrix $R$. Nevertheless, saturation is observed for the out-of-plane component $f_y$. According to the implementation logic (29), in case of saturation, the maximum value of the available “lift” component is applied. When the deviations of the trajectory become small, the saturation disappears, and the calculated control is implemented almost as it is.

Figure 9 shows the evolution of the orbit altitude and atmospheric density according to the GOST model of the atmosphere for the considered example. It can be concluded that during the convergence, the mean value of the orbit height has a faster decay compared to the converged relative motion. The height oscillations are caused by the $J_2$ disturbance influence. The orbit height change $\Delta H$ is the cost of the aerodynamic-based maneuvering. Each maneuver results in a shortening of the mission’s lifetime. For this reason, it is important to study the performance of the proposed control in terms of the resulting $\Delta H$ for achieving the required configuration.

4.2. Monte-Carlo Performance Study

The example above shows that the proposed aerodynamic-based control can establish the required relative trajectories and form images by means of satellites with sunlight reflectors. However, the considered system depends on multiple parameters, and the performance study of the controlled motion is carried out using the Monte Carlo approach. For this study, two characteristics are considered: the change in the orbit height $\Delta H$ of the formation during the controlled motion and the convergence time taken to achieve the required relative trajectories. The formation is assumed to be converged if the mean error of trajectory deviation is less than 50 m during one orbital period. The single-image demonstration scenario without reconfiguration as in the previous section is simulated.

One of the important system parameters is the error of the separation velocity, which is considered to be a random normally distributed variable with standard deviation ${\sigma }_{\delta V_e}$. The results of 20 numerical simulations for each random initial condition set with a defined ${\sigma }_{\delta V_e}$ are analyzed. For each simulation run, the mean value of $\Delta H$ and the convergence time are calculated. In Figure 10, the results of the experiment are presented in boxplots. Each box contains 50% of the results, the red line corresponds to the mean value, above and below the box are 25% of the results, and red crosses are outliers. As one can see from Figure 10, the greater the errors in the launch velocity, the longer the convergence time and the broader the box of the change in the orbit height $\Delta H$. It is interesting that with random errors in launch velocity at some simulations, the resulting $\Delta H$ is less than in the case of low error values of 0.1 m/s. It can be explained that at some conditions the errors lead to less initial relative orbit deviations—that is why even the mean value of $\Delta H$ is 1 km less than $\Delta H\approx 14$ km at ${\sigma }_{\delta V_e}=0.1$ m/s. Nevertheless, the convergence time mean value is evenly increasing with ${\sigma }_{\delta V_e}$.

The lower the orbit height, the higher the atmospheric density and the greater the magnitude of the aerodynamic control force. Other system parameters being the same, the variations in the initial orbit height significantly influence the performance of the formation flying controlled motion. With the fixed value of the launch error ${\sigma }_{\delta V_e}=0.1$ m/s, the initial circular orbit height is varied, and the Monte Carlo simulations are carried out. Figure 11 demonstrates the results of the study in boxplots of the random values of $\Delta H$ and convergence time. As expected, the convergence time is much shorter at lower orbit height. The change in the height $\Delta H\approx 6.5$ km if the initial orbit height is 300 km, although the lifetime of the formation at such low orbits could only be a couple of weeks. The convergence time at an orbit height of 400 km is dramatically longer and amounts to nearly 155 h, although the change in the orbit height is about the same as at an orbit height of 350 km.

The performance study showed that the proposed aerodynamic control for the formation with sunlight reflectors is quite effective at orbit altitudes lower than 400 km, although the change in the orbit height during the required image configuration results in a very limited lifetime of such formation of less than about a month without propulsion. This kind of control can be applied when satellites run out of fuel to extend mission lifetime.

5. Conclusions

The proposed decentralized control using aerodynamic forces allows the satellites to achieve the required relative satellite configuration, as demonstrated for a particular case in this paper. The performance numerical study showed that in the case of 2 × 2 m² sunlight reflectors, used as aerodynamic wings in LEO, the satellites converge to the reference trajectories in about 30 h by changing their attitude relative to the incoming airflow at the initial orbit height of 350 km. On the one hand, the change in the orbit height during the maneuvering leads to accelerated orbit decay, and the lifetime of a mission using only aerodynamic control is limited by about a month. On the other hand, this can be considered as an opportunity to extend a regular mission lifetime after some of its satellites have run out of fuel. As a continuation of the presented study, the controlled motion performance should be investigated, taking into account disturbances caused by the uncertainties in atmospheric density, control implementation errors, and relative navigation features.