Rapid Construction of Aerocapture Attainability Sets Using Sequential Convex Programming

: This paper proposes a novel method to efﬁciently construct the aerocapture attainability set based on convex optimization. Using dimensionality reduction, constructing the attainability set is equivalent to solving a set of discrete points on its boundary. As solving each of the boundary points is a typical nonlinear optimal control problem, a sequential convex programming method is adopted. The efﬁciency and accuracy of the proposed algorithm is demonstrated by high-ﬁdelity numerical results. This is the ﬁrst time that the conﬁguration of the aerocapture attainability set is precisely described by the state variables at atmospheric exit. Since the quantiﬁcation of the set is signiﬁcant for assessing the feasibility of performing an aerocapture maneuver, the proposed method can be used as a reliable tool for systematic design for aerocapture mission.


Introduction
For target planets with an atmosphere, such as Venus, Mars, Titan, and Neptune [1][2][3], aerocapture is an orbital maneuver that captures an interplanetary probe from a hyperbolic orbit or interplanetary approach orbit and places it in a mission orbit around the target planet through planetary atmospheric flight.Compared with traditional propulsive capture, aerocapture can significantly reduce the required fuel consumption [4], which attracts a lot of research attention [5][6][7][8].
The high deceleration efficiency of aerocapture maneuver renders a wide range of potential applications.Typical applications are shown schematically in Figure 1.An aerocapture maneuver can capture a vehicle in a planetary orbit or deliver the vehicle into the periodic orbits around the Sun-planet libration points, such as collinear libration points 1 and 2 (L1 and L2) [9][10][11][12][13][14][15][16].In addition, an aerocapture maneuver can adjust the phase and flight velocity of an interplanetary probe to fly along an ideal trajectory.Such an aerocapture mode is called aerogravity assist (AGA) [17][18][19].Note that in both cases, the target orbit is eventually reached after aerocapture through an intermediate transfer orbit, such as the stable manifold of a halo orbit or an interplanetary transfer orbit.These transfer orbits are determined directly by the terminal state of the atmospheric flight.Therefore, accurately constructing the aerocapture attainability sets is crucial in assessing the vehicle's reachable ability after aerocapture.
An aerocapture attainability set describes the maneuverability of the aerocapture and determines the feasibility of a given aerocapture mission.Constructing the aerocapture attainability set is equivalent to determining the boundary of the attainability set, and numerically solving the aerocapture trajectory for each terminal state on the boundary is a family of optimal control problems.Therefore, a robust and efficient algorithm of aerocapture trajectory planning is key to accurately construct and formulate the complete aerocapture attainability set.To solve the aerocapture trajectory planning problem, many guidance algorithms have been proposed since the 1980s, such as the Apollo entry guidance scheme [20], analytic aerocapture guidance [21][22][23][24], and numerical predictive-corrector aerocapture guidance [25][26][27][28][29].They are reliable but not globally optimal.At the same time, open-loop aerocapture trajectory optimization methods were also investigated [30][31][32][33].In comparison with the guidance algorithms, optimization methods ensure the optimal characteristic of aerocapture trajectories.The aerocapture trajectory can be optimized by an indirect or direct method, such as the direct collocation pseudospectral method [34].However, the open-loop optimization is slow in computation speed and has a worse convergence [35,36].Consequently, existing aerocapture trajectory planning methods, both the guidance algorithms and the open-loop optimization methods, are not suitable for the rapid and accurate construction of an aerocapture attainability set.
This paper develops a robust and efficient algorithm based on convex optimization.Convex optimization is a popular method for rapid planning and has been used widely in aerospace in recent years.Typical applications include planetary powered descent guidance [37,38], space rendezvous and proximity operations [39], planetary entry trajectory planning [40], and path planning for zero-propellant maneuvers [41,42].If a nonlinear programming problem can be formulated to a specific form of convex optimization, such as linear programming, second-order cone programming, and semidefinite programming [43], the optimal solution of the nonlinear problem can be found rapidly and accurately [44].
In Section 2, the aerocapture attainability set is defined in the longitudinal dynamics model, and the problem is formulated as determining the boundary of the attainability set.In order to obtain a complete point set, the boundary of the attainability set is discretized as a series of points by means of dimensionality reduction as in Section 3, and these boundary points are solved by our proposed convex method in Section 4. To adopt convex optimization, the original problem, which is highly nonlinear and nonconvex, is converted into a sequence of convex subproblems.Then, the solution of the original problem is approximated successively.Section 5 presents high-fidelity numerical results to demonstrate the robustness and efficiency of the proposed method.Finally, conclusions are drawn in Section 6.

Formulation of Aerocapture Attainability Set
The main phase of an aerocapture process is the planetary atmospheric flight.For an unpowered atmospheric flight, regardless of the transversal modulation and modeling in the nonrotating planet frame, the vehicle state at any time t during the aerocapture is defined as x(t) = [r, V, γ] T , where r, V, and γ are the radial distance, velocity, and flight-path angle, respectively.The motion of the vehicle is described by the nondimensional dynamic differential equations [28]  To solve the aerocapture trajectory planning problem, many guidance algorithms have been proposed since the 1980s, such as the Apollo entry guidance scheme [20], analytic aerocapture guidance [21][22][23][24], and numerical predictive-corrector aerocapture guidance [25][26][27][28][29].They are reliable but not globally optimal.At the same time, open-loop aerocapture trajectory optimization methods were also investigated [30][31][32][33].In comparison with the guidance algorithms, optimization methods ensure the optimal characteristic of aerocapture trajectories.The aerocapture trajectory can be optimized by an indirect or direct method, such as the direct collocation pseudospectral method [34].However, the open-loop optimization is slow in computation speed and has a worse convergence [35,36].Consequently, existing aerocapture trajectory planning methods, both the guidance algorithms and the open-loop optimization methods, are not suitable for the rapid and accurate construction of an aerocapture attainability set.
This paper develops a robust and efficient algorithm based on convex optimization.Convex optimization is a popular method for rapid planning and has been used widely in aerospace in recent years.Typical applications include planetary powered descent guidance [37,38], space rendezvous and proximity operations [39], planetary entry trajectory planning [40], and path planning for zero-propellant maneuvers [41,42].If a nonlinear programming problem can be formulated to a specific form of convex optimization, such as linear programming, second-order cone programming, and semidefinite programming [43], the optimal solution of the nonlinear problem can be found rapidly and accurately [44].
In Section 2, the aerocapture attainability set is defined in the longitudinal dynamics model, and the problem is formulated as determining the boundary of the attainability set.In order to obtain a complete point set, the boundary of the attainability set is discretized as a series of points by means of dimensionality reduction as in Section 3, and these boundary points are solved by our proposed convex method in Section 4. To adopt convex optimization, the original problem, which is highly nonlinear and nonconvex, is converted into a sequence of convex subproblems.Then, the solution of the original problem is approximated successively.Section 5 presents high-fidelity numerical results to demonstrate the robustness and efficiency of the proposed method.Finally, conclusions are drawn in Section 6.

Formulation of Aerocapture Attainability Set
The main phase of an aerocapture process is the planetary atmospheric flight.For an unpowered atmospheric flight, regardless of the transversal modulation and modeling in the nonrotating planet frame, the vehicle state at any time t during the aerocapture is defined as x(t) = [r, V, γ] T , where r, V, and γ are the radial distance, velocity, and flight-path angle, respectively.The motion of the vehicle is described by the nondimensional dynamic differential equations [28] dr dt = V sin γ (1) The radial distance r is normalized by the radius of planet R M , the relative velocity V is normalized by µ M /R M , and the nondimensional time t is normalized by R 3 M /µ M , where µ M is the gravitational parameter of the planet.The bank angle σ is the only control variable.The lift and drag acceleration L and D are normalized by where m represents the mass of the vehicle, S indicates the dimensional reference area, C L and C D represent the lift coefficient and drag coefficient, respectively, and ρ = ρ 0 exp(−R M h/h s ) is the dimensional atmospheric density, in which ρ 0 is the dimensional atmospheric density at R M , h s is the atmospheric density coefficient, and h is the nondimensional altitude.
Given the initial entry state of the atmospheric flight x(t 0 ) = [r 0 , V 0 , γ 0 ], all the terminal states of the atmospheric flight form the attainability set of the aerocapture, where T is the terminal state of the i-th atmospheric flight trajectory.
Note that as the height of atmosphere h atm is usually a constant for a given planet, both r 0 and r f are constant and equal to the radius of atmosphere r atm , which is the height of atmosphere plus the radius of planet R M .The terminal states are then distinguished only by V f and γ f .Additionally, note that as an infinite small change in bank angle σ leads to a sufficiently small difference in terminal state, the aerocapture attainability set A is a simply connected region in the V f − γ f space.An illustration of an aerocapture attainability set is shown in Figure 2, where the boundary of the attainability set is denoted as O.
Appl.Sci.2022, 12, 6437 3 of 18 The radial distance r is normalized by the radius of planet RM, the relative velocity V is normalized by , and the nondimensional time t is normalized by , where µM is the gravitational parameter of the planet.The bank angle σ is the only control variable.The lift and drag acceleration L and D are normalized by where m represents the mass of the vehicle, S indicates the dimensional reference area, CL and CD represent the lift coefficient and drag coefficient, respectively, and ( ) is the dimensional atmospheric density, in which ρ0 is the dimensional atmospheric density at RM, hs is the atmospheric density coefficient, and h is the nondimensional altitude.
Given the initial entry state of the atmospheric flight , all the terminal states of the atmospheric flight form the attainability set of the aerocapture, where is the terminal state of the i-th atmospheric flight trajectory.
Note that as the height of atmosphere hatm is usually a constant for a given planet, both r0 and rf are constant and equal to the radius of atmosphere ratm, which is the height of atmosphere plus the radius of planet RM.The terminal states are then distinguished only by f V and f γ .Additionally, note that as an infinite small change in bank angle σ leads to a sufficiently small difference in terminal state, the aerocapture attainability set  is a simply connected region in the f f V γ − space.An illustration of an aerocapture attainability set is shown in Figure 2, where the boundary of the attainability set is denoted as  .In this paper, the attainability set consists of only the terminal states of the vehicle exiting the atmosphere after aerocapture.If the vehicle drops into the atmosphere and eventually collides with the planet, the terminal state is excluded from the attainability set because the aerocapture failed.The demarcation of a successful and a failed aerocapture is γ f = 0, which is delineated by the line D-E in Figure 3a.Given that the target orbit after aerocapture is far from the atmosphere, a minimum altitude of periapsis h a,min > h atm is also constrained.Given the terminal state x(t f ) of atmospheric flight, the constrain can be derived from the general energy equation of the elliptical orbit [28], where r a,min = h a,min + R M is the radial distance of periapsis.This constraint imposes a lower boundary O L on A, as represented by the curve B-C in Figure 3b.As a summary, the boundary of the aerocapture attainability set O L is the closed curve A-B-C as in Figure 3b.
Appl.Sci.2022, 12, 6437 4 of 18 In this paper, the attainability set consists of only the terminal states of the vehicle exiting the atmosphere after aerocapture.If the vehicle drops into the atmosphere and eventually collides with the planet, the terminal state is excluded from the attainability set because the aerocapture failed.The demarcation of a successful and a failed aerocapture is 0 f γ = , which is delineated by the line D-E in Figure 3a.Given that the target orbit after aerocapture is far from the atmosphere, a minimum altitude of periapsis where ,min

Algorithm to Construct Aerocapture Attainability Set
Based on the formulation in Section 2, the aerocapture attainability set  can be constructed by determining its boundary  .In this paper, the boundary is calculated by means of dimensionality reduction, as illustrated in Figure 4.The key idea is to determine the reachable bounds of f γ for each fixed f V .For a specific k f V , the line subject to the given initial state x0, where the state vector x = [r, V, γ] T and the control variable σ follow the dynamic differential equations given in Equations ( 2)-( 4).For convenience, the two problems are denoted as Problem 1.
The algorithm flowchart for numerically solving the boundary  is summarized in Figure 5. Firstly, an initial terminal velocity 1 f V is chosen where both the problems of min

Algorithm to Construct Aerocapture Attainability Set
Based on the formulation in Section 2, the aerocapture attainability set A can be constructed by determining its boundary O.In this paper, the boundary is calculated by means of dimensionality reduction, as illustrated in Figure 4.The key idea is to determine the reachable bounds of γ f for each fixed V f .For a specific V k f , the line V = V k f intersects the boundary O at two points, denoted as γ k f ,min and γ k f ,max .These two points are the solutions of the optimal control problems min (7) subject to the given initial state x 0 , where the state vector x = [r, V, γ] T and the control variable σ follow the dynamic differential equations given in Equations ( 2)-( 4).For convenience, the two problems are denoted as Problem 1.
where R a r and R p r are the apoapsis and periapsis radii of the estimated terminal orbit.Then, the fixed terminal velocity is monotonically increased, The algorithm flowchart for numerically solving the boundary O is summarized in Figure 5. Firstly, an initial terminal velocity V 1 f is chosen where both the problems of min γ f and max γ f have solutions.Without loss of generality, V 1 f can be set in accordance with the energy equation [45-47]

With the stopping criterion
where r R a and r R p are the apoapsis and periapsis radii of the estimated terminal orbit.
Appl.Sci.2022, 12, 6437 5 of 18 where R a r and R p r are the apoapsis and periapsis radii of the estimated terminal orbit.Then, the fixed terminal velocity is monotonically increased, With the stopping criterion Then, the fixed terminal velocity is monotonically increased, where ∆γ f > 0 is a constant close to zero and the two reachable bounds of γ f eventually shrink into a common point A.
Next, the terminal velocity is monotonically decreased, f .As the two reachable bounds of γ f do not shrink to a common point, f are adopted as the stopping criterions, where γ B f and γ C f are the flight-path angles at points B and C determined by Equation ( 6), respectively.
Note that the choice of step ∆V k f is crucial for the numerical determination of the boundary O.In order to solve γ k f ,min and γ k f ,max precisely near the point A in the case of increasing V f , and to avoid the vehicle crashing on the planet surface in the case of decreasing V f , the step size ∆V f should be reduced adaptively.

Solving with Sequential Convex Programming
Solving Problem 1 is an essential step in the algorithm, which is marked in Figure 5 with shaded blocks.Key procedures in solving Problem 1 such as convexification, discretization, and successive iteration are introduced in this section.
As the dynamic equations are nonlinear, Problem 1 is a highly nonlinear optimal control problem.To solve it using convex optimization, all the nonconvex items in Problem 1 need to be first converted into convex functions, resulting in a sequence of convex subproblems to approximate the solution of the original Problem 1.Such a method is called sequential convex programming.Note that the subproblems are not equivalent to the original nonlinear problem.As a result, successive approximation is then used until the solution of the convex subproblems is sufficiently close to the solution of the original problem.

Change of the Independent Variable
Because the boundary points are solved for each fixed terminal velocity V f , it is convenient to use a new independent variable directly related to the velocity V to replace the time variable t in the original dynamic differential equations.In this paper, the nondimensional energy e = 1/r − V 2 /2 is chosen as the new independent variable.According to [48], de/dt = DV > 0, indicating that e increases monotonically.Then, the nondimensional velocity can be computed inversely by V = √ 2/r − 2e.During the solving process, it can be approximated by V Z+1 = √ 2/r Z − 2e in the (z + 1)-th iteration of the successive approximation.
Replacing the time t by the new independent variable e, the original dynamic differential Equations ( 1)-( 3) become dr de = sin γ/D (10) Note that using e as the independent variable can also simplify the dynamics.One differential equation is eliminated.The state vector then becomes x = [r, γ] T .

New Control Variable
The bank angle σ is the only control variable which has lower and upper bounds on the magnitude σ min ≤ |σ| ≤ σ max (12) where σ min is greater than 0 deg and σ max must be less than 180 deg, because flying at a bank angle of 0 deg or 180 deg will cause the vehicle to have no crossrange control ability [28].The bank angle is in the form of sine in the dynamic equations, which impedes the convergence of the optimal control problem.To overcome this issue, u = cos σ is defined as the new control variable.Then the constraint becomes where u min = cos σ max and u max = cos σ min .The constraint on u is naturally convex.

Successive Linearization for Dynamic Equations
The dynamics in Equations ( 10) and ( 11) are highly nonlinear.In order to convert Problem 1 into a convex subproblem, the dynamics should be linearized.Firstly, Equations ( 10) and ( 11) are rewritten into a concise form as: .
x is the derivative of x = [r, γ] T with respects to e, the column vector and the coefficient vector of the control variable For simplicity, the argument e is dropped in the remainder of this paper.The nonlinear dynamics in Equation ( 14) can be approximated by a successive small-disturbance-based linearization method [40].Given an existing solution x z from the z-th iteration of the successive approximation process, the dynamic Equation ( 14) is linearized at x z in the (z + 1)-th iteration as .
where the coefficient matrix of the state vector a 11 a 12 a 21 a 22 (18) in which where To ensure the validity of linearization, a trust-region constraint is used, where δ is the constrained boundary and set in accordance with the magnitudes of the elements in x.

Convex Subproblem
Linearizing Equation ( 14) into Equation ( 17), the original optimal control problem can be iteratively solved by a sequence of convex subproblems, subject to r(e f ) = r atm (24) |x − which is denoted as Problem 2. The criterion of convergence is defined as where ε is a constant vector, which represents the accuracy of the solution with respect to the original Problem 1.

Discretization
To solve Problem 2 numerically, the independent variable e is equally discretized into N + 1 points in the interval [e 0 , e f ].The constraints in Equations ( 23) and ( 25) apply for each of these points.The state vector and control variable are denoted as x i and u i for the i-th point, respectively.The variable sets {x i } and {u i } are the variables to be solved.The dynamic Equation ( 22) becomes where i = 0, 1, • • • , N − 1 and ∆e = (e f − e 0 )/N.After rearrangement, Equation ( 27) is then rewritten into where I is a unit matrix with the same dimension as The discretized problem is denoted as Problem 3, which is equivalent to and has the same solution as Problem 2. The performance index is given by Equation (21).The linear equality constraints are given by Equations ( 24) and (28).The linear inequality constraint is given by Equation ( 25).

Procedure of Iteration
The successive procedure to solve Problem 1 is summarized as follows: Step 1: Set z = 0 and the state vectors x 0 i to an initial profile.In accordance with the trajectory characteristic of a normal aerocapture flight, the initial profile of radial distances is given as r 01 r 02 r 03 , as shown in Figure 6a.r 01 is a straight line from r atm to r m = 1 + h atm /2 divided evenly as [N/10] points, r 02 is a straight line from r m to r m having N + 1 − 2[N/10] points, and r 03 is a straight line from r m to r atm having [N/10] points.The symbol [ ] is the floor function mapping a real number onto the largest integer that is less than or equal to the real number.The initial flight-path angle profile is chosen as a constant value γ 0 as shown in Figure 6b.Note that choosing the initial guess is not necessary because the solution of the convex optimization problem is unique.However, a good initial profile can accelerate the convergence of the numerical procedure.Step 2: The solution of the z-th iteration z x is chosen to be the initial profile in the (z + 1)-th iteration to solve x and 1 z u + [49,50].
Step 3: Check the convergence criterion ( ) ( ) If it is satisfied, then proceed to step (4); otherwise, go to step (2) for the next iteration.
Step 4: The solution of Problem 1 is found to be x and It is worth to note that approximating by a sequence of convex optimizations as in Problem 3, although the convergence of sequential approximation to the original Problem 1 remains to be proven, a considerable number of preceding studies have provided sufficient confidence in the convergence [37][38][39][40][41][42][43].In this paper, the numerical simulations of solving the aerocapture attainability sets given in Section 5 showed good performance and reliability of the developed convex method.
In general, aiming at the convergence of the method proposed in this paper, it can mainly be improved in three aspects.First, the dynamic model reduces the state quantity by one dimension through the substitution of independent variables, thereby reducing the complexity of the problem.Secondly, the groove-shaped initial trajectory is given for the solution of boundary points for the first time calculation, which increases the feasibility of the problem.Finally, when solving other boundary points, the initial value is given by the optimal solution of the previous solution.This adjacency idea greatly improves the sensitivity of the initial value.

Numerical Demonstrations
As Mars is a hot destination for deep space exploration, the Mars aerocapture was taken as an example to verify our proposed algorithm.The radius of Mars is RM = 3396 km.The gravitational parameter of Mars is µM = 42,828 km 3 /s 2 .The Orion Multi-Purpose Crew Vehicle (MPCV) [28] was used in our simulations.The mass of the Orion MPCV is 10,387 kg.The spacecraft always flies a trim angle-of-attack profile and the hypersonic lift-drag ratio L/D is 0.27.The lift and drag coefficients are CD = 1.37 and CL = 0.37, respectively.The reference area of the vehicle is 19.86 m 2 .The initial states entering the atmosphere of Mars are listed in Table 1 [3].Step 2: The solution of the z-th iteration x z is chosen to be the initial profile in the (z + 1)-th iteration to solve x z+1 and u z+1 [49,50].
Step 3: Check the convergence criterion max 0≤i≤N If it is satisfied, then proceed to step (4); otherwise, go to step (2) for the next iteration.
Step 4: The solution of Problem 1 is found to be x * = x z+1 and u * = u z+1 .It is worth to note that approximating by a sequence of convex optimizations as in Problem 3, although the convergence of sequential approximation to the original Problem 1 remains to be proven, a considerable number of preceding studies have provided sufficient confidence in the convergence [37][38][39][40][41][42][43].In this paper, the numerical simulations of solving the aerocapture attainability sets given in Section 5 showed good performance and reliability of the developed convex method.
In general, aiming at the convergence of the method proposed in this paper, it can mainly be improved in three aspects.First, the dynamic model reduces the state quantity by one dimension through the substitution of independent variables, thereby reducing the complexity of the problem.Secondly, the groove-shaped initial trajectory is given for the solution of boundary points for the first time calculation, which increases the feasibility of the problem.Finally, when solving other boundary points, the initial value is given by the optimal solution of the previous solution.This adjacency idea greatly improves the sensitivity of the initial value.

Numerical Demonstrations
As Mars is a hot destination for deep space exploration, the Mars aerocapture was taken as an example to verify our proposed algorithm.The radius of Mars is R M = 3396 km.The gravitational parameter of Mars is µ M = 42,828 km 3 /s 2 .The Orion Multi-Purpose Crew Vehicle (MPCV) [28] was used in our simulations.The mass of the Orion MPCV is 10,387 kg.The spacecraft always flies a trim angle-of-attack profile and the hypersonic lift-drag ratio L/D is 0.27.The lift and drag coefficients are C D = 1.37 and C L = 0.37, respectively.The reference area of the vehicle is 19.86 m 2 .The initial states entering the atmosphere of Mars are listed in Table 1 [3].In the construction of the aerocapture attainability set, the lower boundary given in Equation ( 6) was applied with a minimum altitude of apoapsis h a,min = 200 km after the atmospheric flight.The terminal orbit was estimated with the size of 300 × 3000 km to initialize the terminal velocity in Equation (8).The initial adaptive step of terminal velocity was ∆ 1 f = 100 m/s.The stopping criterions were ∆γ f = ∆γ f ,min = γ k f ,max = 0.01 deg.In the solving of optimal control problems, the lower and upper bounds of the bank angle in Equation ( 12) were set to be σ min =15 deg and σ max = 165 deg, respectively.For the case of D ≤ 0.1 g M , the states were integrated directly with zero bank angle because of the weak modulation of the vehicle in the thin atmosphere.The trust-region boundary in Equation (25) was δ = [6000/R M , π/4] T .The stopping criterion of successive approximation in Equation ( 29) was ε = [200/R M , 0.05π/180] T .The problem was discretized by N = 300 and solved by the state-of-the-art primal-dual interior point method (PDIPM), which is efficient and accurate in solving convex programming problems without initial guess [51] using the mature software MOSEK [52].All computations were performed via MATLAB-R2018a on a desktop with Intel i7-4790K CPU @4.00 GHz.
In addition, various simulations with different initial states and parameters of the vehicle were also investigated to determine the impact of different parameters on the constructed aerocapture attainability set.

Geometric Structure of the Aerocapture Attainability Sets
Figure 7a shows the aerocapture attainability set for the initial state given in Table 1, plotted in the state space.The flight-path angle has a wide range for a low terminal velocity.As the terminal velocity increases, the boundary of the attainability set eventually shrinks to point A. Without a lower boundary limit given by Equation ( 6), the attainability set extends infinitely close to γ f = 0 and even to γ f ≤ 0, where the vehicle falls into the atmosphere.The corresponding orbit elements are plotted in Figure 7b.The eccentricity of the terminal orbit is calculated using Keplerian orbital mechanics [28], In the construction of the aerocapture attainability set, the lower boundary given in Equation ( 6) was applied with a minimum altitude of apoapsis ,min In the solving of optimal control problems, the lower and upper bounds of the bank angle in Equation ( 12 The problem was discretized by N = 300 and solved by the state-of-the-art primal-dual interior point method (PDIPM), which is efficient and accurate in solving convex programming problems without initial guess [51] using the mature software MOSEK [52].All computations were performed via MATLAB-R2018a on a desktop with Intel i7-4790K CPU @4.00 GHz.
In addition, various simulations with different initial states and parameters of the vehicle were also investigated to determine the impact of different parameters on the constructed aerocapture attainability set.

Geometric Structure of the Aerocapture Attainability Sets
Figure 7a shows the aerocapture attainability set for the initial state given in Table 1, plotted in the state space.The flight-path angle has a wide range for a low terminal velocity.As the terminal velocity increases, the boundary of the attainability set eventually shrinks to point A. Without a lower boundary limit given by Equation ( 6  Aerocapture attainability set for the initial state given in Table 1. Figure 7. Aerocapture attainability set for the initial state given in Table 1. The nondimensional semimajor axis is a f = 1/(2/r atm − V 2 f ), and the periapsis of the terminal orbit has a nondimensional altitude h p, f = a f (1 − e f ) − 1.As shown in the figure, the larger the terminal flight-path angle, the larger the eccentricity of the terminal orbit.Besides, the left boundary of the attainability set, the curve A-B in Figure 7a, gives a stable upper boundary of the altitude of periapsis, the curve A1-B1 in Figure 7b.
The altitude, velocity, and flight-path angle profiles are shown in Figure 8.The state profiles plotted with black and grey colors have terminal states on the left and right boundaries of the attainability set in Figure 7a, respectively, and the state profiles corresponding to the point A are highlighted with the red color.Figure 8a shows that the atmospheric flights with a terminal state on the right boundary of the attainability set have a lower flight altitude in the atmosphere than those on the left boundary.They require a larger change in flight-path angle, and the vehicle must fly at a lower altitude to obtain enough lift to increase dγ/dt.Figure 8b shows that the time of atmospheric flight for a terminal state on the left boundary of the attainability set is longer than that on the right boundary.Meanwhile, Figure 8c indicates that the flight-path angle for a terminal state on the left boundary changes more smoothly during the atmospheric ascent phase than that on the right boundary.
a lower flight altitude in the atmosphere than those on the left boundary.They require a larger change in flight-path angle, and the vehicle must fly at a lower altitude to obtain enough lift to increase / d dt γ .Figure 8b shows that the time of atmospheric flight for a terminal state on the left boundary of the attainability set is longer than that on the right boundary.Meanwhile, Figure 8c indicates that the flight-path angle for a terminal state on the left boundary changes more smoothly during the atmospheric ascent phase than that on the right boundary.
Note that, in the aforementioned computing environment, it took a total of 7.47 s to calculate the envelope of the attainability set in Figure 7. Specifically, a total of 57 extreme values of points were calculated, and the average calculation time per extreme value was about 0.131 s.That is, it took only about a hundred milliseconds.At the same time, in the process of solving all extreme values, the successive approximation was generally completed between 4~7 steps, and the calculation time of each step was about 0.02~0.03s, which is equivalent to the time of solving a convex subproblem.The statistical results of the calculation time demonstrate the high efficiency of the proposed method, and it also verifies the feasibility of the method to determine the aerocapture attainability set onboard for an aerocapture mission.Note that, in the aforementioned computing environment, it took a total of 7.47 s to calculate the envelope of the attainability set in Figure 7. Specifically, a total of 57 extreme values of points were calculated, and the average calculation time per extreme value was about 0.131 s.That is, it took only about a hundred milliseconds.At the same time, in the process of solving all extreme values, the successive approximation was generally completed between 4~7 steps, and the calculation time of each step was about 0.02~0.03s, which is equivalent to the time of solving a convex subproblem.The statistical results of the calculation time demonstrate the high efficiency of the proposed method, and it also verifies the feasibility of the method to determine the aerocapture attainability set onboard for an aerocapture mission.

Aerocapture Attainability Sets with Different Initial Velocity
Figure 9 shows the aerocapture attainability sets with different initial velocities V 0 = 5.6, 5.9 and 6.2 km/s.A large V 0 increases the range of the attainability set to a higher terminal velocity and flight-path angle.The terminal states in the subset of the attainability set, separated by the cyan line, makes the vehicle escape from the planet.This subset provides the possibility of delivering the vehicle into a periodic orbit around the Sun-planet L1 or L2 points, which renders a fuel-efficient interplanetary transfer.

Aerocapture Attainability Sets with Different Initial Flight-Path Angles
The initial flight-path angle 0 γ is also an important design parameter for aerocap-

Aerocapture Attainability Sets with Different Initial Flight-Path Angles
The initial flight-path angle γ 0 is also an important design parameter for aerocapture.For Mars aerocapture, the permissible range of γ 0 is relatively small.Three different initial flight-path angles γ 0 = −10.

Aerocapture Attainability Sets with Different Hypersonic Lift-Drag Ratios
Although the lift-drag ratios are fixed for the Orion MPCV, the analysis of the aerocapture attainability set with different hypersonic lift-drag ratios can provide some knowledge of the vehicle's structure design.In our simulations, C D is fixed as 1.37, and C L varies with different lift-drag ratios L/D.As shown in Figure 11, the hypersonic lift-drag ratio has a significant impact on the range of the aerocapture attainability set.A higher lift-drag ratio results in a larger attainability set.Additionally, the aerocapture attainability set with a high lift-drag ratio contains the ones with lower lift-drag ratios.In our simulations, A L/D=0.27 ⊂ A L/D=0.54 ⊂ A L/D=0.81 .
aerocapture attainability set with different hypersonic lift-drag ratios can provide some knowledge of the vehicle's structure design.In our simulations, CD is fixed as 1.37, and CL varies with different lift-drag ratios L/D.As shown in Figure 11, the hypersonic lift-drag ratio has a significant impact on the range of the aerocapture attainability set.A higher lift-drag ratio results in a larger attainability set.Additionally, the aerocapture attainability set with a high lift-drag ratio contains the ones with lower lift-drag ratios.In our simulations,

Aerocapture Attainability Sets with Different Upper Bounds of Bank Angle
For the optimal aerocapture process, the optimal bank angle profile has a bang-bang control structure.Hence, the bank angle switches between σmin to σmax.Thus, the bounds of the bank angle are significantly related to the crossrange control ability [28].In our simulation, the lower bound of the bank angle is fixed, and the aerocapture attainability set influence was analyzed with different upper bounds of bank angle σmax.As shown in Figure 12

Aerocapture Attainability Sets with Different Upper Bounds of Bank Angle
For the optimal aerocapture process, the optimal bank angle profile has a bang-bang control structure.Hence, the bank angle switches between σ min to σ max .Thus, the bounds of the bank angle are significantly related to the crossrange control ability [28].In our simulation, the lower bound of the bank angle is fixed, and the aerocapture attainability set influence was analyzed with different upper bounds of bank angle σ max .As shown in knowledge of the vehicle's structure design.In our simulations, CD is fixed as 1.37, and CL varies with different lift-drag ratios L/D.As shown in Figure 11, the hypersonic lift-drag ratio has a significant impact on the range of the aerocapture attainability set.A higher lift-drag ratio results in a larger attainability set.Additionally, the aerocapture attainability set with a high lift-drag ratio contains the ones with lower lift-drag ratios.In our simulations,

Aerocapture Attainability Sets with Different Upper Bounds of Bank Angle
For the optimal aerocapture process, the optimal bank angle profile has a bang-bang control structure.Hence, the bank angle switches between σmin to σmax.Thus, the bounds of the bank angle are significantly related to the crossrange control ability [28].In our simulation, the lower bound of the bank angle is fixed, and the aerocapture attainability set influence was analyzed with different upper bounds of bank angle σmax.As shown in Figure 12     The peak load factor, heating rate, and dynamic pressure together with the atmospheric flight time of aerocapture for each of the points A, B, C and C′ are given in Table 2.The bank angle profiles of the four characteristic points are shown in Figure 14.Despite the case of D ≤ 0.1 M g , wherein the states are integrated directly with a zero bank angle, the bank angle in the optimal profile switches twice, and the max σ always appears in the middle arc.In point C', a singular arc of the bank angle profile exists, that is, the bank angle profile is not exactly the bang-bang control structure used to satisfy the path constraints.The comparison of the path variable profiles between points C and C' is shown in Figure 15, which intuitively reflects the effect of suppressing the path variables.Generally, large peak values of the path variables may threaten the aerocapture vehicle; however, the results in this subsection indicate that the peak values of the path variables are tolerable for Mars aerocapture.The peak load factor, heating rate, and dynamic pressure together with the atmospheric flight time of aerocapture for each of the points A, B, C and C′ are given in Table 2.The bank angle profiles of the four characteristic points are shown in Figure 14.Despite the case of D ≤ 0.1 M g , wherein the states are integrated directly with a zero bank angle, the bank angle in the optimal profile switches twice, and the max σ always appears in the middle arc.In point C', a singular arc of the bank angle profile exists, that is, the bank angle profile is not exactly the bang-bang control structure used to satisfy the path constraints.The comparison of the path variable profiles between points C and C' is shown in Figure 15, which intuitively reflects the effect of suppressing the path variables.Generally, large peak values of the path variables may threaten the aerocapture vehicle; however, the results in this subsection indicate that the peak values of the path variables are tolerable for Mars aerocapture.

Conclusions
In this paper, an algorithm was proposed to determine the exact aerocapture attainability set on the basis of convex optimization.The efficiency and robustness of the proposed algorithm were demonstrated successfully by a variety of examples.To solve the aerocapture attainability set, its boundary was converted into a set of points by means of dimensionality reduction.Finding each element in the point set was an optimal control problem and was achieved by using the convex method.
The simulation showed the aerocapture attainability set was affected by the following parameters: the initial velocity, the initial flight-path angle, the lift-drag ratio, and the maximum magnitude of the bank angle.Moreover, among these parameters, L/D had the greatest effect on the aerocapture attainability set; a larger 0 V can expand the aerocapture attainability and increase the capability of the aerocapture flight, whereas a smaller σmax, such as less than 90 deg, can greatly weaken the longitudinal control ability of the vehicle during the atmospheric flight.Moreover, the influence of path constraints on aerocapture attainability sets was simulated.It revealed that the existence of path constraints reduces significantly the range of aerocapture attainability sets.Numerical simulations showed the proposed method works well during the solving process.Because of the relatively low computational cost, the proposed algorithm can potentially be used onboard to rapidly assess the aerocapture attainability set.
Aiming at the problem of solving the general aerocapture attainability set, in order to further improve the online computing capability, future work should focus on how to analytically characterize the structure of the attainability set and how to switch a pointby-point solution to a characterization parameter calculation.

Conclusions
In this paper, an algorithm was proposed to determine the exact aerocapture attainability set on the basis of convex optimization.The efficiency and robustness of the proposed algorithm were demonstrated successfully by a variety of examples.To solve the aerocapture attainability set, its boundary was converted into a set of points by means of dimensionality reduction.Finding each element in the point set was an optimal control problem and was achieved by using the convex method.
The simulation showed the aerocapture attainability set was affected by the following parameters: the initial velocity, the initial flight-path angle, the lift-drag ratio, and the maximum magnitude of the bank angle.Moreover, among these parameters, L/D had the greatest effect on the aerocapture attainability set; a larger V 0 can expand the aerocapture attainability and increase the capability of the aerocapture flight, whereas a smaller σ max , such as less than 90 deg, can greatly weaken the longitudinal control ability of the vehicle during the atmospheric flight.Moreover, the influence of path constraints on aerocapture attainability sets was simulated.It revealed that the existence of path constraints reduces significantly the range of aerocapture attainability sets.
Numerical simulations showed the proposed method works well during the solving process.Because of the relatively low computational cost, the proposed algorithm can potentially be used onboard to rapidly assess the aerocapture attainability set.
Aiming at the problem of solving the general aerocapture attainability set, in order to further improve the online computing capability, future work should focus on how to analytically characterize the structure of the attainability set and how to switch a point-bypoint solution to a characterization parameter calculation.

Figure 2 .
Figure 2. Illustration of the aerocapture attainability set in state space.Figure 2. Illustration of the aerocapture attainability set in state space.

Figure 2 .
Figure 2. Illustration of the aerocapture attainability set in state space.Figure 2. Illustration of the aerocapture attainability set in state space.
distance of periapsis.This constraint imposes a lower boundary L  on  , as represented by the curve B-C in Figure 3b.As a summary, the boundary of the aerocapture attainability set L  is the closed curve A-B-C as in Figure 3b.

Figure 3 .
Figure 3.The boundary of aerocapture attainability set in state space.
sects the boundary  at two points, denoted as

fγ 1 fV
and max f γ have solutions.Without loss of generality, can be set in accordance with the energy equation[45][46][47]

Figure 3 .
Figure 3.The boundary of aerocapture attainability set in state space.

Figure 4 .
Figure 4. Solving the boundary of aerocapture attainability set by dimensionality reduction.

Figure 5 .
Figure 5.The algorithm flowchart for numerically solving the boundary of aerocapture attainability set.

Figure 4 .
Figure 4. Solving the boundary of aerocapture attainability set by dimensionality reduction.

Figure 4 .
Figure 4. Solving the boundary of aerocapture attainability set by dimensionality reduction.

Figure 5 .
Figure 5.The algorithm flowchart for numerically solving the boundary of aerocapture attainability set.

Figure 5 .
Figure 5.The algorithm flowchart for numerically solving the boundary of aerocapture attainability set.
03  r is a straight line from m r to atm r having [ /10] N points.The symbol [ ] is the floor function mapping a real number onto the largest integer that is less than or equal to the real number.The initial flight-path angle profile is chosen as a constant value 0 γ as shown in Figure6b.Note that choosing the initial guess is not necessary because the solution of the convex optimization problem is unique.However, a good initial profile can accelerate the convergence of the numerical procedure.

Figure 6 .
Figure 6.Illustration of the initial profile of radius distances and flight-path angle.

Figure 6 .
Figure 6.Illustration of the initial profile of radius distances and flight-path angle.

1 fΔ
after the atmospheric flight.The terminal orbit was estimated with the size of 300 × 3000 km to initialize the terminal velocity in Equation(8).The initial adaptive step of terminal velocity was = 100 m/s.The stopping criterions were f ) were set to be min σ =15 deg and max σ = 165 deg, respectively.For the case of D ≤ 0.1 gM, the states were integrated directly with zero bank angle because of the weak modulation of the vehicle in the thin atmosphere.The trust-region boundary in Equation (25) was [6000 / , / 4] Figure7ashows the aerocapture attainability set for the initial state given in Table1, plotted in the state space.The flight-path angle has a wide range for a low terminal velocity.As the terminal velocity increases, the boundary of the attainability set eventually shrinks to point A. Without a lower boundary limit given by Equation (6), the attainability set extends infinitely close to 0 f γ = and even to

Figure 7 .
Figure 7. Aerocapture attainability set for the initial state given in Table1.Figure7.Aerocapture attainability set for the initial state given in Table1.

Figure 8 .
Figure 8. Aerocapture attainability sets described by state variables and orbital elements.Figure 8. Aerocapture attainability sets described by state variables and orbital elements.

Figure 8 .
Figure 8. Aerocapture attainability sets described by state variables and orbital elements.Figure 8. Aerocapture attainability sets described by state variables and orbital elements.

Figure 9
Figure9shows the aerocapture attainability sets with different initial velocities 0 V = 5.6, 5.9 and 6.2 km/s.A large 0 V increases the range of the attainability set to a higher terminal velocity and flight-path angle.The terminal states in the subset of the attainability set, separated by the cyan line, makes the vehicle escape from the planet.This subset provides the possibility of delivering the vehicle into a periodic orbit around the Sunplanet L1 or L2 points, which renders a fuel-efficient interplanetary transfer.
ture.For Mars aerocapture, the permissible range of 0 γ is relatively small.Three differ- ent initial flight-path angles 0 γ = −10.5 deg, −11.0 deg, and −11.5 deg were considered in our simulations.The resulting attainability sets are shown in Figure 10.The larger the initial flight-path angle, the larger the aerocapture attainability set expanded to a higher terminal velocity and flight-path angle.(a) In terminal state space (b) In terminal orbital elements

Figure 10 .
Figure 10.Aerocapture attainability sets with a different initial flight-path angle.
5 deg, −11.0 deg, and −11.5 deg were considered in our simulations.The resulting attainability sets are shown in Figure 10.The larger the initial flight-path angle, the larger the aerocapture attainability set expanded to a higher terminal velocity and flight-path angle.

Figure 9
Figure9shows the aerocapture attainability sets with different initial velocities 0 V = 5.6, 5.9 and 6.2 km/s.A large 0 V increases the range of the attainability set to a higher terminal velocity and flight-path angle.The terminal states in the subset of the attainability set, separated by the cyan line, makes the vehicle escape from the planet.This subset provides the possibility of delivering the vehicle into a periodic orbit around the Sunplanet L1 or L2 points, which renders a fuel-efficient interplanetary transfer.

Figure 9 .
Figure 9. Aerocapture attainability sets with different initial velocities.5.3.Aerocapture Attainability Sets with Different Initial Flight-Path AnglesThe initial flight-path angle 0 γ is also an important design parameter for aerocap- ture.For Mars aerocapture, the permissible range of 0 γ is relatively small.Three differ-

Figure 10 .
Figure 10.Aerocapture attainability sets with a different initial flight-path angle.Figure 10.Aerocapture attainability sets with a different initial flight-path angle.

Figure 10 .
Figure 10.Aerocapture attainability sets with a different initial flight-path angle.Figure 10.Aerocapture attainability sets with a different initial flight-path angle.
In terminal state space (b) In terminal orbital elements
, the attainability sets are similar with σmax = 125 deg, 165 deg, but for σmax = 85 deg, the aerocapture attainability set significantly contracted.It shows that σmax < 90 deg greatly weakens the longitudinal control ability of the vehicle during the atmospheric flight.
(a) In terminal state space (b) In terminal orbital elements

Figure 12 .
Figure 12.Aerocapture attainability sets with different upper bounds of bank angle.
Figure 12, the attainability sets are similar with σ max = 125 deg, 165 deg, but for σ max = 85 deg, the aerocapture attainability set significantly contracted.It shows that σ max < 90 deg greatly weakens the longitudinal control ability of the vehicle during the atmospheric flight.
In terminal state space (b) In terminal orbital elements
, the attainability sets are similar with σmax = 125 deg, 165 deg, but for σmax = 85 deg, the aerocapture attainability set significantly contracted.It shows that σmax < 90 deg greatly weakens the longitudinal control ability of the vehicle during the atmospheric flight.
(a) In terminal state space (b) In terminal orbital elements

Figure 12 .
Figure 12.Aerocapture attainability sets with different upper bounds of bank angle.Figure 12. Aerocapture attainability sets with different upper bounds of bank angle.

Figure 12 .
Figure 12.Aerocapture attainability sets with different upper bounds of bank angle.Figure 12. Aerocapture attainability sets with different upper bounds of bank angle.

Figure 13 .
Figure 13.Aerocapture attainability sets with and without path constraints.

Figure 14 .
Figure 14.Bank angle profiles for each of the points A, B, C, and C′ in Figure 13.

Figure 13 .
Figure 13.Aerocapture attainability sets with and without path constraints.

Figure 14 .
Figure 14.Bank angle profiles for each of the points A, B, C, and C′ in Figure 13.Figure 14.Bank angle profiles for each of the points A, B, C, and C in Figure 13.

Figure 14 .
Figure 14.Bank angle profiles for each of the points A, B, C, and C′ in Figure 13.Figure 14.Bank angle profiles for each of the points A, B, C, and C in Figure 13.

Figure 15 .
Figure 15.Load factor, heating rate, and dynamic pressure profiles corresponding to points C and C′ in Figure 13.

Figure 15 .
Figure 15.Load factor, heating rate, and dynamic pressure profiles corresponding to points C and C in Figure 13.

Table 1 .
Initial conditions of Mars aerocapture.

Table 2 .
The peak load factor, heating rate, and dynamic pressure together with the atmospheric flight time of aerocapture for each of the points A, B, C and C′ in Figure13.

Table 2 .
The peak load factor, heating rate, and dynamic pressure together with the atmospheric flight time of aerocapture for each of the points A, B, C and C in Figure13.

Table 2 .
The peak load factor, heating rate, and dynamic pressure together with the atmospheric flight time of aerocapture for each of the points A, B, C and C′ in Figure13.