Rapid Orbit-to-Orbit Transfer to Asteroid 4660 Nereus Using Solar Electric Propulsion

Abstract

This paper analyzes the rapid rendezvous trajectory of a spacecraft equipped with an advanced solar electric propulsion system towards asteroid 4660 Nereus. In this context, a set of possible minimum-time orbit-to-orbit transfer trajectories is calculated by modeling the propulsion system performance characteristics on those of NASA’s Evolutionary Xenon Thruster-Commercial (NEXT-C). In particular, the actual NEXT-C ion engine throttle table is used to calculate the optimal thrust control law that ensures the flight time is minimized for an assigned value for the spacecraft’s initial mass and the reference (electric) power at the beginning of the transfer. A baseline scenario that considers the actual inertial characteristics of the NASA’s DART spacecraft is analyzed in detail, and a parametric study is proposed to evaluate the transfer performance as a function of the main design parameters as, for example, the spacecraft’s initial mass and the reference power.

1. Introduction

Significant advances in ion thruster development over the past 25 years have contributed significantly to the planning of increasingly complex and ambitious missions, both in terms of achievable scientific goals and the ability to achieve more and better utilization of space resources [1]. So far, several missions using solar electric propulsion (SEP) as their main propulsion system have been carried out or are planned, including Deep Space 1, the first interplanetary mission propelled by SEP [2]; SMART-1, the first European mission to the Moon [3]; Dawn, the first mission to use SEP for a long-duration journey to explore Vesta and Ceres [4], and Psyche, to rendezvous with the largest metal asteroid in the solar system [5], to cite just a few.

The main reason for its success is that SEP makes many space missions attractive, as it can reduce the size of launch vehicles, and reduce the total propellant consumption, when compared with more conventional (i.e., chemical) propulsion systems. All this translates into significant mission cost savings [6]. The development of increasingly reliable SEP thrusters is also having an important impact from an economic point of view, making the concept of space resource exploitation increasingly advantageous, not only from an engineering perspective, but also from that of possible commercial investment. A first step in this direction is the recent missions that have demonstrated the possibility of extracting material from celestial bodies (usually near-Earth asteroids), collecting it, and returning it to Earth. Typical examples of asteroid mining include Hayabusa and its successor Hayabusa 2 and, more recently, OSIRIS-REx, which visited and collected a sample from asteroid 101955 Bennu [7,8,9]. These missions can be thought of as a first step towards a more ambitious goal of collecting large quantities of valuable substances that are hard to find on Earth, or even entire asteroids [10,11].

The starting point for studying the feasibility of reaching an asteroid is to design trajectories and analyze its accessibility considering commercial propulsion systems. Such a study is generally carried out using an optimal control approach, that is, minimizing an appropriate performance index, typically propellant mass, time of flight, or a combination of these two quantities [12]. Because the quality of the solution found is highly dependent on the mathematical model describing the operation of the thruster, it is important to use sufficiently accurate models to obtain reliable estimates of achievable performance. However, examples of studying optimal trajectories with realistic SEP throttle models are rarely addressed in the literature. An interesting approach in this regard is offered by the works of Taheri et al. [13,14], who studied the effects of a refined thrust model on the performance of a space mission.

In most cases the behavior of the SEP is mathematically described as a function of a single control variable, the input power to the power processing unit (PPU), which is assumed to be able to vary continuously over an assigned range of variation. Instead, since a SEP system has a finite number of operation points, the purpose of this paper is to evaluate the performance of an electric thruster with a realistic model, in which each operation point corresponds to given values of thrust, propellant mass flow rate, and PPU input power. In this regard, we updated a mathematical model, previously validated in other mission scenarios [15,16], by including the characteristics of a recent commercial-type propulsion system.

In particular, this paper analyzes minimum-time rendezvous trajectories to asteroid 4660 Nereus for a spacecraft equipped with a next-generation solar electric propulsion system. Asteroid 4660 Nereus, which approached Earth on 11 December 2011 with a minimum distance of about $4\mathrm{Mkm}$, has been considered in the past as a potential target for a robotic exploration mission [17,18], due to both its accessibility (in terms of delta-V) and its physical characteristics, which make it possible to maintain a closed orbit around such a small body in the solar system. It is classified as an E-type asteroid, possibly associated with aubrite meteorites, with a hight albedo (greater than 0.3) and a maximum diameter of $1.3\mathrm{km}$ [17] (see Figure 1).

Its orbit, with an inclination of $1.4\mathrm{deg}$ relative to the ecliptic plane, is illustrated in Figure 2.

Assuming asteroid 4660 Nereus to be the primary target of a deep space exploration mission, a set of possible minimum-time orbit-to-orbit transfer trajectories is calculated by considering the performance characteristics of NASA’s Evolutionary Xenon Thruster-Commercial (NEXT-C) [19,20]. The latter is the propulsion system that was installed onboard the Double Asteroid Redirection Test (DART) spacecraft [21,22], which successfully completed its mission at the end of September 2022. Specifically, the acceleration table of the NEXT-C engine, together with a simplified analytical model of the performance of the power subsystem, is used to calculate the optimal control law that guarantees the minimization of the time of flight for an assigned value of the initial spacecraft mass and reference (electric) power. A baseline scenario that considers the actual inertial characteristics of the DART spacecraft is analyzed in detail, while a parametric study is proposed to evaluate transfer performance as a function of the spacecraft’s initial mass and reference power.

2. Mission Schematization and Mathematical Model Description

This section describes the mathematical model used to analyze the interplanetary mission scenario. Specifically, it illustrates the spacecraft dynamics model, the spacecraft thrust vector scheme, and the trajectory optimization process that was employed to obtain the numerical results. The mathematical approach is largely based on the model proposed by the authors in Ref. [15], where the interested reader can find a detailed discussion of the assumptions made. A simplified version of this model has also been used to study a circle-to-circle transfer to the inner regions of the solar system [16]. In this context, a novelty with respect to the model already presented in the literature [15] is the specific description of the characteristics of the propulsive system (i.e., the detailed table of thrusts), which is explained later in this section.

2.1. Spacecraft Dynamics in the Interplanetary Space

The spacecraft dynamics during the transfer are described by a simplified model, commonly used in the preliminary phase of interplanetary mission design. Specifically, the external forces acting on the spacecraft during flight consist of the gravitational pull of the Sun and the thrust provided by the on-board propulsion system. The latter can be set to zero by turning off the thruster, so that the transfer trajectory may include a coasting arc (i.e., a phase in which the spacecraft motion is purely Keplerian). Since the mass m of the spacecraft varies with time t and using a three-dimensional model for interplanetary transfer, the spacecraft state vector has 7 scalar components, i.e., the instantaneous mass of the vehicle and the 6 elements that define the position and velocity of the spacecraft. Following the common practice for trajectory analysis and optimization [23,24,25], the modified equinoctial orbit element set (MEOE) (An interesting discussion on the use of MEOEs can be found at the following link, retrieved on 10 September 2023: https://spsweb.fltops.jpl.nasa.gov/portaldataops/mpg/MPG_Docs/Source%20Docs/EquinoctalElements-modified.pdf) proposed by Walker et al. [26] was used to describe the spacecraft position and velocity. The relationships between the set of MEOEs $\{p,f,g,h,k,L\}$ and the classical orbital elements are as follows

$$\begin{array}{c}f=ecos(\omega +\Omega ),g=esin(\omega +\Omega ),h=tan(i/2)cos\Omega ,\hfill \\ k=tan(i/2)sin\Omega ,L=\nu +\Omega +\omega ,p=a\left(1-e^2\right)\hfill \end{array}$$

(1)

where a is the semimajor axis, e is the eccentricity, i is the inclination, $\omega$ is the argument of perihelion, $\nu$ is the true anomaly, and $\Omega$ is the right ascension of the ascending node of the spacecraft osculating orbit. Note that p coincides with the osculating orbit semilatus rectum, while the others MEOEs are dimensionless. The set of classical orbit elements can be recovered by inverting Equation (1), as described in Ref. [27], and the result is

$$\begin{array}{c}a=\frac{p}{1-f^2-g^2},e=\sqrt{f^2+g^2},i=2arctan\sqrt{h^2+k^2},sin\omega =gh-fk,\hfill \\ cos\omega =fh+gk,sin\Omega =k,cos\Omega =h,\nu =L-\Omega -\omega \hfill \end{array}$$

(2)

Finally, the Sun-spacecraft distance r can be written as a function of MEOEs as

$$r=\frac{p}{1+fcosL+gsinL}$$

(3)

while the components of the spacecraft position $\mathit{r}$ and velocity $\mathit{v}$ vectors in a heliocentric-ecliptic (Cartesian) reference frame [28] of unit vectors $\{\widehat{\mathit{i}},\widehat{\mathit{j}},\widehat{\mathit{k}}\}$ can be obtained through the following relationships [23,29]

$$\begin{array}{c}\mathit{r}=\frac{r}{1+h^2+k^2}\left\{\left[\left(1+h^2-k^2\right)cosL+2hksinL\right]\widehat{\mathit{i}}+\right.\hfill \\ \left.\left[\left(1-h^2+k^2\right)sinL+2hkcosL\right]\widehat{\mathit{j}}+\left(2hsinL-2kcosL\right)\widehat{\mathit{k}}\right\}\hfill \end{array}$$

(4)

$$\begin{array}{c}\mathit{v}=-\frac{\sqrt{{\mu }_{\odot }/p}}{1+h^2+k^2}\left\{\left[\left(sinL+g\right)\left(h^2-k^2+1\right)-2hk\left(cosL+f\right)\right]\widehat{\mathit{i}}+\right.\hfill \\ \left.\left[\left(cosL+f\right)\left(h^2-k^2-1\right)+2hk\left(sinL+g\right)\right]\widehat{\mathit{j}}+\left[-2h\left(cosL+f\right)-2k\left(sinL+g\right)\right]\widehat{\mathit{k}}\right\}\hfill \end{array}$$

(5)

where ${\mu }_{\odot }$ is the Sun’s gravitational parameter.

Paralleling the procedure illustrated by Betts [29], and according to the model discussed in Ref. [15], the spacecraft state vector $\mathit{x}\in {\mathbb{R}}^{(7\times 1)}$ is defined as

$$\mathit{x}={\left[p,f,g,h,k,L,m\right]}^{\mathrm{T}}$$

(6)

so that the spacecraft dynamics is described by the vectorial equation

$$\dot{\mathit{x}}=\mathbb{A}{\mathit{a}}_p+\mathit{b}$$

(7)

where ${\mathit{a}}_p$ is the electric thruster-induced propulsive acceleration vector, $\mathbb{A}\in {\mathbb{R}}^{(7\times 3)}$ is a matrix whose entries depend on the MEOEs, while $\mathit{b}\in {\mathbb{R}}^{(7\times 1)}$ is a vector with two non-zero components only, which depend on both the propellant mass flow rate ${\dot{m}}_p\ge 0$ and a subset of MEOEs. In particular, according to Refs. [15,29] the expression of $\mathbb{A}$ is

$$\mathbb{A}=\left[\begin{array}{ccc}0& \frac{2p}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}& 0\\ sinL\sqrt{\frac{p}{{\mu }_{\odot }}}& \frac{\left(2+fcosL+gsinL\right)cosL+f}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}& -\frac{g\left(hsinL-kcosL\right)}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ -cosL\sqrt{\frac{p}{{\mu }_{\odot }}}& \frac{\left(2+fcosL+gsinL\right)sinL+g}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}& \frac{f\left(hsinL-kcosL\right)}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ 0& 0& \frac{\left(1+h^2+k^2\right)cosL}{2\left(1+fcosL+gsinL\right)}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ 0& 0& \frac{\left(1+h^2+k^2\right)sinL}{2\left(1+fcosL+gsinL\right)}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ 0& 0& \frac{hsinL-kcosL}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ 0& 0& 0\end{array}\right]$$

(8)

while the expression of $\mathit{b}$ is

$$\mathit{b}=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ \sqrt{{\mu }_{\odot }p}{\left(\frac{1+fcosL+gsinL}{p}\right)}^2\\ -{\dot{m}}_p\end{array}\right]$$

(9)

The two propulsive terms ${\mathit{a}}_p$ and ${\dot{m}}_p$ may be written as a function of the electric thruster characteristics, as described in the next subsection. Note that the components of the propulsive acceleration vector ${\mathit{a}}_p$ in Equation (7) are to be evaluated in a radial-tangential-normal (RTN) rotating reference frame whose unit vectors $\{{\widehat{\mathit{i}}}_{\mathrm{R}},{\widehat{\mathit{i}}}_{\mathrm{T}},{\widehat{\mathit{i}}}_{\mathrm{N}}\}$ are defined as [29]

$${\widehat{\mathit{i}}}_{\mathrm{R}}=\frac{\mathit{r}}{r},{\widehat{\mathit{i}}}_{\mathrm{N}}=\frac{\mathit{r}\times \mathit{v}}{\parallel \mathit{r}\times \mathit{v}\parallel },{\widehat{\mathit{i}}}_{\mathrm{T}}={\widehat{\mathit{i}}}_{\mathrm{N}}\times {\widehat{\mathit{i}}}_{\mathrm{R}}$$

(10)

The RTN reference frame is schematized in Figure 3, which is an adaptation of a similar figure that has recently appeared in Ref. [27]. Note that the plane $({\widehat{\mathit{i}}}_{\mathrm{R}},{\widehat{\mathit{i}}}_{\mathrm{T}})$ coincides with the plane of the spacecraft osculating orbit.

The initial condition required to complete the differential Equation (7) models the state of the spacecraft at the beginning of the interplanetary transfer, that is, at the initial time $t_0=0$. In an orbit-to-orbit transfer, it is assumed that the spacecraft has an assigned initial mass

$$m\left(t_0\right)=m_0$$

(11)

and that it initially travels through a Keplerian (parking) orbit with given characteristics. In this case, such a parking orbit coincides with Earth’s heliocentric orbit, so as to model a typical scenario in which the spacecraft leaves Earth’s sphere of influence with zero hyperbolic excess velocity relative to the planet [28]. The classical orbital elements describing the shape and orientation of the spacecraft osculating orbit at the initial time are obtained through the JPL Horizons system. Specifically, when using data from JPL Horizons (as of 1 August 2023), we have

$$\begin{array}{c}a\left(t_0\right)=1.49458051\times {10}^8\mathrm{km},e\left(t_0\right)=1.75733946\times {10}^{-2},i\left(t_0\right)=3.14715328\times {10}^{-3}\mathrm{deg}\hfill \\ \omega \left(t_0\right)=2.88554336\times {10}^2\mathrm{deg},\Omega \left(t_0\right)=1.75663719\times {10}^2\mathrm{deg}\hfill \end{array}$$

(12)

Note that Equation (12) does not contain the value of the spacecraft’s initial true anomaly $\nu \left(t_0\right)$, because in a typical orbit-to-orbit transfer the angular position of the spacecraft along the parking orbit is left unconstrained and will be computed by solving an optimization problem, as described below. Likewise, the initial value of the true longitude L is left free, while the others elements are calculated by using Equations (1) and (12). In this case, the result is

$$\begin{array}{c}p\left(t_0\right)=1.49411895\times {10}^8\mathrm{km},f\left(t_0\right)=-4.31625224\times {10}^{-3},g\left(t_0\right)=1.70350863\times {10}^{-2}\hfill \\ h\left(t_0\right)=-2.73854762\times {10}^{-5},k\left(t_0\right)=2.07656391\times {10}^{-6}\hfill \end{array}$$

(13)

The latter equation provides five (scalar) boundary constraints that will be used to complete the Two-Point Boundary Value Problem (TPBVP) associated with the optimization process described later in the section. Finally, the spacecraft’s mass $m_0$ at the beginning of the transfer must be considered as an assigned term (i.e., the value of $m_0$ is the sixth constraint of the TPBVP). In particular, the specific value of $m_0$ can be used to obtain the performance of the transfer parametrically.

2.2. Propulsive Characteristics and Thrust Vector Schematization

The propellant mass flow rate ${\dot{m}}_p$ and the propulsive acceleration vector ${\mathit{a}}_p$, i.e., the two “propulsive” terms that appear in the vector equation of motion (7), are usually modeled as a function of the on-board propulsion system. In this paper, the transfer performance is studied by considering an electric thruster modeled on the actual characteristics of the NEXT-C [30,31,32] that was installed aboard the DART spacecraft [22].

Specifically, the NEXT-C throttle table obtained from Refs. [33,34] was used to express the available thrust magnitude T, the nominal propellant mass flow rate ${\dot{m}}_p^n$, and the input power P to the PPU as a function of the throttle level $\mathrm{TL}\in {\mathbb{N}}^{+}$, which is a (positive) integer indicating the possible operating point of the electric thruster. According to [33,34], the NEXT-C has 40 possible throttle levels, which are summarized in

Table 1. Electric thruster throttle table adapted from the actual NEXT-C table.

TL	T  $\left[\mathbf{mN}\right]$	P  $\left[\mathbf{kW}\right]$	${\dot{\mathit{m}}}_{\mathit{p}}^{\mathit{n}}$ $[\mathbf{mg}/\mathbf{s}]$
1	25	0.545	1.85
2	31	0.668	2.05
3	37	0.788	2.05
4	48	1.085	2.05
5	49	1.120	2.05
6	54	1.319	2.05
7	57	1.419	2.05
8	59	1.517	2.05
9	63	1.702	2.05
10	69	1.958	2.05
11	73	2.160	2.05
12	78	2.437	2.05
13	80	2.019	2.60
14	85	2.265	2.60
15	92	2.606	2.60
16	98	2.876	2.60
17	104	3.244	2.60
18	100	2.470	3.16
19	107	2.778	3.16
20	116	3.204	3.16
21	123	3.541	3.16
22	131	4.002	3.16
23	119	2.814	3.92
24	127	3.178	3.92
25	137	3.680	3.92
26	145	4.075	3.92
27	155	4.614	3.92
28	137	3.217	4.46
29	146	3.636	4.46
30	159	4.212	4.46
31	167	4.666	4.46
32	178	5.285	4.46
33	169	4.158	5.12
34	183	4.819	5.12
35	193	5.341	5.12
36	206	6.052	5.12
37	192	4.703	5.76
38	208	5.454	5.76
39	219	6.046	5.76
40	235	6.853	5.76
41	0	0	0

, where $\mathrm{TL}=41$ is a fictitious operating point that was introduced to model the electric thruster shutdown condition (in fact, $\mathrm{TL}=41$ gives a null value of T, P, and ${\dot{m}}_p^n$).

The throttle table of the electric thruster used in the numerical simulations is also shown in Figure 4, where it can be seen that the maximum value of the thrust magnitude (i.e., $235\mathrm{mN}$) is reached when $\mathrm{TL}=40$ and the input power is about $6.8\mathrm{kW}$ with a propellant mass flow rate of about $5.7\mathrm{mg}/\mathrm{s}$.

In a spacecraft based on a solar electric propulsion system, i.e., one that uses solar panels to produce the electricity needed to power all on-board subsystems (including the thruster), the electric power $P_A$ available to drive the electric thruster can be approximated as the difference between the solar array output power $P_S$ and the power required by the other subsystems $P_L$. The latter is considered in this paper to be a fixed (and assigned) value. Specifically, according to [15] the value of $P_A$ can be written as

$$P_A=\left\{\begin{array}{cc}{P}_S-P_L\hfill & \mathrm{if}{P}_S>P_L\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(14)

so that when the solar array output power is less than the value of $P_L$, the electric thruster is shut down by construction, i.e., the level $\mathrm{TL}=41$ is naturally selected; see Table 1. On the other hand, when $P_A>0$ the possible throttle levels that can be used to propel the spacecraft (at a generic point on the interplanetary trajectory) are those that satisfy the condition $P\le P_A$. A conceptual scheme of the electric power selection during the transfer is illustrated in Figure 5.

In Equation (14), the solar array output power $P_S$ is modeled with a simple function of the solar distance r and the reference power $P_0$ is defined as the solar array output power at the beginning of the life cycle and at a Sun-spacecraft distance $r_{\oplus }=1\mathrm{au}$, i.e.,

$$P_S=P_0{\left(\frac{r_{\oplus }}{r}\right)}^2$$

(15)

In particular, the value of $P_0$ is considered in this paper as assigned. Note that more refined mathematical models are presented in the literature to describe the variation of $P_S$ as a function of both solar distance and time (due to the degradation effect [35]), as discussed in Ref. [15]. In this context, an interesting historical approach is proposed in Refs. [6,36,37,38]. However, the simple model of Equation (15) is consistent with a preliminary trajectory design without introducing additional parameters that should be tailored to the specific power subsystem. Note that in this model the solar panels are always perpendicular to the incident sunlight.

The value of the propulsive acceleration vector ${\mathit{a}}_p$ and the (effective) propellant mass flow rate ${\dot{m}}_p$ are obtained, as a function of TL, from the pair $\{T,{\dot{m}}_p^n\}$ by introducing the duty cycle $\eta$ defined, according to [39], as the fraction of time during the deterministic thrust periods when the primary propulsive system is on. When using the duty cycle $\eta$ and the pair $\{T,{\dot{m}}_p^n\}$, we obtain

$${\dot{m}}_p=\eta {\dot{m}}_p^n,{\mathit{a}}_p=\eta \frac{T}{m}{\widehat{\mathit{a}}}_p$$

(16)

where ${\widehat{\mathit{a}}}_p$ is the unit vector of propulsive acceleration indicating the direction of the thrust vector induced by the electric thruster. The components in the rotating RTN frame of the unit vector ${\widehat{\mathit{a}}}_p$ can be expressed as a function of the two control angles $\alpha \in [0,\pi ]\mathrm{rad}$ and $\delta \in [0,2\pi ]\mathrm{rad}$ defined in Figure 6 as

$${\widehat{\mathit{a}}}_p=cos\alpha {\widehat{\mathit{i}}}_{\mathrm{R}}+sin\alpha cos\delta {\widehat{\mathit{i}}}_{\mathrm{T}}+sin\alpha sin\delta {\widehat{\mathit{i}}}_{\mathrm{N}}$$

(17)

When taking into account the data summarized in Table 1 and Equation (16), the propulsive terms $\{{\mathit{a}}_p,{\dot{m}}_p\}$ in the equations of motion (7) depend on the three dimensionless control variables $\{\mathrm{TL},\alpha ,\delta \}$. Accordingly, the throttle level and the two thrust angles are considered in this paper as the control parameters that can be used to design the interplanetary transfer trajectory. In particular, the control law, that is, the time variation of the triplet $\{\mathrm{TL},\alpha ,\delta \}$ is determined by solving the optimization problem described in the next subsection.

2.3. Optimal Trajectory Design

The spacecraft transfer trajectory and the control law are obtained by studying an optimization problem in which the flight time $\Delta t=t_f-t_0\equiv t_f$ is minimized. The spacecraft trajectory is determined so as to maximize the performance index $J=-t_f$, where $t_f$ is the final time when the spacecraft’s heliocentric orbit coincides with that of asteroid 4660 Nereus. This implies that the final value of the spacecraft true anomaly is not constrained. In fact, the final condition to be met is to obtain a final orbit whose orbital elements $\{a,e,i,\omega ,\Omega \}$ coincide with those of the target. By using again the JPL Horizons data at the 1 August 2023 we get

$$\begin{array}{c}a\left(t_f\right)=2.22181926\times {10}^8\mathrm{km},e\left(t_f\right)=3.58678173\times {10}^{-1},i\left(t_f\right)=1.45336936\mathrm{deg}\hfill \\ \hfill \omega \left(t_f\right)=1.59511676\times {10}^2\mathrm{deg},\Omega \left(t_f\right)=3.13127787\times {10}^2\mathrm{deg}\end{array}$$

(18)

which can be rewritten in terms of MEOEs as

$$\begin{array}{c}p\left(t_f\right)=1.93598214\times {10}^8\mathrm{km},f\left(t_f\right)=-1.38066389\times {10}^{-1},g\left(t_f\right)=3.31040336\times {10}^{-1}\hfill \\ h\left(t_f\right)=8.67094405\times {10}^{-3},k\left(t_f\right)=-9.25697001\times {10}^{-3}\hfill \end{array}$$

(19)

The flight time is minimized using an indirect approach [40,41], in which the costate vector $\mathit{\lambda }\in {\mathbb{R}}^{(7\times 1)}$ is defined as

$$\mathit{\lambda }={\left[{\lambda }_p,{\lambda }_f,{\lambda }_g,{\lambda }_h,{\lambda }_k,{\lambda }_L,{\lambda }_m\right]}^{\mathrm{T}}$$

(20)

and the Hamiltonian function $\mathcal{H}$ is obtained from Equations (7) and (20) as

$$\mathcal{H}=\left(\mathbb{A}{\mathit{a}}_p\right)·\mathit{\lambda }+\mathit{b}·\mathit{\lambda }$$

(21)

from which the Euler-Lagrange equations are derived, viz.

$$\dot{\mathit{\lambda }}=-\frac{\partial \mathcal{H}}{\partial \mathit{x}}$$

(22)

The explicit expressions of the seven nonlinear differential equations that can be obtained from Equation (22) are omitted here for brevity.

The analytical form of the Hamiltonian function is used to find the optimal control law $\mathrm{TL}=\mathrm{TL}\left(t\right)$, $\alpha =\alpha \left(t\right)$, and $\delta =\delta \left(t\right)$. Indeed, according to the Pontryagin’s maximum principle [42], the optimal control law is the one that maximizes, at each instant of time t, the Hamiltonian function given by Equation (21). A standard procedure was employed to find the triplet $\{\mathrm{TL},\alpha ,\delta \}$ that maximizes the local value of $\mathcal{H}$. The Hamiltonian function is also useful for writing the transversality condition that gives the following 4 boundary constraints

$${\lambda }_L\left(t_0\right)=0,{\lambda }_L\left(t_f\right),{\lambda }_m\left(t_f\right)=0,\mathcal{H}\left(t_f\right)=1$$

(23)

The TPBVP associated with the optimization process at an unspecified terminal time [40] is given by the 14 scalar differential Equations (7) and (22), with the 15 boundary constraints (recall that the flight time is an output of the process) given by Equations (11), (12), (18), and (23). The solution of the TPBVP, obtained by using the numerical approach used by the authors in [43,44], provides the optimal spacecraft trajectory and mission performance as a function of the duty cycle $\eta$ and the spacecraft design characteristics, i.e., initial mass $m_0$, load power $P_L$, and reference power $P_0$. The numerical method used to solve the associated TPBVP is based on the use of a multiple shooting method [45] and a continuation procedure to extend the solution with for low computational cost. The numerical integration of the differential equations was performed using a standard numerical procedure based on the Adams–Bashforth method [46], with a relative and absolute tolerance of ${10}^{-10}$, whereas the TPBVP solution was obtained with an absolute tolerance of ${10}^{-7}$. The results of the numerical simulations are discussed in the next section.

3. Simulations and Numerical Results

Based on the mathematical model described in the previous section, the generic rapid transfer trajectory depends on the value of the set of design parameters $\{\eta ,m_0,P_0,P_L\}$. This section presents the numerical results of a series of simulations designed to assess the impact of the value of the design parameters on mission performance in terms of minimum flight time $\Delta t$ and required propellant mass $m_p\triangleq m_0-m\left(t_f\right)$.

Inder order to reduce the complexity of the numerical problem, a baseline scenario was assumed, in which the spacecraft characteristics were modeled on those of NASA’s DART spacecraft [21], with the exception of $P_L$. More precisely, by taking into account the data presented in [47], a type of “reference spacecraft” was considered with an initial mass of $m_0=610\mathrm{kg}$, a reference power of $P_0=6.6\mathrm{kW}$, a load power of $P_L=500\mathrm{W}$, and a duty cycle of $\eta =0.9$. In this scenario, the solution of the optimization problem gives a minimum flight time of $\Delta t\simeq 330\mathrm{days}$ with a propellant mass consumption of $m_p\simeq 106\mathrm{kg}$. For comparative purposes, the DART spacecraft carried on board about $50\mathrm{kg}$ of hydrazine propellant onboard for deep space maneuvers (and attitude control) and $60\mathrm{kg}$ of xenon for the NEXT-C engine. As such, the total amount of propellant required to reach asteroid 4660 Nereus in this scenario, using a minimum time (orbit-to-orbit) transfer trajectory, is consistent with the total amount of propellant loaded on board the DART spacecraft. In this scenario, an optimal bi-impulsive orbit-to-orbit transfer trajectory, without ephemeris constraints, requires a minimum total velocity change of about $4.96\mathrm{km}/\mathrm{s}$ with a corresponding flight time of $471\mathrm{days}$. On the other hand, according to the CNEOS website (See the url retrieved in 23 October 2023: https://cneos.jpl.nasa.gov) a case using ephemeris constraints and a departure from a circular LEO with an altitude of $400\mathrm{km}$ requires a velocity change of about $5.67\mathrm{km}/\mathrm{s}$.

The solution of the optimal transfer problem gives the temporal variation of the spacecraft’s mass illustrated in Figure 7. Since the curve has a negative slope everywhere, it is evident that the propulsion system always remains on during the entire transfer. In other words, the optimal transfer trajectory contains no coasting arcs, that is, the (virtual) throttle level $\mathrm{TL}=41$ is never engaged.

This aspect is confirmed by the curves in Figure 8, which show the time variation of the optimal control law in terms of throttle level TL and thrust angles $\{\alpha ,\delta \}$. Note, in the upper part of Figure 8, how the propulsion system engages the $\mathrm{TL}=39$ in the first phase of the transfer and then gradually reduces the throttle level according to the available electric power. The latter, in fact, decreases during the transfer because the solar array output power decreases with time due to the increase in Sun-spacecraft distance; see the spacecraft trajectory plotted in Figure 9. The initial phase with a constant value of TL is consistent with the time variation for the spacecraft’s mass shown in Figure 7, where the constant slope of the curve $m=m\left(t\right)$ during the first $70\mathrm{days}$ of the transfer can be observed.

The optimal orbit-to-orbit transfer in the baseline scenario begins (or ends) when the spacecraft’s true anomaly is $\nu \left(t_0\right)=285\mathrm{deg}$ (or $\nu \left(t_f\right)=152\mathrm{deg}$), as confirmed by Figure 9, which shows the spacecraft transfer trajectory in a Cartesian heliocentric-ecliptic reference frame. In particular, Figure 9 shows the spacecraft trajectory (black line), the heliocentric orbit of the Earth (blue line) and asteroid (red line), the points of departure (black circle) and arrival (black square), the Sun (orange circle), and the perihelion points of the planetary orbits (blue and red stars). In this figure, the scale of the z-axis of the side and isometric views is expanded to better appreciate the three-dimensionality of the interplanetary transfer (recall that the asteroid orbit inclination is about $1.4\mathrm{deg}$). Finally, the time variation of the classical orbital elements $\{a,e,i,\omega ,\Omega \}$ of the spacecraft’s osculating orbit are reported in Figure 10.

Transfer Performance Sensitivity

Transfer performance in the baseline mission scenario was used to evaluate the sensitivity of minimum flight time and propellant mass to a change in the design parameters $\{P_L,P_0,m_0\}$. Specifically, each of the three parameters $\{P_L,P_0,m_0\}$ is designed to vary over an appropriate range in order to quantify the corresponding change in $\Delta t$ and $m_p$.

For example, consider the case of $P_L$ ranging in the interval $[100,1000]\mathrm{W}$. When the values of the other design parameters $\{\eta ,m_0,P_0\}$ are the same as those of the baseline scenario, the optimization process gives the results summarized in Figure 11. According to this figure, both the flight time and the propellant mass increase with $P_L$. This is a rather expected result because a higher value of $P_L$ ($P_0$ being fixed) gives a lower available power, which, in turn, changes the set of possible throttle levels in the engine performance of Table 1.

A thorough analysis of the effects of a variation in $P_L$ on mission performance is illustrated in Figure 12.

Similar trends are observed by varying the spacecraft’s initial mass $m_0$ or the reference electric power $P_0$ while keeping the other design parameters fixed (with a value equal to that used in the baseline scenario); see Figure 13. For example, a spacecraft with an initial mass of $1000\mathrm{kg}$, $P_0=6.6\mathrm{kW}$, and $P_L=500\mathrm{W}$ is capable of completing the interplanetary transfer in less than $500\mathrm{days}$ with a propellant expenditure of about $200\mathrm{kg}$. On the other hand, a spacecraft with an initial mass of about twice the value used in the baseline scenario requires a flight time of about $560\mathrm{days}$, which is nearly $1.7$ times the value obtained in the baseline scenario.

4. Final Remarks and Conclusions

Optimal trajectories to asteroid 4660 Nereus have been analyzed using data from a commercial-type electric propulsion system. By using a realistic dataset, namely the actual NEXT-C ion engine throttle table, it is possible to obtain the control law that minimizes the time of flight for an assigned value of the spacecraft’s initial mass and the electric power which are required at the beginning of the transfer. The simulations show that the transfer is feasible in a reasonable flight time, that is, less than one year and with a propellant mass consumption just over $100\mathrm{kg}$ under the assumption of a load power of $500\mathrm{W}$. This result is compatible with the propellant loaded on board the DART spacecraft, which used the same propellant system. The simulations are complemented by a thorough sensitivity study, in which the transfer performance was parameterized as a function of solar array output power, load power, and initial mass. The obtained results show that it is possible to have an accurate estimate of transfer times using a realistic thruster model.