## 2 Relative Motion Dynamics Model
### Formation Dynamics Modeling
 \tag{11}\]  
where \(s_{\theta}=\sin(o),c_{o}=\cos(o)\), \(s_{\theta}=\sin(o),c_{o}=\cos(o)\), \(\Omega\) is the right ascension of the satellite ascending node, \(\theta\) is the argument of latitude, and \(i\) is the inclination of the formation mass center. Given that \(\mathbf{T}\) represents the orthogonal matrix, hence \(\mathbf{T}^{\text{T}}\mathbf{T}=\mathbf{I}_{3}\), differentiating both sides of this equation yields \(\dot{\mathbf{T}}^{\text{T}}\mathbf{T}+\mathbf{T}^{\text{T}}\dot{\mathbf{T}}=\mathbf{0}\), confirming its skew-symmetric matrix property. Thus, we designate the skew-symmetric matrix function as \(\mathbf{S}=\dot{\mathbf{T}}^{\text{T}}\mathbf{T}\). Subsequently, the relationship between the ECI and Hill frames can be articulated. \[\mathbf{Q}_{a}=\mathbf{T}(\mathbf{q}_{0}+\mathbf{q}_{a}) \tag{12}\]  
As a result, the kinetic energy of the spacecraft within the Hill frame can be delineated as follows:  
\[\mathcal{T}(\mathbf{q}_{a},\dot{\mathbf{q}}_{a})=\frac{1}{2}m_{a}\big{\|}\dot{\mathbf{Q}}_{ a}\big{\|}^{2}=\frac{1}{2}m_{a}\big{\|}(\dot{\mathbf{q}}_{0}+\dot{\mathbf{q}}_{a})+ \mathbf{S}(\mathbf{q}_{0}+\mathbf{q}_{a})\big{\|}^{2} \tag{13}\]  
The gravitational potential energy of the spacecraft can be described as  
\[\mathcal{U}_{\mathrm{g}}(\mathbf{q}_{a})=-\frac{\mu m_{a}}{R_{a}}\bigg{[}1-\frac{ l_{2}}{2}\frac{R_{e}^{2}}{R_{a}^{2}}\bigg{(}3\frac{\left(\mathbf{T}_{3}(\mathbf{q}_{0}+ \mathbf{q}_{a})\right)^{2}}{R_{a}^{2}}-1\bigg{)}\bigg{]} \tag{14}\]  
Putting Equations (5), (13), (14) into Equation (10), the Lagrangian equation can be written as  
\[ \begin{split}\mathcal{L}(\mathbf{q}_{a},\dot{\mathbf{q}}_{a})& =\mathcal{T}-\mathcal{U}_{g}-\mathcal{U}_{e}\\ &=\frac{1}{2}m_{a}\big{\|}(\dot{\mathbf{q}}_{0}+\dot{\mathbf{q}}_{a})+\bm {S}(\mathbf{q}_{0}+\mathbf{q}_{a})\big{\|}^{2}\\ &+\frac{\mu m_{a}}{R_{a}}\bigg{[}1-\frac{l_{2}}{2}\frac{R_{e}^{2} }{R_{a}^{2}}\bigg{(}3\frac{\left(\mathbf{T}_{3}(\mathbf{q}_{0}+\mathbf{q}_{a})\right)^{2} }{R_{a}^{2}}-1\bigg{)}\bigg{]}\\ &+\frac{\mu_{0}}{4\pi}\mathbf{\mu}_{a}^{\mathsf{T}}\bigg{(}\frac{3 \mathbf{\mu}_{b}^{\mathsf{T}}(1+k)\mathbf{q}_{a}}{\big{\|}(1+k)\mathbf{q}_{a}\big{\|}^{5}} (1+k)\mathbf{q}_{a}-\frac{\mathbf{\mu}_{b}}{\big{\|}(1+k)\mathbf{q}_{a}\big{\|}^{3}}\bigg{)} \end{split}  \tag{15}\]  
Considering the focus on the motion of the mass center and follower spacecraft within the Hill frame, \(\mathbf{q}_{a}\) is selected as the generalized coordinates. Following the Legendre transformation, the corresponding generalized momentum is derived. \[\mathbf{p}_{a}=\frac{\partial\mathcal{L}}{\partial\mathbf{q}_{a}}=m_{a}[\dot{\mathbf{q}}_{ 0}+\dot{\mathbf{q}}_{a}+\mathbf{S}(\dot{\mathbf{q}}_{0}+\dot{\mathbf{q}}_{a})] \tag{16}\]  
Consequently, the dynamics of relative motion are determined by the Euler-Lagrange equation. \[ \begin{split}&\frac{\mathrm{d}}{\mathrm{d}t}\bigg{(}\frac{ \partial\mathcal{L}}{\partial\dot{\mathbf{q}}_{a}}\bigg{)}-\frac{\partial\mathcal{ L}}{\partial\mathbf{q}_{a}}=\mathbf{u}\\ &\Rightarrow m_{a}[\ddot{\mathbf{q}}_{0}+\ddot{\mathbf{q}}_{a}+\hat{\mathbf{S }}(\mathbf{q}_{0}+\mathbf{q}_{a})+\mathbf{S}(\dot{\mathbf{q}}_{0}+\dot{\mathbf{q}}_{a})]\\ &-m_{a}\mathbf{S}^{\mathsf{T}}[(\dot{\mathbf{q}}_{0}+\dot{\mathbf{q}}_{a})- \mathbf{S}(\mathbf{q}_{0}+\mathbf{q}_{a})]+\
abla_{\mathbf{q}_{a}}\mathcal{U}_{g}+\
abla_{\mathbf{ q}_{a}}\mathcal{U}_{e}=\mathbf{u}\end{split}  \tag{17}\]  
Rectifying the equation yields  
\[\begin{split}\bigg{[} \begin{array}{c}\dot{\mathbf{q}}_{a}\\ \dot{\mathbf{p}}_{a}\end{array} \bigg{]}=\bigg{[} \begin{array}{c}\frac{1}{m_{a}} \mathbf{p}_{a}-\dot{\mathbf{q}}_{0}-\mathbf{S}(\mathbf{q}_{0}+\mathbf{q}_{a})\\ m_{a}\mathbf{S}^{\mathsf{T}}[(\dot{\mathbf{q}}_{0}+\dot{\mathbf{q}}_{a})-\mathbf{S}(\mathbf{q}_{0}+ \mathbf{q}_{a})]-\
abla_{\mathbf{q}_{a}}\mathcal{U}_{g}-\
abla_{\mathbf{q}_{a}}\mathcal{U} _{e}\end{array} \bigg{]}+\bigg{[} \begin{array}{c}\mathbf{0}_{3}\\ \mathbf{I}_{3}\end{array} \bigg{]}\mathbf{u}\\