# Low-Thrust Nonlinear Orbit Control for Very Low Lunar Orbits
## 4 Numerical Simulations  
The spacecraft is assumed to be equipped with a low-thrust propulsion system, characterized by the effective exhaust velocity \(c\) and the maximum thrust acceleration \(u_{Tmax}\),  
\[c=30\ \mathrm{km}/\mathrm{s},\ \ \ u_{Tmax}=5\cdot 10^{-5}g_{0}, \tag{75}\]  
where \(g_{0}=9.81\ \mathrm{m}/\mathrm{s}^{2}\) is the gravity acceleration at sea level. Canonical units are employed for the numerical simulations; the distance unit (DU) is equal to the Moon radius, while the time unit (TU) is selected in order to obtain a unit gravitational parameter,  
\[DU=1738.1\ \mathrm{km}, \tag{76}\] \[TU=\sqrt{\frac{1}{\mu_{m}}DU^{3}}. \tag{77}\]  
The following weighting coefficients are employed:  
\[k_{1}=10^{5},\ \ \ k_{2}=10^{5},\ \ \ k_{3}=10^{7}. \tag{78}\]  
Two target orbits are selected, with a different inclination of the orbital plane. Table 1 collects the initial conditions for these two orbits. The initial argument of latitude is set to \(10^{\circ}\).  
\begin{table}
\begin{tabular}{c c c} \hline \hline
**Orbital Element** & \(i_{d}\)**[deg]** & \(\Omega_{d,0}\)**[deg]** \\ \hline Target orbit 1 & 60 & 300 \\ Target orbit 2 & 90 & 300 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Target orbits.  
Instead, the initial state of the spacecraft is specified by assigning a random value of the initial radii at periselenium (\(r_{p}^{0}\)) and aposelenium (\(r_{A}^{0}\)), selected using a stochastic uniform distribution over the following intervals:  
\[r_{p}^{0}\in[50,100]\ \text{km}, \tag{79}\] \[r_{A}^{0}\in[100,500]\ \text{km}. \tag{80}\]  
Furthermore, the initial inclination and RAAN are chosen according to a normal distribution with mean value, respectively, \(i_{d}\) and \(\Omega_{d}^{0}\), and standard deviation is set to 5 deg.  
For each target orbit, two different operational scenarios are analyzed: (i) final circular orbit at an altitude of 100 km; (ii) final elliptic orbit with assigned minimum periselenium radius (\(r_{Pmin}=R_{m}+50\) km) and maximum aposelenium radius (\(r_{Amax}=R_{m}+150\) km). The latter operational scenario can be specified through the following relations:  
\[r_{P}=\frac{p}{1+e}\geq r_{Pmin}\to p\geq r_{Pmin}(1+e), \tag{81}\] \[r_{A}=\frac{p}{1-e}\leq r_{Amax}\to p\leq r_{Amax}(1-e). \tag{82}\]  
These two inequalities identify a triangular region in the \((p,e)\) plane, where the motion of the spacecraft is constrained. In particular, after orbit injection, the propulsive effort is employed in order to release the vehicle inside this specific region and the thrust is turned off when the inequalities (81) and (82) are satisfied, together with  
\[|\Psi_{3}|<10^{-6}. \tag{83}\]  
For each Monte Carlo campaign, 100 simulations are run. The initial reference epoch is set to 1 January 2025 and the total simulation time is about 1 year (370 days, 2 h, 51 min). For each simulation, the mean values of the semimajor axis (\(\overline{a}\)), eccentricity (\(\overline{e}\)), inclination (\(\overline{i}\)), and RAAN deviation (\(\Delta\overline{\Omega}\)) are computed as  
\[\overline{a}=\int_{t_{s}}^{t_{f}}a\ dt, \tag{84}\] \[\overline{e}=\int_{t_{s}}^{t_{f}}e\ dt,\] (85) \[\overline{i}=\int_{t_{s}}^{t_{f}}i\ dt,\] (86) \[\Delta\overline{\Omega}=\int_{t_{s}}^{t_{f}}(\Omega-\Omega_{d}) \ dt, \tag{87}\]  
where \(t_{f}\) is the final time, while \(t_{s}\) is conventionally set to 1 month, when the acquisition phase is completed. For each campaign, the mean values and standard deviations of each mean variable are provided. The run time associated with a single 1-year simulation is about 215 s on an Apple _M2 Pro_ processor, 32 GB RAM, 10 cores (4 high-efficiency + 6 high-performance).