## 3 Nonlinear Estimation of the Perturbed Keplerian Parameters  
Let us consider a scenario in which \(m\) satellites in MEO are all equipped with two terminals for ISLs. The constellation is maintained fully connected via a link scheduling that accounts for the satellite positions, orientations and terminal field of regards; see details \begin{table}
\begin{tabular}{l l} \hline \hline \multicolumn{2}{c}{**(a) Ephemeris data definition.**} \\ \hline  & **Parameter Description** \\ \hline \(M_{0}\) & Mean anomaly \\ \(\Delta n\) & Mean motion correction \\ \(e\) & Eccentricity \\ \(\sqrt{a}\) & Square root of semi-major axis \\ \(\Omega_{0}\) & Longitude of ascending node \\ \(\Omega\) & Rate of change in longitude of ascending node \\ \(i_{0}\) & Inclination angle \\ \(i\) & Rate of change in inclination angle \\ \(\omega\) & Argument of perigee \\ \(C_{uc}\) & Amplitude of cosine correction to the argument of latitude \\ \(C_{us}\) & Amplitude of sine correction to the argument of latitude \\ \(C_{rc}\) & Amplitude of cosine correction to the orbit radius \\ \(C_{rs}\) & Amplitude of sine correction to the orbit radius \\ \(C_{ic}\) & Amplitude of cosine correction to the inclination \\ \(C_{is}\) & Amplitude of sine correction to the inclination \\ \hline \multicolumn{2}{c}{**(b) The user algorithm.**} \\ \hline \hline \multicolumn{2}{c}{**Computation**} \\ \hline \(a=(\sqrt{a})^{2}\) \\ \(n=\sqrt{\dfrac{\mu}{a^{3}}}+\Delta n\) \\ \(\Delta t=t-t_{0}\epsilon\) \\ \(M=M_{0}+n\Delta t\) \\ \(M=E-e\sin E\) \\ \(\sin\
u=\dfrac{\sqrt{1-e^{2}}\sin E}{1-\epsilon\cos E}\) \\ \(\cos\
u=\dfrac{\cos E-e}{1-\epsilon\cos E}\) \\ \(\phi=\
u+\omega\) \\ \(\Omega=\Omega_{0}+\left(\hat{\Omega}-\omega_{\mathrm{c}}\right)\Delta t-\omega_ {\mathrm{c}}t_{0\mathrm{c}}\) \\ \(u=\phi+C_{us}\sin(2\phi)+C_{uc}\cos(2\phi)\) \\ \(r_{\mathrm{s}}=a\cdot(1-e\cos E)+C_{rs}\sin(2\phi)+C_{rc}\cos(2\phi)\) \\ \(i=i_{0}+i\cdot\Delta t+C_{is}\sin(2\phi)+C_{ic}\cos(2\phi)\) \\ \(\left\{s_{\mathrm{r}}=r_{\mathrm{s}}\cos\iota\cos\Omega-r_{s}\sin\iota\cos \iota\sin\Omega\right\}\) \\ \(\left\{y_{\mathrm{s}}=r_{\mathrm{s}}\cos\iota\sin\Omega+r_{s}\sin\iota\cos \iota\cos\Omega\right\}\) \\ \(z_{\mathrm{s}}=r_{\mathrm{s}}\sin\iota\sin i\) \\ \hline \(\mu=3.986004418\times 10^{14}\) m\({}^{3}\)/s\({}^{2}\) \\ \(\omega_{\mathrm{c}}=7.2921151467\times 10^{-5}\)rad /s \\ \(c=299792458\) m/s \\ \hline \hline \end{tabular}
\end{table} 
Table 1: The 15 F/NAV orbital parameters and the user algorithm. in [7; 8]. Moreover, GSR measurements are obtained from \(n\) ground stations. We assume that each ground station provides GSR observables to just one MEO satellite at a time, and that every satellite is only observed by one ground station at a time. In closed-loop link topology, all satellites are assumed linked in a 'closed-ring' manner, with \((m+n)\) range observations (ISR and GSR) available at each epoch, and \(15m\) unknowns. The ephemeris state vector \(\mathbf{p}\) can be solved by using a non-linear least squares adjustment to minimize objective function  
\[\sum_{k=1}^{K}\Big{\{}\big{(}\mathbf{y}_{t_{k}}-\mathbf{f}_{t_{k}}\big{(} \mathbf{p}(t_{oe})\big{)}^{\mathsf{T}}\mathbf{Q}_{\mathbf{y}_{k}}^{-1}\big{(} \mathbf{y}_{t_{k}}-\mathbf{f}_{t_{k}}\big{(}\mathbf{p}(t_{oe})\big{)}\big{)} \Big{\}}, \tag{3}\]  
where \(K\) is the number of epochs used, \(\mathbf{y}_{t_{k}}\) is the vector containing the \((m+n)\) observed ranges at time \(t_{k}\), \(\mathbf{f}_{t_{k}}\big{(}\mathbf{p}(t_{oe})\big{)}\) is the vector of nonlinear functions mapping the \(15m\) broadcast parameters to the range at time \(t_{k}\) (given in Table 2b), and \(\mathbf{Q}_{\mathbf{y}_{t_{k}}}^{-1}\) is the inverse of the covariance matrix of the observables.