Prerpints.org logo
Preprint
Article

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

Altmetrics

Downloads

143

Views

77

Comments

0

Submitted:

14 January 2024

Posted:

15 January 2024

You are already at the latest version

Alerts
Abstract
Recent advancements have revolutionized interplanetary spacecraft missions by combining low-thrust propulsion systems with planetary Gravity-Assist (GA) maneuvers. This study investigates/improved a novel dual-step efficient-robust homotopy algorithm, grounded in indirect optimal control principles, to analyze fuel-optimal bang-bang control problems. The methodology integrates efficient techniques, including the homotopic method, compatible normalization, and improved switching function detection, aimed at enhancing optimal convergence in solving low-thrust gravity-assist spacecraft missions. To illustrate practical application, an Earth-to-Jupiter interplanetary scenario utilizing Mars-GA is examined. Results demonstrate the algorithm's superior efficacy in achieving fuel efficiency, robust convergence, and remarkable accuracy compared to existing solutions.
Keywords: 
Subject: Engineering  -   Aerospace Engineering

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:
x ˙ = { r ˙ = v , v ˙ = − μ r 3 r + T max u m α , m ˙ = − T max I s p g 0 u ,
here, r and v denote position and velocity vectors of the spacecraft. The spacecraft's initial and final conditions are given as ψ t 0 = r 0 v 0 m 0 T and ψ t f = r f v f T , respectively, with the final mass m f left as a free parameter. Also, T m a x and u ∈ [ 0 , 1 ] signify the maximum thrust and engine throttle magnitudes. α represents the unit vector in the control direction, shaping the thrust vector as T = u T m a x α . m refers to the spacecraft’s instantaneous mass, while I s 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 ∞ 0 = 0 . According to Figure 1, while v ∞ 0 > 0 , the spacecraft initial mass is given by Tsiolkovsky’s rocket equation as follows:
m 0 ( v ∞ 0 ) = m 0 exp ( − Δ v v e x h a u s t ) ,
v ∞ 0 = v ∞ 0 u ∞ 0   ,   u ∞ 0 = ± v E ( t 0 ) ‖ v E ( t 0 ) ‖ , Δ v = v LEO ( 2 + ( v ∞ 0 v LEO ) 2 − 1 ) ,
here, v 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 r t G − = r t G + = r p l t G , 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 v t G − = v p l t G + v ∞ − and v t G + = v p l t G + v ∞ + , where v ∞ − = v ∞ + . Here, v ∞ ± and v p l t G represent GA hyperbolic excess velocity vectors and the velocity vector of the GA planet, respectively. Meanwhile, δ denotes the GA turn angle, and the hyperbolic periapsis radius of the GA orbit is parameterized as:
δ = arccos ( v ∞ − . v ∞ + v ∞ 2 )   ,   δ ∈ [ 0 , δ max ] ,
r p = μ p l v ∞ − v ∞ + ( 1 sin ( δ / 2 ) − 1 ) ,
The minimum acceptable GA radius ( r m i n ) limits the maximum turn angle ( δ m a x ⁡ ) as follows:
δ max = 2 arcsin ( 1 1 + r min μ p l v ∞ 2 ) ,
Hence, v ∞ −   a n d   v ∞ + possess identical magnitudes but differ only in the angle δ . Also, the gravity-assist framework can be written as:
i ^ = v ∞ − v ∞ ,   k ^ = v p l ( t G ) × v ∞ − ‖ v p l ( t G ) × v ∞ − ‖ , j ^ = k ^ × i ^ ,  
In Figure 3, α and δ are utilized to parameterize the velocity vector. Following this parameterization, the GA impulse is obtained from the velocity increment in the following manner:
Δ v G = v ( t G + ) − v ( t G − ) = v ∞ ( cos δ   i ^ + sin δ sin α   j ^ + sin δ cos α   k ^ ) ,

3. Trajectory optimization

The performance index aimed at minimizing propellant consumption is expressed as follows:
J = ∫ t 0 t f T max I S P g 0 u   d t ,
The homotopy performance index is employed to address the challenge of the optimal control problem [21,30]:
J a = λ 0 ∫ t 0 t f T max I S P g 0 [ u − ε u ( 1 − u ) ]   d t ,
In this context, ε = 1 corresponds to the energy-optimal problem, whereas ε = 0 corresponds to the fuel-optimal problem. The constant auxiliary multiplier   λ 0 belongs to the interval ( 0,1 ] . Since Pontryagin’s Maximum Principle (PMP) [49],   λ 0 , and the co-state vector   λ = λ 0 , λ r T , λ v T , λ m T are used to build the Hamiltonian as follows:
H = λ r T . v + λ v T . ( − μ r 3 r + T max u m α ) − λ m T max I S P g 0 u + λ 0 T max I S P g 0 [ u − ε u ( 1 − u ) ] ,
The Primer Vector [50], α = − λ v / λ v , is used to determine the optimal thrust direction for minimizing the Hamiltonian. In this way, the Euler–Lagrange conditions, λ ˙ = − ∂ H / ∂ x T , and state dynamics system (1) are used to establish the ordinary differential equations (ODEs):
x ˙ = { r ˙ = v , v ˙ = − μ r 3 r + T max u m α , m ˙ = − T max I s p g 0 u , λ ˙ r = μ r 3 λ v − 3 μ λ v .   r r 5 r , λ ˙ v = − λ r , λ ˙ m = − T max u m 2 ‖ λ v ‖ ,
Where the switching function (SF) and the extremal control magnitude are formed by ρ , which is also referred to as the switching function:
ρ = ( 1 − ‖ λ v ‖ I s p g 0 λ 0 m − λ m λ 0 ) ,
u = { 0         ρ > ε , 1         ρ < − ε , u = 1 2 − ρ 2 ε   | ρ | ≤ ε ,
Eqs. (15) should be satisfied as gravity-assist transversal conditions:
{ λ r ( t G + ) = λ r ( t G − ) − Z 1 ∼ 3 , λ v ( t G + ) = Z 4 u + + 1 r min κ B , λ v ( t G − ) = Z 4 u − − 1 r min κ A ,
where
u ± = v ∞ ± v ∞ , A = r p v ∞ − [ 1 4 sin 2 δ / 2 ( 1 − sin δ / 2 ) ( u + − cos δ u − ) − u − ] , B = r p v ∞ + [ 1 4 sin 2 δ / 2 ( 1 − sin δ / 2 ) ( u − − cos δ u + ) − u + ] ,
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 ( t G − ) − H ( t G + ) − Z 1 ~ 3 ⋅ v p l ( t G ) + Z 4 ( u + − u − ) ⋅ a p l ( t G ) − 1 r min κ C = 0
where a p l ( t G ) denotes the acceleration of the GA planet at the t G and C are written as:
a p l ( t G ) = − μ r p l 3 ( t G ) r p l ( t G ) C = r p { − 1 4 sin 2 δ / 2 ( 1 − sin δ / 2 ) [ 1 v ∞ + ( u − − cos δ u + ) + 1 v ∞ − ( u + − cos δ u − ) ] } . a p l ( t G )

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 r E and   v E ii) gravity-assist maneuver conditions are equal to the planet’s state vectors r p l and   v p l iii) the spacecraft's terminal state vectors are equal to the target body’s state vectors r f and   v f . Additionally, the spacecraft’s final mass is treated as a free parameter. Moreover, the terminal transversality condition is set to zero:
λ m ( t f , ε ) = 0 ,
Hence, according to Eqs. (2-3), we update the boundary conditions as follows:
ψ 0 = [ r ( t 0 ) − r E v ( t 0 ) − v E − v ∞ 0 u ∞ 0 m ( t 0 ) − m 0 ( v ∞ 0 ) ] = 0 ,
ψ f = [ r ( t f ) − r f v ( t f ) − v f λ m ( t f ) ] = 0   ,   m ( t f ) = m f ,
The normalization technique regarding references [21,30] is used as:
λ ( t , ε ) = λ ( t , ε ) ‖ λ ( t 0 , ε ) ‖ ;   λ ( t 0 , ε ) = [ λ 0   λ r ( t 0 )   λ v ( t 0 )   λ m ] ,
Moreover, Eqs. (23) define the intermediate gravity-assist constraints:
ψ ( t G ) = [ r ( t G ) − r p l ( t G ) ‖ v ∞ − ‖ − ‖ v ∞ + ‖ ] = 0 ,
σ G = 1 − r p r min ≤ 0 ,
Besides, rigid conditions need to be accounted for due to the inequality constraint, Eq. (24), as follows:
κ . σ G = 0 ,   κ ≥ 0 ,
Ultimately, Eq. (26) satisfies the normalization condition for the gravity-assist:
‖ [ λ ( t 0 , ε ) ;   Z 1 ~ 4 ;   κ ] ‖ = 1 ,

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   λ v t G − 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 Φ ) should be satisfied as:
Φ ( λ ( t 0 , ε ) ) = [ r ( t f ) − r f v ( t f ) − v f λ m ( t f ) r ( t G ) − r p l ( t G ) v ∞ − − v ∞ + λ v ( t G − ) − Z 4 u − + 1 r min κ A H ( t G − ) − H ( t G + ) − Z 1 ~ 3 ⋅ v p l ( t G ) ... + Z 4 ( u + − u − ) ⋅ a p l ( t G ) − 1 r min κ C κ . σ G ,   κ ≥ 0 ‖ [ λ ( t 0 , ε ) ;   Z 1 ~ 4 ;   κ ] ‖ − 1 ] = 0 ,
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 ( X ( 1 , k ) ) = ‖ Φ ( λ ( t 0 , 1 ) ) ‖ + J a ,
Since, the 12-D angle variables φ , considering the fixed values of t 0 ,   t G , and t f , are constructed utilizing the 12-D search-variables X ε , k in both steps by
φ 1 , 2 , 3 = π 2 X 1 , 2 , 3 φ 4 , 5 = π ( X 4 , 5 - 1 2 ) φ 6 , 7 = 2 π X 6 , 7 φ 8 = π ( X 8 - 1 2 ) φ 9 = π 2 X 9 φ 10 = 2 π X 10 φ 11 , 12 = π 2 X 11 , 12 ∈ [ 0 , π 2 ] , ∈ [ − π 2 , π 2 ] , ∈ [ 0 , 2 π ] , ∈ [ − π 2 , π 2 ] , ∈ [ 0 , π 2 ] , ∈ [ 0 , 2 π ] , ∈ [ 0 , π 2 ] ,
So, the 8-D co-states are constrained through a new approach during the computational process, restricted within the domain [-1,1] by
λ 0 = sin φ 1 , λ m ( t 0 , ε ) = cos φ 1 sin φ 2 cos φ 12 , λ r ( t 0 , ε ) = cos φ 1 cos φ 2 cos φ 3 cos φ 12 [ cos φ 4 cos φ 6 cos φ 4 sin φ 6 sin φ 4 ] T , λ v ( t 0 , ε ) = cos φ 1 cos φ 2 sin φ 3 cos φ 12 [ cos φ 5 cos φ 7 cos φ 5 sin φ 7 sin φ 5 ] T ,
Furthermore, the 5-D numerical multipliers related to the intermediate GA maneuver are normalized within the domain [-1,1] by
Z 1 = cos φ 1 sin φ 8 sin φ 12 , Z 2 ~ 4 = cos φ 1 cos φ 8 cos φ 9 sin φ 12 [ cos φ 10 cos φ 11 cos φ 10 sin φ 11 sin φ 10 ] T ,
κ = cos φ 1 cos φ 8 sin φ 9 sin φ 12 ,
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   X ∗ 1 , k 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   − 1,1 . The Runge-Kutta method is employed as the integrator, and   λ 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,   Z 1 ~ 4   κ , and the 3-D GA velocity vector, v t G + = v p l t G + v ∞ + . 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, λ v t G − , 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   Δ v ⃑ G Eq. (8), instead of the GA velocity vector v t G + . It is essential to note that in Eq. (15), λ r t G + and   λ v t G + 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 ' h ' 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:
ρ i + 1 = ρ i + ρ ˙ i h ,
The updated SF detection method includes five steps within the interval   t ∈ [ t 0 , t f ] :
Step (1): Initiated according to Flowchart Figure 5, detects   ρ i x :
Step (2): As long as   ρ i + 1 is greater than threshold ε , 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 ρ i + 1 falls below negative ε , the thrust magnitude is set to u = 1 . The SF detection procedure for this scenario is shown in Figure 7:
Step (4): When ρ i + 1 lies within the range − ε , ε , the thrust magnitude is set as u = 1 / 2 − ρ / 2 ε . 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 x i + 1 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 ∞ . 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 s p 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 ∞ 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.
Δ v G = X 15 μ p l r min [ cos δ cos α sin δ cos α sin α ] T , { α = π ( X 13 − 0.5 ) ∈ [ − π 2 , π 2 ] , δ = 2 π X 14 ∈ [ 0 , 2 π ] ,
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:
X X = [ v ∞ 0 ;   λ 0 ;   λ r ( t 0 ) ;   λ v ( t 0 ) ;   λ m ;   Z 1 ~ 4 ;   κ ;   Δ v G ] ,
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 ∞ 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, v ∞ ± , measures 3.654 km/s, along with a minimum periapsis radius r m i n of 500 km. The turn angles δ for fuel-optimal GA maneuvers are valued at 53.74°, occurring at the periapsis radius equal to r p = r m i n .
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 v ∞ ± , the GA impulsive velocity increment ∆ v G , and the 4-dimensional GA numerical multipliers Z 1 ~ 4 and the 1-dimensional GA numerical multiplier κ , 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 ρ ( t ,   0 ) attains a value of zero at six distinct points within the solution. Moreover, when ρ ( t ,   0 ) identifies a negative value, the resultant thrust profile escalates to maximum magnitude. Conversely, the thrust magnitude drops to zero whenever ρ ( 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 × 10 − 5 km and ∆ v f = 2.631 × 10 − 14 km/s.
In addition, Figure 11 shows that   λ m ∗ t f , 0 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   λ 0 = 0.74369984308156 between the sequence solutions from the energy- to fuel-optimal problems.
Instantaneous mass ratio, m ( t ) / 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 λ ∗ ( t , 0 ) = λ 0 , λ r T , λ v T , λ m T remains consistently confined within the   − 1,1 domain, in sharp contrast to broader solution spaces ranging   − 100000,100000 , where global convergence poses significant computational challenges for interpolating energy- and fuel-optimal solutions. Moreover, Figure 11 illustrates that   λ m ∗ t f , 0 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   λ 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.

References

  1. Wagner, A.S. ; Automated trajectory design for impulsive and low thrust interplanetary mission analysis. Ph.D. dissertation, Iowa State University, 2014.
  2. Bernelli-Zazzera, F.; Lavagna, M.; Armellin, R.; Di Lizia, P.; Morselli, A.; Olympio, J.; Izzo, D.; Summerer, L. Trajectory optimisation under uncertainties. ESA/ACT, Ariadna Final Report, id., 2012, 10, 4101. [Google Scholar]
  3. Williams, S.N.; Coverstone, V. Benefits of solar electric propulsion for the next generation of planetary exploration missions. Journal of the Astronautical Sciences, 1997, 45, 143–159. [Google Scholar] [CrossRef]
  4. Miele, A.; Wang, T. Multiple-subarc gradient-restoration algorithm, part 1: algorithm structure. Journal of Optimization Theory and Applications 2023, 116, 1–17. [Google Scholar] [CrossRef]
  5. Carter, T.; Humi, M. Planar optimal two-impulse transfers with closed-form solutions of the transverse transfers. Journal of Optimization Theory and Applications, 2016, 169, 262–279. [Google Scholar] [CrossRef]
  6. Pontani, M.; Miele, A. Theorem of optimal image trajectories in the restricted problem of three bodies. Journal of Optimization Theory and Applications, 2016, 168, 992–1013. [Google Scholar] [CrossRef]
  7. Gong, S.; Macdonald, M. Review on solar sail technology. Astrodynamics, 2019, 3, 93–125. [Google Scholar] [CrossRef]
  8. Mereta, A.; Izzo, D. Target selection for a small low-thrust mission to near-Earth asteroids. Astrodynamics, 2018, 2, 249–263. [Google Scholar] [CrossRef]
  9. Rayman, M.D.; Varghese, P.; Lehman, D.H.; Livesay, L.L. Results from the deep space 1 technology validation mission. Acta Astronautica, 2000, 47, 475–487. [Google Scholar] [CrossRef]
  10. Kuninaka, H.; Nishiyama, K.; Funakai, I.; Tetsuya, K.; Shimizu, Y.; Kawaguchi, J.I. Asteroid rendezvous of HAYABUSA explorer using microwave discharge ion engines. In 29th International Electric Propulsion Conference, 2005, number IEPC-2005-10.
  11. Kawaguchi, J.; Fujiwara, A.; Uesugi, T.K. The ion engines cruise operation and the earth swingby of “Hayabusa” (Muses-C). Space Technology, 2005, 25, 105–115. [Google Scholar]
  12. Rayman, M.D.; Fraschetti, T.C.; Raymond, C.A.; Russell, C.T. Dawn: A mission in development for exploration of main belt asteroids Vesta and Ceres. Acta Astronautica, 2006, 58, 605–616. [Google Scholar] [CrossRef]
  13. Kugelberg, J.; Bodin, P.; Persson, S.; Rathsman, P. Accommodating electric propulsion on SMART-1. Acta Astronautica, 2004, 55, 121–130. [Google Scholar] [CrossRef]
  14. Betts, J. Survey of Numerical Methods for Trajectory Optimization. Journal of Guidance, Control, and Dynamics, 1998, 21, 193–207. [Google Scholar] [CrossRef]
  15. Conway, B.A. A Survey of Methods Available for the Numerical Optimization of Continuous Dynamic Systems. Journal of Optimization Theory and Applications, 2011, 152, 271–306. [Google Scholar] [CrossRef]
  16. Wang, J.; Fečkan, M.; Debbouche, A. Time Optimal Control of a System Governed by Non-instantaneous Impulsive Differential Equations. Journal of Optimization Theory and Applications, 2019, 182, 573–587. [Google Scholar] [CrossRef]
  17. Betts, J.T.; Frank, P.D. A sparse nonlinear optimization algorithm. Journal of Optimization Theory and Applications, 1994, 82, 519–541. [Google Scholar] [CrossRef]
  18. Ranieri, C.L.; Ocampo, C.A. Indirect Optimization of Three-Dimensional Finite-Burning Interplanetary Transfers Including Spiral Dynamics. Journal of Guidance, Control, and Dynamics 2009, 32, 444–454. [Google Scholar] [CrossRef]
  19. Murray, D.; Yakowitz, S. Differential dynamic programming and Newton’s method for discrete optimal control problems. Journal of Optimization Theory and Applications, 1984, 43, 395–414. [Google Scholar] [CrossRef]
  20. Saghamanesh, M.; Novinzadeh, A. Optimal guidance of spacecraft rendezvous via free initial velocity vector. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2014, 228, 1020–1032. [Google Scholar] [CrossRef]
  21. Saghamanesh, M.; Baoyin, H. A robust homotopic approach for continuous variable low-thrust trajectory optimization. Advances in Space Research, 2018, 62, 3095–3113. [Google Scholar] [CrossRef]
  22. Tang, G.; Hauser, K. A data-driven indirect method for nonlinear optimal control. Astrodynamics. 2019, 3, 345–59. [Google Scholar] [CrossRef]
  23. Enright, P.J.; Conway, B.A. Discrete Approximations to Optimal Trajectories Using Direct Transcription and Nonlinear Programming. Journal of Guidance, Control, and Dynamics, 1992, 15, 994–1002. [Google Scholar] [CrossRef]
  24. Hargraves, C.R.; Paris, S.W. Direct trajectory optimization using nonlinear programming and collocation. Journal of Guidance, Control, and Dynamics, 1987, 10, 338–342. [Google Scholar] [CrossRef]
  25. Han, S.P. A globally convergent method for nonlinear programming. Journal of Optimization Theory and Applications, 2012, 22, 297–309. [Google Scholar] [CrossRef]
  26. Lantoine, G.; Russell, R.P. A hybrid differential dynamic programming algorithm for constrained optimal control problems. Part 1: Theory. Journal of Optimization Theory and Applications, 2012, 154, 382–417. [Google Scholar] [CrossRef]
  27. Lantoine, G.; Russell, R.P. A hybrid differential dynamic programming algorithm for constrained optimal control problems. part 2. Application. Journal of Optimization Theory and Applications, 2012, 154, 418–442. [Google Scholar] [CrossRef]
  28. Cheng, L.; Wang, Z.; Jiang, F. Real-time control for fuel-optimal moon landing based on an interactive deep reinforcement learning algorithm. Astrodynamics, 2019, 3, 375–86. [Google Scholar] [CrossRef]
  29. Bertrand, R.; Epenoy, R. New smoothing techniques for solving bang–bang optimal control problems-numerical results and statistical interpretation. Optimal Control Applications and Methods, 2012, 23, 171–197. [Google Scholar] [CrossRef]
  30. Jiang, F.; Baoyin, H.; Li, J. Practical techniques for low-thrust trajectory optimization with homotopic approach. Journal of Guidance, Control, and Dynamics, 2012, 35, 245–258. [Google Scholar] [CrossRef]
  31. Hans, S.; Renjith, R.K. Finite difference scheme for automatic costate calculation. Journal of Guidance, Control, and Dynamics, 1996, 19, 231–239. [Google Scholar]
  32. Martell, C.A.; Lawton, J.A. Adjoint variable solutions via an auxiliary optimization problem. Journal of Guidance, Control, and Dynamics, 1995, 18, 1267–1272. [Google Scholar] [CrossRef]
  33. Saghamanesh, M.; Baoyin, H. Practical Mission Design for Mars Trajectory Optimization Based on the Ephemeris Model and Full Perturbation System. Advances in the Astronautical Sciences, 2018, 165, 1293–1313. [Google Scholar]
  34. Saghamanesh, M.; Taheri, E.; Baoyin, H. Systematic Low-Thrust Trajectory Design to Mars Based on a Full Ephemeris Modelling. Advances in Space Research, 2019, 64, 2356–78. [Google Scholar] [CrossRef]
  35. Saghamanesh, M.; Taheri, E.; Baoyin, H. Interplanetary Gravity-Assist Fuel-Optimal Trajectory Optimization with Planetary and Solar Radiation Pressure Perturbations. Celestial Mechanics and Dynamical Astronomy, 2019, 132, 1–21. [Google Scholar] [CrossRef]
  36. Saghamanesh, M.; Baoyin, H. Optimal Gravity-Assist Low-Thrust Trajectory Using a Robust Homotopic Approach. IEEE/CSAA Guidance, Navigation and Control 2018, 1–8. [Google Scholar]
  37. Hansen, C.J.; Bolton, S.J.; Matson, D.L.; Spilker, L.J.; Lebreton, J.P. The Cassini/ Huygens flyby of Jupiter. Icarus 2004, 172, 1–8. [Google Scholar] [CrossRef]
  38. McNutt Jr, R.L.; Solomon, S.C.; Grard, R.; Novara, M.; Mukai, T. An international program for Mercury exploration: synergy of MESSENGER and BepiColombo. Advances in Space Research, 2004, 33, 2126–2132. [Google Scholar] [CrossRef]
  39. McNutt Jr, R.L.; Solomon, S.C.; Gold, R.E.; Leary, J.C. MESSENGER Team: The MESSENGER mission to Mercury: Development history and early mission status. Advances in Space Research, 2006, 38, 564–571. [Google Scholar] [CrossRef]
  40. Hechler, M. Rosetta mission design. Advances in Space Research 1997, 19, 127–136. [Google Scholar] [CrossRef]
  41. Guo, Y.; Farquhar, R.W. New Horizons mission design. Space science reviews 2008, 140, 49–74. [Google Scholar] [CrossRef]
  42. D’Amario, L.A.; Bright, L.E.; Wolf, A.A. Galileo trajectory design. Space Science Reviews, 1992, 60, 23–78. [Google Scholar]
  43. Curtis, H. Orbital mechanics for engineering students. Butterworth-Heinemann, 2013.
  44. Strange, N.J.; Longuski, J.M. Graphical method for gravityassist trajectory design. Journal of Spacecraft and Rockets, 2002, 30, 9–16. [Google Scholar] [CrossRef]
  45. Prado, A.F.B. Powered swingby. Journal of Guidance, Control, and Dynamic, 1996, 19, 1142–1147. [Google Scholar]
  46. Felipe, G.; Prado, A.F.B. Trajectory selection for a spacecraft performing a two-dimensional swing-by. Advances in Space Research 2004, 34, 2256–2261. [Google Scholar]
  47. Broucke, R.A. The celestial mechanics of gravity assist. AIAA/AAS Astrodynamics Conference, 1988, 15-18.
  48. Yam, C.H.; McConaghy, T.T.; Chen, K.J.; Longuski, J.M. Preliminary design of nuclear electric propulsion missions to the outer planets. AIAA/AAS Astrodynamics Specialist Conference, 2004, 5393, 16–19. [Google Scholar]
  49. Trélat, E. Optimal control and applications to aerospace: some results and challenges. Journal of Optimization Theory and Applications, 2012, 154, 713–758. [Google Scholar] [CrossRef]
  50. Russell, R.P. Primer Vector Theory Applied to Global Low-Thrust Trade Studies. Journal of Guidance, Control, and Dynamics, 2007, 30, 460–472. [Google Scholar] [CrossRef]
1
https://ssd.jpl.nasa.gov/horizons.cgi [Retrieved on: Feb. 20, 2019]
Preprints 96321 g001
Figure 1. Departure of the spacecraft from Earth parking orbit.
Figure 1. Departure of the spacecraft from Earth parking orbit.
Preprints 96321 g001
Preprints 96321 g002
Figure 2. GA maneuver in a two-body system.
Figure 2. GA maneuver in a two-body system.
Preprints 96321 g002
Preprints 96321 g003
Figure 3. GA impulsive model.
Figure 3. GA impulsive model.
Preprints 96321 g003
Preprints 96321 g004
Figure 4. Computational flowchart of the dual-step optimization algorithm.
Figure 4. Computational flowchart of the dual-step optimization algorithm.
Preprints 96321 g004
Preprints 96321 g005
Figure 5. Initial step of SF detection.
Figure 5. Initial step of SF detection.
Preprints 96321 g005
Preprints 96321 g006
Figure 6. Flowchart illustrating the SF detection procedure when ρ i exceeds ε .
Figure 6. Flowchart illustrating the SF detection procedure when ρ i exceeds ε .
Preprints 96321 g006
Preprints 96321 g007
Figure 7. Flowchart illustrating the SF detection procedure when ρ i is less than − ε .
Figure 7. Flowchart illustrating the SF detection procedure when ρ i is less than − ε .
Preprints 96321 g007
Preprints 96321 g008
Figure 8. SF detection procedure for when ρ i is within the range − ε , ε .
Figure 8. SF detection procedure for when ρ i is within the range − ε , ε .
Preprints 96321 g008
Preprints 96321 g009
Figure 9. Earth-Mars-Jupiter energy- and fuel-optimal thrust profiles with switching function detection for bang-bang optimal-control.
Figure 9. Earth-Mars-Jupiter energy- and fuel-optimal thrust profiles with switching function detection for bang-bang optimal-control.
Preprints 96321 g009
Preprints 96321 g010
Figure 10. Earth-Mars-Jupiter GA fuel-optimal low-thrust spacecraft trajectory.
Figure 10. Earth-Mars-Jupiter GA fuel-optimal low-thrust spacecraft trajectory.
Preprints 96321 g010
Preprints 96321 g011
Figure 11. Fuel-optimal consumption with co-states vectors.
Figure 11. Fuel-optimal consumption with co-states vectors.
Preprints 96321 g011
Table 1. Boundary values and parameters pertinent to the Earth-Mars-Jupiter mission.
Table 1. Boundary values and parameters pertinent to the Earth-Mars-Jupiter mission.
Parameter Value
Initial Position, km [-1.25326142025922E+08, 7.83493051133693E+07, -3.62097454475238E+03]
Initial Velocity, km/s [-1.62697077142250E+01, -2.53616129608122E+01, 1.47459860098919E-03]
Final Position, km [-6.20180296761522E+08, 5.04570725212524E+08, 1.17797453566354E+07]
Final Velocity, km/s [-8.40428050226635E+00, -9.53464925506379E+00, 2.27651689048206E-01]
Initial time, TDB Feb 17, 2022
Time of Flight, day 1809
GA time, TDB March 18, 2024
Final time, TDB Feb 01, 2027
m 0 , kg; ( v ∞ 0 = 0 , k m / s ) 20000
m 0 , kg; ( v ∞ 0 = 0.8 , k m / s ) 19819
I s p , s 6000
T m a x , N 2.26
Table 2. Optimal co-states for E-M-J mission.
Table 2. Optimal co-states for E-M-J mission.
λ ∗ t 0 , ε = [ λ 0 ; λ r ; λ v ; λ m ]
ε = 1 ε = 0
[0.743699843081560;
0.352463239633363, -0.277598991926813, -0.079343331137383;
0.023718570580293, 0.049511669434066, 0.003327795085571;
0.457631920898898]		 [0.743699843081560;
0.192035090839654, -0.149833850785203, -0.063080789937559;
0.011833279537100, 0.026101523172379, 0.002858568363144;
0.370611693389334]
Table 3. The Comparison Results.
Table 3. The Comparison Results.
Our results GALLOP
GA time, t G , TDB March 18, 2024 March 18, 2024
TOF, day 1,809 1,809
GA altitude, r m i n , km 500 500
v ∞ 0 , km/s 0.8 1.04
GA v ∞ ± , km/s 3.654 --
δ ° 53.74 --
σ G 1.2256e-16 --
Final mass, m f , kg 15,068.17 14,879
GA time, t G , TDB March 18, 2024 March 18, 2024
Table 4. GA Numerical Results of E-M-J Mission.
Table 4. GA Numerical Results of E-M-J Mission.
Parameter Value
r M a r s ( t G ) , km [1.16465658375768E+08, -1.73955158925876E+08, -6.50236091435127E+06]
v M a r s ( t G ) , km/s [2.10481455623537E+01, 1.55555848401697E+01, -1.90278457157558E-01]
GA v ∞ − , AU/yr [0.524178202715554, -0.565039214276219, 0.008672119576309]
GA v ∞ + , AU/yr [0.765632779887656, 0.088191564273440, -0.011661984497357]
GA ∆ v G , AU/yr [0.241454577172102, 0.653230778549659, -0.020334104073666]
GA   Z 1 ~ 4 ; κ [0.008360758766825, -0.062264141998319, 0.078717645223935, -0.053488824816944; 0.030663561324637]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated