Framework for the Multi-Objective Design Optimization of Aerocapture Missions

1. Introduction

A critical step in planetary research is capturing a spacecraft from its hyperbolic trajectory to a closed orbit around the target planet. The traditional insertion burn method requires significant fuel allocation, reducing payload mass fraction. The aerocapture maneuver offers a solution that bypasses the need for large fuel fractions, potentially increasing payload capacity. Despite extensive study and optimization efforts, its implementation in interplanetary missions is yet to be realized due to spacecraft architecture difficulties and mission outcome sensitivity to minor, unpredictable uncertainties.

Recent optimization efforts in this field [1,2,3,4] have centered on the pursuit of efficient trans-atmospheric trajectories through meticulous manipulation of drag and lift forces to achieve an optimal aerocapture maneuver. Recent years have witnessed a push towards enhancing aerocapture technology through drag modulation (DM) techniques. These methods enable control over the spacecraft’s ballistic coefficient ($\beta$) by manipulating either the reference area or the drag coefficient ($C_D$) through continuous or discrete strategies, simplifying avionics and obviating the need for complex attitude control systems. However, certain technical challenges remain [5]. In contrast, with a history of successful missions, lift modulation relies on complex avionics algorithms and control systems for continuous bank angle modulation. Despite the complexities associated with this technique, it remains better understood and established than DM techniques. Optimal aerocapture problems have often overlooked the complex engineering intricacies inherent in effecting control over these parameters. Furthermore, the influence of aerodynamic control on the resulting captured orbit is limited, and its nature is highly dependent upon the spacecraft’s state at the atmospheric interface $\left(AI\right)$. Consequently, optimizing the atmospheric trajectory yields only incremental benefits to the overall mission, primarily impacting atmospheric phenomena such as heat loads.

Current optimization strategies face challenges due to highly nonlinear physics, requiring intricate nonlinear programming methods to address potential non-convexity. While Liu et al. [2] explored optimizing nonlinear, heavily constrained re-entry trajectories by reframing the nonlinear problem as a convex one, their approach compromises generalizability, limiting its applicability. Hongwei et al. [6] tackle the inherent non-linearity of the problem by reformulating it into a convex sub-problem. This reformulation enables the efficient resolution of the original problem through a sequence of convex optimization tasks. Nevertheless, this approach is also applied to the optimal guidance of the spacecraft during re-entry. Additionally, current strategies often employ multi-objective optimization (MOO) [1,2,3,4] for optimal re-entry trajectory problems, presenting numerical and decision-making challenges. MOO involves navigating a multi-dimensional objective space, resulting in trade-offs between objectives. The exploration of vast potential solutions and the non-convex and irregular Pareto front further complicates optimal solution searches. Selecting a single solution from a multitude of trade-offs involves considering subjective preferences, priorities, and uncertainty, requiring informed judgment, which can be challenging for algorithms due to ambiguous decision-making criteria. Moreover, the AMAT framework, developed by Girija et al. [7] based on previous work by Girija [8], offers a Python-based tool to assist in the design of re-entry vehicles. This framework is specifically geared towards facilitating trade-off analysis rather than optimizing the design process.

In light of these multifaceted challenges, we propose a novel framework, based around a new algorithm for the Determination of Aerocapture Successful Trajectories and Robust Optimization (D-ASTRO) algorithm. Developed within the MATLAB environment, D-ASTRO is an efficient and robust tool for rapidly analyzing and optimizing generic aerocapture missions. D-ASTRO’s unique capability lies in its ability to rapidly optimize non-thrusting, fixed angle-of-attack aerocapture-capable spacecraft and trajectories. By requiring only minimal aeroshell information and mission parameters, D-ASTRO employs a nonlinear, single-objective constrained optimization strategy to identify optimal values for the flight path angle ($\gamma$) at the AI and the fixed ballistic coefficient ($\beta$) of the aeroshell.

This paper’s structure proceeds as follows: In Section 2, we outline the equations of motion and physics models employed by D-ASTRO. Section 3 introduces the concept of the aerocapture corridor. The optimization problem, along with novel objective functions and performance metrics, is discussed in Section 4. The novel methodology used to rapidly compute the aerocapture corridor rapidly is covered in Section 5 while Section 6 focuses on the Numerical implementation and pipeline of the algorithm. Finally, Section 7 presents the algorithm’s numerical stability and convergence properties applied to a typical Martian aerocapture mission test case, with computational considerations and performance of D-ASTRO explored in Section .

2. Equations of Motion and Physics Model

2.1. Atmospheric Flight Dynamics

This framework supports the preliminary design and optimization of aerocapture-capable spacecraft. At this stage, key geometric parameters such as the inertia matrix, center of mass, and center of pressure are undefined, limiting the analysis of rotational kinematics during re-entry. Consequently, D-ASTRO employs a three-degree-of-freedom (3-DOF) model, following the equations outlined by Kumar and Tewari [9]:

$$\begin{array}{c}\dot{\delta }={\displaystyle \frac{Vcos\left(\gamma \right)cos\left(\chi \right)}{rcos\left(\lambda \right)}}-\Omega \\ \dot{\lambda }={\displaystyle \frac{Vcos\left(\gamma \right)sin\left(\chi \right)}{r}}\\ \dot{r}=Vsin\left(\gamma \right)\\ \dot{V}=-{\displaystyle \frac{\mu sin\left(\gamma \right)}{r^2}}-{\displaystyle \frac{D}{m}}+{\Omega }^{2}rcos\left(\lambda \right)[sin\left(\gamma \right)cos\left(\lambda \right)-cos\left(\gamma \right)sin\left(\chi \right)sin\left(\lambda \right)]\\ \begin{array}{c}\hfill \dot{\gamma }=\left({\displaystyle \frac{V}{r}}-{\displaystyle \frac{\mu }{V{r}^2}}\right)cos\left(\gamma \right)+{\displaystyle \frac{L}{Vm}}+2\Omega cos\left(\chi \right)cos\left(\lambda \right)+\\ \hfill {\displaystyle \frac{{\Omega }^{2}rcos\left(\lambda \right)}{V}}[cos\left(\gamma \right)cos\left(\lambda \right)+sin\left(\gamma \right)sin\left(\chi \right)sin\left(\lambda \right)]\end{array}\\ \dot{\chi }={\displaystyle \frac{-V}{r}}cos\left(\gamma \right)cos\left(\chi \right)tan\left(\lambda \right)+2\Omega \left[tan\left(\gamma \right)sin\left(\chi \right)cos\left(\lambda \right)-sin\left(\lambda \right)\right]\end{array}$$

(1)

where r is the radial distance from the planet’s center to the vehicle, V is the planetocentric velocity, $\gamma$ is the flight path angle, $\delta$ and $\lambda$ are the longitude and latitude, respectively, and $\chi$ is the heading angle. Additionally, m denotes the spacecraft’s mass, L and D are the lift and drag forces, respectively, and $\mu$ and $\Omega$ represent the planet’s gravitational parameter and rotational period, respectively.

Let the trajectory state in Equation (1) be $\mathbf{x}={\left[r,V,\gamma ,\delta ,\lambda ,\chi \right]}^T$. The initial condition of $\mathbf{x}$ is assumed to be user-specified and coincident with the atmospheric interface (AI).

2.2. Aerodynamics Model

The aerodynamic model assumes the Mach independence principle, simplifying the determination of aerodynamic drag for hypersonic continuum flows, resulting in a single $C_D$ and $L/D$ that characterize the entire high-altitude trajectory. However, phenomena such as gas rarefaction at high altitudes can cause deviations from the Mach independence principle. To account for this, D-ASTRO allows users to input $C_D$ and $C_L$ as a function of Mach number, M, Reynolds number, $Re$, and Knudsen number, $Kn$ for greater accuracy, as described by Equations (2) and ().

$$\begin{array}{c}D={\displaystyle \frac{1}{2}}\rho V^2{S}_{ref}{C}_D\left(Kn,M,Re\right)\end{array}$$

(2)

$$\begin{array}{c}L={\displaystyle \frac{1}{2}}\rho V^2{S}_{ref}{C}_L\left(Kn,M,Re\right)\mathrm{or}D\left({\displaystyle \frac{L}{D}}\right)\end{array}$$

(3)

where $\rho$ is the freestream density and $S_{ref}$ is the reference area of the aeroshell

2.3. Aerothermal Model

To compute the heat flux, ${\dot{q}}_{max}$ and heat load, Q, of the trajectory is a simplistic manner, the thermal model focuses on stagnation point aerothermal heating, as this is the expected peak heat transfer region of the aeroshell. Thus, the Sutton-Graves relation [10], Equation (4), approximates the convective stagnation point heating. Furthermore, D-ASTRO incorporates the Tauber-Sutton relations [11], Equation (), to account for radiative heat transfer, ${\dot{q}}_{rad}$, in cases where re-entry speeds are moderately high.

$$\begin{array}{c}{\dot{q}}_{convective}=k{\left({\displaystyle \frac{\rho }{R_N}}\right)}^{0.5}{V}^3\end{array}$$

(4)

$$\begin{array}{c}{\dot{q}}_{radiative}=C{R_N}^a{\rho }^{m}f\left(V\right)\end{array}$$

(5)

where $R_N$ is the nose radius, k, C, a and m are atmospheric-specific constants, and $f\left(V\right)$ is an empirically determined function that is also planet-specific.

In cases where wall temperature is a concern, radiative equilibrium approximates the equilibrium wall temperature, $T_W$. Assuming an energy balance between total heating, ${\dot{q}}_{tot}={\dot{q}}_{convective}+{\dot{q}}_{radiative}$, and re-radiative cooling, ${\dot{q}}_{re-rad}$, calculated using the Stefan-Boltzmann law, Equation (7) $T_W$ is approximated by (Equation 8).

$${\dot{q}}_{re-radiation}=ϵ\sigma {\left(T_W\right)}^4$$

(7)

$$T_W=\sqrt[4]{{\displaystyle \frac{{\dot{q}}_{tot}}{ϵ\sigma }}}$$

(8)

where $ϵ$ is the emissivity and $\sigma$ is the Boltzmann constant.

2.4. Atmospheric Model

A reliable atmospheric model is essential for accurately calculating the vehicle’s aerodynamic forces and aerothermal heating during the aerocapture maneuver. This report investigates Terrestrial and Martian atmospheres, highlighting their specific models below.

2.4.1. Earth Atmospheric Model

Given the high altitude characteristic of aerocapture trajectories, the U.S. standard atmospheric model [12] was deemed inadequate. Instead, the Earth atmosphere model proposed by NASA [13] (Equation (8)) will be utilized as an initial approximation. As discussed in Section 3, the robustness of an aerocapture corridor hinges on uncertainties in atmospheric properties, particularly density. More reliable atmospheric models exist that incorporate seasonal and global positioning variations in atmospheric quantities, like NASA’s Earth Global Reference Atmospheric Model (Earth-GRAM) [14] and COSPAR International Reference Atmosphere (CIRA) model [15].

To account for atmospheric perturbations, the research by Leslie and Justus [14] is used, indicating a maximum density fluctuation of 25% for ±3 standard deviations in the 0-80$km$ altitude range. As Terrestrial aerocapture maneuvers are not the main focus of this study, ±30% fluctuations in density corresponding to 3 standard deviation perturbations are assumed to approximate uncertainties in density.

$$\begin{array}{cc}\hfill T& =141.94+0.00299h,\hfill \\ \hfill p& =2.488{\left[{\displaystyle \frac{T}{216.7}}\right]}^{-11.388},\hfill \end{array}$$

(9)

which is applicable for $h>25\mathrm{km}$.

2.4.2. Mars Atmospheric Model

Mars-GRAM [16] was employed to generate atmospheric data for Martian aerocapture maneuvers. A simplified table was generated from the extensive data available in the Mars-GRAM database, encompassing atmospheric variations concerning seasons, latitude and longitude positions, and numerous other variables. This table is based on a yearly average of equatorial properties at zero longitude. Additionally, a statistical analysis was conducted to ascertain an accurate variation of atmospheric properties for $\pm 3\sigma$, which encompasses 99.73% of cases.

2.5. Orbital Element Determination

Determining orbital elements is crucial for characterizing a spacecraft’s orbit after an atmospheric pass. Identifying the nature of this orbit is essential to assess the success of the aerocapture maneuver. Additionally, these elements play a vital role in the proposed optimization strategy, as they will dictate the required fuel for corrective maneuvers. D-ASTRO employs an orbital determination algorithm based on the method proposed by Curtis [17], summarized in the pseudo-code shown in Algorithm , to compute the orbital elements.

3. Aerocapture Corridor

Without employing lift or drag modulation, the flight path angle $\gamma$ at AI and $\beta$ are the primary factors governing the vehicle dynamics during the atmospheric pass, thus determining the success of the aerocapture maneuver. The aerocapture corridor defines the range of $\gamma$ and $\beta$ values that result in a successful maneuver. The lower limit occurs when excessive drag causes the spacecraft to collide with the surface. In contrast, the upper limit occurs when insufficient drag causes the vehicle to remain in a hyperbolic orbit. The lower limit is easily computed when the altitude reaches zero at any point in the trajectory. The spacecraft’s post-pass orbital eccentricity e determines the upper limit, where $e<1$ indicates a closed orbit.

Figure 1. Example of Martian aerocapture corridors with post-atmospheric pass orbital eccentricity contours: (a) Ballistic re-entry profile, $L/D=0.0$. (b) Lifted re-entry profile, $L/D=0.2$

Figure 1. Example of Martian aerocapture corridors with post-atmospheric pass orbital eccentricity contours: (a) Ballistic re-entry profile, $L/D=0.0$. (b) Lifted re-entry profile, $L/D=0.2$

Figure compares the aerocapture corridors for ballistic and lifting re-entry profiles of a spacecraft approaching Mars with a hyperbolic excess velocity, $V_e$, of 3.5 ${\mathrm{kms}}^{-1}$. The maneuver’s success is represented by green hashing in Figure . These contours are depicted in the figures. These contours illustrate the intricate relationship among $\gamma$, $\beta$, and the lift-to-drag ratio necessary for successful maneuvers. A lifting re-entry profile provides a more favorable corridor, reducing sensitivity to precise insertion attitudes. While D-ASTRO supports various re-entry profiles, this study focuses on non-zero $L/D$ profiles for Martian aerocapture, utilizing a fixed value of 0.2, consistent with established norms for Martian re-entry vehicles [18].

Uncertainties in atmospheric properties significantly influence the width of the aerocapture corridor. Variations in atmospheric density lead to changes in drag forces, causing the spacecraft to dissipate energy differently than expected. Near the aerocapture boundaries, even minor density fluctuations can result in maneuver failure, leading to either a hyperbolic orbit or surface collision. Two corridors are computed to establish a robust corridor: one with reduced atmospheric density and another with increased density. The `robust’ aerocapture corridor represents the intersection of these two, resulting in a narrowed corridor. These corridors are illustrated in Figure 2, depicting Terrestrial and Martian corridors for various inbound $V_e$, where the green and yellow regions represent the robust and non-robust corridors.

4. Optimal Aerocapture Maneuver via Parameter Optimization

Multiple approaches can be taken to optimize an aerocapture trajectory. However, we focus on the optimization of design parameters $\gamma$ and $\beta$ which yield a successful aerocapture trajectory. Determining these two parameters allows mission planners to identify the most suitable interplanetary orbit and structural design of the heat shield to ensure a successful and optimal trajectory.

4.1. Formulation of Optimization Problem

The optimal aerocapture corridor problem is a constrained nonlinear optimization problem of the form:

where f represents the objective function, and $g_u$ and $g_l$ are functions of $\beta$ that define the upper and lower boundaries of the aerocapture corridor. The vector $\mathbf{\theta }$ comprises selected performance indices.

Insertion inaccuracies during aerocapture missions can arise from accumulating minor guidance errors during the interplanetary trajectory. This leads to significant deviations from the intended atmospheric insertion coordinates. Deviating from the intended ${\gamma }_{AI}$ leads to different atmospheric trajectories, impacting the aerocapture maneuver’s optimality and sometimes the success. To account for this, an expected insertion inaccuracy is incorporated as $\Delta \gamma$, which narrows the aerocapture corridor to ensure mission success is robust to $\pm \Delta \gamma$ insertion inaccuracies.

4.2. Performance Indices

This section introduces a series of typically used performance metrics for re-entry trajectories.

4.2.1. Thermal Environment

The corresponding aerocapture trajectory is computed for every candidate point, and the ${\dot{q}}_{tot}$ profile is available. Hence, ${\dot{q}}_{tot,max}$ is computed by locating the maximum value in the profile and Q by using Equation (10).

$$Q={\int }_{t_0}^{t_t}{\dot{q}}_{tot}dt$$

(10)

where $t_t$ is the total duration of the atmospheric flight.

4.2.2. Volumetric Constraints

To address fairing constraints and current technological limitations, this optimization process aims to minimize the volume of the aeroshell by reducing the “equivalent axisymmetric aeroshell radius” , $R_B$, computed by rearranging Equation (11). This metric is preferred because it can represent any aeroshell configuration, from axisymmetric to complex high lift-to-drag designs, as long as $C_D$ and $S_{ref}$ are known.

$$\beta ={\displaystyle \frac{m}{C_{D}S}}={\displaystyle \frac{m}{C_D\pi R_B^2}}$$

(11)

4.2.3. Orbital Corrective Maneuvers

$\Delta V$ represents the fuel needed to adjust the post-atmospheric pass orbit to the predetermined target orbit. In reality, each candidate trajectory will have its unique optimal orbital corrective maneuver composed of complex transfers. To simplify comparisons between all candidate trajectories, trivial orbit transfers are used. These include two-impulse Hohmann transfers to raise the orbit’s periapsis and match the desired target orbital eccentricity. In addition, a plane change burn at the apoapsis is assumed to match the target orbital inclination. The $\Delta V$ required can be computed using Equation (12).

$$\begin{array}{c}\hfill \Delta V_{TOTAL}=\sqrt{\mu }\left[\sqrt{{\displaystyle \frac{1}{a_t}}}+\sqrt{{\displaystyle \frac{2}{a_i\left(e_i+1\right)}}-{\displaystyle \frac{2}{a_i\left(e_i+1\right)+a_t}}}-\sqrt{{\displaystyle \frac{2}{a_t}}-{\displaystyle \frac{2}{\left(a_i\left(e_i+1\right)+a_t\right)}}}\right.\\ \hfill -\left.\sqrt{{\displaystyle \frac{2}{a_i\left(e_i+1\right)}}-{\displaystyle \frac{1}{a_i}}}\right]+2\sqrt{\mu \left({\displaystyle \frac{2}{a_i\left(e_i+1\right)}}-{\displaystyle \frac{1}{a_i}}\right)}sin\left({\displaystyle \frac{i_t-i_i}{2}}\right)\end{array}$$

(12)

where a is the semi-major axis of the orbit and subscripts `i’ and `t’ , denote the quantity $\left(·\right)$ evaluated at the post-atmospheric pass orbit (determined at AE) and the target operational orbit, respectively.

4.3. Single-Objective Function Optimization

This paper considers two formulations of optimization Problem Section 4.1. The first strategy is the most widely employed multi-objective optimization (MOO), where $i\ne 1$, and a novel single objective optimization formulation is presented, where $i=1$.

MOO is the preferred approach for re-entry trajectory optimization, addressing problems with multiple conflicting criteria that cannot be combined into a single objective. However, MOO methods face challenges, including non-unique optimal solutions and subjective decision-making. Additionally, achieving both convergence and diversity—how close and spread the candidate solutions are relative to the Pareto front— is complex. Consequently, a new optimization problem may need to be formulated to address these issues and improve the performance of the MOO process.

This paper introduces a novel optimization problem that combines all the discussed performance indices into a single objective function. This enables the use of traditional gradient descent (GD) methods to solve the problem. The proposed objective function utilizes the weighted least squares (WLS) approach and is presented below.

$$\begin{array}{c}\hfill (\mathbf{\theta };\mathbf{x})=\sum _i\mathbf{w}\left(i\right){\left[\tilde{\mathbf{\theta }}\left(i\right)-{\tilde{\mathbf{\theta }}}_{target}\left(i\right)\right]}^2\end{array}$$

(13)

where w is the weight vector that includes the biasses of each performance index, and the definition of $\tilde{\mathbf{\theta }}$ follows.

Unlike MOO, the proposed Single-Objective Optimization (SOO) problem requires feature scaling due to significant differences in magnitudes among the performance indices. Specifically, Q is typically 5 to 7 orders of magnitude higher than $\Delta V$ and $R_B$. To avoid biasing the results towards minimizing the performance index with the most significant magnitude, feature scaling is employed. The normalization technique rescales the performance indices between 0 and 1, as shown in Equation (15).

$$\tilde{\mathbf{\theta }}\left(i\right)={\displaystyle \frac{\mathbf{\theta }\left(i\right)-\underline{\mathbf{\theta }}\left(i\right)}{\overline{\mathbf{\theta }}\left(i\right)-\underline{\mathbf{\theta }}\left(i\right)}}$$

(15)

where $\overline{\left(·\right)}$ and $\underline{\left(·\right)}$ represent the upper and lower bounds of property $\left(·\right)$. To perform normalization accurately, it is essential to know the range of possible values that each feature can take. However, computing this range is not trivial, and an approximate procedure can be employed, which will be explained in Section 6.

5. Approximating Aerocapture Corridor Boundaries

The primary challenge in aerocapture optimization lies in handling the nonlinear constraints introduced by $g_u$ and $g_l$ in Problem P1. Current optimization techniques necessitate computing $g_u$ and $g_l$ at each internal candidate point, significantly increasing execution time despite efficient search algorithms. This is quantitatively characterized in Section .

Multiple aerocapture corridors for Martian and Terrestrial missions with various hyperbolic approach orbits were analyzed. Figure 2 displays three examined corridors, with the upper boundary exhibiting an exponential trend across all cases, while the lower boundary differs from this trend in some scenarios. This discrepancy is imparted by the differing atmospheric structures of the planets at low altitudes. Several candidate fundamental functions were considered to rapidly map $\beta$ to $g_u$ and $g_l$, including exponential and polynomial curves.

Initial testing revealed limitations of a primary exponential curve ($n=1$ in Equation (16), leading to the adoption of a second-degree exponential fit ($n=2$), significantly enhancing predictive accuracy with error bounded by $\pm 4\%$ and $\pm 0.8\%$ for Terrestrial and Martian aerocaptures respectively. Additionally, sixth-degree polynomial curves were examined due to the potential flat nature of some boundaries. As illustrated by Figure 3, the polynomial curves resulted in excellent predictive capabilities, exceeding the former’s accuracy.

$$f\left(x\right)=\sum _{i=1}^n{a}_i{e}^{b_{i}x}$$

(16)

Exploration into the minimum number of points required for accurate fitting demonstrated that the accuracy of the fits was very sensitive to the selected fitting points. Randomly selecting seven points within the interest $\beta$ range demonstrated poor fit accuracy, as illustrated by Figure A1. However, when the fitting points were evenly distributed across the range, both fits exhibited increased accuracy with fewer fitting points. Specifically, fitting the sixth-degree polynomial with seven linearly spaced points proved to capture boundaries accurately, showcasing superior capabilities in approximating flat boundaries, resulting in approximations within $\pm 2.6\%$ and $\pm 0.6\%$ for Terrestrial and Martian maneuvers, respectively. Figure 3 illustrates the error between the actual boundary and the proposed fits, comparing the second-degree exponential fit with the sixth-degree polynomial fit, both fitted with 7 linearly spaced points.

6. The D-ASTRO Framework

The novelty of D-ASTRO lies in its pre-optimization phase, where complex physics are efficiently analyzed to enhance the robustness and speed of the optimization. The key aspects of this stage are detailed in the following sections.

6.1. Computation of Normalization Quantities

As discussed in Section 4.3, computing the extremal values of the performance metrics in $\mathbf{\theta }$ is crucial for casting the optimization problem into SOO form. We investigate which $\gamma$ and $\beta$ values may achieve these targets.

For volumetric constraints, the bounds are easily computed based on $\beta$, with maximum and minimum volumes corresponding to the lower and upper limits of $\beta$, respectively. Thermal metrics, such as ${\dot{q}}_{max}$, increase with deeper atmospheric entries, so ${\dot{q}}_{max}$ reaches its maximum at the lower corridor boundary and upper $\beta$ limit, while the minimum occurs at the upper corridor boundary and lower $\beta$ limit. The same applies to Q, which is maximized for deeper entries. Determining the bounds of $\Delta V$ is more complex, as post-atmospheric trajectory eccentricities can span from 0 to 1 within the aerocapture corridor. If circularization and plane change maneuvers are neglected, $\underline{\Delta V}\approx 0$. Hence, we assign $\underline{\Delta V}=0$, while $\overline{\Delta V}$ depends on the operational orbit and requires assessment across the $\beta$.

This analysis focuses on identifying the maximum and minimum values within the corridor. While the optimizer can explore $\beta$ and $\gamma$ values beyond those used for normalization, this discrepancy does not significantly impact the process. Normalization ensures that values remain on a similar scale, with values exceeding the estimated maximum normalized above 1 and smaller values becoming negative. However, the non-linear constraints $g_l$ and $g_u$ enforce candidate points within the corridor, mitigating this issue.

6.2. Computation of Polynomial Fits

A key step is generating polynomial fits to approximate the corridor boundaries. As described in Section 5, this requires computing $g_l$ and $g_u$ at seven evenly spaced $\beta$ values to capture boundary dynamics accurately. This allows the normalization metric search to integrate with the boundary value computations for curve fitting. The complete D-ASTRO framework is detailed in the next section.

6.3. The Framework

A user-specified $\beta$ range is sampled with seven linearly spaced points, including the upper and lower limits. From this analysis, values for $\overline{{\dot{q}}_{max}},\underline{{\dot{q}}_{max}},\overline{Q},\underline{Q},\overline{R_B},\underline{R_B}$ are computed, while the remaining points focus on $\overline{\Delta V}$ and corridor boundaries $g_u$ and $g_l$, with $\underline{\Delta V}$ set to zero. A binary search determines $\gamma$’s lower and upper boundary values at each $\beta$, which are then fitted with a sixth-order polynomial. After computing normalization quantities and polynomial fits, MATLAB’s fmincon optimizer is used, employing sequential quadratic programming with an active set strategy, as described by Gill et al. [19]. The robust aerocapture corridor and polynomial fits are shown in Figure 4, with the normalization values listed in Table 1.

7. Results and Analysis

A Martian aerocapture mission is considered to assess the robustness and convergence properties of the proposed optimization strategy. The Martian atmosphere is modeled using data from Mars Global Reference Atmospheric Model 2010 (Mars GRAM 2010) [16] (an engineering-level atmospheric tool developed for mission design applications). D-ASTRO requires multiple input parameters, including spacecraft preliminary properties, interplanetary insertion trajectory, and target operational orbit. These input quantities are provided in Table 2 and Table 3.

7.1. Dependencies on Solver Initial Conditions

The optimization algorithm was subjected to various initial conditions to assess its convergence properties. A bias vector $\mathbf{w}$ assigning equal importance to all performance indices (denoted as $\mathbf{w}=\left[1,1,1,1\right]$) was employed for this examination. The range of considered initial conditions comprised extreme and random values both within and outside the aerocapture corridor, with some having limited physical meaning. Table 4 provides detailed information on the initial conditions used (${\beta }_0$ and ${\gamma }_0$).

In Figure 5, level lines of the cost function are depicted within and around the aerocapture corridor. The region delineated by magenta solid lines represents the aerocapture corridor, defined by the novel polynomial fit approach discussed in Section 6. The figure illustrates that the novel objective function formulation is well-suited for gradient descent (GD) algorithms, as it has no visible sharp discontinuities. Notably, areas where aerocapture fails, such as surface collisions in the lower left corner of the figure, incur high costs, expediting solution convergence.

Based on Table 4 and Figure 5, the problem’s minimum is near $\left[13.20,-9.{76}^{\circ }\right]$, though further investigation is needed to resolve discrepancies in the suggested optima and validate trajectory feasibility. While Case 2 indicates a slightly different optimal point (explored further in Section 7.2, the algorithm consistently converges to the global minimum, even with initial conditions lacking physical significance (e.g., Cases 1, 3, and 4). The strong convergence properties of D-ASTRO are further demonstrated by the similarity in trajectory profiles shown in Figure 6.

7.2. Minimum Fuel Mission Design

The adoption of an aerocapture maneuver is often motivated by a desire to minimize fuel requirements to achieve a desired operational orbit. The minimum fuel trajectory for the proposed case study can be identified by selecting $\mathbf{w}=[0,1,0,0]$. The same initial conditions as the previous section are studied to assess convergence properties. Additionally, various target semi-major axis values, specifically $3500km$, $4621km$, and $5621km$, are analyzed to explore the behavior of the optimality valley introduced earlier.

As illustrated by the $\Delta V$ (${\mathrm{kms}}^{-1}$) cost contours in Figure 7 and the information in Table 5, this optimization problem yields optima that are scattered over a broader range. However, it is crucial not to interpret this as a sign of poor convergence of D-ASTRO’s optimizer. Instead, understanding the underlying physics and formulation of the optimization problem is essential to justify this apparent discrepancy.

Figure 7 illustrates the presence of a low-cost valley in the $\Delta V$ cost function, directing the global minimum along a specific direction. This flat valley is the reason behind the discrepancies in this scenario, with some minor discrepancies in the overall optimal trajectory. As shown in Figure 5 and Table 4, the optimum for Case 2 slightly deviates from other test cases but results in identical cost (accurate to four decimal places) as Cases 1 and 5. Although the low $\Delta V$ valley affects the minimum fuel mission design more, its impact on the optimal design with multiple metrics is minimal, ensuring good local agreement in the suggested optima.

The valley-like structure of the $\Delta V$ cost function characterizes a typical aerocapture corridor, with the trough’s vertical position determined by the target a. Smaller values require deeper atmospheric entries, shifting the valley downward, while larger values lead to shallower entries, moving it upward. The valley’s convex structure arises because one side requires increasing orbital energy due to excessive energy depletion, while the other side requires decreasing energy due to insufficient depletion. The valley in $\Delta V$ does not span the entire $\beta$ range; the increase in cost at low $\beta$ is linked to plane change maneuvers. Figure 8 illustrated the $\Delta V$ contours without considering these maneuvers, revealing the anticipated valley that spans the entire $\beta$ range. The evolution of plane changes $\Delta V$$(\Delta V_{pc})$ result from the forces acting on the spacecraft during atmospheric entry, which do not align with the orbital plane, thereby altering the inclination of the orbit’s plane.

To enhance convergence and mitigate the impact on other cases, the cost function for $\Delta V$ must be revised. Comparing $\Delta V$ contours in Figure 8 and Figure 7a shows that including $\Delta V_{pc}$ significantly increases fuel requirement. Applying feature standardization to the total $\Delta V=\Delta Vpc+\Delta V_{Hohmann}$ flattens the valley as $\Delta V_{pc}>>\Delta V_{Hohmann}$, causing $\Delta V_{pc}$ to dominate the problem. To address this, a potential solution is to decompose the $\Delta V$ into two components and apply separate standardizations.

Figure 9 depicts the optimal trajectory profiles. The extensive variation in both $\beta$ and $\gamma$ results in a multitude of trajectories. It is crucial to emphasize that despite discrepancies in the h-v profiles shown in Figure 9a, they all converge to a similar velocity (V) at the atmospheric exit, thereby leading to approximately comparable values of $\Delta V$. As anticipated, thermal quantities exhibit significant differences across all trajectories, reflecting the lack of assigned importance to these features in the optimization process.

7.3. Minimum Heat Load Mission Design

Ensuring mission success during re-entry requires equipping the spacecraft with a suitably robust TPS. By minimizing the heat load experienced during the atmospheric trajectory, the TPS can be made more lightweight, allowing for an increased payload mass fraction.

As outlined by Equation (10), calculating stagnation point heat loads requires complete trajectory knowledge. In practice, researchers often use Equation (17) for simplified heat load estimations during Entry, Descent, and Landing (EDL), based only on parameters like $\beta$ and $\gamma$ at AI. Utilizing the Sutton-Graves model for ${\dot{q}}_{conv}$ and an exponential atmospheric model, Allen and Eggers [20] formulated the complete expression, suggesting that heat load can be reduced by minimizing $\beta$ and maximizing $sin\gamma$. This approach is later explored as a substitute for the current Q cost function for aerocapture missions.

$$Q\propto \sqrt{{\displaystyle \frac{\beta }{sin\gamma }}}$$

(17)

Figure 10 illustrates the exploration of multiple initial conditions aimed at determining the optimal trajectory with the minimum heat load, revealing partially anticipated trends. Table 6 details the outcomes of the explored cases. All cases suggest $\beta$ values near the lower limit of feasible $\beta$$\left(3kg{m}^{-2}\right)$, consistent with Equation (17). However, only Case 1 aligns with the minimum of Q predicted by Equation (17), with $\gamma$ at the lower boundary of feasibility, while other cases show $\gamma$ near the upper boundary.

To investigate the disparity between the anticipated minimum of Equation (17) and the complete trajectory analysis, both heat loads are normalized with respect to their corresponding minima and maxima and illustrated in Figure 11. This figure reveals differences in the contours of Q, particularly in the slopes with respect to $\gamma$$(\partial Q/\partial \gamma$). The approximate expression in Equation (17) incorrectly suggests that increments on Q with $\gamma$$(\partial Q/\partial \gamma >0)$, whereas the complete trajectory analysis shows the opposite trend $(\partial Q/\partial \gamma <0)$. This discrepancy stems from the assumption of an EDL trajectory made by Allen and Eggers [20].

Optimizing for minimum Q highlighted the importance of the optimizer’s step tolerance in D-ASTRO robustness. Initially, fmincon’s default convergence criteria yielded reliable and consistent results, as demonstrated in Figure 5and Figure 7, and later shown in Figure 13. However, minimizing Q initially showed significant variations in ${\gamma }_f$ due to the flat nature of the normalized Q ($\tilde{Q}$) in the $\gamma$ direction, with $\tilde{Q}$ varying only slightly at the leftmost region of the corridor, of the order of $\mathcal{O}\left({10}^{-4}\right)$. This was mitigated by reducing the function and step tolerance of the solver, causing D-ASTRO to identify the global minimum consistently, and in Case 1, a local minimum coincident with the feasibility boundary (see Figure A4). However, local minima can be avoided by selecting a reasonable initial condition. The decrement tolerance only benefits the convergence properties of the minimal Q trajectory, hindering the other problems with no changes in optimal trajectory for the other test cases. The resulting optimal trajectory profiles are shown in Figure 12, where the difference between the global and local minima trajectories can be seen.

7.4. Minimum Peak Heat transfer trajectory

Similar to the parameter Q discussed earlier, the peak heat transfer rate ${\dot{q}}_{max}$ plays a crucial role in determining the design specifications of the TPS. This metric undergoes an independent optimization process to validate the underlying physics principles and the cost function employed.

The minimization of peak heat transfer demonstrates strong convergence properties, as shown in Figure 13 and Table 7, where multiple initial points converge to the same optima. The variation in convergence quality between Q and ${\dot{q}}_{max}$ arises from the latter’s greater sensitivity to $\gamma$, illustrated by curved contours in Figure 13 versus the quasi-vertical contours of Figure 10. Figure 14 shows indistinguishable trajectory profiles, highlighting the optimizer’s excellent convergence

The results presented in Table 7 align with established expectations that smaller values of $\beta$ lead to reduced ${\dot{q}}_{max}$. However, the practical implementation of this solution is limited, as $\beta =3.00kg{m}^{-2}$ corresponds to an unfeasibly small value for a realistic mission as unprecedented aeroshell surface areas (A) would be required to achieve this $\beta$. However, by introducing volumetric considerations in the optimization (non-zero first entry of $\mathbf{w}$) will lead to more realistic $\beta$. The resulting optimal points of combining thermal loads and volumetric considerations are illustrated in Table 9, where $\mathbf{w}=\left[1,0,1,0\right]$ and $\left[5,0,1,0\right]$ are considered.

7.5. Direct Comparison with MOO

This section compares D-ASTRO with a basic MOO algorithm using MATLAB’s fgoalattain optimizer. Multiple Initial conditions are again tested to evaluate the solver’s convergence dynamics and stability. Figure 15 shows the convergence dynamics of both approaches, revealing significant disagreement in optimal solutions. Table 8 presents additional comparison, incorporating data from Table 4.

The subjective decision-making bias of MOO strategies is disregarded as the primary cause behind the disagreement between D-ASTRO and MOO results. A wide range of initial conditions consistently converge to one of two possible optima, corresponding to a local and global minimum. This is evident in Table 8, which shows reductions in $\Delta V$, Q, and ${\dot{q}}_{max}$ without compromising other performance metrics for Cases 4 and 5 when compared to other, similar to the local optima in the Q cost function (Figure A4). It is also worth noting that the optimization problem may not be ill-conditioned despite the absence of metric regularization in the MOO strategy, as the MOO optimization conducted with and without D-ASTRO’s normalization technique yielded identical results.

Further analysis revealed that the MOO strategy prioritized the optimization of Q and ${\dot{q}}_{max}$ over $\Delta V$ and $R_B$. As the MOO strategy suggested identical results for normalized and raw metrics, the discrepancy does not stem from dissimilar metric magnitudes but from the individual convergence dynamics of the metrics. As illustrated by Figure 10 and Figure 13, the heat load and peak heating rate cost functions exhibit a strong convex behavior in the $\beta$ direction, facilitating rapid convergence of GD algorithms to the leftmost regions. In contrast, the volumetric consideration cost function experiences an inverse proportionality to $\beta$, leading to rapid convergence to the rightmost regions. Finally, the $\Delta V$ cost function experiences a convex-like behavior in the $\gamma$ direction, causing the optimality valley explored in Section 7.2.

The dynamics of these metrics indicate three convergence behaviors: one toward the left region of the corridor, another toward the rightmost regions, and a third toward the corridor’s center. We conclude that the MOO optimizer is attempting to globally optimize Q and ${\dot{q}}_{max}$ at the expense of $R_B$. To later find the optimal $\Delta V$ in regions where Q and ${\dot{q}}_{max}$ are quasi-globally optimal, as $\Delta V$ has relatively low cost over the entire $\beta$ range.

For completeness, all optimal trajectory profiles presented in this paper are illustrated in Figure 17 to aid the comparison between the different optima.

8. Computational Performance and Considerations of D-ASTRO

Section 7 demonstrated the superior convergence properties of D-ASTRO over conventional MOO strategies. This section examines the computational performance and considerations within the D-ASTRO framework. The analysis is divided into two parts: the first addresses the computational costs and considerations of the binary search algorithm used to compute the robust aerocapture corridor; the second evaluates the performance of the optimization process.

8.1. Dependency of Corridor fits on Binary Search Tolerance

To balance accuracy and computational efficiency, a study is conducted to analyze how the accuracy of the computed corridor boundaries varies with different tolerance levels. To evaluate this, we compute the robust coarse aerocapture corridor depicted in Figure 4 using varying binary search tolerances. The overall accuracy of the resulting polynomial fits is assessed by calculating the Euclidean L-2 norm between actual corridor boundaries, computed from a finely discretized aerocapture corridor (570 points), and those computed using the fits. Figure A5 illustrates the evolution of polynomial fit accuracy for both lower and upper corridor boundaries over the considered $\beta$ range. Additionally, Figure 18 presents the cumulative accuracy via the Euclidean L-2 norm approach and the execution time of the corridor, examining Martian and Terrestrial aerocapture maneuvers. This figure shows that the Euclidean norm reaches a plateau at a tolerance of approximately $2\times {10}^{-5}$ in both scenarios, indicating that more accurate fits are not achievable with a further decrement in tolerance. Hence, for expedited solutions of D-ASTRO, a binary search tolerance of $2\times {10}^{-5}$ is selected as the optimal balance between computational efficiency and accuracy. For the selected tolerance level, 31 trajectory simulations are required to compute $g_l$ and $g_u$, with the number of simulations being unaffected by the chosen $\beta$.

8.2. Execution Time of Optimizers

To evaluate the optimizers’ performance, 10 iterations of the optimization process were conducted using a consistent initial guess (${\beta }_0=40.60kg{m}^{-2}$ and ${\gamma }_0=-10.113\circ$) to minimize execution timing outliers. Various weightings ($\mathbf{w}$) were also considered to further explore the convergence properties of MOO and D-ASTRO, highlighting the challenges of the multi-objective approach in optimizing these problems. Table 9 presents the execution time, number of function calls, time per function call, and the suggested optima for both methods. This table indicates a significant dependency of execution time on the weight vector and optimization strategy employed. However, the time required per function call is roughly consistent across all cases, with the MOO function requiring marginally less time than D-ASTRO’s. This difference arises from the distinct optimizers used, as fgoalattain and fmincon are fundamentally different algorithms.

IIn general, D-ASTRO converges more quickly and with fewer iterations to the global optimum compared to the MOO strategy, except in the case where $\mathbf{w}=\left[1,1,1,1\right]$ case. Under this condition, D-ASTRO requires more time but yields a solution that reflects a more balanced trade-off among all metrics, in contrast to the MOO strategy, which prioritizes Q and ${\dot{q}}_{max}$.

Comparing rows 3 and 4 reveals that different weightings do not impact the solution provided by MOO, whereas they clearly affect D-ASTRO’s solution. This behavior was consistent across multiple initial guesses and weightings. The only modification that enabled the MOO solution to depend on the weight vector was performing feature regularization. These findings underscore the superior capabilities of the D-ASTRO strategy over the multi-objective approach.

Table 9. Computational Performance and Suggested Optima of D-ASTRO and MOO strategies with

${\beta }_0=40.60kg{m}^{-2}$ and ${\gamma }_0=-10.113\circ$ for Martian test case.

weights $\left(\mathit{w}\right)$	Strategy	Execution time (s)	Function counts	Time per function call (ms)	${\mathit{\beta }}_{\mathit{f}}$	${\mathit{\gamma }}_{\mathit{f}}$${(}^{\circ })$
$\left[1,1,1,1\right]$	D-ASTRO	1.889	1820	10.381	13.23	-9.765
	MOO	0.381	360	9.979	3.03	-9.638
$\left[1,1,0,0\right]$	D-ASTRO	2.512	2010	12.496	43.97	-10.434
	MOO	3.066	3020	10.151	40.61	-10.113
$\left[1,0,1,0\right]$	D-ASTRO	0.564	500	11.286	15.58	-9.653
	MOO	2.788	3010	9.263	12.39	-9.853
$\left[5,0,1,0\right]$	D-ASTRO	0.506	470	10.773	37.84	-10.101
	MOO	2.927	3010	9.723	12.39	-9.853
$\left[10,5,1,0\right]$	D-ASTRO	2.033	1630	12.473	41.49	-10.401
	MOO	2.767	3010	9.192	12.40	-9.163

Table 9. Computational Performance and Suggested Optima of D-ASTRO and MOO strategies with

${\beta }_0=40.60kg{m}^{-2}$ and ${\gamma }_0=-10.113\circ$ for Martian test case.

weights $\left(\mathit{w}\right)$	Strategy	Execution time (s)	Function counts	Time per function call (ms)	${\mathit{\beta }}_{\mathit{f}}$	${\mathit{\gamma }}_{\mathit{f}}$${(}^{\circ })$
$\left[1,1,1,1\right]$	D-ASTRO	1.889	1820	10.381	13.23	-9.765
	MOO	0.381	360	9.979	3.03	-9.638
$\left[1,1,0,0\right]$	D-ASTRO	2.512	2010	12.496	43.97	-10.434
	MOO	3.066	3020	10.151	40.61	-10.113
$\left[1,0,1,0\right]$	D-ASTRO	0.564	500	11.286	15.58	-9.653
	MOO	2.788	3010	9.263	12.39	-9.853
$\left[5,0,1,0\right]$	D-ASTRO	0.506	470	10.773	37.84	-10.101
	MOO	2.927	3010	9.723	12.39	-9.853
$\left[10,5,1,0\right]$	D-ASTRO	2.033	1630	12.473	41.49	-10.401
	MOO	2.767	3010	9.192	12.40	-9.163

As described in Section 6, determining the robust aerocapture corridor is essential to compute the scaling quantities required for normalizing the performance metrics and performing the single-objective optimization. Therefore, the time required to compute the robust corridor must be added to D-ASTRO’s optimizer execution times listed in Table 9 to obtain the total execution time. From Figure 18a and recalling that a binary search tolerance of $2\times {10}^{-5}$ is used, the overhead time required for computing the robust corridor is approximately 2.2 seconds.

To quantify the computational benefits of using polynomial fits over a binary search for computing point feasibility, D-ASTRO and MOO strategies are executed with a binary search algorithm instead of polynomial fits to evaluate $g_l$ and $g_u$. Table 10 performs the same optimization as Table 9, including this modification. The significant difference in execution times is evident, with 1 to 2 orders of magnitude increments. Although there is a slight disagreement in the optima suggested by D-ASTRO when using polynomial fits versus a binary search, this discrepancy results in minimal changes to performance metrics. Table 11 shows the percentage difference for the first three rows. This slight difference is likely introduced by the accuracy of the constraint gradients employed by fmincon, as slight variations in this quantity can lead to different optimization paths.

9. Conclusions

This study presents a novel framework for the Determinations of Aerocapture Successful Trajectories and Robust Optimization (D-ASTRO) for the rapid optimization of aerocapture-capable spacecraft geometry and trans-atmospheric trajectory. D-ASTRO robustly identifies globally optimal trajectories by optimizing the flight path angle ($\gamma$) at the atmospheric interface and ballistic coefficient ($\beta$) for specific target orbits and design requirements. This approach aims to assist mission designers in streamlining the architecture of aerocapture-capable spacecraft, facilitating the integration of aerocapture maneuvers into deep space exploration missions.

The proposed algorithm employs sixth-order polynomials to define the aerocapture corridor boundaries precisely, obviating the need for a double binary search at each candidate point to assess feasibility. This method substantially enhances computational efficiency and reduces execution time. While this approach necessitates a higher number of function evaluations to reach convergence, the execution time per function call is markedly reduced, typically around $10ms$ per evaluation. In contrast, implementing a double binary search at each candidate point results in an execution time per function call ranging from 300 to $400ms$.

Departing from prevalent multi-objective optimization (MOO) paradigms, the D-ASTRO algorithm introduces an innovative single-objective constrained optimization approach, demonstrating clear superiority through direct comparison. By recognizing that performance metrics are bounded within aerocapture constraints, D-ASTRO incorporates feature normalization, effectively integrating knowledge of metric ranges into the optimization framework. This method accelerates gradient descent optimizers and enhances convergence to global optima. Furthermore, this modification allows for a representative trade-off analysis, where the weight vectors appropriately impose the corresponding impact of the design driver on cost.

A Martian aerocapture test mission showcases D-ASTRO’s capabilities and validates the feasibility of the proposed solutions. D-ASTRO successfully determines the optimal geometries and trajectories for various mission criteria. Optimal points attained through D-ASTRO exhibit lower costs relative to multi-objective optimal trajectories, primarily driven by reduced fuel requirements and volumetric considerations. However, this is achieved at the expense of augmented heat loads and peak heat transfers for the overall optimal trajectory. The algorithm consistently converges towards global minima, demonstrating enhanced convergence properties and numerical stability even under challenging initial conditions. This underscores its effectiveness in addressing complex optimization challenges within aerocapture trajectory design.

While the analysis conducted in this study has provided valuable insights, it is essential to acknowledge certain limitations and outline potential directions for future research. One significant limitation lies in the assumption of a fixed lift-to-drag ratio (L/D) of 0.2, tailored for Martian re-entry missions. To enhance the algorithm’s applicability, future iterations should explore the optimization of L/D as a variable, thereby broadening its utility. Additionally, the assumed knowledge of the spacecraft’s state at orbital insertion presents a challenge, as these parameters are intricately linked to the design of the interplanetary insertion trajectory. Future research endeavors may involve coupling D-ASTRO with an interplanetary trajectory design tool to address this. This coupling would enable iterative refinement, ensuring coherence between the optimization of aerocapture trajectories and the design of the broader mission architecture. Through such enhancements, the algorithm can evolve to better meet the demands of complex aerocapture missions in the dynamic realm of space exploration.