## 4 Numerical Simulation Results
### 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,  
\begin{table}
\begin{tabular}{l c} \hline \hline
**Parameter** & **Value** \\ \hline 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 \(\times\) 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\({}^{\circ}\) \\ Simulation time & 30 + 30 h \\ Integration step & 60 s \\ \hline \hline \end{tabular}
\end{table}
Table 1: Main parameters of numerical simulations.  
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 \(\mathbf{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 6: Two sets of demonstrated letters formed by satellites with sunlight reflectors after 30 h of simulation from initial conditions (**a**) and the relative position after 30 h after the reconfiguration (**b**).  
Figure 7: Average relative trajectory error for all the satellites (**a**) and for the first satellite (**b**).  
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.