Low-Thrust Trajectory Optimization for Earth-Mars-Jupiter Gravity-Assist Spacecraft Mission

1. Introduction

Space missions face significant challenges, including enhancing payloads, cost reduction, and mitigating uncertainties [1]. To overcome these challenges, it is necessary to use modern space propulsion systems and innovative optimization methods. The feasibility, cost, and uncertainties of an interplanetary mission are significantly dependent on the trajectory optimization design [2,3,4,5,6,7,8].

Recent developments in low-thrust propulsion systems have introduced operational interplanetary missions, exemplified by missions like Deep Space-1 [9], Hayabusa [10,11], Dawn [12], and SMART-1 [13]. The numerical optimization algorithms are widely adopted in low-thrust trajectory optimization due to the lack of analytical solutions [14,15,16,17]. Numerical methods are basically categorized into indirect [18,19,20,21,22] and direct methods [23,24,25], although some are implemented by hybrid solutions [22,23,24,25,26,27,28].

To address challenges inherent in indirect methods, the homotopy/homotopic method introduced by Bertrand and Epenoy aimed to mitigate initial guess sensitivity issues [29]. The homotopy method involves a continuous process of random variable searching, which can reduce efficiency. Some alternative methods tackle this challenge by initializing co-states through particle swarm optimization (PSO) [21,30], optimal auxiliary unconstrained solutions [31], and the adjoint variable method [32]. Subsequent refinements by Jiang et al. [30] and Saghamanesh and Baoyin [21,33,34,35,36] have further improved this approach using various effective-practical techniques.

Low-thrust spacecraft propulsion systems effectively integrate with planetary gravity-assist (GA) maneuvers in order to reduce fuel-consumption. Notable missions utilizing GA include Cassini-Huygens [37], Messenger [38,39], Rosetta [40], New Horizons [41], and Jupiter's Galileo [42]. GA allows spacecraft to capture a target planet with reduced fuel consumption.

Most GA research traditionally relies on the patched-conics method [43] within the heliocentric two-body (2B) model, approximating the actual ephemeris model [33,34,35]. However, limitations arise in performing trajectory constraints in practical space environments. Recent studies explore alternatives like focusing on energy considerations [44], combining GA and impulsive scenarios [45], implementing two-dimensional swing-bys [46], and applying the patched-conics 2B model.

This paper introduces a practical methodology for deriving the state/co-state dynamics in solving the GA fuel-optimal problem based on the 2B model in Cartesian coordinates. Introducing this GA fuel-optimal problem adds further nonlinearity to the differential equations of motion. Firstly, solving the continuous GA-MPBVPs (Multiple Point Boundary Value Problem), specifically in Cartesian coordinates, becomes notably more complex than in other coordinate systems. Secondly, this paper devises an effective set of shooting function equation constraints, encompassing trajectory, GA inner, and GA intermediate constraints. These effectively manage the highly nonlinear and sensitive state/co-state dynamics inherent in continuous GA-MPBVPs, ensuring accurate control over the fuel-optimal GA computational process. Thirdly, a robust homotopic approach, serving as a hybrid solution, is systematically refined to address the challenges of solving the GA-MPBVPs. The proposed algorithm effectively addresses convergence challenges in solving GA-MPBVPs by integrating homotopy methods with normalized co-state vectors. Fourthly, a modified switching detection algorithm is introduced, comprising five iterative steps and applied in a practical application of the methodology. Fifthly, A novel parameterization of the search space is introduced, providing a comprehensive set of relations. Additionally, the methodology incorporates efficient techniques to reduce computational complexity, ensuring a robust understanding of unknowns and accurate satisfaction of penalty conditions. The objective functions in GA energy-to-fuel optimal problems are formed using new constraints as penalty conditions and are satisfied with high accuracy.

The proposed algorithm is applied to solve fuel-optimal low-thrust missions, including scenarios with and without intermediate Mars-GA maneuvers for Jupiter rendezvous. Challenges arise in the E-M-J scenario due to vast distances and extended flight times, especially within the GA maneuver dynamics model. Comparison of the solution with the GALLOP space mission [48] shows the improvement of the spacecraft mass and the overall performance of the solution.

2. System dynamics model

We focus on the dynamic system design, encompassing a low-thrust interplanetary mission from Earth to Jupiter, considering both direct transfer and the Mars gravity-assist scenario. The spacecraft encounters the gravitational force of the Sun, while its low-thrust propulsion system induces acceleration. The state dynamics system is modeled as follows:

$$\dot{\mathit{x}}=\{\begin{array}{l}\dot{\mathit{r}}=\mathit{v},\hfill \\ \dot{\mathit{v}}=-\frac{\mu }{r^3}\mathit{r}+\frac{T_{\max }u}{m}\mathit{\alpha },\hfill \\ \dot{m}=-\frac{T_{\max }}{I_{sp}{g}_0}u,\hfill \end{array}$$

(1)

here, $\mathbf{r}$ and $\mathbf{v}$ denote position and velocity vectors of the spacecraft. The spacecraft's initial and final conditions are given as $\mathsf{\psi }\left({\mathrm{t}}_0\right)={\left[{\mathbf{r}}_0{\mathbf{v}}_0{\mathrm{m}}_0\right]}^{\mathrm{T}}$and $\mathsf{\psi }\left({\mathrm{t}}_{\mathrm{f}}\right)={\left[{\mathbf{r}}_{\mathbf{f}}{\mathbf{v}}_{\mathbf{f}}\right]}^{\mathrm{T}}$, respectively, with the final mass ${\mathrm{m}}_{\mathrm{f}}$ left as a free parameter. Also, ${\mathrm{T}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and $\mathrm{u}\in [0,1]$ signify the maximum thrust and engine throttle magnitudes. $\mathsf{\alpha }$ represents the unit vector in the control direction, shaping the thrust vector as $\mathbf{T}=\mathrm{u}{\mathrm{T}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\mathsf{\alpha }$. $\mathrm{m}$ refers to the spacecraft’s instantaneous mass, while ${\mathrm{I}}_{\mathrm{s}\mathrm{p}}$ stands for the thruster's specific impulse.

It is assumed that the upper stage of the launcher includes a Heavy Lifter with a chemical propulsion system. The Heavy Lifter is able to inject the spacecraft with a maximum initial mass of $m_0$ (kg) to the interplanetary parabolic trajectory, while the launch hyperbolic excess velocity is ${v_{\mathrm{\infty }}}_0=0$. According to Figure 1, while ${v_{\mathrm{\infty }}}_0>0$, the spacecraft initial mass is given by Tsiolkovsky’s rocket equation as follows:

$$m_0(v_{{\infty }_0})=m_0\exp (-\frac{\Delta v}{v_{exhaust}}),$$

(2)

$$\begin{array}{l}{\mathit{v}}_{{\infty }_0}=v_{{\infty }_0}{\mathit{u}}_{{\infty }_0},{\mathit{u}}_{{\infty }_0}=\pm \frac{{\mathit{v}}_{\mathrm{E}}(t_0)}{\Vert {\mathit{v}}_{\mathrm{E}}(t_0)\Vert },\\ \Delta v=v_{\mathrm{LEO}}\left(\sqrt{2+{(\frac{v_{{\infty }_0}}{v_{\mathrm{LEO}}})}^2}-1\right),\end{array}$$

(3)

here, ${\mathit{v}}_{\mathrm{E}}$ represents the Earth's velocity vector. In this scenario, the objective here is to minimize fuel consumption. It's noted that the spacecraft's position vector coincides with the planet's position vector at the instantaneous GA time, $t_G$. Additionally, the time spent within the GA mission is disregarded, simplifying to $t_G^{-}=t_G^{+}=t_G$, and denoting that $\mathit{r}\left(t_G^{-}\right)=\mathit{r}\left(t_G^{+}\right)={\mathit{r}}_{\mathit{p}\mathit{l}}\left(t_G\right)$, where – and + signify the time immediately before and after the GA scenario in HERF.

Based on Figure 2, the spacecraft's velocity vectors are expressed as $\mathit{v}\left(t_G^{-}\right)={\mathit{v}}_{pl}\left(t_G\right)+{\mathit{v}}_{\infty }^{-}$ and $\mathit{v}\left(t_G^{+}\right)={\mathit{v}}_{pl}\left(t_G\right)+{\mathit{v}}_{\infty }^{+}$, where $‖{\mathit{v}}_{\infty }^{-}‖=‖{\mathit{v}}_{\infty }^{+}‖$. Here, ${\mathit{v}}_{\mathrm{\infty }}^{\pm }$ and ${\mathit{v}}_{pl}\left(t_G\right)$ represent GA hyperbolic excess velocity vectors and the velocity vector of the GA planet, respectively. Meanwhile, $\delta$ denotes the GA turn angle, and the hyperbolic periapsis radius of the GA orbit is parameterized as:

$$\delta =\arccos \left(\frac{{\mathit{v}}_{\infty }^{-}.v_{\infty }^{+}}{{\mathit{v}}_{\infty }^2}\right),\delta \in [0,{\delta }_{\max }],$$

(4)

$$r_p=\frac{{\mu }_{pl}}{v_{\infty }^{-}{v}_{\infty }^{+}}\left(\frac{1}{\sin (\delta /2)}-1\right),$$

(5)

The minimum acceptable GA radius ($r_{min})$ limits the maximum turn angle ${(\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}})$ as follows:

$${\delta }_{\max }=2\arcsin \left(\frac{1}{1+\frac{r_{\min }}{{\mu }_{pl}}{v}_{\infty }^2}\right),$$

(6)

Hence, ${\mathit{v}}_{\infty }^{-}and{\mathit{v}}_{\infty }^{+}$ possess identical magnitudes but differ only in the angle $\delta$. Also, the gravity-assist framework can be written as:

$$\begin{array}{l}\widehat{i}=\frac{{\mathit{v}}_{\infty }^{-}}{v_{\infty }},\\ \widehat{k}=\frac{{\mathit{v}}_{pl}(t_G)\times {\mathit{v}}_{\infty }^{-}}{\Vert {\mathit{v}}_{pl}(t_G)\times {\mathit{v}}_{\infty }^{-}\Vert },\\ \widehat{j}=\widehat{k}\times \widehat{i},\end{array}$$

(7)

In Figure 3, $\alpha$ and $\delta$ are utilized to parameterize the velocity vector. Following this parameterization, the GA impulse is obtained from the velocity increment in the following manner:

$$\Delta {\mathit{v}}_G=v(t_G^{+})-\mathit{v}(t_G^{-})={\mathit{v}}_{\infty }(\cos \delta \widehat{i}+\sin \delta \sin \alpha \widehat{j}+\sin \delta \cos \alpha \widehat{k}),$$

(8)

3. Trajectory optimization

The performance index aimed at minimizing propellant consumption is expressed as follows:

$$J={\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{T_{\max }}{I_{SP}{g}_0}udt},$$

(9)

The homotopy performance index is employed to address the challenge of the optimal control problem [21,30]:

$$J_a={\lambda }_0{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{T_{\max }}{I_{SP}{g}_0}[u-\epsilon u(1-u)]dt},$$

(10)

In this context, $\epsilon =1$ corresponds to the energy-optimal problem, whereas $\epsilon =0$ corresponds to the fuel-optimal problem. The constant auxiliary multiplier ${\lambda }_0$ belongs to the interval $\left(\mathrm{0,1}\right]$. Since Pontryagin’s Maximum Principle (PMP) [49],${\lambda }_0$, and the co-state vector${\mathit{\lambda }=\left[{\lambda }_0,{\mathit{\lambda }}_{\mathit{r}}^T,{\mathit{\lambda }}_{\mathit{v}}^T,{\lambda }_m\right]}^T$ are used to build the Hamiltonian as follows:

$$H={\mathit{\lambda }}_r^T.\mathit{v}+{\mathit{\lambda }}_{\mathit{v}}^T.\left(-\frac{\mu }{r^3}\mathit{r}+\frac{T_{\max }u}{m}\mathit{\alpha }\right)-{\lambda }_m\frac{T_{\max }}{I_{SP}{g}_0}u+{\lambda }_0\frac{T_{\max }}{I_{SP}{g}_0}[u-\epsilon u(1-u)],$$

(11)

The Primer Vector [50], $\mathit{\alpha }=-{\mathit{\lambda }}_{\mathit{v}}/‖{\mathit{\lambda }}_{\mathit{v}}‖$, is used to determine the optimal thrust direction for minimizing the Hamiltonian. In this way, the Euler–Lagrange conditions,$\dot{\mathit{\lambda }}=-{\left[\partial H/\partial \mathit{x}\right]}^T$, and state dynamics system (1) are used to establish the ordinary differential equations (ODEs):

$$\dot{\mathit{x}}=\{\begin{array}{l}\dot{\mathit{r}}=\mathit{v},\hfill \\ \dot{\mathit{v}}=-\frac{\mu }{r^3}\mathit{r}+\frac{T_{\max }u}{m}\mathit{\alpha },\hfill \\ \dot{m}=-\frac{T_{\max }}{I_{sp}{g}_0}u,\hfill \\ {\dot{\mathit{\lambda }}}_r=\frac{\mu }{r^3}{\mathit{\lambda }}_{\mathit{v}}-\frac{3\mu {\mathit{\lambda }}_{\mathit{v}}.r}{r^5}\mathit{r},\hfill \\ {\dot{\mathit{\lambda }}}_v=-{\mathit{\lambda }}_r,\hfill \\ {\dot{\lambda }}_m=-\frac{T_{\max }u}{m^2}\Vert {\mathit{\lambda }}_{\mathit{v}}\Vert ,\hfill \end{array}$$

(12)

Where the switching function (SF) and the extremal control magnitude are formed by $\rho$, which is also referred to as the switching function:

$$\rho =\left(1-\frac{\Vert {\mathit{\lambda }}_{\mathit{v}}\Vert I_{sp}{g}_0}{{\lambda }_{0}m}-\frac{{\lambda }_m}{{\lambda }_0}\right),$$

(13)

$$u=\{\begin{array}{l}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\rho >\epsilon ,\hfill \\ 1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\rho <-\epsilon ,\hfill \\ u=\frac{1}{2}-\frac{\rho }{2\epsilon } \left|\rho \right|\le \epsilon ,\hfill \end{array}$$

(14)

Eqs. (15) should be satisfied as gravity-assist transversal conditions:

$$\{\begin{array}{l}{\mathit{\lambda }}_r\left(t_G^{+}\right)={\mathit{\lambda }}_r\left(t_G^{-}\right)-Z_{1\sim 3},\\ {\mathit{\lambda }}_v\left(t_G^{+}\right)=Z_4{\mathit{u}}^{+}+\frac{1}{r_{\min }}\kappa \mathit{B},\\ {\mathit{\lambda }}_{\mathit{v}}\left(t_G^{-}\right)=Z_4{\mathit{u}}^{-}-\frac{1}{r_{\min }}\kappa \mathit{A},\end{array}$$

(15)

where

$$\begin{array}{l}{\mathit{u}}^{\pm }=\frac{{\mathit{v}}_{\infty }^{\pm }}{v_{\infty }},\\ \mathit{A}=\frac{r_p}{v_{\infty }^{-}}\left[\frac{1}{4{\sin }^2\delta /2\left(1-\sin \delta /2\right)}\left({\mathit{u}}^{+}-\cos \delta {\mathit{u}}^{-}\right)-{\mathit{u}}^{-}\right],\\ \mathit{B}=\frac{r_p}{v_{\infty }^{+}}\left[\frac{1}{4{\sin }^2\delta /2\left(1-\sin \delta /2\right)}\left({\mathit{u}}^{-}-\cos \delta {\mathit{u}}^{+}\right)-{\mathit{u}}^{+}\right],\end{array}$$

(16)

The velocity and co-state vectors are different at the GA maneuver times, $t_G^{-}$and$t_G^{+}$. Therefore, Eqs. (17) must be satisfied as the stationarity conditions of gravity-assist:

$$H\left(t_G^{-}\right)-H\left(t_G^{+}\right)-Z_{1~3}\cdot {\mathit{v}}_{pl}\left(t_G\right)+Z_4\left({\mathit{u}}^{+}-{\mathit{u}}^{-}\right)\cdot {\mathit{a}}_{pl}\left(t_G\right)-\frac{1}{r_{\min }}\kappa \mathit{C}=0$$

(17)

where ${\mathit{a}}_{pl}\left(t_G\right)$ denotes the acceleration of the GA planet at the $t_G$ and $C$ are written as:

$$\begin{array}{l}{\mathit{a}}_{pl}\left(t_G\right)=-\frac{\mu }{{\mathit{r}}_{pl}^3(t_G)}{r}_{pl}(t_G)\\ \mathit{C}=r_p\{-\frac{1}{4{\sin }^2\delta /2\left(1-\sin \delta /2\right)}[\frac{1}{v_{\infty }^{+}}\left({\mathit{u}}^{-}-\cos \delta {\mathit{u}}^{+}\right)+\frac{1}{v_{\infty }^{-}}\left({\mathit{u}}^{+}-\cos \delta {\mathit{u}}^{-}\right)]\}.{\mathit{a}}_{pl}\left(t_G\right)\end{array}$$

(18)

3.1. Constraints introduced by the GA scenario

We consider the fuel-optimal bang-bang control trajectory using a Mars-GA to showcase the distinguishing attributes of the algorithm. Predetermined values include the initial time$(t_0$), GA time ($t_G$), and terminal time ($t_f$). Consequently, the conditions governing the optimal trajectory design are divided into three parts: i) initial transfer conditions are set to match Earth’s state vectors ${\mathit{r}}_E$ and${\mathit{v}}_E$ ii) gravity-assist maneuver conditions are equal to the planet’s state vectors ${\mathit{r}}_{pl}$ and${\mathit{v}}_{pl}$ iii) the spacecraft's terminal state vectors are equal to the target body’s state vectors ${\mathit{r}}_f$ and${\mathit{v}}_f$. Additionally, the spacecraft’s final mass is treated as a free parameter. Moreover, the terminal transversality condition is set to zero:

$${\lambda }_m(t_f,\epsilon )=0,$$

(19)

Hence, according to Eqs. (2-3), we update the boundary conditions as follows:

$${\mathit{\psi }}_0=\left[\begin{array}{c}\mathit{r}(t_0)-{\mathit{r}}_E\\ \begin{array}{l}\mathit{v}(t_0)-{\mathit{v}}_E-v_{{\infty }_0}{\mathit{u}}_{{\infty }_0}\\ m(t_0)-m_0(v_{{\infty }_0})\end{array}\end{array}\right]=0,$$

(20)

$${\mathit{\psi }}_f=\left[\begin{array}{c}\mathit{r}(t_f)-{\mathit{r}}_f\\ \begin{array}{c}\mathit{v}(t_f)-{\mathit{v}}_f\\ {\lambda }_m(t_f)\end{array}\end{array}\right]=\mathbf{0},m(t_f)=m_f,$$

(21)

The normalization technique regarding references [21,30] is used as:

$$\mathit{\lambda }(t,\epsilon )=\frac{\mathit{\lambda }(t,\epsilon )}{\Vert \mathit{\lambda }(t_0,\epsilon )\Vert };\mathit{\lambda }(t_0,\epsilon )=\left[{\lambda }_0{\mathit{\lambda }}_r(t_0){\mathit{\lambda }}_{\mathit{v}}(t_0){\lambda }_m\right],$$

(22)

Moreover, Eqs. (23) define the intermediate gravity-assist constraints:

$$\mathit{\psi }(t_G)=\left[\begin{array}{c}\mathit{r}(t_G)-{\mathit{r}}_{pl}(t_G)\\ \Vert {\mathit{v}}_{\infty }^{-}\Vert -\Vert {\mathit{v}}_{\infty }^{+}\Vert \end{array}\right]=\mathbf{0},$$

23

$${\sigma }_G=1-\frac{r_p}{r_{\min }}\le 0,$$

24

Besides, rigid conditions need to be accounted for due to the inequality constraint, Eq. (24), as follows:

$$\kappa .{\sigma }_G=0,\kappa \ge 0,$$

(25)

Ultimately, Eq. (26) satisfies the normalization condition for the gravity-assist:

$$\Vert \left[\mathit{\lambda }(t_0,\epsilon );{\mathrm{Z}}_{1~4};\kappa \right]\Vert =1,$$

(26)

4. Summary of the optimization algorithm

In this discussion, let's summarize the robust dual-step homotopic algorithm. The trajectory design problem undergoes transformation into a multi-point boundary value problem (MPBVP) by integrating various elements including the initial co-state vectors from Eq. (22), the boundary conditions outlined in Eqs. (20-21), the intermediate GA constraints represented by Eq. (23), the GA transversal condition ${\mathit{\lambda }}_{\mathit{v}}\left(t_G^{-}\right)$ derived from Eqs. (15), the GA stationarity condition from Eq. (17), the rigid conditions indicated in Eq. (25), and the GA normalization condition expressed in Eq. (26). Ultimately, achieving satisfaction of the shooting function equation $\mathsf{\Phi })$ should be satisfied as:

$$\Phi \left(\mathit{\lambda }(t_0,\epsilon )\right)=\left[\begin{array}{c}\begin{array}{c}\mathit{r}(t_f)-{\mathit{r}}_f\\ \mathit{v}(t_f)-{\mathit{v}}_f\\ {\lambda }_m\left(t_f\right)\\ \mathit{r}(t_G)-{\mathit{r}}_{pl}(t_G)\\ v_{\infty }^{-}-v_{\infty }^{+}\end{array}\\ \begin{array}{c}{\mathit{\lambda }}_v\left(t_G^{-}\right)-{\mathrm{Z}}_4{\mathit{u}}^{-}+\frac{1}{r_{\min }}\kappa \mathit{A}\\ H\left(t_G^{-}\right)-H\left(t_G^{+}\right)-Z_{1~3}\cdot v_{pl}\left(t_G\right)\mathrm{...}\\ +Z_4\left({\mathit{u}}^{+}-{\mathit{u}}^{-}\right)\cdot {\mathit{a}}_{pl}\left(t_G\right)-\frac{1}{r_{\min }}\kappa \mathit{C}\\ \kappa .{\sigma }_G,\kappa \ge 0\end{array}\\ \\ \Vert \left[\mathit{\lambda }(t_0,\epsilon );{\mathrm{Z}}_{1~4};\kappa \right]\Vert -1\end{array}\right]=0,$$

(27)

Certainly, by employing normalized co-states consistently within the algorithm, the constraint presented in Eq. (26) is satisfied throughout the optimization process. Therefore, it can be disregarded from the shooting function mentioned in Eq. (27).

The algorithm process is outlined in Figure 4, comprising two distinct steps. The initial step encompasses the Energy-Optimal step, divided into two segments: the PSO segment, functioning as a direct method with an objective function represented by Eq. (28), and the homotopy segment, an indirect method with the inclusion of the penalty factor outlined in Eq. (27). Employing the homotopy performance index, $J_a$, the PSO segment aids in establishing an approximate solution utilized within the homotopy segment to resolve the energy-optimal problem.

$$I\left(\mathit{X}(1,k)\right)=\Vert \Phi \left(\mathit{\lambda }(t_0,1)\right)\Vert +J_a,$$

(28)

Since, the 12-D angle variables $\phi$, considering the fixed values of $t_0$,$t_G$, and $t_f$, are constructed utilizing the 12-D search-variables $X\left(\epsilon ,k\right)$ in both steps by

$$\begin{array}{cc}\begin{array}{l}{\phi }_{1,2,3}=\frac{\pi }{2}{X}_{1,2,3}\\ {\phi }_{4,5}=\pi (X_{4,5}-\frac{1}{2})\\ {\phi }_{6,7}=2\pi X_{6,7}\\ {\phi }_8=\pi (X_8-\frac{1}{2})\\ {\phi }_9=\frac{\pi }{2}{X}_9\\ {\phi }_{10}=2\pi X_{10}\\ {\phi }_{11,12}=\frac{\pi }{2}{X}_{11,12}\end{array}& \begin{array}{l}\in \left[0,\frac{\pi }{2}\right],\\ \in \left[-\frac{\pi }{2},\frac{\pi }{2}\right],\\ \in \left[0,2\pi \right],\\ \in \left[-\frac{\pi }{2},\frac{\pi }{2}\right],\\ \in \left[0,\frac{\pi }{2}\right],\\ \in \left[0,2\pi \right],\\ \in \left[0,\frac{\pi }{2}\right],\end{array}\end{array}$$

(29)

So, the 8-D co-states are constrained through a new approach during the computational process, restricted within the domain [-1,1] by

$$\begin{array}{l}{\lambda }_0=\sin {\phi }_1,\\ {\lambda }_m(t_0,\epsilon )=\cos {\phi }_1\sin {\phi }_2\cos {\phi }_{12},\\ {\mathit{\lambda }}_r(t_0,\epsilon )=\cos {\phi }_1\cos {\phi }_2\cos {\phi }_3\cos {\phi }_{12}{[\cos {\phi }_4\cos {\phi }_6\cos {\phi }_4\sin {\phi }_6\sin {\phi }_4]}^T,\\ {\mathit{\lambda }}_v(t_0,\epsilon )=\cos {\phi }_1\cos {\phi }_2\sin {\phi }_3\cos {\phi }_{12}{[\cos {\phi }_5\cos {\phi }_7\cos {\phi }_5\sin {\phi }_7\sin {\phi }_5]}^T,\end{array}$$

(30)

Furthermore, the 5-D numerical multipliers related to the intermediate GA maneuver are normalized within the domain [-1,1] by

$$\begin{array}{l}{Z}_1=\cos {\phi }_1\sin {\phi }_8\sin {\phi }_{12},\\ Z_{2~4}=\cos {\phi }_1\cos {\phi }_8\cos {\phi }_9\sin {\phi }_{12}{[\cos {\phi }_{10}\cos {\phi }_{11}\cos {\phi }_{10}\sin {\phi }_{11}\sin {\phi }_{10}]}^T,\end{array}$$

(31)

$$\kappa =\cos {\phi }_1\cos {\phi }_8\sin {\phi }_9\sin {\phi }_{12},$$

(32)

The residual errors of both the objective function (Eq. 28) and the penalty factor (Eq. 27) oversee both segments of the algorithm. Additionally, the Runge-Kutta method is employed as the integrator.

The Fuel-Optimal step, the second step, determines a series solution from the energy-optimal to the fuel-optimal problem using the homotopy algorithm. It is important to note that the optimal search-variables${\mathit{X}}^{\ast }\left(1,k\right)$ obtained in the first step are used as the starting point.

Throughout the algorithm, the co-states are normalized throughout the algorithm within the domain$\left[-\mathrm{1,1}\right]$. The Runge-Kutta method is employed as the integrator, and ${\lambda }_0$ maintains a constant value from the initial step. The residual error of the penalty factor Eq. (27) oversees this step. Additionally, it is crucial to note specific details concerning the Intermediate GA maneuver:

1)

It involves eight unknowns, including the 5-D numerical multipliers, $\left[Z_{1~4}\kappa \right]$ , and the 3-D GA velocity vector, $\mathit{v}\left(t_G^{+}\right)={\mathit{v}}_{pl}\left(t_G\right)+{\mathit{v}}_{\infty }^{+}$. The GA date (t_G) remains fixed throughout the solution.

2)

Additionally, nine identical equations must be satisfied: the 4-D intermediate GA constraints in Eq. (23), the 1-D rigid conditions in Eq. (25) utilizing the inequality condition from Eq. (24), the 3-D GA transversal condition, ${\mathit{\lambda }}_{\mathit{v}}\left(t_G^{-}\right)$, derived from Eqs. (15), and the 1-D stationarity condition specified in Eq. (17)

3)

According to Figure 3, the gravity-assist generates the velocity increment$\mathsf{\Delta }{\stackrel{⃑}{v}}_G$ Eq. (8), instead of the GA velocity vector $\mathit{v}\left(t_G^{+}\right)$. It is essential to note that in Eq. (15), ${\mathit{\lambda }}_{\mathit{r}}\left(t_G^{+}\right)$ and ${\mathit{\lambda }}_{\mathit{v}}\left(t_G^{+}\right)$ are exclusively utilized to revise the position and velocity co-state vectors at the date $t_G^{+}$ immediately after the GA, and they are omitted from the penalty factor Eq. (27).

4)

Furthermore, the mass, its corresponding co-state, and the position vector remain continuous immediately at $t_G^{-}$and$t_G^{+}$ during the GA date, and hence are fixed within the GA scenario.

5. Switching function detection

This section presents a robust solution to the intrinsic complexities faced in optimal control problems. Although the magnitude of thrust is modelled as continuous, the switching function is susceptible to obtaining critical values at specific intermediate and terminal points. Nevertheless, by adopting a step size $\text{'}h\text{'}$ of suitable smallness, continuity of the switching function’s execution across the designated intervals in Eq. (14) is maintained. Innovatively, we constrain the admissible trajectories for the detection of the switching function—bridging energy and fuel optimization—within a tighter domain. This is achieved through the systematic integration of a normalized co-state approach across all stages of the algorithm. Consequently, this refinement ensures the attainment of global convergence with notable accuracy, thereby expediting the computational process towards an optimal solution with greater efficiency.

When employing the normalization technique throughout the computational solution, the subsequent iteration's switching function (SF) value can be approximated as:

$${\rho }_{i+1}={\rho }_i+{\dot{\rho }}_{i}h,$$

(33)

The updated SF detection method includes five steps within the interval$t\in [t_0,t_f]$:

Step (1): Initiated according to Flowchart Figure 5, detects${\rho }_i\left(x\right)$:

Step (2): As long as ${\rho }_{i+1}$ is greater than threshold $\epsilon$, the thrust magnitude is assigned a zero value $(u=0)$. The procedure for SF detection is then outlined in the flowchart depicted in Figure 6:

Step (3): If ${\rho }_{i+1}$ falls below negative $\epsilon$, the thrust magnitude is set to $u=1$. The SF detection procedure for this scenario is shown in Figure 7:

Step (4): When ${\rho }_{i+1}$ lies within the range $\left[-\epsilon ,\epsilon \right]$, the thrust magnitude is set as $u=1/2-\rho /2\epsilon$. The corresponding SF detection procedure is captured in flowchart Figure 8:

Step (5): This is the decision-making step where the normalized thrust magnitude $u_{i+1}=u\left({\mathit{x}}_{i+1}\right)$ is computed. Following this, the detection procedure switches to the next iteration loop.

6. Numerical analysis

The Earth-Mars-Jupiter Gravity Assist (E-M-J) Gravity-Assis scenario, employing a deep-space orbit and nuclear electric propulsion, is analyzed to demonstrate the proposed algorithm's efficiency and automation. The homotopy algorithm is utilized to transition from energy-optimal to fuel-optimal problem-solving sequences. The complexity of achieving a fuel-efficient rendezvous, given the extensive Earth-Jupiter transfer, lengthy flight duration, and intricate GA maneuvers, necessitates robust convergence mechanisms aimed at high-accuracy. The mission implements an optimal escape velocity $v_{\mathrm{\infty }}.$ This mission scenario has been previously computed using the GALLOP program [48], engaging numerous optimization parameters and parallel computational processes. Reference [48] asserts the utilization of a Heavy Lifter, incorporating a chemical propulsion system with a specific impulse $I_{sp}$ of 404 seconds, to propel the spacecraft, weighing an initial mass of 20,000 kg, onto an interplanetary parabolic trajectory with a commencement hyperbolic excess velocity of ${v_{\mathrm{\infty }}}_0=0$.

The fixed departure date $t_0$ (February 17, 2022), the GA date$t_G$ (March 18, 2024), and the capture date $t_f$ (February 01, 2027) are synchronized to the GALLOP program's schedule. The parameters and conditions for orbital transfers are tabulated in Table 1, identical to those specified in the reference material. The mission trajectory is designed such that the spacecraft escapes Earth, utilizes a Mars GA, and captures Jupiter, aligning with the heliocentric position and velocity vectors of Earth, Mars, and Jupiter, respectively. The planets' orbital elements are derived from the JPL/HORIZON1 database.

In line with Section 4, the initial approximation of the solution is obtained using the Particle Swarm Optimization (PSO) algorithm. This estimate then serves as a basis for computing the energy-optimal route via an indirect method. As illustrated in Figure 3 and articulated in Eq. (8), the increment in GA velocity required for the numerical analysis is determined.

$$\begin{array}{l}\Delta {\mathit{v}}_G=X_{15}\sqrt{\frac{{\mu }_{pl}}{r_{\min }}}{\left[\cos \delta \cos \alpha \sin \delta \cos \alpha \sin \alpha \right]}^T,\\ \{\begin{array}{l}\begin{array}{cc}\alpha =\pi \left(X_{13}-0.5\right)& \in \left[-\frac{\pi }{2},\frac{\pi }{2}\right],\end{array}\\ \begin{array}{cc}\delta =2\pi X_{14}& \in \left[0,2\pi \right],\end{array}\end{array}\end{array}$$

(34)

As previously stated, the constraint embodied in Eq. (26) is eliminated from the shooting function delineated in Eq. (27). Taking into account the intermediate GA maneuver, the resulting Multiple Point Boundary Value Problem (MPBVP) encompasses 17 unknowns:

$$XX=\left[v_{{\infty }_0};{\lambda }_0;{\mathit{\lambda }}_r(t_0);{\mathit{\lambda }}_{\mathit{v}}(t_0);{\lambda }_m;Z_{1~4};\kappa ;\Delta {\mathit{v}}_G\right],$$

(35)

The 8-dimensional optimal co-state values are listed in Table 2, elucidating results for both ε=1 (energy-optimal) and ε=0 (fuel-optimal) scenarios. The findings are further detailed in Table 3, which provides a comparison with those of the GALLOP program. According to Table 3, Figure 1, and Eqs. (2-3), the optimal launch hyperbolic excess velocity, denoted by ${v_{\mathrm{\infty }}}_0$, is 0.8 km/s, correlating to an initial mass of $m_0$ of 19,819 kg. Additionally, the gravity-assist's hyperbolic excess velocity, $‖{\mathit{v}}_{\mathrm{\infty }}^{\pm }‖$, measures 3.654 km/s, along with a minimum periapsis radius $r_{min}$ of 500 km. The turn angles $\delta$ for fuel-optimal GA maneuvers are valued at 53.74°, occurring at the periapsis radius equal to $r_p=r_{min}$.

Moreover, the suggested solution demonstrates a 'bang-bang' optimal-control trajectory, which outperforms existing solutions in terms of reduced fuel consumption. It delivers an additional 189.17 kg of mass (a 1% net improvement, albeit modest). The numerical illustrations of the GA scenario are summarized in Table 4, which includes the GA timing, GA hyperbolic excess velocity vectors ${\mathit{v}}_{\mathrm{\infty }}^{\pm }$, the GA impulsive velocity increment ${∆\mathit{v}}_G$, and the 4-dimensional GA numerical multipliers $\left[Z_{1~4}\right]$ and the 1-dimensional GA numerical multiplier $\kappa$, which confined within the domains [-1,1] and [0,1], respectively.

Figure 9 charts the energy-optimal and fuel-optimal solutions throughout the Mars-GA maneuver. This illustrates the SF (Switching Function) detection associated with the bang-bang optimal-control strategy. Specifically, the function $\rho (t,0)$ attains a value of zero at six distinct points within the solution. Moreover, when $\rho (t,0)$ identifies a negative value, the resultant thrust profile escalates to maximum magnitude. Conversely, the thrust magnitude drops to zero whenever $\rho (t,0)$ recognizes a positive value. The employment of sophisticated techniques confines the SF within a narrow region, thereby facilitating a smoother thrust profile.

Figure 10 delineates the fuel-optimal trajectory connecting Earth to Jupiter, incorporating the Mars-GA. Featuring three zero-thrust and three maximum-thrust phases, Figure 10 differentiates the pre- and post-Mars-GA optimal trajectories via solid and dash-dot lines, respectively. The results emphasize remarkable accuracy against the shooting function Equation (27) constraints, evidenced by the negligible final position and velocity errors $∆r_f=7.215\times {10}^{-5}$ km and $∆v_f=2.631\times {10}^{-14}$km/s.

In addition, Figure 11 shows that ${\lambda }_m^{\ast }\left(t_f,0\right)$ achieves zero value to satisfy terminal transversality condition, Eq. (19). Moreover, Fig 11 shows that the variations of the mass and its corresponding co-state are continuous just at $t_G^{-}$and$t_G^{+}$ of the GA date. Therefore, they are not considered for an update at the GA scenario. The corresponding auxiliary multiplier is${\lambda }_0=0.74369984308156$ between the sequence solutions from the energy- to fuel-optimal problems.

Instantaneous mass ratio, $m\left(t\right)/m_0$, is plotted in Figure 11, showcasing a final spacecraft mass that preserves 76.029% of its initial value, thereby indicating a fuel consumption of 23.971%. In the same Figure, the fuel-optimal co-states vector ${{\mathit{\lambda }}^{\mathit{\ast }}(t,0)=\left[{\lambda }_0,{\mathit{\lambda }}_{\mathit{r}}^T,{\mathit{\lambda }}_{\mathit{v}}^T,{\lambda }_m\right]}^T$ remains consistently confined within the $\left[-\mathrm{1,1}\right]$ domain, in sharp contrast to broader solution spaces ranging$\left[-\mathrm{100000,100000}\right]$, where global convergence poses significant computational challenges for interpolating energy- and fuel-optimal solutions. Moreover, Figure 11 illustrates that${\lambda }_m^{\ast }\left(t_f,0\right)$ settles at zero to satisfy the terminal transversality condition, Eq. 19, and it also indicates continuous variations of mass and its associated co-state at $t_G^{-}$and$t_G^{+}$, exempting them from updates within the GA framework. The inter-solution auxiliary multiplier is ${\lambda }_0=0.74369984308156$, bridging energy- and fuel-optimal solution spaces.

7. Conclusions

This study presented a systematic design framework with computational efficiency at its core, leveraging a robust homotopic method in conjunction with an indirect optimization approach. A gradated sequence of solutions, transitioning from energy optimization to fuel optimization, was effectively actualized through a dual-step hybrid optimization algorithm, integrating Particle Swarm Optimization (PSO) with Homotopy. Initially, PSO facilitated the derivation of an approximate starting solution, forming the foundation for subsequent refinements. The research was particularly concentrated on the formulation of a fuel-optimal, low-thrust gravity-assist trajectory. The application of advanced techniques enhanced both the efficiency and robustness of the optimization process, culminating in optimal solutions. The numerical assessments carried out demonstrated that objective functions and constraints were precisely upheld throughout the computational execution. Unknowns and terminal values were adeptly ascertained, which bolstered the guarantee of global convergence when addressing the bang-bang optimal-control challenge. The performance of this method was further corroborated by a case study involving a Mars gravity assist with a final approach to Jupiter. Ultimately, the outcomes indicated improvements in the terminal spacecraft mass and the robustness of the solution, standing out favorably against existing methodologies.