## 2 Problem Statement and Equations of Motion
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{\varphi}_{ij}\,+\,2\omega\dot{\rho}_{ij}/a_{0} \,=\,w_{\varphi}^{ij},\\ \ddot{\theta}_{ij}\,+\,\omega^{2}\theta_{ij}\,=\,w_{\theta}^{ij},\\ \ddot{\rho}_{ij}\,-\,2\omega a_{0}\dot{\varphi}_{ij}\,-\,3\omega^{2}\rho_{ij} \,=\,w_{\rho}^{ij}.\end{array} \right. \tag{4}\]  
Orbital plane of the \(i\)-th satellite  
Figure 2: Curvilinear coordinates definition. The solution to the free motion Equation (4) without control is similar to Equation (2) and is given by:  
\[\left\{ \begin{array}{l}\varphi_{ij}(t)=-3\widetilde{C}_{1}^{ij}\omega t+2 \widetilde{C}_{2}^{ij}\cos(\omega t)-2\widetilde{C}_{3}^{ij}\sin(\omega t)+ \widetilde{C}_{4}^{ij},\\ \theta_{ij}(t)=\widetilde{C}_{5}^{ij}\sin(\omega t)+\widetilde{C}_{6}^{ij} \cos(\omega t),\\ \rho_{ij}(t)=2\widetilde{C}_{1}^{ij}a_{0}+\widetilde{C}_{2}^{ij}a_{0}\sin( \omega t)+\widetilde{C}_{3}^{ij}a_{0}\cos(\omega t),\end{array} \right. \tag{5}\]  
where constants \(\widetilde{C}_{1}^{ij}-\widetilde{C}_{6}^{ij}\) depend on the initial conditions at time \(t_{0}=0\) similarly to Equation (3) as follows:  
\[ \begin{array}{l}\widetilde{C}_{1}^{ij}=\frac{\dot{\varphi}_{ij}(0)}{\omega} +2\frac{\rho_{ij}(0)}{a_{0}},\ \widetilde{C}_{2}^{ij}=\frac{1}{\omega}\frac{\dot{\rho}_{ij}(0)}{a_{0}},\ \widetilde{C}_{3}^{ij}=-2\frac{\dot{\varphi}_{ij}(0)}{\omega}-3\frac{\dot{ \rho}_{ij}(0)}{\theta_{0}},\\ \widetilde{C}_{4}^{ij}=-\frac{2}{\omega}\frac{\dot{\rho}_{ij}(0)}{\theta_{0}} +\varphi_{ij}(0),\ \widetilde{C}_{5}^{ij}=\frac{\dot{\theta}_{ij}(0)}{\omega},\ \widetilde{C}_{6}^{ij}=\theta_{ij}(0).\end{array}  \tag{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{\mathbf{R}}_{i}=-\frac{\mu_{g}\mathbf{R}_{i}}{R_{i}^{3}}+\mathbf{f}_{i }^{\text{grav}}+\mathbf{f}_{i}^{\text{aero}}, \tag{7}\]  
where \(\mu_{g}=3.986\cdot 10^{14}\text{m}^{3}/\text{s}^{2}\) is the Earth gravitational parameter, \(\mathbf{R}_{i}=[X_{i},Y_{i},Z_{i}]^{T}\) is the radius-vector of the \(i\)-th satellite in the Earth-centered inertial reference frame (ECI), and \(\mathbf{f}_{i}^{\text{grav}}\) is the vector of acceleration caused by the second zonal harmonic of Earth gravitational potential:  
\[\mathbf{f}_{i}^{\text{grav}}= \begin{pmatrix}\frac{\delta}{R_{i}^{3}}X_{i} \Big{(}\frac{5Z_{i}^{2}}{R_{i}^{2}}-1\Big{)}\\ \frac{\delta}{R_{i}^{3}}Y_{i}\Big{(}\frac{5Z_{i}^{2}}{R_{i}^{2}}-1\Big{)}\\ \frac{\delta}{R_{i}^{3}}Z_{i}\Big{(}\frac{5Z_{i}^{2}}{R_{i}^{2}}-3\Big{)}\end{pmatrix}