## 2 Modeling the Orbital Dynamics  
In the geocentric equatorial inertial coordinate systems, the orbit motion of a satellite can be reduced to a perturbed two-body problem, and the motion equation is a set of ordinary differential equations [23]  
\[\begin{cases}\ddot{r}=F_{0}(r)+F_{k}(r,\dot{r},t;\varepsilon),\\ r(t_{0})=r_{0},\ \dot{r}(t_{0})=\dot{r}_{0},\end{cases} \tag{1}\]  
where \(r\) and \(\dot{r}\) are the position vectors and velocity vectors of the satellite, \(t\) is the propagation time, \(F_{0}\) is the two-body motion acceleration of the central body, \(F_{k}\) is the perturbation acceleration, and \(\varepsilon\) is a small parameter. Utilizing asymptotic expansion and other approaches, the power series form solution of small parameter \(\varepsilon\) can be constructed [24], such as  
\[\kappa(t)=\kappa^{(0)}(t)+\varepsilon\beta(t)+\varepsilon^{2}\beta^{2}(t)+\dots+ \varepsilon^{l}\beta^{l}(t), \tag{2}\]  
where \(\kappa=(a,e,i,\Omega,\omega,M)^{\mathrm{T}}\) are the classical orbital elements: semimajor axis, eccentricity, inclination, right ascension of the ascending node, argument of perigee, and mean anomaly at epoch, respectively. By substituting the formal solution into the original perturbed motion equation and comparing the coefficients of the same power \(\varepsilon^{l}\) at both ends, the perturbation solutions of each order can be obtained. By the average approach, the perturbation solutions can be divided into secular terms, long-period terms, and short-period terms according to forms of presentation.  
\[\begin{cases}\kappa(t)=f(\overline{\kappa}_{0},\,t)=\overline{\kappa}(t)+ \kappa_{\mathrm{lp}}(t)+\kappa_{\mathrm{sp}}(t),\\ \overline{\kappa}(t)=\overline{\kappa}_{0}+(\delta n+\kappa_{\mathrm{sec}}^{(1 )}+\kappa_{\mathrm{sec}}^{(2)}+\dots)(t-t_{0}).\end{cases} \tag{3}\]  
Among them, \(t_{0}\) is the epoch time, \(\kappa(t)\) are the osculating elements at time \(\kappa(t)\), \(\overline{\kappa}_{0}\) is the initial mean elements, \(\overline{\kappa}(t)\) is the mean elements at time \(t\), \(\delta=[0,0,0,0,0,1]^{T}\), \(n\) is the mean motion of satellite in two body, \(\kappa_{\mathrm{sec}}(t-t_{0})\) makes up the secular terms, \(\kappa_{\mathrm{lp}}(t)\) consists of the long period terms, and \(\kappa_{\mathrm{sp}}(t)\) consists of the short-period terms. Besides this, Equation (3) describes the identity mapping between mean elements and the osculating elements, which can be solved by numerical algorithm such as Newton iteration [25].  
LEO constellations are mainly affected by atmospheric drag and non-spherical Earth gravity. The first-order and second-order secular term of \(J_{2}\), the first-order long-period term of \(J_{2}\), the first-order short-period term of \(J_{2}\), and the secular term of atmospheric drag are involved in the semi-analytical dynamic model. It should be noted that in order to avoid the abnormal matrix caused by the large numerical difference among the Keplerian elements, distance and time variables need to be normalized by \([R_{e}]\) and \([(Re^{3}/\mu)^{1/2}]\), where \(R_{e}\) is the Earth's equatorial radius and \(\mu\) is the Earth's gravitational constant. Thus, the mean motion \(n\) is expressed by \(a^{-3/2}\). Detailed semi-analytical solutions are recorded in [26].