Trajectory Optimization and Control Applied to Landing Maneuvers on Phobos from Mars-Phobos Distant Retrograde Orbits

1. Introduction

In order to facilitate the exploration of the Moon and Mars, intermediate staging locations have been proposed. The Lunar Gateway represents an intermediate orbital platform capable of linking, organizing, and supporting missions between the Earth and the Moon [1,2,3]. However, a lesser amount of literature regarding the use of similar platforms in the Martian environment exists. Such Martian platform would be referred to as Mars Base Camp [4]. As well as the Lunar Gateway, Mars Base Camp is a proposed platform that would enable and aid exploration of Mars and its surroundings, including its moons, Phobos and Deimos [5,6,7]. Although a final location for Mars Base Camp has not been chosen, it has been proposed that it could be located in the vicinity of Phobos or Deimos [8], such as a Mars-Phobos Distant Retrograde Orbit (MP DRO) [9]. This would avoid locating Mars Base Camp deep in Mars’ gravity well while keeping the Martian moons an accessible target. Some studies have shown the effectiveness of intermediate staging locations such as the aforementioned Mars Base Camp to facilitate missions in the vicinity of Mars and its moons [10,11]. Thus, the establishment of an infrastructure capable of performing In-Situ Resource Utilization (ISRU) of the material found in the regolith and under the surface of Phobos would prove to increase human and robotic exploration capabilities on Mars and its vicinity. For the purpose of this research, we assumed the existence of a Phobos Base capable of resupplying Mars Base Camp located in a Distant Retrograde Orbit (DRO). We analyzed the necessary trajectories required to transport material from the surface of Phobos to Mars Base Camp. We applied trajectory design, optimization, and control to an orbital transfer in the Circular Restricted Three-Body Problem (CR3BP), specifically from a Mars-Phobos Distant Retrograde Orbit to the surface of Phobos [12]. Given an MP DRO and a landing location on Phobos, the goal is to find the trajectory for which total $\Delta v$ for a direct 2-burn maneuver is minimum [8,13]. Here, the total $\Delta v$ corresponds to the sum of $\Delta v$’s required to initiate the transfer in the MP DRO and complete the transfer, i.e. land on the surface of Phobos. Considerations on Time-Of-Flight (TOF) were done in order to ensure that any necessary phasing and/or repositioning maneuvers are taken into account. A 2-D analysis is sufficiently accurate since we considered the Mars-Phobos orbit plane as the reference frame, and motion will mostly happen in such plane. However, further analysis will take into account a 3-dimensional space problem in which the gravitational potential of Phobos is also considered [14]. Furthermore, the initial position among the departing DRO will be implemented as the third optimization parameter, given that the minimum $\Delta v$ does not necessarily occur at a fixed position along the orbit for a given landing location. In this paper, a Particle Swarm Optimization (PSO) is used to determine the initial MP DRO given two parameters, $[A_x,V_y]$. Once the initial DRO is found, a trajectory optimization through a PSO is performed in order to determine the lowest $\Delta v$ trajectory to land a specific location on Phobos’ surface, given three parameters to optimize $[\alpha ,\varphi ,\Delta v]$. Here, we also discuss which PSO parameters are used for the analysis and which ones lead to the most efficient (i.e., the fastest) computations, that can be used for other optimization problems in the CR3BP. For further details on such parameters, refer to the nomenclature session at the beginning of the paper or to the detailed Section 2.3 and Section 3.

Phobos is the largest of two moons orbiting Mars, with an approximate size of $13.4$ km $\times 11.2$ km $\times 9.2$ km [15,16]. Phobos has an unusual shape (somewhat similar to a potato) which creates numerous difficulties for spacecrafts to orbit and therefore land on the moon itself. In the literature, the gravitational potential and the density of Phobos have been of high interest, since it becomes essential for a precise and complete analysis to acknowledge such characteristics, especially for long-duration stability of orbits in the vicinity of the Moon [17,18,19].

In this paper, the following approximations are made. Phobos orbits Mars on an elliptical orbit with an eccentricity of $e=0.0151$. For this study, we modeled Mars and Phobos using the CR3BP, which assumes that the eccentricity is 0. The semi-major axis of such orbit (i.e. average distance between Mars and Phobos) is 9376 km. Finally, the orbit’s inclination with respect to Mars’ equatorial plane is around 1°. Since Phobos’ Sphere of Influence (SOI) is below its physical surface, it becomes impossible to orbit the moon in the classical Keplerian sense. In the CR3BP, considering Phobos and Mars as the two primary masses, some periodic orbits exist, such as Halo Orbits or Lyapunov Orbits, as shown in Figure 1. However, such periodic orbits are either highly unstable or come dangerously close to the surface [20]. That being considered, the family of Distant Retrograde Orbits is of use in the following analysis. Such family of orbits proved to be “far” enough from Phobos’ surface and to be sufficiently stable as well, up to multiple hundreds of years, when considering a full-force model [20]. The word “far” is here used loosely, since these orbits are within tens of hundreds of km from Phobos.

Mission analysis for Earth to Mars-Phobos DROs is of high interest when associated with human and robotic exploration [4,10,11,16]. Phobos could prove itself extremely useful for future Mars missions as a popular location for future ISRU plants. It is widely recognized that having such in-orbit and in-situ infrastructures would make space exploration much more affordable, sustainable, and effective [16,21,22]. Several advantages of such stations include the ability not to bring all the necessary material, including propellant, for deep space missions [23]. Such projects are already in development (such as NASA’s Artemis [24]) for what concerns Earth’s Moon: building a refueling in-orbit station along with several on-the-surface infrastructures can easily be imagined as an intermediate step to what would be a deep space exploration mission.

In the following analysis, optimization algorithms are developed in order to generate a specific DRO given a desired distance from Phobos and optimize a landing maneuver to a specific landing site on the surface of Phobos with respect to the $\Delta v$ needed, i.e. minimize the necessary impulsive maneuver applied for landing.

2. Background Theory

For the purpose of giving a better insight of the problem, a brief overview of the basic aspects regarding Distant Retrograde Orbits and the Particle Swarm Optimization method are given.

2.1. Distant Retrograde Orbits

Distant Retrograde Orbits (DROs) are large orbits around the smaller primary body in the 3BP. Perfectly periodic DROs exist in the CR3BP and in-plane velocity perturbations create quasi-periodic orbits. Such orbits are favorable when looking for quarantine orbits since their stability demonstrates a noteworthy capability to resist perturbations [22]. Indeed, DROs can stay quasi-periodic when larger perturbations than other families of three-body orbits can withstand are applied. DROs have been recently proposed to be the ideal family of orbits when it comes to locating in-space infrastructures, and this is due to their outstanding stability and ease of access in terms of a gravity well [22]. DROs are usually characterized by the so-called x-amplitude ($A_x$), which represents the larger distance from $m_2$ (in a CR3BP) in the x-axis direction, given a $xyz$ rotating reference frame orbiting the primaries, as shown in Figure 2. In such analysis, it is typical to use a frame of coordinates with its center located at the barycenter of the primary masses [16]. According to this, several essential parameters can be determined in order to define a DRO. Such parameters are the radius at periapsis, and its corresponding velocity, required at the intersection of the $xz$-plane (i.e. where $y=0$) in order to create a periodic orbit. Since the radius has only one component along the x-axis, which is $A_x$-dependent, we can say that the essential parameters necessary to define a DRO are $[A_x,V_y]$, where $V_y$ indicates the y-component of the velocity vector at such location (which is in facts the only component, the x-component and z-component of the velocity vector at the considered starting location are zero). DROs have been identified and studied in the literature, for both Earth-Moon DROs and Mars-Phobos DROs [16,25], as shown in Figure 3. However, in this study a PSO algorithm will be used to determine a DRO given its x-amplitude, which will be the starting point for the landing trajectory optimization.

2.2. Particle Swarm Optimization

The Particle Swarm Optimization method is used in this study as an optimization algorithm to generate periodic DROs and compute $\Delta v$-optimal landing trajectories. PSO is a heuristic method used several times in previous astrodynamics problems [26,27,28]. It is based on the unpredictable motion of bird flocks while searching for food, taking advantage of the mechanism of information sharing that affects the overall behavior of a swarm [29,30]. The initial population that composes a swarm is randomly generated at the first iteration of the algorithm. Each particle is associated with a position vector and a velocity vector. In the PSO terminology, the words “position” and “velocity” are not to be interpreted in the classical sense (i.e. position and velocity as intended in classical physics and mechanics). Instead, they indicate the position vector containing the unknown parameters to optimize and the velocity vector determining the position update. Each particle represents a possible solution to the problem. At the end of the iterations, the particle corresponding to the optimal solution is selected, i.e. the particle that minimizes the parameters of interest, such as $\Delta v$. At each iteration, the position vector’s elements move in the velocity vector’s corresponding direction, to a new position, which represents a new possible solution. The updating of position and velocity is based on the objective function evaluation at the end of every iteration. The objective function (henceforth referred to as “cost function”) is an expression that needs to be minimized (or maximized). Once this condition is satisfied or the maximum number of iterations set $N_{max}$ is reached, the algorithm will have obtained an optimal solution. At every iteration, the best cost function is used to determine how to update both velocity and position. Indeed, with the approach formulated above, the algorithm will converge on the best solution in every iteration. For each particle, the formula for velocity update includes three components with stochastic weights. In the PSO terminology, these are known as the inertial ($c_I$), cognitive ($c_C$), and social ($c_S$) components. The inertial parameter is usually chosen randomly, but it has appeared to be proportional to each particle’s velocity in the previous iteration for some applications [28]; the cognitive parameter is based on the best position experienced by the particle; finally, the social component is direct toward the personal best position (i.e. the best location yet located by any particle in the swarm). The algorithm terminates when the maximum number of iterations (user decision) is reached. For this study, the following values are used in the algorithm. The function rand stands for a single uniformly distributed random number in the interval (0,1).

$$\begin{array}{cc}\hfill c_I& =\frac{(1+rand)}{2}\hfill \\ \hfill c_C& =1.49445·rand\hfill \\ \hfill c_S& =1.49445·rand\hfill \end{array}$$

(1)

The constants in Equation (1) are derived from the literature and optimized for this problem [31,32]. The implementation of the PSO was run in MatLab where the function unifrnd was used to assign the initial values at each particle element for each particle. The function unifrnd creates a uniformly distributed array of random values included in a specific range as given by the boundary values in Table (Table 2). Once that is implemented, the initial position is assigned

$${\overrightarrow{r}}_0=[A_x+(1-\lambda )]\widehat{i}+0\widehat{j}$$

(2)

where

$$\lambda =\frac{m_{Phobos}}{(m_{Mars}+m_{Phobos})}$$

(3)

represents the mass ratio. Note that, as already mentioned, the following analysis will be conducted assuming 2-D motion. Therefore the $i^{th}$ particle (i.e. the particle corresponding to the $j^{th}$ iteration) is

$$P_i=[V_y\left(i\right);T\left(i\right)]$$

(4)

so that

$${\overrightarrow{v}}_0=0\widehat{i}+P_i\left(1\right)\widehat{j}$$

(5)

The PSO algorithm here used is summarized in the flowchart in Figure 4. The coordinate system used is represented in Figure 5. Note that these two parameters are picked for the MP DRO analysis described in the Mars-Phobos DRO section below and are here used to better describe the general problem and give the reader a full comprehension of the general approach. Once the Initial Conditions (ICs) are determined (${\overrightarrow{r}}_0$ and ${\overrightarrow{v}}_0$), numerical integration is performed in order to compute the resulting trajectory using the CR3BP equations:

$$\begin{array}{cc}\hfill \ddot{x}-2\dot{y}-x& =-\frac{(1-\lambda )(x+\lambda )}{r_1^3}-\frac{\lambda (x-1+\lambda )}{r_2^3}\hfill \\ \hfill \ddot{y}+2\dot{x}-y& =-\frac{(1-\lambda )y}{r_1^3}-\frac{\lambda y}{r_2^3}\hfill \\ \hfill \ddot{z}& =-\frac{(1-\lambda )z}{r_1^3}-\frac{\lambda z}{r_2^3}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill r_1& =\sqrt{{(x+\lambda )}^2+y^2+z^2}\hfill \\ \hfill r_2& =\sqrt{{(x+\lambda -1)}^2+y^2+z^2}\hfill \end{array}$$

(7)

Once the integration is done, if the final position and velocity of the resulting trajectory are within a small tolerance (here ${10}^{-10}$) to the initial conditions ${\overrightarrow{r}}_0$ and ${\overrightarrow{v}}_0$, then a periodic orbit is achieved and the optimization algorithm has reached a final solution.

Therefore, the cost function, i.e. function to minimize, is

$$\begin{array}{cc}\hfill J& =\mid {\overrightarrow{r}}_{final}-{\overrightarrow{r}}_0\mid +\mid {\overrightarrow{v}}_{final}-{\overrightarrow{v}}_0\mid \hfill \\ \hfill & =\mid r_{final}\widehat{i}-[A_x+(1-\lambda )]\mid +\mid v_{final}\widehat{j}-V_y\left(i\right)\mid \hfill \end{array}$$

(8)

Satisfying the condition in Equation (8) results in having a periodic orbit, i.e. $\overrightarrow{r}\left(0\right)=\overrightarrow{r}\left(t_f\right)$ and $\overrightarrow{v}\left(0\right)=\overrightarrow{v}\left(t_f\right)$. Once the cost function is evaluated, each particle velocity (intended in the PSO terminology) is updated as follows

$$\begin{array}{c}\hfill V_i=c_I{V}_i+c_C(P_{best}-P_i)+c_S(G_{best}-P_i)\end{array}$$

(9)

where $P_{best}$ is the personal best of each particle, and $G_{best}$ is the overall best position between all particles and all iterations as of the $i^{th}$ iteration. The coefficients $c_S,c_C,c_I$ are defined in Equation (1). Note that the subscript i in the Equations (4, 8, 9) indicates the $i^{th}$ element at the $j^{th}$ iteration of the algorithm. Naturally, the greater number of particles initialized, the more probability the algorithm has to achieve an optimal result. Similarly, a higher number of maximum iterations has more probability to converge to a solution. Nevertheless, a compromise between computational complexity and the number of iterations and particles needs to be made. In fact, increasing either the number of particles or the maximum iterations allowed can make the computational time increase significantly. However, due to the random nature of the algorithm, one cannot guarantee that an absolute optimal solution can be reached. On the other hand, a `good’ initial guess is not necessary in order to initiate the algorithm. Furthermore, unlike gradient-based methods, heuristic optimization methods are capable of finding absolute maxima/minima regardless of ICs [31].

2.3. Mars-Phobos DRO

Several DROs have been defined in the literature [12,16] for given $A_x$ values. Here a PSO algorithm has been implemented to generate DROs. The optimization was implemented starting from the two parameters [T,$V_y$], where $V_y$ was defined in Section 2 and T represents the orbital period of the DRO. The initial data are summarized in Table 1.

Starting from data taken from [16] the value of $V_y$ for $A_x=100$ km is known to be $-0.045620256764708$ km/s, that identifies a DRO with an orbital period $T=2.731044880670166\times {10}^4$ s [33]. Knowing such parameters helps the implementation of boundary values for the two parameters used in this analysis (i.e. the maximum and minimum values that the parameters can assume in the particles). In particular, the initial conditions for the optimization are listed in Table 2. Similarly, any other DRO of interest can be derived, although only one is here chosen for the results shown.

Table 1. Orbital parameters for Mars-Phobos and initial values for DRO Optimization.

Variable	Value	Description
${\mu }_{Phobos}$ (km3 s−2)	$7.11358812096305·{10}^{-4}$	Gravitational parameter for Phobos
${\mu }_{Mars}$ (km3 s−2)	$42828.375214$	Gravitational parameter for Mars
$\lambda$	$1.660952106463386·{10}^{-8}$	Mass ratio
R (km)	9376	Average distance Mars-Phobos
$A_x$ (km)	125	DRO Amplitude

Table 1. Orbital parameters for Mars-Phobos and initial values for DRO Optimization.

Variable	Value	Description
${\mu }_{Phobos}$ (km3 s−2)	$7.11358812096305·{10}^{-4}$	Gravitational parameter for Phobos
${\mu }_{Mars}$ (km3 s−2)	$42828.375214$	Gravitational parameter for Mars
$\lambda$	$1.660952106463386·{10}^{-8}$	Mass ratio
R (km)	9376	Average distance Mars-Phobos
$A_x$ (km)	125	DRO Amplitude

Table 2. Boundary Values for DRO Optimization.

Variable	Value	Description
$B_{lo}$($V_y$)	$\frac{3}{2}{V}_y$	Lower limit for $V_y$ (km/s)
$B_{up}$($V_y$)	$V_y$	Upper limit for $V_y$ (km/s)
$B_{lo}$($\alpha$)	T	Lower limit for T (s)
$B_{up}$($\alpha$)	$2\pi$	Upper limit for T (s)
$N_{max}$	80	Maximum number of iterations
$N_{particles}$	40	Number of particles

Table 2. Boundary Values for DRO Optimization.

Variable	Value	Description
$B_{lo}$($V_y$)	$\frac{3}{2}{V}_y$	Lower limit for $V_y$ (km/s)
$B_{up}$($V_y$)	$V_y$	Upper limit for $V_y$ (km/s)
$B_{lo}$($\alpha$)	T	Lower limit for T (s)
$B_{up}$($\alpha$)	$2\pi$	Upper limit for T (s)
$N_{max}$	80	Maximum number of iterations
$N_{particles}$	40	Number of particles

The parameters to optimize are the initial velocity $v_0$ (Equation (5)) at position $r_0$ (Equation (2)) and the orbital period T (Equation (4)). In the literature, different DROs have been identified for $A_x\in [15,300]$ km [16]. DROs with $A_x$ values below 15 km result in a certain impact with Phobos’ surface and the ones with higher $A_x$ values than 300 km are considered too far away to maintain the characteristics in which we are interested (vicinity to Phobos and low-$\Delta v$ for landing). Given the initial orbital parameters for Mars-Phobos and the initial values for the DRO, specified in Table 1, the resulting DRO is shown in Figure 6, and its reference values for position and velocity are summarized in Equation (10).

$$\begin{array}{cc}\hfill \overrightarrow{r}& =[9500.920971466154,-0.006046323363597]km\hfill \\ \hfill \overrightarrow{v}& =[0.0007100917133235669,-0.056800560873732]km/s\hfill \\ \hfill T& =27310.98049729574s\hfill \end{array}$$

(10)

where T represents the total orbital period, and v is the velocity vector at position r. Given the $N_{max}$ and $N_{particles}$ specified in Table 2, the result was achieved with a minimum cost function J in the order of ${10}^{-4}$ approximately, as shown in Figure 7, so the final values are such that $\mid \overrightarrow{r_f}-\overrightarrow{r_0}\mid \approx$0.1 mm and $\mid \overrightarrow{v_f}-\overrightarrow{v_0}\mid \approx$0.1 mm/s . The total average computation time was ∼20 seconds. Several other approaches were tried, in terms of the number of particles and the number of maximum iterations, but these parameters proved to give a sufficiently accurate result in a reasonable amount of time. If, for instance, we doubled the $N_{max}$ up to 160, the computation time increases up to ∼40 seconds, generating a minimum cost function J in the order of ${10}^{-5}$. So doubling the computational time will generally give a result that is just one order of magnitude more accurate than the initial one. As already outlined, one has to compromise between computational efficiency, computational speed, and results’ accuracy.

Many computational issues may occur while implementing the algorithm. However, parallel programming and computing is a suitable solution to speed up the computation time. In fact, every possible trajectory computed at each iteration is independent with respect to all the other trajectories, until every particle is updated. Note that for this study a 2016 MacBook Pro is used. Hardware and Software details are summarized in Table 3.

3. Landing Trajectory Optimization

In this section, the trajectory optimization for multiple landing maneuvers is analyzed through the use of PSO. Given an initial stable parking DRO around Phobos like the example computed in Section 2.3, and a desired landing location, the analysis aims to find the minimum $\Delta v$ necessary to encounter that specific location. The problem is computationally more complex than the one described in the MP DRO Section and it presents several more difficulties in the formulation. First of all, the choice of the parameters to optimize is non-trivial and there are different ways to approach the optimization. Here, the following parameters were chosen, minimizing the dimensionality of the trade space to three:

$$P_i=[\alpha \left(i\right);\Delta v\left(i\right);\varphi \left(i\right)]$$

(11)

where the element $\Delta v\left(i\right)$ represents the $\Delta v$ to apply, and $\alpha \left(i\right)$ indicates the angle with respect to the x-axis at which the $\Delta v$ is applied. The parameter $\varphi \left(i\right)$ represents the angle with respect to the x-axis which consequently identifies a location on the DRO $r\left(\varphi \right)$, as shown in Figure 8. Note that in this analysis, the index i represents the $i^{th}$ particle at the $j^{th}$ iteration. This analysis requires more computational time than the DRO optimization because there are 3 parameters to optimize instead of 2, and we do not have specific analytic expressions to relate them to each other. The starting DRO is equal to the one defined in the optimization described in Section 2.3. Therefore, as obtained in Equation (10), the initial values and the particle boundary values are summarized in Table 4 and Table 5 respectively.

In this particular optimization, a penalty scaling factor $\gamma$ was introduced. The addition of such a factor is necessary in order to prioritize the terms of the cost function [34]. In fact, the essential condition for a valid outcome of the optimization is to land at the desired location on Phobos. Indeed, a trajectory is considered valid only if it intersects the surface of Phobos. Secondarily, given that a valid trajectory to accomplish such a condition has been found, the total $\Delta v$ needs to be minimized. The expression for $\gamma$ was derived in the literature [34], and the following proved to be an effective value for this factor

$$\gamma \left(i\right)=\beta \left[{\gamma }_{initial}-\frac{{\gamma }_{initial}-{\gamma }_{final}}{N_{max}}\sqrt{(2N_{max}-i)i}\right]$$

(12)

where i is the current iteration number, $N_{max}$ is the maximum number of iterations and ${\gamma }_{initial}$ and ${\gamma }_{final}$ represent the initial and the final weighting values. In this analysis, they are equal to 1 and 0.1 respectively. The factor $\beta$ is an order scaling factor, set equal to 100 for this optimization, which is used to make sure that the location checking term has a higher order than the $\Delta v$ term. Therefore, the cost function becomes

$$J=\left\{\begin{array}{cc}\gamma \mid {\overrightarrow{r}}_{final}-{\overrightarrow{r}}_{landing}\mid +\Delta v_{total},\hfill & if\mid {\overrightarrow{r}}_{final}-{\overrightarrow{r}}_{landing}\mid >tol\hfill \\ \Delta v_{total},\hfill & if\mid {\overrightarrow{r}}_{final}-{\overrightarrow{r}}_{landing}\mid <tol\hfill \end{array}\right.$$

(13)

For this optimization, the tolerance (tol) is set to 0.1 km (100 meters). It is clear how the scaling factor prioritizes a specific term in Equation (13), given that the $\Delta v$ term is a relatively small value.

For a matter of simplicity, Phobos was ideally divided into 12 sectors, and a landing location within each sector was picked as a reference location for such area. This means that the $\Delta v$ required to achieve landing in such sector can be considered quantitatively similar to the reference location’s $\Delta v$. Results for the landing locations in each sector are summarized in Table 7. Computation times for each example in Table 7 are indicated in Table 6. As expected, the average computation time is slightly higher than the one for the DRO optimization described in the MP DRO analysis.

The landing trajectories for each sector are represented in Figure 10 and a heat map of all the landing locations is shown in Figure 9. To have a better insight into the results, a fully loaded SpaceX Starship is considered. SpaceX’s Starship will have an empty mass of ∼120 mt (260,000 lb) and a gross mass of ∼1,320 mt (2,910,000 lb) with a propellant capacity of ∼1,200 mt (2,600,000 lb)[35]. Recall the Tsiolkovsky rocket equation

$$\Delta v=v_{e}ln\left(\frac{m_0}{m_f}\right)=g_0{I}_{sp}ln\left(\frac{m_0}{m_f}\right)$$

(14)

where $m_0$ is the total mass of the spacecraft, including the propellant, $m_f$ is the final total mass assuming all propellant to be consumed (also known as dry mass), $I_{sp}$ is the spacecraft’s specific impulse and $g_0$ is the standard gravity at Earth’s sea level. We can now assume to land on Phobos’ surface at sector 7, which latitude corresponds to the Stickney Crater on the surface of Phobos and which corresponds to a total $\Delta v$ of 27.41 m/s (see Table 7). From Equation (14) we can estimate the necessary propellant mass to land a full-load Starship on Phobos.

We consider the spacecraft to empty its tanks at landing. Given a specific impulse of 380 s in a vacuum and a payload of 100 mt, Starship will achieve landing with a total propellant mass equal to

$$\begin{array}{cc}\hfill m_p& =m_f(e^{\frac{\Delta v}{g_0{I}_{sp}}}-1)\hfill \\ \hfill & =(120+100)·(e^{\frac{27.41}{9.81·380}}-1)\hfill \\ \hfill & \approx 1.62mt\hfill \end{array}$$

(15)

$$\frac{m_p}{m_f}=0.7\%$$

(16)

where $m_p$ is the total propellant mass used to achieve the indicated $\Delta v$. Since the landing $\Delta v$ is quantitatively equal to the $\Delta v$ required to go from the surface of Phobos to the initial DRO (the trajectory will just look different from a qualitative point of view), and we assumed Phobos to have a re-fueling station on its surface, the total propellant mass necessary to leave and arrive the DRO is simply two times the mass indicated in Equation (15). Therefore, assuming one travel per terrestrial day for a year (365 days) back and forth the Stickney Crater, the total propellant mass per year necessary to achieve such transfers is ∼1,185.22 mt of propellant. If we consider a 1,200 mt propellant capacity, Starship can achieve a back-and-forth transfer ∼369 times before refueling. Although it is challenging to precisely estimate the propellant price per tonne, SpaceX claims that refueling Starship would cost approximately $900,000. This means that a complete transfer would cost around $2,435, transporting a 100 mt payload.

4. Conclusions

This study led to valuable results for both DRO implementations and landing trajectory optimization. Several DROs have been defined in the literature for given $A_x$ values, but not much exists about landing trajectories from DROs to the surface of Phobos. Regarding the landing trajectory optimization, which is the main purpose of this study, it is clear from Table 7 how the $\Delta v$ is on average equal to ∼30 m/s. This value is reasonably low and makes this orbit feasible for an in-orbit refueling station or in general for a stationary orbit around Phobos which one can easily reach and leave. Several landing locations were reached, and the impulsive $\Delta v$ was always found to be in approximately the same range. Therefore, nearly every point on Phobos’ surface can be reached with a $\Delta v$ of approximately ∼30 m/s. Furthermore, as summarized in Table 6, the average computation time for such PSO algorithm is reasonable and makes the optimization feasible for future mission designs.

Mars’ moons exploration has been of high interest within the scientific community in the past years. The existence of periodic orbits around Phobos, such as DROs, makes the future design of in-orbit refueling and supporting stations more than feasible. Moreover, the ease of departing and arriving to such orbits with low $\Delta v$’s gives the possibility of exploring Phobos and establishing a Mars Base Camp, in order to explore Mars and its surroundings with great versatility and with lower flight times between transfers from and to Mars itself. Having stations or spacecrafts orbiting around Mars (in Low-Mars Orbits for example) would require continuous trajectory adjustments because of Mars’ gravity potential acting on the spacecraft/station, and will then result in significant $\Delta v$ used over the course of a long-term mission. Besides, given the distance between Earth and Mars, it would not be convenient to perform frequent trajectory adjustments or remote control operations, given an average transmission delay of ∼20 minutes between Mars and Earth. Hence, having a long-term stability solution is crucial, and DROs appear to be feasible for such an objective. The method here presented can be used for multiple Mars-Phobos DROs optimizations at different amplitudes. Plus, a Mars-Deimos DRO can utilize this very same approach.