# Lunar Satellite Constellations in Frozen Low Orbits
## 3 Design of Frozen Orbits
This global optimization solver augmented with the sequential domain reduction tool allows speeding up the convergence process by adaptively squeezing the search domain. Another source of computational efficiency is due to translating such low-level routines as numerical integration to Fortran and then calling them by the F2 PY wrapper from the NumPy library. Now it is time to give the technical details of the optimization problem posed. The objective function is defined as the squared Euclidean norm of the difference in the eccentricity vector over the propagation interval:  
\[(\Delta e_{x})^{2}+(\Delta e_{y})^{2}\rightarrow\min\]  
where \(\Delta e_{x}=e_{x,f}-e_{x,0}\) and \(\Delta e_{y}=e_{y,f}-e_{y,0}\) stand for the difference between the initial and final values of the eccentricity vector components. To search for Walker-type frozen constellations of a given altitude \(h=h_{\text{ref}}\), we impose the box constraints:  
\[a_{0}\in[a_{\min},a_{\max}],\ \ \ e_{x,0}\in[-e_{\max},e_{\max}],\ \ \ e_{y,0}\in[-e_{\max},e_{\max}]\]  
on the optimized vector \(\boldsymbol{\xi}=(a_{0},\,e_{x,0},\,e_{y,0})\), with the bounds being set to \(e_{\max}=0.02\) for the eccentricity and:  
\[a_{\min} =R_{M}+h_{\text{ref}}\cdot 0.95\] \[a_{\max} =R_{M}+h_{\text{ref}}\cdot 1.05\]  
for the semimajor axis. Here, \(R_{M}=1737.4\) km is the Moon's mean radius. The initial values of the inclination and the right ascension of the ascending node can be set fixed: we found that they do not affect the iterative procedure convergence. The choice of RAAN values, as well as of the reference altitude, comprises the issue of preliminary constellation design, which is covered in the next section. What concerns the inclination, it is assumed to be \(i_{0}=84^{\circ}\) for all the orbital planes in any design considered. Such a value is close to one of the four inclinations at which frozen LLOs can be near-circular ([PERSON], 2017). Suitable for observing lunar polar regions, the \(i_{0}=84^{\circ}\) inclination is at the same time more safe: strictly polar (\(i=90^{\circ}\)) constellations are more susceptible to the danger of collisions between the constellation satellites at high latitudes. For definiteness, a spacecraft is assumed to start orbiting from the ascending node and the initial moment of time \(t_{0}\) corresponds to the midnight of 1 January 2022. The propagation interval can be as long as one year (further, we considered \(t_{f}=t_{0}+365\) days). The obvious initial guess \(a_{0}=R_{M}+h_{\text{ref}}\), \(e_{x,0}=e_{y,0}=0\) appeared to be adequate in all the cases. The sequential domain reduction option was applied as follows. Starting with the 20 th iteration, the search box size is multiplied by 0.95. The total number of Bayesian algorithm iterations is manually set to 130, which is enough for deep capture in the basin of attraction of a correct local minimum. After that, to further decrease the objective function value to as low as \(10^{-8}\) and less, the Nelder-Mead non-gradient solver from the SciPy library is advised to be started from the point the Bayesian solver stopped with the same box constraints. The stopping condition of the Nelder-Mead solver was defined in terms of the optimized vector tolerance: \(|\Delta\boldsymbol{\xi}|<\texttt{xatol}=10^{-5}\). To validate the suggested two-stage procedure of frozen orbit design, we performed a lot of numerical experiments for various orbital plane orientations in space and different reference altitudes. Particularly, we have reproduced two repeat ground track orbits (RGT) from the paper of [PERSON] and [PERSON] (1998). The eccentricity vector evolution for two RGTs, medium (\(h\approx 3308\) km) and low (\(h\approx 124\) km), drawn in that paper (reprinted for convenience in Figures 2(a) and 3(a)), very closely matches--almost coincides--with what is observed in our simulations (see Figures 2(b) and Figure 3(b), respectively).