# Lunar Satellite Constellations in Frozen Low Orbits
## 3 Design of Frozen Orbits  
The concept of a frozen orbit, an orbit with constant average eccentricity and argument of periapsis values, was first introduced for near-Earth orbits at the dawn of the space era and has soon become associated with an equilibrium of the averaged or doubly averaged system [12]. Unlike the Earth or an Earth-like oblate planet, the Moon has an irregular gravity field with no single dominant harmonic. Asymptotic expansions of the lunar orbiter perturbation theory, mostly derived in simplified models with low-degree harmonics and third-body attraction [13; 14; 15], fail to represent an accurate approximation to low-altitude frozen orbits, whereas more realistic models require tedious and lengthy symbolic computations [16; 17]. Not surprisingly, mission designers prefer numerical techniques since they are much easier to understand and implement. Along with the brute-force search in the high-precision dynamical model [18; 19], smarter approaches are also developed that are based on the numerical continuation of periodic orbit families via different predictor-corrector schemes [20; 21]. Recent advances include a clustering-based approach for analyzing numerically generated trajectories, which successfully identified 15 candidates for low lunar frozen or quasi-frozen orbits [22]. At the same time, it is worth noticing that the condition of periodicity in the synodic frame, which is equivalent to having a repeating ground track [21], is formally a stronger condition than the frozenness because the latter assumes periodic behavior for the argument of periapsis and the eccentricity only, not imposing the constraint of synodic resonance on the altitude. It is, therefore, essential to build a more flexible technique of lunar frozen orbit design, capable of searching for long-lived orbits under user-specified constraints on their orbital elements. The equal spacing of circular near-polar constellation orbits around the Moon is a relevant example of such a constraint. The major obstacle towards the robust design of frozen LLOs is very high sensitivity when an initial-guess orbit is propagated for a long time interval. It deteriorates the local correction procedure by a gradient method. So, we came up with a natural idea of using non-gradient methods. The best performance, among both gradient and non-gradient optimization solvers, is demonstrated by the Bayesian algorithm, one of the most powerful non-gradient techniques of optimizing a function that is computationally costly to evaluate [23]. In an earlier paper [9], MATLAB's bayesopt function was exploited. \begin{table}
\begin{tabular}{c c c c c c} \hline \hline \(h=260\) km & \(|\Delta a|_{t}\), km & \(|\Delta e|\) & \(|\Delta i|_{t}\) deg & \(|\Delta\Omega|_{t}\) deg & \(|\Delta\Lambda|_{t}\) deg \\ \hline min & \(5.03\cdot 10^{-3}\) & \(3.20\cdot 10^{-5}\) & \(9.91\cdot 10^{-5}\) & \(3.86\cdot 10^{-4}\) & \(1.39\cdot 10^{-1}\) \\ med & \(1.27\cdot 10^{-2}\) & \(5.57\cdot 10^{-5}\) & \(1.47\cdot 10^{-4}\) & \(9.06\cdot 10^{-4}\) & \(3.34\cdot 10^{-1}\) \\ max & \(2.31\cdot 10^{-2}\) & \(1.09\cdot 10^{-4}\) & \(2.32\cdot 10^{-4}\) & \(2.10\cdot 10^{-3}\) & \(6.11\cdot 10^{-1}\) \\ \hline \hline \(h=520\) km & \(|\Delta a|_{t}\), km & \(|\Delta e|\) & \(|\Delta i|_{t}\) deg & \(|\Delta\Omega|_{t}\) deg & \(|\Delta\Lambda|_{t}\) deg \\ \hline min & \(5.79\cdot 10^{-3}\) & \(3.31\cdot 10^{-5}\) & \(5.46\cdot 10^{-5}\) & \(9.93\cdot 10^{-4}\) & \(1.91\cdot 10^{-1}\) \\ med & \(9.99\cdot 10^{-3}\) & \(4.14\cdot 10^{-5}\) & \(8.46\cdot 10^{-5}\) & \(1.49\cdot 10^{-3}\) & \(3.54\cdot 10^{-1}\) \\ max & \(1.45\cdot 10^{-2}\) & \(7.35\cdot 10^{-5}\) & \(1.10\cdot 10^{-4}\) & \(2.11\cdot 10^{-3}\) & \(5.36\cdot 10^{-1}\) \\ \hline \hline \end{tabular}
\end{table} 
Table 1: One-year propagation errors of the truncated dynamic model. Currently, we use the improved Python version of the algorithm based on the open-source code developed by [PERSON] (2018).