Exploiting Lunar Navigation Constellation for GNC Enhancement in Landing Missions

Abstract

To support the increasing number of planned lunar missions, a collaborative international initiative is underway to conceptualise and establish a lunar satellite constellation for communication and navigation. In this context, the goal of the current paper is to analyse what the obtainable performance is for a lunar lander that executes state estimation employing one-way ranging signals from such a Lunar Navigation Service (LNS). In particular, a small-sized optimised navigation constellation is considered as the main source of measurements, which, coupled with an accelerometer and an altimeter, is used to estimate the lander absolute trajectory during the main braking phase. The guidance is extracted on board by interpolation of a ground-optimised trajectory, followed by a reference-tracking regulator. Two alternative control tuning cases are presented, one targeting high performance, the other targeting low propellant mass. Nominal performance and associated sensitivity analyses assessed the feasibility of supporting such a critical phase with a reduced LNS constellation, reaching final control errors below 500$\mathrm{m}$, with the better performing one going down to 56$\mathrm{m}$. Among the two proposed alternatives, the one targeting low fuel expenditure has proven, however, to also be more robust against time and state uncertainty, providing much larger success rates.

1. Introduction

The International Space Exploration Coordination Group (ISECG) recognises lunar exploration as a pivotal milestone for advancing crucial technologies necessary for sustained human presence in space and deep-space exploration. As a result, lunar missions have regained much interest globally, attracting involvement from both space agencies and private companies. While these missions pursue different objectives, they all demand precise knowledge of spacecraft positioning and trajectory, along with a robust communication infrastructure for efficient data transmission back to Earth.

Up to now, lunar missions have heavily relied on Direct-to-Earth (DTE) communications and ranging radiometric measurements from the ground. However, as the number of missions grows, several space agencies have proposed dedicated lunar communication and navigation infrastructures to ease the Earth ground-segment load and establish more efficient and reliable communication links [1,2,3]. In this context, ESA initiated the Moonlight project, aimed at offering an affordable and high-performance Lunar Communication and Navigation Service to support upcoming institutional and commercial lunar missions. Such infrastructure will maintain continuous Earth contact, even during DTE link downtime (e.g., on the Moon’s far side), and provide high-precision navigation data. Furthermore, a dedicated Global Navigation Satellite System (GNSS) constellation could enable onboard autonomous lunar navigation systems for rovers, landers, and spacecraft, meeting the stringent ISECG requirements for safe Lunar surface operations, such as landing within a 90-metre 3$\sigma$ uncertainty from the intended location.

Looking instead specifically at the LNS constellation design, the research community is currently quite active, focusing on different procedures for the definition of the orbital configuration. When looking for a global coverage, the well-known principle of Walker constellation is a possible alternative [4,5]. Such a solution, however, may not be optimal in terms of cost effectiveness, in particular if an incremental approach is forsaken for the constellation deployment. A more advanced design proposed in two consecutive publications [6,7] considered multi-objective optimisation to provide high navigation performance availability to the lunar global surface, while still also considering other system-level objectives such as costs, station-keeping expenditure, and robustness to failure. Other studies considered instead constellations deployed in Lagrange Point Orbits (LPO), which provide many geometrical and dynamical features that classical two-body orbits do not possess [8,9]. A hybrid approach to the design of the architecture was instead followed in [10], where both Keplerian and LPOs are considered in the optimisation, having performance and coverage as maximisation indexes.

In any case, looking at an incremental constellation approach, during the initial phases a primary issue will be the limited availability of navigation signals due to the substantially fewer LNS servicing satellites compared to consolidated Earth GNSS systems. The smaller constellation size implies the impossibility of having a continuous view of at least four satellites from any location on the lunar surface at the same time. For this reason, the initial LNS architecture will be designed to prioritise enhancing the navigation performance of surface and Low Lunar Orbital (LLO) users around the South Pole region (which is a subject of high scientific interest) at the cost of a reduced LNS availability in the northern hemisphere. Nonetheless, even for these South Pole users, custom navigation algorithms will need to be developed to meet the mission navigation requirements and cope with the limited LNS-visibility windows.

To counteract such problems, recent years have seen numerous studies investigating achievable performance levels for various lunar missions utilising potential LNS constellations. Notably, for the state estimation, ref. [11] utilised a sequential Extended Kalman Filter (EKF) to highlight the system complexity benefits of implementing the LNS infrastructure over traditional visual-based navigation sensors and ground-tracking techniques. Their work revealed that a formal horizontal dilution of precision below 30 metres can be achieved with a minimum of three satellites in view. Another case study is provided by [12], where the LNS system is used coupled with a high-fidelity Digital Elevation Model (DEM) to navigate a rover under various service availability conditions. Conversely, ref. [13] employed a Batch filter to fuse LNS signals with IMU and altimeter readings, resulting in navigation errors below 30 metres. In both cases, the kinematic approach used for the estimation implies that navigation errors would rapidly diverge during signal outages. The advantages of a dynamic estimation process were discussed in [14,15]. The attained performance significantly depends on the accuracy of the onboard dynamical model, necessitating precise Sun and Earth ephemerides and harmonic coefficients for the spherical expansion of the lunar gravity field. However, these demands might conflict with low-cost hardware and Commercial-Off-The-Shelf (COTS) components, especially when the filter must operate at high frequencies, which is why a proper navigation architecture is fundamental (see [16]). Comprehensive introductions to the dynamic state estimation theory, together with a description of the main algorithms and comparisons between sequential and batch filtering techniques, can be found in [17,18,19].

The current paper is intended as a continuation of a series of successive articles in the framework of lunar constellations. Indeed, the first two studies proposed a method to design an LNS constellation to provide communication and navigation services to different users, optimising the orbital configuration considering multiple objectives. In particular, in [20] we focused on surface users only, and multi-objective optimisation was used to service different Moon regions. Then, in [14], we focussed our attention on orbital users in polar inclinations. Also in this case, an optimisation combining different clashing objectives was performed in order to derive optimal constellation architectures. The deriving navigation performance for the orbital user was assessed in accordance with a proposed filtering architecture, which has been tested under different conditions also in [16]. From the set of optimal configurations extracted in [14], the best performing one was selected to be used for further assessments, such as in this current work. In our previous analyses and in the investigated literature, however, only the navigation task was taken into consideration. Even the works concerning landers or rovers considered the dynamics of the user as defined independently from the navigation tasks, thus excluding a complete GNC scheme. The goal of this paper is to evaluate the behaviour of the LNS architecture when inserted into a more dynamically challenging scenario. Indeed, testing the final performance of the GNC in closed loop provides a testbed to understand if, even with a reduced constellation, the critical tasks of a landing operation can be completed successfully. In particular, the main braking phase of the landing is considered here, as the last phase in a landing which could be handled with an absolute navigation strategy, before handover to a relative one for the final (or terminal) descent.

The results in a nominal scenario are presented, after the description of the added guidance and control algorithms are provided. After that, a Monte Carlo analysis with dispersion in the initial state is performed and the overall results collected.

After this introduction, the complete modelling environment used for spacecraft dynamics and sensor simulation is presented in Section 2. The description of the core of the employed GNC algorithms is provided in Section 3. Section 4 shows the obtained results for both the nominal scenario and the sensitivity analysis, and the conclusions of the manuscript are provided in Section 5.

2. Simulation Environment

The simulation environment used to performed the analyses of this study is described in the following section. In particular, we first present the dynamic environment in which the servicers and the lander are propagated, together with the associated assumptions. Then, we provide an overview of the sensor models used in the proposed navigation strategy.

2.1. Ground-Truth Dynamics

To accurately simulate the ground-truth dynamics of the landing user and LNS servicers, the major gravitational and non-gravitational forces acting in the cislunar environment have been modelled. The former include the gravitational attraction of the Moon, the Earth, and the Sun. The remaining solar system bodies have not been included because their gravitational influence (${10}^{-12}\mathrm{m}/{\mathrm{s}}^2$ for Jupiter) is negligible over the short-term. The irregularities of the lunar gravity field are considered by modelling the lunar gravitational potential, $\mathcal{U}$, in a body-fixed reference frame with a Spherical Harmonic Expansion (SHE) [21], as in Equation (1):

$$\begin{array}{cc}\hfill \mathcal{U}& =\frac{\mu }{r}+\frac{\mu }{r}\sum _{n=2}^N\sum _{m=0}^n{\left(\frac{R_0}{r}\right)}^n\left[{\overline{C}}_{nm}cos\left(m\lambda \right)+{\overline{S}}_{nm}sin\left(m\lambda \right)\right]{\overline{\mathcal{P}}}_{nm}(sin\varphi )\hfill \end{array}$$

(1)

where $\mu$ and $R_0$ are the gravitational parameter and radius of the Moon, r is the radial distance from the center of mass, $\lambda$ is the east longitude, $\varphi$ is the latitude, and ${\overline{\mathcal{P}}}_{nm}$ are the fully normalised associated Legendre polynomials. The normalised gravity coefficients (${\overline{C}}_{nm}$ and ${\overline{S}}_{nm}$) are based on the LP165P gravity model [22]. The harmonic expansion is truncated at order and degree 60 and is such that the error introduced by neglecting the higher-order terms at these altitudes is below the sensitivity of traditional on-board accelerometers.

For the Sun and Earth contributions, only the point-mass gravitational acceleration is relevant. In particular, Earth’s $J_2$ term produces an acceleration in the order of only ${10}^{-11}\mathrm{m}/{\mathrm{s}}^2$ for a lunar orbiter in a 10 km LLO. According to the Restricted N-Body Problem formulation [21], the perturbing acceleration induced on the spacecraft, ${\mathcal{P}}_i$, by the point-mass gravitational sources, ${\mathcal{P}}_j$, as seen by the central body, ${\mathcal{P}}_k$, is modelled as per Equation (2):

$${\ddot{\mathbf{r}}}_{ki}=\sum _{\begin{array}{c}j=1\\ j\ne i,k\end{array}}^{\mathrm{N}}{\mu }_j\left(\frac{{\mathbf{r}}_{ij}}{r_{ij}^3}-\frac{{\mathbf{r}}_{kj}}{r_{kj}^3}\right)$$

(2)

where the generic vector ${\mathbf{r}}_{ij}$ is the position of the mass ${\mathcal{P}}_i$ relative to ${\mathcal{P}}_j$ and is computed as:

$${\mathbf{r}}_{ji}={\mathbf{r}}_{kj}-{\mathbf{r}}_{ki}$$

(3)

In this work, NAIF’s SPICE library (see [23]) is exploited to retrieve the exact Earth and Sun positions from the DE440 (https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/, accessed on 20 July 2023) ephemerides. All the remaining field forces, such as general relativistic effects (∼${10}^{-10}\mathrm{m}/{\mathrm{s}}^2$ for a 10 km LLO), can be safely neglected for short-term propagations.

On the other hand, the remaining forces acting on the spacecraft are associated with non-gravitational perturbations and must be accounted for in a high-fidelity lunar dynamical model. These include contributions such as the Solar Radiation Pressure (SRP), the spacecraft thermal emissions, the lunar albedo and thermal infrared emissions, and the navigation antenna thrust for the LNS servicers. In this analysis, only the SRP contribution is taken into account since it is by far the largest non-gravitational acceleration. A standard cannonball model, which assumes that the satellite is a homogeneous sphere of given radius [24], is used to express the magnitude of the SRP acceleration, as shown in Equation (4):

$$a_{SRP}=\frac{S_{⨀}}{c}\frac{{\left(1\mathrm{AU}\right)}^2}{d_{⨀}^2}{c}_R\frac{A_{⨀}}{m}$$

(4)

where $S_{⨀}=1367 \mathrm{W}{\mathrm{m}}^{-2}$ is the Sun mean flux at 1AU, c = 299,792,458 m/s⁻¹ is the speed of light, $d_{⨀}$ is the current Sun–spacecraft distance, $c_R$ is the reflectivity, and $A_{⨀}$ is the cross-sectional area of the satellite exposed to the radiation. Although a traditional box-wing model provides a more accurate representation of the SRP acceleration [25], by modelling the spacecraft as a sphere, the cannonball model allows one to decouple the user trajectory from the geometry and actual orientation of the spacecraft. Additionally, we assume that throughout the landing mission, the spacecraft is in constant Sun illumination.

2.2. User Navigation Equipment

Dealing with a navigation filter, it is fundamental to correctly analyse the measurements that are exploited by the estimation algorithm in order to simulate the behaviour of the sensors providing such observations and define the measurement functions used within the filter itself.

2.2.1. Radiofrequency-Based Measurements

In this work, the user LNS receiver terminal retrieves from any ith visible element of the constellation the range $\rho$, the range-rate, $\dot{\rho }$, and the servicer ephemeris under the form of the state vector ${\mathbf{x}}_{s,i}$. When dealing with such measurements, the observable data is usually defined as pseudorange $\tilde{\rho }$ and pseudorange-rate $\tilde{\dot{\rho }}$, in order to indicate that the two measurements are affected by errors. In particular, looking at the range measurement, the geometric range obtained from the ith servicer can be defined as per Equation (5):

$${\rho }_i=c\Delta t_i=∥{\mathbf{r}}_{s,i}-\mathbf{r}∥$$

(5)

where c is the speed of light and $\Delta t_i$ is the time required for the signal to travel from the servicer position, ${\mathbf{r}}_{s,i}$, to the user one, $\mathbf{r}$. To simplify the analysis, the availability of the signal from the i-th servicer to the user is constrained only by point-to-point and Field-Of-View (FOV) considerations. In one-way ranging, this time difference is obtained with standard Time-of-Arrival (ToA) techniques by subtracting the servicer clock-time at signal emission from the user clock-time at reception, i.e., $\Delta t_i=t_u-t_{s,i}$. However, this measure is affected by errors since both the receiver and servicer clock present some offsets with respect to the LNS system time. As such, the measured time difference from the user terminal will provide the pseudorange, as per Equation (6):

$${\tilde{\rho }}_i=c(t_u+\delta t_u-t_{s,i}-\delta t_{s,i})={\rho }_i+c\delta t_u+{\epsilon }_{\delta t_{s,i}}$$

(6)

where $t_u$ is the system time at which the signal would have reached the user in the absence of errors, $\delta t_u$ is the receiver clock bias, $t_{s,i}$ is the system time at which the signal left the ith servicer, and $\delta t_{s,i}$ is the offset of the ith servicer clock from the system time. From this equation, it is clear that a bias in the measurements is present. Additionally, there are other possible sources of errors in the pseudorange readings, such as interference effects, receiver noise, instrumental delays, multipath losses, and relativistic effects. With this in mind, the final pseudorange formulation is reported in Equation (7), where all these error contributions (including ${\epsilon }_{\delta t_{s,i}}$) are collapsed in a single component, ${\epsilon }_{{\rho }_i}$.

$${\tilde{\rho }}_i={\rho }_i+c\delta t_u+{\epsilon }_{{\rho }_i}$$

(7)

A similar reasoning can be performed for the pseudorange-rate observable, which can be derived directly from the observed frequency Doppler shift, $\Delta f_D$, but again is inherently affected by the same effects as in the previous case, summed up in a ${\epsilon }_{{\dot{\rho }}_i}$ term in addition to the time bias derivative, ${\dot{\delta t}}_u$, as shown by Equation (8).

$${\tilde{\dot{\rho }}}_i=c\frac{\Delta f_D}{f_{s,i}}+c{\dot{\delta t}}_u+{\epsilon }_{{\dot{\rho }}_i}$$

(8)

In this work, to simulate these effects and generate the measurements fed to the filter, we implemented the model in Equations (9) and (10).

$$\begin{array}{cc}\hfill {\tilde{\rho }}_i& =∥{\tilde{\mathbf{r}}}_{s,i}-\mathbf{r}∥+b_c+{\epsilon }_{{\rho }_i}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\tilde{\dot{\rho }}}_i& =({\tilde{\mathbf{v}}}_{s,i}-\mathbf{v})·({\tilde{\mathbf{r}}}_{s,i}-\mathbf{r})/\rho +d_c+{\epsilon }_{{\dot{\rho }}_i}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{x}}}_{s,i}& ={\left[\begin{array}{c}{\tilde{\mathbf{r}}}_{s,i}^{\top },{\tilde{\mathbf{v}}}_{s,i}^{\top }\end{array}\right]}^{\top }={\mathbf{x}}_{s,i}+{\left[\begin{array}{c}{\epsilon }_r^{1\times 3},{\epsilon }_v^{1\times 3}\end{array}\right]}^{\top }\hfill \end{array}$$

(11)

For both the observables, all the error effects are directly collected in single additive Gaussian noise terms, except for the receiver clock bias, $b_c$, and drift, $d_c$, contributions, which are treated as an additional parameter to be estimated. The standard deviations associated with ${\epsilon }_{\rho }$ and ${\epsilon }_{\dot{\rho }}$ are ${\sigma }_{\rho }=10 \mathrm{m}$ and ${\sigma }_{\dot{\rho }}=0.1 \mathrm{m}/\mathrm{s}$, values assumed in accordance with the ongoing studies [11,13].

The receiver clock bias and drift, which represent $c\delta t_u$ and $c{\dot{\delta t}}_u$, respectively, are modelled using the two-state clock model presented in [26], in which the frequency deviation (i.e., the clock drift) originates from two types of noise, a White Frequency Modulation (WFM) and a Random Walk Frequency Modulation (RWFM). The resulting clock bias will then be represented by a Wiener noise plus an integrated Wiener noise. The dynamical system that simulates the evolution of these quantities is summarised in Equations (12) and (13), in which ${\epsilon }_b$ and ${\epsilon }_d$ are Gaussian white noises, whose standard deviation can be properly tuned to match the Allan variance of the desired type of receiver clock.

$$\begin{array}{cc}\hfill {\dot{b}}_c& =d_c+{\epsilon }_{b_c}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\dot{d}}_c& ={\epsilon }_{d_c}\hfill \end{array}$$

(13)

The servicers’ ephemerides are used inside the navigation filter to predict the pseudorange and pseudorange-rate quantities. In this work, the residual errors of the LNS predictions coming from the Orbit Determination and Time Synchronisation (ODTS) process performed on the ground, also known as Signal In Space Error (SISE), are modelled as simple three-dimensional position and velocity additive Gaussian white noises, with standard deviations of ${\sigma }_r$ and ${\sigma }_v$, respectively, here assumed as 15 $\mathrm{m}$ and $0.15$ $\mathrm{m}$$\mathrm{s}$, respectively, in line with the Galileo Service Definition Document [27] and the expected Moonlight performance [11].

A good reference to further investigate the details of GNSS systems can be found in [28].

2.2.2. Accelerometer

Concerning the accelerometers contained in the Inertial Measurement Unit (IMU), it is fundamental to recall that such sensors are insensitive to volume forces, such as all the gravitational effects included in the environment of our simulation.

For this reason, they are particularly useful for obtaining an estimate of the real control acceleration provided by the thrusters, as well as the perturbations induced by the SRP.

A high-fidelity model for an accelerometer typically includes misalignment and cross-axis sensitivity errors, sensor biases, and scale factor errors [29]. For the purpose of this work, the model has been simplified and all the errors have been collected into an additive Gaussian white noise, with a standard deviation that is representative of commercial high-fidelity accelerometers. As a consequence, the simulated sensor reading is defined as per Equation (14).

$${\tilde{\mathbf{a}}}_{IMU,I}={\mathbf{a}}_{SRP,I}+{\mathbf{a}}_{thr,I}+{\epsilon }_{\mathbf{IMU}}$$

(14)

It is also relevant to recall that the accelerations recorded by the IMU are usually expressed in the sensor frame, S. However, since in this framework the attitude dynamics and kinematics are not taken into consideration, the spacecraft attitude is assumed perfectly known and the accelerations of the IMU are directly expressed in the inertial frame, I.

A consolidated practice to exploit such sensors in the orbital filters is the so-called Dynamic Compensation Mode [17]. The basic principle is not to use the accelerometer readings among the measurement vector, but to leverage them within the propagation step of the filter. This concept, also exploited with gyroscopes, is backed by the high accuracies that such sensors show. Moreover, some computation savings are thus achievable, avoiding the evaluation of the non-volume forces, which generally have quite complex mathematical formulations, e.g., a high-fidelity SRP or propulsive acceleration representation.

2.2.3. Altimeter

In landing applications, altimeters are very relevant sensors since they provide a direct estimate of the vertical position in a user-centred East–North–Up (ENU) reference frame. Additionally for this case, there is the possibility to employ fancy models of the sensor, considering the exploited technology and the considered topography of the Moon. Indeed, regarding this latter, there is a wide range of possibilities, from assuming a perfect sphere or using detailed DEM. In the current work, using the simpler spherical Moon model, the real altitude is obtained as in Equation (15) employing Moon average radius, $R_{☽}$.

$$\zeta =∥\mathbf{r}∥-R_{☽}$$

(15)

The error effects introduced by the sensor are provided by a zero-mean additive white Gaussian noise, which, to reflect the behaviour of laser altimeter technology, considers a standard deviation of 1% of the current real height [30]. The resulting measurement function is thus defined as in Equation (16):

$$\tilde{\zeta }=\zeta +{\epsilon }_{\zeta }\left(\zeta \right)$$

(16)

with ${\epsilon }_{\zeta }\left(\zeta \right)=\mathcal{N}(0,1)\zeta /100$ defined as the noisy component.

3. GNC Algorithms

The Guidance Navigation and Control strategy exploited in this landing scenario relies on the architecture shown in Figure 1.

From the real world simulation, the dynamics of the lander and of the servicers is propagated in a high-fidelity environment, and the sensor readings are generated. These are fed to the navigation block, whose output is the absolute state estimation. The guidance and control block employs this state, together with the targeted trajectory, to generate a control action, ${\mathbf{u}}_k$.

3.1. Navigation

Similarly to our previous work in [14], the navigation is performed through an Extended Kalman Filter, which employs, as in the previous study, the GNSS/INS strategy. In the current scenario, however, an altimeter is considered as well, deemed as fundamental for both state estimation enhancement and safety reason. The work we performed in [16] already analysed the improvements that can be achieved by combining these sensors for orbital users. Specifically, it highlighted that by providing a direct estimate of the spacecraft vertical components, the altimeter allows the reduction by one of the actual number of visible LNS satellites that are requested to reconstruct the user position, and thus becomes particularly beneficial for small-sized LNS constellations. The resulting navigation architecture is schematised in Figure 2.

The presence of the IMU, used in dynamic model replacement, outputs the acceleration, ${\tilde{\mathbf{a}}}_{IMU}$, which is added to the on-board dynamics for the state prediction step in the filter, providing the a priori estimate, ${\widehat{\mathbf{x}}}_{k+1}^{-}$, at time step $k+1$. Such an approach provides the possibility to directly retrieve a measurement of control acceleration, ${\tilde{\mathbf{u}}}_k$, to be fed in the EKF prediction step, avoiding using the acceleration predicted by the controller ${\mathbf{u}}_k$. Indeed, the measurement error introduced by the accelerometer is much lower than the error between the commanded and executed actuator action. The accelerometer is employed, however, only during the thrusted phase, while in the preceding one it is disregarded. Indeed, during natural motion, the error introduced by the accelerometer is larger than the unmodelled dynamics of SRP and higher-order irregular mass distribution.

The observables of the LNS are composed of pseudorange, ${\tilde{\rho }}_i$, pseudorange-rate, ${\tilde{\dot{\rho }}}_i$, and the servicers’ ephemerides. Such measurements are passed to the filter, together with the altimeter reading, $\tilde{\zeta }$, in order to provide, through the Kalman update step, the a posteriori corrected state, ${\widehat{\mathbf{x}}}_{k+1}^{+}$. Such GNSS/INS formulation, known as tightly coupled, takes a different approach than the loosely coupled alternative (see [31]), since the GNSS (or LNS) receiver is not treated as a mere black box that outputs the complete state. Instead, this approach offers distinct advantages by exploiting LNS signals, even when fewer than four satellites are visible, a condition which is a likely scenario in our proposed context due to the limited constellation elements. In scenarios where LNS signals are absent, our proposed architecture seamlessly continues the onboard propagation of the user’s position, using only the information provided by the IMU.

Opting to incorporate the clock bias, $b_c$, and drift, $d_c$, as part of the parameter estimation, we augment the spacecraft state, ${\mathbf{x}}_{S/C}$, by including these components. This is obtained by considering a stacked filter state, $\mathbf{x}={[{\mathbf{x}}_{S/C}^{\top },b_c,d_c]}^{\top }$.

The prediction step of the EKF is characterised by a simplified spacecraft dynamical model, ${\mathbf{f}}_{S/C}({\mathbf{x}}_{S/C},{\tilde{\mathbf{a}}}_{IMU})$, and a clock dynamical model governing the propagation and measurement function, $\mathbf{h}(\mathbf{x},{\tilde{\mathbf{x}}}_{s,i})$, outlined respectively in Equations (17) and (18).

$${\mathbf{f}}_{S/C}({\mathbf{x}}_{S/C},{\tilde{\mathbf{a}}}_{IMU})={\mathbf{f}}_{S/C}{\left({\mathbf{x}}_{S/C}\right)}_{☽}+{\mathbf{f}}_{S/C}{\left({\mathbf{x}}_{S/C}\right)}_{♁}+{\tilde{\mathbf{a}}}_{IMU}$$

(17)

$$\mathbf{h}(\mathbf{x},{\tilde{\mathbf{x}}}_{s,i})={\left[{\rho }_1,{\dot{\rho }}_1,{\rho }_2,{\dot{\rho }}_2,\cdots ,{\rho }_n,{\dot{\rho }}_n,\zeta \right]}^{\top }$$

(18)

Equation (17) shows the strategy of dynamic model replacement to leverage the accelerometer within the navigation filter, encompassing both Solar Radiation Pressure (SRP) and the thrusters’ acceleration contributions. It is also possible to see that only gravitational terms of Moon and Earth point mass are introduced. The Moon irregular mass distribution is then simply modelled by adding the oblateness, $J_2$, contribution, in order to reduce the computational effort. The spacecraft state is estimated in the Lunar Mean Equator at J2000 (LME2000) frame, whose origin is the Moon’s centre of mass, and its X–Y plane is defined by the Moon equator at J2000. In particular, the +Z axis is parallel to the Moon’s rotation axis, pointing toward the north side of the invariant plane. The +X axis is defined by the intersection of the Moon’s equator at the J2000 with the Earth Mean Equator of J2000, and the +Y axis completes the right-handed triad.

Moreover, despite the fact that the broadcast ephemeris of traditional GNSS constellations are provided in the Earth-Centred–Earth-Fixed (ECEF) frame [32], in this scenario it is further assumed that the retrieved states of the LNS satellites are also directly expressed in the LME2000 frame, as a proper definition of an official Moon-Centred–Moon-Fixed (MCMF) frame is still not available.

Additionally, it is worth underlining that the measurement function of Equation (18) integrates the estimated clock bias and drift to refine the pseudorange and range-rate a priori predictions, guided by Equations (9) and (10). The propagation model for these parameters is expressed by ${\mathbf{f}}_c$, as defined in Equation (19).

$${\mathbf{f}}_c=\left\{\begin{array}{cc}{\dot{b}}_c\hfill & =d_c\hfill \\ {\dot{d}}_c\hfill & =0\hfill \end{array}\right.$$

(19)

Combining Equations (17) and (19) results in the comprehensive dynamical model for the filter, denoted as $\mathbf{f}={[{\mathbf{f}}_{S/C}^{\top },{\mathbf{f}}_c^{\top }]}^{\top }$.

Initialisation of the filter entails setting the initial conditions for the state, ${\widehat{\mathbf{x}}}_0$, and the state covariance matrix, ${\mathbf{P}}_0$, at the beginning of the simulation. For ${\mathbf{P}}_0$, we considered the values of ${\sigma }_{r,0}=1 \mathrm{k}\mathrm{m}$ and ${\sigma }_{v,0}=1 \mathrm{m}\mathrm{s}$, which also serve as the standard deviations for ${\widehat{\mathbf{x}}}_0$. $\mathbf{Q}$ and $\mathbf{R}$ represent the process and measurement noise covariance matrices. If the former can be regarded as a tuning parameter, the latter is dictated by the accuracy of the sensors.

The evaluation of the innovation vector defined by ${\tilde{\mathbf{y}}}_{k+1}-h({\widehat{\mathbf{x}}}_{k+1}^{-},{\tilde{\mathbf{x}}}_{s,i,k+1})$ is influenced by two distinct errors: the pseudorange and range-rate errors, resulting from the physics of signal transmission and RF terminals, along with the errors arising in the ephemerides of the servicers, an outcome of their individual navigation budgets.

Given the complex and more dynamically challenging scenario with respect to the natural motion studied in [14], the filter update frequency has been set to 10 $\mathrm{Hz}$.

3.2. Guidance and Control

The guidance and control strategy employed for this paper is to provide an offline optimised landing trajectory that is interpolated on board and then fed to a reference tracking controller, based on a PID regulator. Such a strategy is not completely autonomous, since the guidance optimisation is performed on the ground and uploaded to the lander just before starting this critical phase. It is, however, performed completely autonomously during the actual thrusting period.

The trajectory used for the guidance is generated offline through an optimisation similar to the one presented in [33,34], with the goal of minimising fuel consumption. The procedure considers an indirect optimisation step to solve a two-point boundary value problem with prescribed initial and final conditions in a simplified 2D dynamics on the polar orbital plane. A constant thrust value is considered, and the minimisation of propellant is obtained by constructing a cost function considering the integral of the control acceleration, thus including the evolving mass state variable. The solution to the indirect optimisation consists of finding the unknown initial conditions of the adjoint states by exploiting a shooting method. The result of this process is a history of the control action and the associated 2D dynamics, which can be remapped to a 3D trajectory by knowing the mission scenario geometry.

Such a trajectory is imported only as a sequence of six states, discretised in the sampling times output of the optimisation. To exploit this sequence, the on-board guidance uses the offline optimised six states as grid points to perform the interpolation and to retrieve the target states at any given query point. In such a way, the controller is fed at any instant with a control error computed between the obtained reference trajectory and the estimated current state, coming from the navigation. The control action of the optimal guidance is instead not exploited, since such an action is recomputed using the PID controller described in Figure 3.

The input of this block, $\mathbf{\Delta }{\widehat{\mathbf{x}}}_{\mathbf{k}}$, is computed as the error between the set-point trajectory state at time $t_k$ and the current estimation of the state, ${\widehat{\mathbf{x}}}_{\mathbf{k}}$. This state error is decomposed into its position and velocity parts, $\mathbf{\Delta }{\widehat{\mathbf{r}}}_{\mathbf{k}}$ and $\mathbf{\Delta }{\widehat{\mathbf{v}}}_{\mathbf{k}}$. The former is used for the proportional and the integral parts of the control, while the latter is used for the derivative one. The sum of the three components then passes through a saturation block, in order to impose the maximum thrust constraints. To prevent uncontrollability due to the actuator saturation caused by the integral term, a simple anti-windup logic [35] is introduced. The computed control action before passing the saturation step is compared in absolute value to the output of the limiter. A non-null difference means that the saturation is active, so to avoid controllability issues, the input to the integral term is driven to zero. In such a way, the regulator is able to avoid the windup condition, for which a saturated control action would continue increasing the integral of the error. The anti-windup logic remains active as long as the control acts to reduce the error and consequently re-enter the saturation limit.

The only available tuning parameters for such a controller are the three PID coefficients, $k_P$, $k_I$, and $k_D$, representing proportional, integral, and derivative coefficients.

The control action computed by the controller is perturbed before entering the lander dynamics, in order to simulate both pointing and magnitude zero mean additive Gaussian errors. The former presents a standard deviation of 1° on each of the two angles defining the thrusting direction, while the latter is defined with standard deviation equal to 0.1% of the maximum thrust.

4. Results

After presenting in Section 4.1 the adopted LNS constellation architecture, the considered GNC strategy is applied to a specific scenario, described in the following Section 4.2. The results obtained in the nominal landing scenario are reported in Section 4.3. Finally, the sensitivity of the approach is tested through a Monte Carlo analysis in Section 4.4.

4.1. LNS Constellation Architecture

The LNS constellation considered for this assessment has been obtained in our previous paper through a multi-objective optimisation, based on maximising the Dilution Of Precision (DOP) performance for a set of polar LLO users. The details of how we set up the objectives and the genotype can be found in [14], together with the related navigation performance for a naturally evolving LLO user. The constellation is composed of six servicers posed in Elliptic Lunar Frozen Orbits (ELFO, details provided in [36]), considered as the best candidates for an LNS constellation due to their relatively high stability in the long-term, thus substantially reducing the associated station-keeping $\Delta v$. The selected orbital architecture (tagged llo_sgl6A in the original work) is presented in Figure 4.

The six orbits share the same value of semi-major axis, eccentricity, inclination, and argument of pericentre. The former is imposed in order to have a period of 24 h, which is a prominent feature for operative reasons for such an asset. In addition, the argument of pericentre is imposed a priori, with a value of 90° in order to have the aposelene above the South Pole and prioritise the visibility of such a scientifically relevant location. The eccentricity of 0.489 is instead an output of the optimisation procedure. With the previous parameters fixed, the ELFO constraint leads to an inclination of 47.5°. The right ascension of the ascending nodes and the true anomalies are instead the other optimised parameters, different per each element of the constellation. As visible in Figure 4, this architecture presents a clear clustering of the servicers, with only two planes considered, with right ascension values of approximately ±90°. The true anomalies instead are almost equally distributed.

As thoroughly described in [14], this constellation is able to provide around-metre navigation accuracy when more than three servicers are visible. Such a feature is fundamental for a landing scenario, in which satisfying the accuracy requirements is often extremely challenging. Indeed, in this case, a complete closed-loop GNC strategy is considered, adding propulsion to the orbital user, which introduces more complex navigation needs.

4.2. Landing Mission Overview

The case of a South Pole lander is analysed here considering only the main braking phase of the landing. Starting from the apocentre of an elliptical orbit ($100\times 15 \mathrm{k}\mathrm{m}$ altitude), the lander follows a natural path up to the periselene, at half of its orbital period, i.e., after roughly 3400 $\mathrm{s}$. Here it starts the thrusted arcs from an altitude of 15 $\mathrm{k}$$\mathrm{m}$, which ends at 1 $\mathrm{k}$$\mathrm{m}$, before relative navigation hand-over. The terminal landing phase is not included in this analysis, since it cannot be separated from the exploitation of a fully relative navigation strategy.

The trajectory followed by the lander can be seen in Figure 5, superposed to the landing guidance trajectory. A zoomed trajectory of the propelled arc is instead shown in Figure 6, with the view taken directly on the landing trajectory plane. The lander is assumed to have a mass of 1600 kg, a ballistic coefficient of 150 $\mathrm{k}\mathrm{g}$/$\mathrm{m}$${}^2$, a maximum thrust of $2.5$ $\mathrm{k}$$\mathrm{N}$, and a specific impulse of 320 $\mathrm{s}$.

4.3. Nominal Scenario Performance

The current analysis entails and compares two different controllers in the presented scenario. In particular, these two controllers have been obtained through two separated tuning procedures: a minimisation of the final trajectory control error (case A) and a minimisation of the propellant expenditure (case B), whilst maintaining a final error below the threshold of 500 $\mathrm{m}$. Such a threshold is considered enough for the main braking phase, which is targeted by this study, before handling the final descend where a precision below 90 $\mathrm{m}$ is targeted.

The results from a GNC point of view are displayed for the two cases in Figure 7 and Figure 8, where both the estimation and control errors are plotted over the simulated time on a logarithmic scale. The total number of LNS satellites in view by the lander is plotted in Figure 9. Note that the control error plot starts slightly before $t=1 \mathrm{h}$, i.e., when the propulsive leg begins.

The navigation errors depend on the number of satellites visible (see Figure 9), for both cases, evolving practically identically for the non-propelled phase. Indeed, for the first part of the trajectory the error remains quite high, until a new servicer enters the visibility of the lander, just before $t=0.4 \mathrm{h}$. At that instant, the navigation error drops to the order of ∼ 100 $\mathrm{m}$ and then reduces even more as soon as the number of satellites in view increases up to a peak of five servicers. The reached estimation accuracy stabilises in the order of a few metres in this part. When the control action starts, the initial oscillatory behaviour of the thrust acceleration is reflected on both the navigation and control errors. The former recovers the accuracy of a few metres after stabilisation. It is useful to remark that the filter is in a stable and converged state just before the start of the control phase, which is reflected by the low values in the covariance matrix, $\mathbf{P}$. This effect gives a high stiffness to the filter, which prevents its convergence to the new correct status at the beginning of the thrusting phase. In order to overcome this issue, the covariance matrix is re-initialised to 10% of its initial condition, ${\mathbf{P}}_0$, at the beginning of the thrusted arc. In such a way, the filter is made more flexible and driven to convergence to the new state.

The control error instead behaves differently in the two cases, as seen in both Figure 7, Figure 8 and Figure 10, where such errors are plotted in a Local Vertical Local Horizontal (LVLH) frame relative to the set-point trajectory alongside the propellant mass expenditure.

Case A, tuned for best control error performance, after the initial oscillating behaviour, is able to reach a final error of $55.6$ $\mathrm{m}$, which is kept stable during the whole phase, at the expense of $13.3$ $\mathrm{k}\mathrm{g}$ of propellant mass. On the other hand, case B reaches a coarser value of $380.4$ $\mathrm{m}$, using only less than half of the propellant mass, with a value of $6.6$ $\mathrm{k}\mathrm{g}$. The performance of the two controllers is summarised in

Table 1. Collection of performance for the two proposed controllers in the nominal landing scenario. Final navigation and control errors are reported, with total propellant mass and the three PID tuning coefficients.

CTRL	$\mathit{\delta }{\mathit{r}}_{\mathit{N}\mathit{A}\mathit{V}}$ (m)	$\mathit{\delta }{\mathit{r}}_{\mathit{C}\mathit{T}\mathit{R}\mathit{L}}$ (m)	$\mathit{\delta }{\mathit{x}}_{\mathit{C}\mathit{T}\mathit{R}\mathit{L}}$ (m)	$\mathit{\delta }{\mathit{y}}_{\mathit{C}\mathit{T}\mathit{R}\mathit{L}}$ (m)	$\mathit{\delta }{\mathit{z}}_{\mathit{C}\mathit{T}\mathit{R}\mathit{L}}$ (m)
A	1.4	55.6	48.6	27.0	0.2
B	1.3	380.4	331.1	186.7	12.2
	${\mathit{m}}_{\mathit{P}\mathit{R}\mathit{O}\mathit{P}}$(kg)	${\mathit{k}}_{\mathit{P}}$	${\mathit{k}}_{\mathit{I}}$	${\mathit{k}}_{\mathit{D}}$
A	13.3	250	70	1000
B	6.6	70	7	693.7

, where the values of the PID coefficients are also displayed.

It is possible to see, in particular by the vertical component of the error in Figure 10, how case A presents higher frequency oscillations, reflecting the more aggressive behaviour of this tuning condition. This high responsiveness of the controller impacts largely the propellant mass expenditure, where it is possible to see that the difference in propellant mass between the two cases is reached in the first minute of controlled action.

The navigation error during the thrusted phase is characterised by a Root Mean Squared Error (RMSE) of around $2.5$ $\mathrm{m}$ in the vertical component, $9.7$ $\mathrm{m}$ and $2.7$ $\mathrm{m}$ in the along-track and cross-track components respectively. The associated control errors show $122.6$ $\mathrm{m}$, $1273.5$ $\mathrm{m}$, and $259.5$ $\mathrm{m}$ for case A and $172.1$ $\mathrm{m}$, $3049.9$ $\mathrm{m}$, and $37.6$ $\mathrm{m}$ for case B, in the vertical, along-, and cross-track components, respectively. It is possible to highlight that, in both cases, most of the error is present in the along-track component, due to the larger velocities involved by the orbital motion in this direction.

Overall, we may say that the GNC performance is satisfactory for the nominal scenario under study, giving with both tuning cases the possibility to follow the prescribed guidance, with errors below 500 $\mathrm{m}$, a threshold considered reasonable for the end of the main braking phase, before the terminal descent.

4.4. Monte Carlo Analysis

To consolidate the landing scenario, the sensitivity of the filter performance and the robustness of the two controllers are assessed with a Monte Carlo analysis on a set of 500 trajectories. These sample trajectories are generated by perturbing both the initial condition and the start time of the landing phase, $t_S$. In particular, by changing $t_S$, the visibility condition of the LNS constellation is slightly changed as well, modifying the observability level accordingly. On the other hand, the error in the initial condition is added to perturb the nominal scenario and test the robustness of the GNC chain.

The landing starting time, $t_S$, is sampled from a uniform distribution in the range 3000$÷$5000 $\mathrm{s}$, recalling that the nominal scenario sees a $t_S$ equal to 3400 $\mathrm{s}$. In practice, this is imposed by shifting the initial true anomaly of the orbit backward or forward in time, such that, in the nominal case, at $t=t_S$, the lander is at the periselene of its orbit when the thrusting begins. The error dispersion in the initial condition is instead introduced as additive zero-mean white Gaussian terms with a standard deviation of 100 $\mathrm{m}$ and $0.1$ $\mathrm{m}$/$\mathrm{s}$ in magnitude per each position and velocity component, respectively. This small perturbation of the state at the beginning of the simulation becomes quite high when propagated forward up to $t=t_S$. This directly impacts the control error with respect to the prescribed guidance that the PID regulator will face at the beginning of the thrusted arc, hence making convergence more and more difficult.

A summary of the simulation settings is reported in

Table 2. Summary of the navigation simulation settings.

Parameter	Value
LCNS SISE ($1\sigma$)	Position ${\epsilon }_r$ (x, y, z): 15 m
	Velocity ${\epsilon }_v$ (x, y, z): 0.15 m/s
	Clock bias ${\epsilon }_{\rho }$: 10 m
	Clock drift ${\epsilon }_{\dot{\rho }}$: 0.1 m/s
User clock Allan variance	$h_0$: 2 × 10−25
	$h_{-2}$: 6 × 10−25
Altimeter noise ($1\sigma$)	1% of $\zeta$
IMU noise ($1\sigma$)	1 $\mathsf{\mu }$g
Initial filter uncertainty ($1\sigma$)	Position ${\sigma }_r$ (x, y, z): 1 km
	Velocity ${\sigma }_v$ (x, y, z): 10 m/s
	Clock bias ${\sigma }_{b_c}$: 100 m
	Clock drift ${\sigma }_{d_c}$: 1 m/s
Filter rate	10 Hz
Landing start time $t_S$	Uniform distribution in (3000 s, 5000 s)
Initial state perturbation ($1\sigma$)	Position (x, y, z): 100 m
	Velocity (x, y, z): 0.1 m/s

, where both stochastic variables and Monte Carlo dispersed parameters are summarised.

When the 1 $\mathrm{k}$$\mathrm{m}$ altitude threshold is reached, the simulations are stopped and the final position errors with respect to the prescribed guidance trajectory are recorded. In order to evaluate the quality of the provided GNC architecture and compare the two controllers, each sample of the Monte Carlo analysis is labelled as a successful or failed sample in accordance with its final error. Considering a threshold of 500 $\mathrm{m}$ of control error, all samples with a larger final control error are classified as failed samples, while the remaining ones as successful.

Figure 11 and Figure 12 provide the logarithmic evolution of the navigation (upper plot) and control (lower plot) position errors for the successful samples during the thrusted arc only, for case A and B, respectively. The colour of the lines reflects the value of $t_S$, from the lowest in dark violet to the highest in yellow. The cyan lines represent the average (solid) and the $3\sigma$ upper bound (dashed), considering the represented lines.

Table 3. Collection of performance for the two proposed controllers in the Monte Carlo analysis. Statistical distribution of navigation and control errors are reported for successful samples, alongside the success rate.

CTRL	Success Rate	${\mathit{\mu }}_{\mathit{\delta }{\mathit{r}}_{\mathit{N}\mathit{A}\mathit{V}}}$ (m)	$3{\mathit{\sigma }}_{\mathit{\delta }{\mathit{r}}_{\mathit{N}\mathit{A}\mathit{V}}}$ (m)	${\mathit{\mu }}_{\mathit{\delta }{\mathit{r}}_{\mathit{C}\mathit{T}\mathit{R}\mathit{L}}}$ (m)	$3{\mathit{\sigma }}_{\mathit{\delta }{\mathit{r}}_{\mathit{C}\mathit{T}\mathit{R}\mathit{L}}}$ (m)
A	51.2%	1.8	2.8	56.5	20.1
B	86.0%	1.7	2.8	388.9	22.8

presents the aggregated results for navigation and control errors in the two cases.

As we can see by looking first at the line colour distribution, in both cases there is not a specific correlation between the landing start time, $t_S$, and the navigation and control errors, $\delta r_{\mathrm{NAV}}$ and $\delta r_{\mathrm{CTRL}}$. Overall, the control error in particular is quite uniform once convergence is reached.

$\delta r_{\mathrm{NAV}}$ is, as expected, comparable for the two cases, providing final errors below 5 $\mathrm{m}$, also in the worst cases. The control error instead varies in accordance with the nominal performance. As one can see in Table 3, the average values are very close to the nominal values for both case A and case B. Both cases also show a similar value for the standard deviation, around ∼$20\mathrm{m}$.

Regarding case A (Figure 11), a quite unsatisfactory success rate of $51.2\%$ is reached, while a much larger value of $86.0\%$ is obtained with the milder controller of case B. This reflects how the larger gains of the more aggressive tuning are less reliable and robust when a dispersion from the nominal scenario is introduced.

An interesting visualisation of the successful trajectories is obtained by looking at the final points observed on the horizontal plane of the Moon-fixed IAU_MOON frame at 1 $\mathrm{k}$$\mathrm{m}$ altitude. This is what the scatter plots in Figure 13 represent, together with two sets of covariance ellipses.

The upper plot reports the distributions of the two cases together, with the aim point of the South Pole visible. The two lower plots instead present the distribution of each case individually. In all representations, the grey ellipses represent the $3\sigma$ covariance of the estimated trajectory in the nominal scenario, while the remaining three ellipses represent the three sigma levels of the successful Monte Carlo samples. The aspect ratio and the orientation of nominal scenario covariance ellipses (the grey ones) are comparable with those of the samples, thus providing a sensible expectation of the landing mission dispersion ratios among the various components. The additional increment of the area is associated with the introduced control error.

The $3\sigma$ uncertainty ellipses of these sample is characterised by a semi-major and semi-minor axes of $11.7$ $\mathrm{m}$ and $4.2$ $\mathrm{m}$ for case A, and $14.0$ $\mathrm{m}$ and $1.9$ $\mathrm{m}$ for case B. In both cases, the semi-major axis is mostly aligned with the X_{IAU_MOON} axis.

This elongation is consistent with the nominal scenario errors presented in the previous section, where most of the errors and dispersion were in the along-track component.

Such effects, reflected in the navigation error as well, are also caused by the geometry involving the user and some of the servicers that are visible in this specific time-frame. Their LoS with respect to the South Pole (where the lander is practically placed) during the few minutes of thrusted arc have a relatively large variation along the Y_{IAU_MOON} axis and a practically null one along the X_{IAU_MOON} axis, causing the y component to be much more observable, thus with a reduced error.

Analysing instead the different failed samples, it is possible to say that the vast majority of them failed by crossing the 1 $\mathrm{k}$$\mathrm{m}$ altitude prematurely. Even though such an outcome is not desirable, it is possible to say that the failure is due to very high state discrepancies between the desired and actual conditions at the beginning of the thrusted phase. Indeed, by looking at Figure 14, it is possible to see the distribution of such an error at $t_S$.

In particular, the two plots, the upper for case A and the lower for case B, present histograms with differentiation between successful and failed samples with respect to the control error at the start of the landing phase. Such errors are the consequence of the imposed initial state dispersion to each Monte Carlo sample, which, as previously highlighted, build up to high values quite easily. The outcome of this visualisation is able to assess that, for the more aggressive controller, the basin of convergence is limited to an initial error between 50 $\mathrm{k}$$\mathrm{m}$, while for the milder controller of case B, such a boundary is much larger, around 130 $\mathrm{k}$$\mathrm{m}$. Both boundary values are very large, hinting that the proposed approach is in any case feasible to provide reliable navigation to landing missions, even though further optimisations in the control algorithms are still possible.

The histogram presented in Figure 15 reports instead the propellant mass consumption of the successful samples for case A and case B, stressing again how the tuning parameters of case B provide less expensive manoeuvres, with an upper boundary of approximately ∼$30 \mathrm{k}\mathrm{g}$, in contrast to the ∼$50 \mathrm{k}\mathrm{g}$ of case A.

With such Monte Carlo analysis, we gained clear insights on the statistical behaviour of the landing GNC strategy for both cases, which provide reasonable alternatives in the trade-off between landing precision and robustness. In order to widen the robustness analysis of such a GNC strategy, the Monte Carlo samples could be generated to include other additional dispersed parameters. In particular, the presented analysis, with the goal of being representative of a safe landing operation, considers only scenarios with very good LNS visibility, with a total number of servicers reaching five before starting the propelled phase. Such conditions are obviously forsaken when operationally such a critical stage is planned; however, looking at the behaviour of this GNC architecture also at different epochs, thus with a worse visibility, could help in understanding the limitations of such an approach and how to improve it.

5. Conclusions

The current paper analysed the performance obtained by a GNC scheme employing GNSS-like signals from a Lunar navigation constellation to assist a landing mission. The current work is presented as a continuation of previous manuscripts analysing the topic of Lunar navigation constellations [14,16,20], where the specific case under study considers a landing mission to the Lunar South Pole, a target of many programmed and identified space missions.

In particular, the lunar constellation optimised in [14] was exploited as a baseline service for the navigation task of a South Pole landing mission. The presented GNC chain comprises an offline optimised guidance to be used as a set-point for a reference tracking PID controller. Two different controller tunings have been considered, one targeting high landing precision (case A) and the other low propellant expenditure (case B). The nominal results provided in both cases demonstrate very good performance, with final navigation errors below 2 $\mathrm{m}$ and control errors around 56 $\mathrm{m}$ and 380 $\mathrm{m}$ for case A and case B, respectively. Naturally, given the much larger control errors with respect to the navigation ones, in both cases further improvements in the final landing precision can be obtained with alternative and more refined control laws. A Monte Carlo analysis was performed to assess the robustness of the proposed solutions and verify their success rate, considering a 500 $\mathrm{m}$ accuracy threshold in the horizontal plane at 1 $\mathrm{k}$$\mathrm{m}$ altitude. The high-performing case A provided a relatively low success rate of 51%, against the 86% of case B, which also resulted in the lowest fuel mass expenditure. For this reason, it is overall advisable to exploit the lower-performing tuning parameters.

Future works should look in the direction of improving the guidance and control strategy, with the goal of combining the good pinpoint landing performance of case A with the low propellant mass expenditure and the high robustness of case B.