# Multiphase Trajectory Optimization of a Lunar Return Mission to an LEO Space Station

Authors: Luyi Yang; Haiyang Li; Jin Zhang; Yazhong Luo
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

Lunar exploration architecture can be made more flexible and reliable with the support of a low-Earth orbit (LEO) space station. This study therefore evaluated a proposed hybrid optimization scheme to design the entire trajectory of a reusable spacecraft starting from trans-Earth injection (EI) at the perilune and ending at an LEO space station. As such a trajectory has multiple constraints and multiple dynamical models, it is divided into the trans-Earth phase, aerocapture phase, and postatmospheric phase. The optimization scheme is performed at two levels: sublevel and top level. At the sublevel, two novel pseudo rules are proposed to optimize the trans-Earth trajectory so that it satisfies the coplanar constraints of the space station. Then, in the aerocapture phase, the bank angle is optimized to satisfy the mission constraints, and in the atmospheric phase, the one-impulsive maneuver is performed and optimized to insert the spacecraft into the target space station orbit. The multiple phases are connected to each other by boundary conditions where the terminal state of the previous phase is transformed into the initial state of the following phase. At the top level, the vacuum perigee height is selected as the mission design variable based on problem characteristics analysis and a hybrid optimization scheme is conducted to minimize the total velocity increment. The simulation results demonstrate that the proposed hybrid optimization method is effective for the design of an entire trajectory with acceptable velocity cost which is less than that in the previous study. The coplanar constraints of the space station and other mission constraints in each phase are also satisfied. Furthermore, the proposed trajectory design method is shown to be applicable to a reusable spacecraft returning to an LEO space station parked in any arbitrary orbital plane.


---

## Body

## 1. Introduction

In January 2019, the Chang’e-4 lunar probe of China became the first spacecraft to land on the far side of the Moon, marking a milestone in the history of space exploration [1]. The Moon is the nearest celestial body to the Earth, and it is an important destination for human visitation because of its abundant resources and its potential as the foundation of future deep space exploration. In recent years, lunar exploration missions have been scheduled by ISRO [2], NASA [3], and other space agencies [4] around the world. Indeed, the next step in the third phase of China’s lunar exploration program is to return samples from the Moon to the Earth [5], providing a good foundation for manned lunar missions in the future. Thus, it is important to develop a suitable mission planning method for the lunar returning mission.The flight mode of the lunar exploration architecture plays a key role in the successful implementation of the lunar mission. Wang et al. [6] proposed a mission architecture matrix to assist in the analysis and selection of the optimal flight mode. With the support of the space stations in different orbits, Peng and Yang [7] found that lunar exploration architecture would be more flexible and reliable with the support of an LEO space station. Particularly for the lunar returning mission, Murtazin [8] pointed out that instead of returning to the Earth ground, the lunar returning spacecraft could perform an aerocapture maneuver near the Earth and then dock with the space station in LEO. The resulting reusage could help to reduce the mission cost and improve the utilization of the space station. The reusable spacecraft could also be used for future missions after propellant supplement. Therefore, it is meaningful to develop a trajectory design method for travel from the Moon to an LEO space station.One of the key technical issues in reusing a lunar return spacecraft by employing an LEO space station is the Moon-to-Earth trajectory design method in the cislunar space. Li et al. [9] proposed a modified patched conic model to generate trans-Earth trajectories from high-latitude lunar landing sites. Gao et al. [10] studied two-impulse transfer trajectories from a lunar orbital station to the Earth. An iterative algorithm was combined with the patched conic approximation, and the orbital windows were also analyzed. Yang et al. [11] improved the multiconic method to design the Moon-to-Earth with good convergence and found that the velocity increment decreased with the transfer time. Another type of Moon-to-Earth trajectory is the point return orbit, in which the spacecraft directly enters the Earth’s atmosphere to land at a desired destination. Shen et al. [12] proposed a serial design method for fixed-point return orbits, beginning with an initial value design using patched conic equations and ending with a precision solution obtained using numerical perturbation models. Zheng and Zhou [13] presented an efficient iterative method to determine trans-Earth trajectories for direct reentry into the Earth’s atmosphere. Most of the previously published research work above focused on how to return a manned or unmanned spacecraft safely to Earth. In contrast, few studies proposed appropriate trajectory optimization strategies to return a reusable spacecraft to an LEO space station. It is challenging to design such trajectories because not only the reentry corridor conditions but also the coplanar constraints of LEO space station rendezvous have to be considered.In order to realize the reuse of the lunar returning spacecraft by an LEO space station, aerocapture should be employed to slow the spacecraft while reducing propellant consumption. Specifically, aerocapture refers to a maneuver in which an interplanetary spacecraft flies into a planet’s atmosphere in order to reduce its entry velocity and reach a specified orbit around the planet [14]. This kind of atmospheric trajectory optimization of the reentry spacecraft has been widely studied. Lu et al. [15] developed a two-phase numerical predictor–corrector guidance algorithm to place a spacecraft in a given-altitude orbit by employing aerocapture. Chai et al. [16] proposed a violation learning differential evolution algorithm that provided a good initial guess to the optimization atmospheric trajectory. A new three-layer-hybrid optimal control solver was designed to optimize the constrained trajectory of space maneuver vehicles [17], and a chance-constrained optimal control mode is established considering both deterministic and probabilistic constraints [18]. Sagliano [19] combined the pseudospectral and convex optimization methods to solve the optimal control problem in the descent phase of the NASA Mars Science Laboratory. It shows good advantage of high-accuracy and real-time control compared with the standard convex methods. Sagliano [20] further proposed a generalized hp method that provided a good trade-off between accuracy and computation efficiency. Chen et al. [21] established an aerocapture maneuver model to allow a reusable human spacecraft returning from a lunar mission to reach a given space station orbit. The trajectory was designed using a modified particle swarm optimization algorithm, and the vehicle mass architecture was also analyzed. The previous studies presented above primarily focused on the control optimization algorithm but neglected the coupling effect of the aerocapture phase and the trans-Earth phase in the cislunar space. So far, to the best of the author’s knowledge, there is no adequate work that has been reported to investigate the multiphase trajectory optimization including the trans-Earth phase, aerocapture phase, and postatmospheric phase.Compared with the traditional lunar return trajectory, the optimization of this multiphase trajectory has the following challenges: (1) The trajectory is composed of multiple phases with different dynamics and coupling constraints [21–23], and it is difficult to apply the standard algorithms to optimize the entire trajectory. The connection conditions between the phases must be considered, and an adaptive trajectory optimization method should be developed to design the entire trajectory. (2) The motion of the spacecraft is influenced by different dominant perturbative forces in different phases [24]. The flight time in the cislunar space is usually 3 to 5 days. This is quite long compared with the duration of the atmospheric trajectory near the Earth of approximately hundreds of seconds, which amplifies the nonlinear characteristics of the problem. (3) The mission design variables are complex, and the trajectory sensitive parameters vary with time which could cause convergence difficulty. To overcome the above challenges, in this study a corresponding trajectory optimization method is proposed to optimize each segment and to analyze the connection conditions between phases that are linked to each other. Then, based on problem characteristic analysis, a hybrid trajectory optimization method is established to optimize the entire trajectory of the reusable spacecraft returning to an LEO space station with sufficient robustness and good convergence. The main contributions of the presented work are as follows:
(1)
Analytical methods are proposed to estimate the orbital transfer window and to provide good initial values for optimizing the trans-Earth trajectory returning to an LEO space station, which satisfies the LEO coplanar constraints and reentry corridor conditions(2)
A hybrid optimization scheme is proposed for designing the multiphase trajectory of the reusable spacecraft. The connection nodes between different phases are carefully analyzed, and the top-level design variable is selected according to the problem characteristicsThe remainder is organized as follows. Section2 restates the problem and divides the entire trajectory into three segments. Section 3 establishes the corresponding optimization models for each phase according to their dynamic characteristics. Section 4 presents a hybrid optimization method used to solve the entire trajectory. Section 5 provides the simulation results and verifies the effectiveness of the proposed method. Finally, Section 6 states the major conclusions of this study.


## 2. Problem Statement

The manned lunar mission architecture employing an LEO space station is shown in Figure1. The Lunar Module (LM), Crew Module (CM), consisting of a Reentry Module and a Propulsion Module (PM), and Orbit Transfer Module (OTM) are delivered into LEO where they dock with the space station for in-orbit assembly. Then, the OTM performs an impulsive maneuver to transfer the LM and CM into the Earth-to-Moon transfer orbit. At the perilune, the LM and CM perform orbital injection into LLO. After the lunar surface tasks have been completed, the PM performs an impulsive maneuver to take the CM into the Earth transfer orbit (ETO). Upon arriving at the Earth’s atmosphere, the PM undocks from the Reentry Module, which directly returns the crew to the Earth, touching down on land. The PM then travels through the atmosphere to reduce its reentry speed and attain a specific postatmospheric orbit that enables the PM to enter the space station orbit with an acceptable velocity increment and dock with the space station. This procedure allows the PM to be reused after propellant supplement, reducing the rocket launch costs in future missions. This paper accordingly concentrates on how to design the trajectory of the PM from the Moon to the LEO space station, as described in Section 2.2.Figure 1
The lunar mission architecture employing a space station in LEO.The return trajectory of the PM is a multiphase process that begins at trans-Earth injection (TEI) and ends at the LEO space station. As shown in Figure2, the entire trajectory of this mission is separated into three phases: trans-Earth phase, aerocapture phase, and postatmospheric phase.Figure 2
Entire trajectory consisting of three phases.


### 2.1. Trans-Earth Phase

After completing its lunar surface task, the LM docks with the CM parked in a circular orbit around the Moon. At the appropriate timetprl, an impulse vprl parallel to the velocity of the parking orbit is executed by the PM, transferring the spacecraft to trans-Earth orbit.


### 2.2. Aerocapture Phase

The spacecraft enters the Earth’s atmosphere at a geocentric heighth, using the aerodynamic lift and drag forces to reduce its reentry speed and achieve the desired postatmospheric LEO station orbit, requiring only a small amount of propellant. Note that during the aerocapture phase, the thermal projection and overload constraints should be considered.


### 2.3. Postatmospheric Phase

The spacecraft exits the atmosphere and enters the specified postatmospheric orbit. In this orbit, the spacecraft transfers to the target space station LEO using a single- or double-impulse maneuver. The objective of this phase is to minimize the insertion velocity increment.The entire trajectory of the reusable spacecraft is a multiphase process with different dynamic models and nonlinear features. For the trans-Earth phase, the trajectory of the reusable PM is mainly influenced by the Moon gravity and the Earth gravity. The Sun also has a perturbative effect on PM dynamics. However, in the aerocapture phase and postatmospheric phase near the Earth, the motion of the PM is mainly determined by the Earth’s gravity. Particularly for the aerocapture phase, the aerodynamic forces could not be neglected.An underlying idea is to establish a hybrid optimization methodology at the subproblem level and top-problem level. At the subproblem level, each trajectory is solved in sequence by employing the appropriate optimization algorithm according to the dynamic features and mission constraints. After the sublevel problems are solved, the entire trajectory is optimized at the top level using the system design parameters. Note that the total mission goal can be obtained after each iteration at the sublevel. Therefore, in the next section, the specific optimization algorithm for each subproblem is introduced first and then a hybrid optimization framework is presented to select the top-level design variables and optimize the entire trajectory of the PM.


## 2.1. Trans-Earth Phase

After completing its lunar surface task, the LM docks with the CM parked in a circular orbit around the Moon. At the appropriate timetprl, an impulse vprl parallel to the velocity of the parking orbit is executed by the PM, transferring the spacecraft to trans-Earth orbit.


## 2.2. Aerocapture Phase

The spacecraft enters the Earth’s atmosphere at a geocentric heighth, using the aerodynamic lift and drag forces to reduce its reentry speed and achieve the desired postatmospheric LEO station orbit, requiring only a small amount of propellant. Note that during the aerocapture phase, the thermal projection and overload constraints should be considered.


## 2.3. Postatmospheric Phase

The spacecraft exits the atmosphere and enters the specified postatmospheric orbit. In this orbit, the spacecraft transfers to the target space station LEO using a single- or double-impulse maneuver. The objective of this phase is to minimize the insertion velocity increment.The entire trajectory of the reusable spacecraft is a multiphase process with different dynamic models and nonlinear features. For the trans-Earth phase, the trajectory of the reusable PM is mainly influenced by the Moon gravity and the Earth gravity. The Sun also has a perturbative effect on PM dynamics. However, in the aerocapture phase and postatmospheric phase near the Earth, the motion of the PM is mainly determined by the Earth’s gravity. Particularly for the aerocapture phase, the aerodynamic forces could not be neglected.An underlying idea is to establish a hybrid optimization methodology at the subproblem level and top-problem level. At the subproblem level, each trajectory is solved in sequence by employing the appropriate optimization algorithm according to the dynamic features and mission constraints. After the sublevel problems are solved, the entire trajectory is optimized at the top level using the system design parameters. Note that the total mission goal can be obtained after each iteration at the sublevel. Therefore, in the next section, the specific optimization algorithm for each subproblem is introduced first and then a hybrid optimization framework is presented to select the top-level design variables and optimize the entire trajectory of the PM.


## 3. Optimization Model for Trajectory in Each Phase

There are different dynamic models and constraints in each different phase of the entire trajectory of the PM. In this section, the optimal models and connection node boundary conditions are presented for the three considered phases.


### 3.1. Trans-Earth Phase Optimization in Cislunar Space

When the spacecraft arrives at the Earth’s atmosphere, the orbital inclination and right ascension of the ascending node (RAAN) is determined according to the ETO. The orbital plane of the ETO should be nearly the same as the plane of the station LEO owing to the high cost of orbital plane changes near the Earth. Therefore, a novel trans-Earth trajectory design method is proposed to satisfy the coplanar constraints between ETO and LEO.


#### 3.1.1. Moon-to-Earth Trajectory Optimization Method

The station LEO should be near-circular with a geocentric altitude ofhS, orbital inclination of iS, and RAAN of ΩS. The vacuum perigee (VCP) is defined as the perigee without considering atmospheric perturbations, as shown in Figure 3. Supposing that ivcp=iS and Ωvcp=ΩS, the coplanar constraint can be satisfied. Therefore, the design variables of the ETO are selected as
(1)x=vvcp,uvcp,tvcp,Δtperi‐vcp,where vvcp is the velocity, uvcp is the argument of latitude, Δtperi‐vcp is the transfer time, and tvcp is the time at VCP. Given these four parameters and the VCP height hvcp, the position vector rvcp and the velocity vector vvcp at VCP can be obtained by
(2)rvcp=cosΩScosuvcp−sinΩScosiSsinuvcpsinΩScosuvcp+cosΩScosiSsinuvcpsiniSsinuvcphvcp+RE,(3)vvcp=−cosΩSsinuvcp−sinΩScosiScosuvcp−sinΩScosuvcp+cosΩScosiScosuvcpsiniScosuvcpvvcp,where RE is the radius of the Earth. The rvcp and vvcp determine the ETO through the backward numerical integral method using the high-fidelity model. At TEI, the radius and inclination constraints are considered as follows:
(4)rprl=hLLO∗+RM,iprl=iLLO∗,where rLLO∗ and iLLO∗ are the selenocentric radius and orbital inclination of LLO, respectively, rprl and iprl are the perilune altitude and inclination of ETO, respectively, and RM is the radius of the Moon.Figure 3
Moon-to-Earth trajectory and vacuum perigee.In the trans-Earth phase, the objective function is to minimize TEI velocity increment as follows:(5)J1=minΔvprl.The sequential quadratic programming (SQP) algorithm has a good advantage of computation efficiency [25], and it is employed to perform this optimization. However, the performance of the algorithm depends closely on the initial guesses for the design variables. Thus, in Section 3.1.2, the Earth-centered Moon-to-Earth plane coordinate (ECMEPC) system is proposed to help estimate the orbital returning window and provide appropriate initial guesses for the design variables by two proposed pseudo rules.


#### 3.1.2. Trans-Earth Window Estimation by Analytical Methods

The geocentric relationship between the Earth, the Moon, ETO, and LEO is illustrated using the ECMEPC in Figure4. The Earth center is located at point O, TEI at point B, and VCP at point A. The ETO intersects the Moon-to-Earth plane at point E and the Earth equatorial plane at point D. At TEI time tprl, the angle between OE and the Moon position OM is ∠MOE=λprl, which is called the pseudo perilune longitude. For the returning orbit in cislunar space, the spacecraft is primarily influenced by the Earth’s gravity for most of the flight time. The orbital plane of the ETO and the station LEO could be considered to be approximately coplanar, yielding
(6)λprl=180∘.Figure 4
Sketch of ETO in the Earth-centered Moon-to-Earth plane coordinate system.Similarly, for the preliminary analysis, the position of the VCP on the Moon-to-Earth plane can be estimated by(7)AE¯=0∘,which is called the pseudo argument of latitude relative to the Moon-to-Earth plane.Then, the initial values of the design variables can be computed based on Equations (6) and (7) as follows. First, define the frame O‐xyz, where the x-axis increases from O to C, the z-axis increases along the direction of the Moon’s angular momentum, and the y-axis is arranged to complete a right-handed Cartesian coordinate system with the other two axes. In this coordinate system, the RAAN and inclination of the space station LEO are ΩL and iL, respectively, while in the Earth-centered J2000 inertial coordinate system, the RAAN and inclination of the Moon orbit are ΩM and iM, respectively.In the spherical triangleΔCDE, CE¯=ΩM−ΩS, ∠CDE=π−iS, and ∠CED=iL. The inclination iL can then be obtained using the cosine formula as follows:
(8)cosiL=cosiMcosiS+siniMsiniScosΩM−ΩS.Applying the sine formula to the spherical triangleΔCDE yields
(9)sinΩLsinπ−iS=sinΩM−ΩSsiniL=sinDE¯siniM.By substitutingiL into Equation (9), the ΩL can be obtained by
(10)ΩL=arcsinsiniSsinΩM−ΩSsiniL.In the ECMEPC, the longitude of the Moon at TEI is(11)Ωprl=π−ΩL.Therefore, the perilune timetprl can be computed by interpolation using Ωprl and the VCP point time tvcp as follows:
(12)tvcp=tprl+Δtprl‐vcp.In Figure4, uA is expressed by DA¯ and uL is expressed by EA¯, which satisfiesDA¯+AE¯=DE¯. Combining Equations (9) and (10), the argument of latitude at VCP is
(13)uvcp=DE¯=−arcsinsiniMsinΩLsiniS.The velocityvvcp is mainly decided by the Moon-to-Earth distance dM and can be estimated by
(14)vvcp=2μE1rvcp−1dM.In the above analysis, the initial values oftvcp and uvcp can be estimated based on Equations (12) and (13), which constitute the pseudo perilune longitude rule and the pseudo argument of latitude rule, respectively. Once the parameters in Equation (1) are optimized by the algorithm, the returning orbit from the Moon to the Earth is determined by backward integration. The orbital elements at the Earth entry surface (EI) can then be obtained to serve as the initial condition of the aerocapture phase, as presented in Section 3.2.


### 3.2. Aerocapture Phase Optimization near the Earth

Before arriving at VCP, the PM enters the Earth’s atmosphere, where the atmospheric perturbation is an important factor influencing its dynamics. Therefore, a different dynamic equation and optimization method that is employed for the trans-Earth trajectory optimization is proposed to solve this subproblem.


#### 3.2.1. Dynamic Equation in the Atmosphere

When the PM enters the Earth’s atmosphere, the six-degree-of-freedom equations of motion are given as [26]
(15)drdt=Vsinγ,dθdt=Vcosγsinψrcosφ,dφdt=Vcosγcosψr,dVdt=−Dm−gsinγ+ω2rcosφsinγcosφ−cosγsinφcosψ,dγdt=1VLcosσm+V2r−gcosγ+2ωVcosφsinψ+ω2rcosφcosγcosφ+sinγsinφcosψ,dψdt=1VLsinσmcosγ+V2rcosγsinψtanφ−2ωVcosφtanγcosψ−sinφ+ω2rcosγsinφcosφsinψ,where r is the radial distance from the Earth’s center to the PM, θ and φ are the longitude and geocentric latitude, respectively, V is the relative velocity of the Earth, γ is the flight-path angle of the relative velocity vector, ψ is the heading angle of the relative velocity vector measured clockwise from the north in the local horizontal plane, g is the gravitational acceleration, and m is the PM mass. Geometric descriptions of these six parameters are shown in Figure 5. Note that the term σ in Equation (15) denotes the bank angle, which is the rotation angle about the relative velocity vector. By appropriately controlling σ, the spacecraft can take advantage of the aerodynamic lift force L and drag force D in order to decrease its velocity. In most aerocapture missions, the bank angle is used as the control variable and optimization problems have been widely studied accordingly. The state and control variables are discretized at a series of Legendre-Gauss points, and then, the Gauss pseudospectral optimization method is employed to optimize the aerocapture trajectory.Figure 5
Parameters of three-dimensional dynamics in the Earth’s atmosphere.TheL and D can be obtained by
(16)L=12ρV2CLS,D=12ρV2CDS,where CD is the drag coefficient, CL is the lift coefficient, S is the reference area, and ρ is the atmospheric density.Several typical trajectory constraints must be considered during the flight in the Earth’s atmosphere. The maxima of the dynamic pressureq, load factor n, and heat rate Q̇ should not exceed the corresponding extreme values defined for the safe entry of the PM. Furthermore, the minimum geocentric altitude h should be greater than the extreme value hmin. All these constraints are formulized as
(17)q=12ρv2≤qmax,n=L2+D2mg≤nmax,Q̇=C1Rdρρ00.5vvc3.15≤Q̇max,h≥hmin,where C1 is a constant coefficient based on the PM characteristics, Rd is the nose radius of the PM, ρ0 is the density of the atmosphere at sea level, and vc is the reference velocity of a circular orbit with an altitude equal to the Earth’s radius. In a lunar exploration mission, qmax, nmax, Q̇max, and hmin are usually set to 30 kPa, 4.5 g, 2.6 MW/m2, and 40 km, respectively [17].In order to prevent the spacecraft’s heating shield system from being burned, the integrated heatQS should be minimized as follows:
(18)J2=minQS=min∫t2,0t2,fQ·dt,where t2,0 and t2,f are the starting and ending times, respectively, of the aerocapture phase.


#### 3.2.2. Boundary Conditions of the Aerocapture Phase

The boundary conditions in the aerocapture problem comprise two factors: the initial states and the terminal states. For the initial states, the boundary conditions are given by the reentry states of the trans-Earth trajectory as follows:(19)r2,i=rreen,θ2,i=θreen,φ2,i=φreen,v2,i=vreen,γ2,i=γreen,ψ2,i=ψreen,where the subscript “reen” denotes the reentry orbital states of the trans-Earth trajectory. The initial states are specified as follows.In the trans-Earth trajectory, the true anomalyfreen at the atmospheric reentry point is obtained by
(20)freen=arccosrvcp/rreen1+e−1e.The corresponding eccentric anomalyEreen is
(21)Ereen=2arctan1−e1+etanfreen2.The flight time between VCP and the reentry injection,Δtreen‐vcp, can be obtained by
(22)Δtreen‐vcp=areen3μEEreen−esinEreen,where areen is the semimajor axis of ETO and e the eccentricity. The reentry time treen is
(23)treen=tvcp−Δtreen‐vcp.By backward integratingrvcp,vvcp in Equations (2) and (3), rreen,vreen can be obtained and should then be transformed into the initial states of the aerocapture phase in Equation (19).The terminal boundary is given by(24)r2,f=r3,i,θ2,f=θ3,i,φ2,f=φ3,i,v2,f=v3,i,γ2,f=γ3,i,ψ2,f=ψ3,i,where the subscript “2,f” denotes the terminal states of the atmospheric trajectory and the subscript “3,i” denotes the initial orbital states of the postatmospheric trajectory, which is discussed in Section 3.3.


### 3.3. Postatmospheric Phase

In the postatmospheric phase, the optimization objective is to minimize the velocity increment that the spacecraft must accommodate when transferring from postatmospheric orbit into space station orbit. This objective function is formulized as(25)J3=minvpost‐atm.At the atmospheric exit point, letr2,f and v2,f be the position and velocity vectors, respectively, of the aerocapture trajectory. It should be noted that these two vectors are relative to the rotating Earth, and given t2,f, they can be transformed into the inertial vectors rexit and vexit in the Earth-centered inertial coordinate system.The semimajor axis of the postatmospheric orbit is(26)aexit=μE2μE/rexit−vexit2,and the apoapsis and periapsis radius is
(27)ra=aexit1+1−vexit2rexit2cos2γμaexit,rp=aexit1−1−vexit2rexit2cos2γμaexit.To reach the space station LEO, two engine burns are planned that are similar to a Hohmann transfer. The first tangential impulse is applied atra to raise rp to rp∗=ra+rS/2. The second impulse is applied at rp∗ to circularize the orbit and bring the apoapsis radius to rS.The total two-velocity increment is(28)vpost‐atm=μE2ra−2ra+rp−μE2ra−1aexit+μE2rs−2ra+rp−μErs.From Equation (28), it can be found that vpost‐atm is minimized when ra=rS, indicating that the apoapsis altitude after the aerocapture phase is equal to the orbital radius of the space station, which is the objective of the apoapsis-targeting problem. Therefore, the velocity increment in Equation (28) can be rewritten as
(29)vpost‐atm=μE2ra−2ra+rp−μE2ra−1aexit.Another set of boundary conditions in atmospheric flight is given by(30)rar3,i,θ3,i,φ3,i,v3,i,γ3,i,ψ3,i=rS.Substituting Equation (29) into Equation (25), the minimum velocity increment of the objective function can be specified.


## 3.1. Trans-Earth Phase Optimization in Cislunar Space

When the spacecraft arrives at the Earth’s atmosphere, the orbital inclination and right ascension of the ascending node (RAAN) is determined according to the ETO. The orbital plane of the ETO should be nearly the same as the plane of the station LEO owing to the high cost of orbital plane changes near the Earth. Therefore, a novel trans-Earth trajectory design method is proposed to satisfy the coplanar constraints between ETO and LEO.


### 3.1.1. Moon-to-Earth Trajectory Optimization Method

The station LEO should be near-circular with a geocentric altitude ofhS, orbital inclination of iS, and RAAN of ΩS. The vacuum perigee (VCP) is defined as the perigee without considering atmospheric perturbations, as shown in Figure 3. Supposing that ivcp=iS and Ωvcp=ΩS, the coplanar constraint can be satisfied. Therefore, the design variables of the ETO are selected as
(1)x=vvcp,uvcp,tvcp,Δtperi‐vcp,where vvcp is the velocity, uvcp is the argument of latitude, Δtperi‐vcp is the transfer time, and tvcp is the time at VCP. Given these four parameters and the VCP height hvcp, the position vector rvcp and the velocity vector vvcp at VCP can be obtained by
(2)rvcp=cosΩScosuvcp−sinΩScosiSsinuvcpsinΩScosuvcp+cosΩScosiSsinuvcpsiniSsinuvcphvcp+RE,(3)vvcp=−cosΩSsinuvcp−sinΩScosiScosuvcp−sinΩScosuvcp+cosΩScosiScosuvcpsiniScosuvcpvvcp,where RE is the radius of the Earth. The rvcp and vvcp determine the ETO through the backward numerical integral method using the high-fidelity model. At TEI, the radius and inclination constraints are considered as follows:
(4)rprl=hLLO∗+RM,iprl=iLLO∗,where rLLO∗ and iLLO∗ are the selenocentric radius and orbital inclination of LLO, respectively, rprl and iprl are the perilune altitude and inclination of ETO, respectively, and RM is the radius of the Moon.Figure 3
Moon-to-Earth trajectory and vacuum perigee.In the trans-Earth phase, the objective function is to minimize TEI velocity increment as follows:(5)J1=minΔvprl.The sequential quadratic programming (SQP) algorithm has a good advantage of computation efficiency [25], and it is employed to perform this optimization. However, the performance of the algorithm depends closely on the initial guesses for the design variables. Thus, in Section 3.1.2, the Earth-centered Moon-to-Earth plane coordinate (ECMEPC) system is proposed to help estimate the orbital returning window and provide appropriate initial guesses for the design variables by two proposed pseudo rules.


### 3.1.2. Trans-Earth Window Estimation by Analytical Methods

The geocentric relationship between the Earth, the Moon, ETO, and LEO is illustrated using the ECMEPC in Figure4. The Earth center is located at point O, TEI at point B, and VCP at point A. The ETO intersects the Moon-to-Earth plane at point E and the Earth equatorial plane at point D. At TEI time tprl, the angle between OE and the Moon position OM is ∠MOE=λprl, which is called the pseudo perilune longitude. For the returning orbit in cislunar space, the spacecraft is primarily influenced by the Earth’s gravity for most of the flight time. The orbital plane of the ETO and the station LEO could be considered to be approximately coplanar, yielding
(6)λprl=180∘.Figure 4
Sketch of ETO in the Earth-centered Moon-to-Earth plane coordinate system.Similarly, for the preliminary analysis, the position of the VCP on the Moon-to-Earth plane can be estimated by(7)AE¯=0∘,which is called the pseudo argument of latitude relative to the Moon-to-Earth plane.Then, the initial values of the design variables can be computed based on Equations (6) and (7) as follows. First, define the frame O‐xyz, where the x-axis increases from O to C, the z-axis increases along the direction of the Moon’s angular momentum, and the y-axis is arranged to complete a right-handed Cartesian coordinate system with the other two axes. In this coordinate system, the RAAN and inclination of the space station LEO are ΩL and iL, respectively, while in the Earth-centered J2000 inertial coordinate system, the RAAN and inclination of the Moon orbit are ΩM and iM, respectively.In the spherical triangleΔCDE, CE¯=ΩM−ΩS, ∠CDE=π−iS, and ∠CED=iL. The inclination iL can then be obtained using the cosine formula as follows:
(8)cosiL=cosiMcosiS+siniMsiniScosΩM−ΩS.Applying the sine formula to the spherical triangleΔCDE yields
(9)sinΩLsinπ−iS=sinΩM−ΩSsiniL=sinDE¯siniM.By substitutingiL into Equation (9), the ΩL can be obtained by
(10)ΩL=arcsinsiniSsinΩM−ΩSsiniL.In the ECMEPC, the longitude of the Moon at TEI is(11)Ωprl=π−ΩL.Therefore, the perilune timetprl can be computed by interpolation using Ωprl and the VCP point time tvcp as follows:
(12)tvcp=tprl+Δtprl‐vcp.In Figure4, uA is expressed by DA¯ and uL is expressed by EA¯, which satisfiesDA¯+AE¯=DE¯. Combining Equations (9) and (10), the argument of latitude at VCP is
(13)uvcp=DE¯=−arcsinsiniMsinΩLsiniS.The velocityvvcp is mainly decided by the Moon-to-Earth distance dM and can be estimated by
(14)vvcp=2μE1rvcp−1dM.In the above analysis, the initial values oftvcp and uvcp can be estimated based on Equations (12) and (13), which constitute the pseudo perilune longitude rule and the pseudo argument of latitude rule, respectively. Once the parameters in Equation (1) are optimized by the algorithm, the returning orbit from the Moon to the Earth is determined by backward integration. The orbital elements at the Earth entry surface (EI) can then be obtained to serve as the initial condition of the aerocapture phase, as presented in Section 3.2.


## 3.1.1. Moon-to-Earth Trajectory Optimization Method

The station LEO should be near-circular with a geocentric altitude ofhS, orbital inclination of iS, and RAAN of ΩS. The vacuum perigee (VCP) is defined as the perigee without considering atmospheric perturbations, as shown in Figure 3. Supposing that ivcp=iS and Ωvcp=ΩS, the coplanar constraint can be satisfied. Therefore, the design variables of the ETO are selected as
(1)x=vvcp,uvcp,tvcp,Δtperi‐vcp,where vvcp is the velocity, uvcp is the argument of latitude, Δtperi‐vcp is the transfer time, and tvcp is the time at VCP. Given these four parameters and the VCP height hvcp, the position vector rvcp and the velocity vector vvcp at VCP can be obtained by
(2)rvcp=cosΩScosuvcp−sinΩScosiSsinuvcpsinΩScosuvcp+cosΩScosiSsinuvcpsiniSsinuvcphvcp+RE,(3)vvcp=−cosΩSsinuvcp−sinΩScosiScosuvcp−sinΩScosuvcp+cosΩScosiScosuvcpsiniScosuvcpvvcp,where RE is the radius of the Earth. The rvcp and vvcp determine the ETO through the backward numerical integral method using the high-fidelity model. At TEI, the radius and inclination constraints are considered as follows:
(4)rprl=hLLO∗+RM,iprl=iLLO∗,where rLLO∗ and iLLO∗ are the selenocentric radius and orbital inclination of LLO, respectively, rprl and iprl are the perilune altitude and inclination of ETO, respectively, and RM is the radius of the Moon.Figure 3
Moon-to-Earth trajectory and vacuum perigee.In the trans-Earth phase, the objective function is to minimize TEI velocity increment as follows:(5)J1=minΔvprl.The sequential quadratic programming (SQP) algorithm has a good advantage of computation efficiency [25], and it is employed to perform this optimization. However, the performance of the algorithm depends closely on the initial guesses for the design variables. Thus, in Section 3.1.2, the Earth-centered Moon-to-Earth plane coordinate (ECMEPC) system is proposed to help estimate the orbital returning window and provide appropriate initial guesses for the design variables by two proposed pseudo rules.


## 3.1.2. Trans-Earth Window Estimation by Analytical Methods

The geocentric relationship between the Earth, the Moon, ETO, and LEO is illustrated using the ECMEPC in Figure4. The Earth center is located at point O, TEI at point B, and VCP at point A. The ETO intersects the Moon-to-Earth plane at point E and the Earth equatorial plane at point D. At TEI time tprl, the angle between OE and the Moon position OM is ∠MOE=λprl, which is called the pseudo perilune longitude. For the returning orbit in cislunar space, the spacecraft is primarily influenced by the Earth’s gravity for most of the flight time. The orbital plane of the ETO and the station LEO could be considered to be approximately coplanar, yielding
(6)λprl=180∘.Figure 4
Sketch of ETO in the Earth-centered Moon-to-Earth plane coordinate system.Similarly, for the preliminary analysis, the position of the VCP on the Moon-to-Earth plane can be estimated by(7)AE¯=0∘,which is called the pseudo argument of latitude relative to the Moon-to-Earth plane.Then, the initial values of the design variables can be computed based on Equations (6) and (7) as follows. First, define the frame O‐xyz, where the x-axis increases from O to C, the z-axis increases along the direction of the Moon’s angular momentum, and the y-axis is arranged to complete a right-handed Cartesian coordinate system with the other two axes. In this coordinate system, the RAAN and inclination of the space station LEO are ΩL and iL, respectively, while in the Earth-centered J2000 inertial coordinate system, the RAAN and inclination of the Moon orbit are ΩM and iM, respectively.In the spherical triangleΔCDE, CE¯=ΩM−ΩS, ∠CDE=π−iS, and ∠CED=iL. The inclination iL can then be obtained using the cosine formula as follows:
(8)cosiL=cosiMcosiS+siniMsiniScosΩM−ΩS.Applying the sine formula to the spherical triangleΔCDE yields
(9)sinΩLsinπ−iS=sinΩM−ΩSsiniL=sinDE¯siniM.By substitutingiL into Equation (9), the ΩL can be obtained by
(10)ΩL=arcsinsiniSsinΩM−ΩSsiniL.In the ECMEPC, the longitude of the Moon at TEI is(11)Ωprl=π−ΩL.Therefore, the perilune timetprl can be computed by interpolation using Ωprl and the VCP point time tvcp as follows:
(12)tvcp=tprl+Δtprl‐vcp.In Figure4, uA is expressed by DA¯ and uL is expressed by EA¯, which satisfiesDA¯+AE¯=DE¯. Combining Equations (9) and (10), the argument of latitude at VCP is
(13)uvcp=DE¯=−arcsinsiniMsinΩLsiniS.The velocityvvcp is mainly decided by the Moon-to-Earth distance dM and can be estimated by
(14)vvcp=2μE1rvcp−1dM.In the above analysis, the initial values oftvcp and uvcp can be estimated based on Equations (12) and (13), which constitute the pseudo perilune longitude rule and the pseudo argument of latitude rule, respectively. Once the parameters in Equation (1) are optimized by the algorithm, the returning orbit from the Moon to the Earth is determined by backward integration. The orbital elements at the Earth entry surface (EI) can then be obtained to serve as the initial condition of the aerocapture phase, as presented in Section 3.2.


## 3.2. Aerocapture Phase Optimization near the Earth

Before arriving at VCP, the PM enters the Earth’s atmosphere, where the atmospheric perturbation is an important factor influencing its dynamics. Therefore, a different dynamic equation and optimization method that is employed for the trans-Earth trajectory optimization is proposed to solve this subproblem.


### 3.2.1. Dynamic Equation in the Atmosphere

When the PM enters the Earth’s atmosphere, the six-degree-of-freedom equations of motion are given as [26]
(15)drdt=Vsinγ,dθdt=Vcosγsinψrcosφ,dφdt=Vcosγcosψr,dVdt=−Dm−gsinγ+ω2rcosφsinγcosφ−cosγsinφcosψ,dγdt=1VLcosσm+V2r−gcosγ+2ωVcosφsinψ+ω2rcosφcosγcosφ+sinγsinφcosψ,dψdt=1VLsinσmcosγ+V2rcosγsinψtanφ−2ωVcosφtanγcosψ−sinφ+ω2rcosγsinφcosφsinψ,where r is the radial distance from the Earth’s center to the PM, θ and φ are the longitude and geocentric latitude, respectively, V is the relative velocity of the Earth, γ is the flight-path angle of the relative velocity vector, ψ is the heading angle of the relative velocity vector measured clockwise from the north in the local horizontal plane, g is the gravitational acceleration, and m is the PM mass. Geometric descriptions of these six parameters are shown in Figure 5. Note that the term σ in Equation (15) denotes the bank angle, which is the rotation angle about the relative velocity vector. By appropriately controlling σ, the spacecraft can take advantage of the aerodynamic lift force L and drag force D in order to decrease its velocity. In most aerocapture missions, the bank angle is used as the control variable and optimization problems have been widely studied accordingly. The state and control variables are discretized at a series of Legendre-Gauss points, and then, the Gauss pseudospectral optimization method is employed to optimize the aerocapture trajectory.Figure 5
Parameters of three-dimensional dynamics in the Earth’s atmosphere.TheL and D can be obtained by
(16)L=12ρV2CLS,D=12ρV2CDS,where CD is the drag coefficient, CL is the lift coefficient, S is the reference area, and ρ is the atmospheric density.Several typical trajectory constraints must be considered during the flight in the Earth’s atmosphere. The maxima of the dynamic pressureq, load factor n, and heat rate Q̇ should not exceed the corresponding extreme values defined for the safe entry of the PM. Furthermore, the minimum geocentric altitude h should be greater than the extreme value hmin. All these constraints are formulized as
(17)q=12ρv2≤qmax,n=L2+D2mg≤nmax,Q̇=C1Rdρρ00.5vvc3.15≤Q̇max,h≥hmin,where C1 is a constant coefficient based on the PM characteristics, Rd is the nose radius of the PM, ρ0 is the density of the atmosphere at sea level, and vc is the reference velocity of a circular orbit with an altitude equal to the Earth’s radius. In a lunar exploration mission, qmax, nmax, Q̇max, and hmin are usually set to 30 kPa, 4.5 g, 2.6 MW/m2, and 40 km, respectively [17].In order to prevent the spacecraft’s heating shield system from being burned, the integrated heatQS should be minimized as follows:
(18)J2=minQS=min∫t2,0t2,fQ·dt,where t2,0 and t2,f are the starting and ending times, respectively, of the aerocapture phase.


### 3.2.2. Boundary Conditions of the Aerocapture Phase

The boundary conditions in the aerocapture problem comprise two factors: the initial states and the terminal states. For the initial states, the boundary conditions are given by the reentry states of the trans-Earth trajectory as follows:(19)r2,i=rreen,θ2,i=θreen,φ2,i=φreen,v2,i=vreen,γ2,i=γreen,ψ2,i=ψreen,where the subscript “reen” denotes the reentry orbital states of the trans-Earth trajectory. The initial states are specified as follows.In the trans-Earth trajectory, the true anomalyfreen at the atmospheric reentry point is obtained by
(20)freen=arccosrvcp/rreen1+e−1e.The corresponding eccentric anomalyEreen is
(21)Ereen=2arctan1−e1+etanfreen2.The flight time between VCP and the reentry injection,Δtreen‐vcp, can be obtained by
(22)Δtreen‐vcp=areen3μEEreen−esinEreen,where areen is the semimajor axis of ETO and e the eccentricity. The reentry time treen is
(23)treen=tvcp−Δtreen‐vcp.By backward integratingrvcp,vvcp in Equations (2) and (3), rreen,vreen can be obtained and should then be transformed into the initial states of the aerocapture phase in Equation (19).The terminal boundary is given by(24)r2,f=r3,i,θ2,f=θ3,i,φ2,f=φ3,i,v2,f=v3,i,γ2,f=γ3,i,ψ2,f=ψ3,i,where the subscript “2,f” denotes the terminal states of the atmospheric trajectory and the subscript “3,i” denotes the initial orbital states of the postatmospheric trajectory, which is discussed in Section 3.3.


## 3.2.1. Dynamic Equation in the Atmosphere

When the PM enters the Earth’s atmosphere, the six-degree-of-freedom equations of motion are given as [26]
(15)drdt=Vsinγ,dθdt=Vcosγsinψrcosφ,dφdt=Vcosγcosψr,dVdt=−Dm−gsinγ+ω2rcosφsinγcosφ−cosγsinφcosψ,dγdt=1VLcosσm+V2r−gcosγ+2ωVcosφsinψ+ω2rcosφcosγcosφ+sinγsinφcosψ,dψdt=1VLsinσmcosγ+V2rcosγsinψtanφ−2ωVcosφtanγcosψ−sinφ+ω2rcosγsinφcosφsinψ,where r is the radial distance from the Earth’s center to the PM, θ and φ are the longitude and geocentric latitude, respectively, V is the relative velocity of the Earth, γ is the flight-path angle of the relative velocity vector, ψ is the heading angle of the relative velocity vector measured clockwise from the north in the local horizontal plane, g is the gravitational acceleration, and m is the PM mass. Geometric descriptions of these six parameters are shown in Figure 5. Note that the term σ in Equation (15) denotes the bank angle, which is the rotation angle about the relative velocity vector. By appropriately controlling σ, the spacecraft can take advantage of the aerodynamic lift force L and drag force D in order to decrease its velocity. In most aerocapture missions, the bank angle is used as the control variable and optimization problems have been widely studied accordingly. The state and control variables are discretized at a series of Legendre-Gauss points, and then, the Gauss pseudospectral optimization method is employed to optimize the aerocapture trajectory.Figure 5
Parameters of three-dimensional dynamics in the Earth’s atmosphere.TheL and D can be obtained by
(16)L=12ρV2CLS,D=12ρV2CDS,where CD is the drag coefficient, CL is the lift coefficient, S is the reference area, and ρ is the atmospheric density.Several typical trajectory constraints must be considered during the flight in the Earth’s atmosphere. The maxima of the dynamic pressureq, load factor n, and heat rate Q̇ should not exceed the corresponding extreme values defined for the safe entry of the PM. Furthermore, the minimum geocentric altitude h should be greater than the extreme value hmin. All these constraints are formulized as
(17)q=12ρv2≤qmax,n=L2+D2mg≤nmax,Q̇=C1Rdρρ00.5vvc3.15≤Q̇max,h≥hmin,where C1 is a constant coefficient based on the PM characteristics, Rd is the nose radius of the PM, ρ0 is the density of the atmosphere at sea level, and vc is the reference velocity of a circular orbit with an altitude equal to the Earth’s radius. In a lunar exploration mission, qmax, nmax, Q̇max, and hmin are usually set to 30 kPa, 4.5 g, 2.6 MW/m2, and 40 km, respectively [17].In order to prevent the spacecraft’s heating shield system from being burned, the integrated heatQS should be minimized as follows:
(18)J2=minQS=min∫t2,0t2,fQ·dt,where t2,0 and t2,f are the starting and ending times, respectively, of the aerocapture phase.


## 3.2.2. Boundary Conditions of the Aerocapture Phase

The boundary conditions in the aerocapture problem comprise two factors: the initial states and the terminal states. For the initial states, the boundary conditions are given by the reentry states of the trans-Earth trajectory as follows:(19)r2,i=rreen,θ2,i=θreen,φ2,i=φreen,v2,i=vreen,γ2,i=γreen,ψ2,i=ψreen,where the subscript “reen” denotes the reentry orbital states of the trans-Earth trajectory. The initial states are specified as follows.In the trans-Earth trajectory, the true anomalyfreen at the atmospheric reentry point is obtained by
(20)freen=arccosrvcp/rreen1+e−1e.The corresponding eccentric anomalyEreen is
(21)Ereen=2arctan1−e1+etanfreen2.The flight time between VCP and the reentry injection,Δtreen‐vcp, can be obtained by
(22)Δtreen‐vcp=areen3μEEreen−esinEreen,where areen is the semimajor axis of ETO and e the eccentricity. The reentry time treen is
(23)treen=tvcp−Δtreen‐vcp.By backward integratingrvcp,vvcp in Equations (2) and (3), rreen,vreen can be obtained and should then be transformed into the initial states of the aerocapture phase in Equation (19).The terminal boundary is given by(24)r2,f=r3,i,θ2,f=θ3,i,φ2,f=φ3,i,v2,f=v3,i,γ2,f=γ3,i,ψ2,f=ψ3,i,where the subscript “2,f” denotes the terminal states of the atmospheric trajectory and the subscript “3,i” denotes the initial orbital states of the postatmospheric trajectory, which is discussed in Section 3.3.


## 3.3. Postatmospheric Phase

In the postatmospheric phase, the optimization objective is to minimize the velocity increment that the spacecraft must accommodate when transferring from postatmospheric orbit into space station orbit. This objective function is formulized as(25)J3=minvpost‐atm.At the atmospheric exit point, letr2,f and v2,f be the position and velocity vectors, respectively, of the aerocapture trajectory. It should be noted that these two vectors are relative to the rotating Earth, and given t2,f, they can be transformed into the inertial vectors rexit and vexit in the Earth-centered inertial coordinate system.The semimajor axis of the postatmospheric orbit is(26)aexit=μE2μE/rexit−vexit2,and the apoapsis and periapsis radius is
(27)ra=aexit1+1−vexit2rexit2cos2γμaexit,rp=aexit1−1−vexit2rexit2cos2γμaexit.To reach the space station LEO, two engine burns are planned that are similar to a Hohmann transfer. The first tangential impulse is applied atra to raise rp to rp∗=ra+rS/2. The second impulse is applied at rp∗ to circularize the orbit and bring the apoapsis radius to rS.The total two-velocity increment is(28)vpost‐atm=μE2ra−2ra+rp−μE2ra−1aexit+μE2rs−2ra+rp−μErs.From Equation (28), it can be found that vpost‐atm is minimized when ra=rS, indicating that the apoapsis altitude after the aerocapture phase is equal to the orbital radius of the space station, which is the objective of the apoapsis-targeting problem. Therefore, the velocity increment in Equation (28) can be rewritten as
(29)vpost‐atm=μE2ra−2ra+rp−μE2ra−1aexit.Another set of boundary conditions in atmospheric flight is given by(30)rar3,i,θ3,i,φ3,i,v3,i,γ3,i,ψ3,i=rS.Substituting Equation (29) into Equation (25), the minimum velocity increment of the objective function can be specified.


## 4. Algorithm Design for Multiphase Trajectory Optimization

The entire trajectory optimization scheme proposed in this paper is actually a typical multiphase optimization problem. It is difficult to apply the standard algorithms to optimize the entire trajectory directly because of the different dynamics and coupling constraints in each phase. In the trans-Earth phase, the spacecraft is mainly influenced by solar, lunar, and Earth gravities, whereas in the atmospheric phase, the atmospheric forces must be considered. Furthermore, the flight time in the cislunar space is usually several days and it is quite long compared with the duration of the atmospheric phase near the Earth, which amplifies the nonlinearity of the problem. An underlying idea is to introduce a hybrid optimization scheme at the sublevel and top level to solve such a problem [23, 24].Letxi denote the optimization variable in phase i. Each subproblem is then solved in sequence as follows:
(31)Ji=minfiX,xi,xi·=giX,xi,s.t.eqiX,xi=0,ieqiX,xi≤0,where Ji is the goal function, gi is the dynamic equation, eqi is the equality constraint, and ieqi is the inequality constraint in phase i.The terminal states of each previous phase are used as the initial states of each following phase and are referred to as the phase connection conditions. The linkage constraint of such phase connection conditions is formulated as follows:(32)rk+1,i=rk,f,vk+1,i=vk,f,where the subscript “k+1,i” denotes the initial states in the k+1 phase, the subscript “k,f” denotes the final states in the k phase, and r and v are the position and velocity vectors of the PM, respectively.The top optimization level of the proposed method deals with the phase connection nodes and combines the phases with each other. The mission design variableX is defined to minimize the total mission objective function J as follows:
(33)J=minFX,x=minFf1X,x1,⋯,fnX,xn, where n represents the number of phases.By splitting the entire trajectory into the trans-Earth phase, aerocapture phase, and postatmospheric phase, each of them can be regarded as a subproblem and optimized separately in sequence as described in Section3. For the trans-Earth phase trajectory, two pseudo rules are proposed to estimate the orbital transfer window and the SQP algorithm is used to optimize the orbital design parameters to satisfy the coplanar constraints of the space station orbit. For the aerocapture phase trajectory, the bank angle is used as the control variable and is discretized at a series of Legendre-Gauss points. This is an optimal control problem, and the Gauss pseudospectral optimization is employed to solve it through conversion into a nonlinear programming problem. For the postatmospheric phase trajectory, the impulse Δvpost‐atm that inserts the spacecraft into the target station orbit is optimized which is determined by the terminal orbital states of the aerocapture phase and the phase is similar to a Hohmann transfer.The general objective of the top-level optimization of the entire trajectory is to minimize the total velocity increment by(34)J=min∑i=1NΔvi,where N is the number of maneuver times and Δvi is the velocity increment in each phase. The objective function J includes the impulse Δvperi at the perilune in the first phase and Δvpost‐atm at the apoapsis in the third phase.At the top level, there are many mission design variables that can be selected to optimize the entire trajectory. For the multiphase problem addressed in this paper, the VCP heighthvcp is selected as the design parameter for the top-level optimization according to problem characteristic analysis. Firstly, it strongly couples the first and second subproblems. Indeed, the VCP height hvcp not only serves to determine the trans-Earth phase trajectory in Equation (2) and to influence the size of the first velocity impulse at the perilune but also influences the initial states of the aerocapture phase. The relationship between hvcp and γ can be estimated as
(35)tanγ=−rreen21+e+2rvcprreen−rreen21−erreen1+e.Secondly, it has been found that the reentry angleγ also influences the aerocapture phase and can result in an unacceptable velocity increment [15]. Therefore, the VCP height is selected as the mission design variable in the top level and this kind of problem could be optimized by the golden section method. This algorithm is very suitable for searching without derivative for the extrema of unimodal objective functions [27, 28].The entire optimization scheme is illustrated in Figure6. The detailed implementation process is described as follows.Figure 6
Flow chart of the hybrid optimization algorithm.Step 1. According to theiS,ΩS of the space station, the trans-Earth window is estimated based on the pseudo perilune longitude rule in Equation (6) and the pseudo argument of latitude rule in Equation (7).Step 2. Set the initial value ofhvcp randomly between 0 and 100 km.Step 3. Estimate the initial values of the design variablesvvcp,fvcp,tvcp,Δtperi‐vcp, employing the SQP algorithm to optimize the trans-Earth phase in cislunar space.Step 4. Determine the Earth entry interface and the initial statesr,θ,φ,v,γ,ψ2,i of the aerocapture phase.Step 5. In the aerocapture phase, the Gauss pseudospectral optimization method is employed to optimize the bank angleσ and to minimize the integrated total heat flux in Equation (18).Step 6. In the postatmospheric phase, calculate the initial statesr,θ,φ,v,γ,ψ3,i and optimize the insertion velocity increment vpost‐atm to the target space station LEO which is similar to a Hohmann transfer.Step 7. Apply the golden section method to optimizehvcp according to the total velocity increment J=minΔvprl+vpost‐atm and repeat Steps 3–7 until the mission constraints are all satisfied and the algorithm at the top-level is converged.


## 5. Simulation

### 5.1. Optimization Results for the Entire Trajectory

#### 5.1.1. Trans-Earth Trajectory in Cislunar Space

The TEI departure timetpri is set between 15 April 2028 08:00:00 and 15 May 2028 08:00:00. The orbital elements of the LEO space station are assumed as given in Table 1. Using the pseudo perilune longitude rule, the initial value of tpri is determined to be 19 Apr 2028 15:29:1252. The altitude and inclination constraints of LLO are set to hLLO∗=200km and iLLO∗=90∘. A desktop PC with a 3.6 GHz Intel Core i7-4790 and 4 GB of RAM was used to run this simulation.Table 1
Orbital elements of an LEO space station.
hSkmeiSoΩS∘343.0042.430.0The four design parameters are constrained by the following inequality:(36)10500m/s<vvcp<11500m/s,−10∘<Δuvcp<10∘,−10h<Δtvcp<10h,2.8day<Δtperi‐vcp<3.2day.For the top level in the optimization scheme, the optimal VCP height is determined to behvcp=76.08km with a total velocity increment of 914.8666 m/s. It consists of 837.20 m/s at the perilune and 77.66 m/s in the postatmospheric insertion. The total computation time is 261.68 s. The other numerical results of the trans-Earth subproblem are listed in Tables 2 and 3. A three-dimensional demonstration of the trans-Earth trajectory in cislunar space is given in Figure 7.Table 2
Solutions for ETOs.
Design parameterOptimal valuevvcpm/s11030.9310uvcp∘39.2871uL∘-0.6291tpri19 April 2028 19:8:31.3362λprl∘181.96Δtperi‐vcpday3.20Table 3
Orbital parameters at TEI.
Design parameterOptimal valueiLLO∘89.96ΩLLO∘354.10hLLOkm199585.42uLLO∘50.31vBm/s837.20Figure 7
Three-dimensional demonstration of ETO in cislunar space.It can be observed in the optimization results in Table2 that the pseudo perilune longitude is λprl=181.96∘ and the pseudo argument of latitude is uL=−0.63∘, which verifies the assumption of the initial values in Equations (6) and (7). The flight time Δtpri‐peri is determined to be 3.2 days, reaching the upper boundary of the time constraint because the velocity cost at TEI decreases with flight time for a Moon-returning orbit [11].


#### 5.1.2. Aerocapture Trajectory near the Earth

The initial states of the aerocapture phase at the atmospheric entry injection point are listed in Table4. Taking the Apollo spacecraft as the reference [29], the drag coefficient is set to 1.2891 and the lift coefficient to 0.3877. The radius of the reference area is 3.0 m. The other spacecraft parameters are listed in Table 4. The initial states of the aerocapture process are shown in Table 5. The time history of the altitude, speed, longitude, latitude, heading angle, and flight path angle is shown in Figures 8–10.Table 4
Aerocapture parameters for the PM.
ParameterOptimal valueMassm (kg)6000Reference area radiusS (m)3.0Lift coefficientCL0.3877Drag coefficientCD1.2891Reentry altitudeh (km)122Nose radiusRd (m)2.5Table 5
Initial states of the aerocapture phase.
h2,imθ2,i∘φ2,i∘v2,im/sγ2,i∘ψ2,i∘122241.8236-52.310518.958010639.7154-5.337949.9770Figure 8
Time history of (a) geocentric altitudeh and (b) velocity v.
(a)(b)Figure 9
Time history of (a) geocentric longitudeθ and (b) latitude φ.
(a)(b)Figure 10
Time history of (a) heading angleψ and (b) flight-path angle γ.
(a)(b)It can be observed in Figure8 that the altitude of the PM decreases sharply to h=67.48 km at t=100.2 s, satisfying the altitude constraint in Equation (17), and then increases modestly to h=122 km. The velocity v increases by approximately 50 m/s at the beginning of the airbrake but then consistently decreases to a minimum of 7.888 km/s before the PM exits the atmosphere. The total deceleration is approximately 2.7517 km/s. The entire flight time of the aerocapture phase is about 782.8 s.The geocentric longitudeλ and latitude φ are shown in Figure 9 to both increase with time. The longitude increases from 18.96∘ to 43.52∘ following a nearly logarithmic trend whereas the latitude increases almost linearly from −52.31∘ to 12.51∘.In Figure10, the heading angle ψ increases steadily from 50.00∘ at t=0 s to 87.22∘ at t=100.2 s. The flight-path angle γ first increases sharply from −5.33∘ to 0.57∘ prior to t=107.3 s, preventing the PM from hitting the Earth, then varies modestly from 0.57∘ to 0.85∘ at the end of the flight to achieve the postatmosphere apoapsis.The optimization result of the bank angle profileσ is shown in Figure 11. The PM moves with a bank angle of σ=0∘ during the early part of its flight. At t=108.9 s, the angle switches from 0∘ to 90∘ before immediately switching from 90∘ to −90∘ following a “bang-bang” profile. Note that it was theoretically proven in [15] that an optimal solution with such a “bang-bang” profile is reflective of the aerocapture problem.Figure 11
Bank angleσ with time t.The time histories of the constraint parameters including the load factorn, heat rate Q·, and dynamic pressure q are shown in Figure 11, in which it can be observed that all three constraints initially increase rapidly but then slowly decrease with time. The maximum values of all three constraints are within the allowable range, with a load factor of 2.5024 g, heat rate of 708.2966 kw/m2, and dynamic pressure of 5566.9255 kPa as demonstrated in Figure 12. Furthermore, the integrated heat minimized in Equation (18) by the optimization process is 123.4263 MJ/m2.Figure 12
(a) Load, (b) heat rate, and (c) dynamic pressure vs. timet.
(a)(b)(c)


#### 5.1.3. Postatmospheric Trajectory

At its atmospheric exit point, the spacecraft reaches a postatmospheric orbit with an apoapsis altitudeha=343.00km and a periapsis altitude hp=78.95km. After a flight time of 1966.88 s, the spacecraft transfers from its atmospheric exit point to its apoapsis. Using Equation (29), the velocity increment that inserts the spacecraft into the station LEO is about vpost‐atm=77.66m/s, which is less than the result vpost‐atm=95.9623m/s obtained in Ref. [21]. The resulting three-dimensional trajectory of the postatmospheric orbit is shown in Figure 13.Figure 13
Three-dimensional demonstration of the three-phase trajectory near the Earth.After the transfer of the trans-Earth phase, aerocapture phase, and postatmospheric phase, the mission constraints are recalculated to determine whether the spacecraft could return to the target LEO space station. The errors of the coplanar constraint and orbital height areΔiS=0.0020∘, ΔΩS=0.0022∘, and ΔhS=0.51m, respectively, which are acceptable for the engineering requirements. It can be concluded that the proposed algorithm is effective and satisfies the nonlinear constraints in the whole mission.


### 5.2. Performance of the Hybrid Optimization Algorithms

To validate the trajectory design method based on the “perilune” and “argument” pseudo estimation rules shown in Equations (6) and (7), respectively, 1000 space station LEOs with random inclinations and RAANs were generated and optimized. The results of the subsequent convergence analysis are shown in Table 6. It can be observed that all 1000 generated orbits optimized using the proposed orbit design method based on the combined perilune and argument pseudo estimation rules converge to corresponding optimal solutions, a convergence rate of 100.0%. However, note that if no good initial value is provided to the algorithm, none of these orbits will converge. It can also be observed that the pseudo perilune rule plays a more important role than the pseudo argument of latitude rule, as they exhibit convergence rates of 87.5% and 21.9%, respectively. The results in Table 6 thus show that the proposed orbital design method provides good convergence and can effectively solve the ETO to a space station in LEO with an arbitrary orbital plane.Table 6
Convergence analysis.
Initial estimated valueOrbitConvergenceAgreementPerilune+argument10001000100.0%Perilune100087587.5%Argument100021921.9%None100000.00%The variation in the velocity incrementvpost‐atm with the reentry angle γreen can be observed in Figure 14 to be uniformly dispersed between −6.5∘ and −4.3∘. It is interesting that vpost‐atm decreases dramatically from 900 m/s to 100 m/s for γreen<−6.0∘ and that vpost‐atm is greater than 100 m/s when γreen>−4.5∘. According to the relationship between the VCP height and reentry angle shown in Figure 15, it can be determined that the preliminary VCP height corresponding to a velocity increment of less than 100 m/s is 50–80 km, which is acceptable from the perspective of the propellant consumption of the reusable PM.Figure 14
Variation in velocity increment with reentry angle.Figure 15
Relationship between VCP height and reentry angle.


## 5.1. Optimization Results for the Entire Trajectory

### 5.1.1. Trans-Earth Trajectory in Cislunar Space

The TEI departure timetpri is set between 15 April 2028 08:00:00 and 15 May 2028 08:00:00. The orbital elements of the LEO space station are assumed as given in Table 1. Using the pseudo perilune longitude rule, the initial value of tpri is determined to be 19 Apr 2028 15:29:1252. The altitude and inclination constraints of LLO are set to hLLO∗=200km and iLLO∗=90∘. A desktop PC with a 3.6 GHz Intel Core i7-4790 and 4 GB of RAM was used to run this simulation.Table 1
Orbital elements of an LEO space station.
hSkmeiSoΩS∘343.0042.430.0The four design parameters are constrained by the following inequality:(36)10500m/s<vvcp<11500m/s,−10∘<Δuvcp<10∘,−10h<Δtvcp<10h,2.8day<Δtperi‐vcp<3.2day.For the top level in the optimization scheme, the optimal VCP height is determined to behvcp=76.08km with a total velocity increment of 914.8666 m/s. It consists of 837.20 m/s at the perilune and 77.66 m/s in the postatmospheric insertion. The total computation time is 261.68 s. The other numerical results of the trans-Earth subproblem are listed in Tables 2 and 3. A three-dimensional demonstration of the trans-Earth trajectory in cislunar space is given in Figure 7.Table 2
Solutions for ETOs.
Design parameterOptimal valuevvcpm/s11030.9310uvcp∘39.2871uL∘-0.6291tpri19 April 2028 19:8:31.3362λprl∘181.96Δtperi‐vcpday3.20Table 3
Orbital parameters at TEI.
Design parameterOptimal valueiLLO∘89.96ΩLLO∘354.10hLLOkm199585.42uLLO∘50.31vBm/s837.20Figure 7
Three-dimensional demonstration of ETO in cislunar space.It can be observed in the optimization results in Table2 that the pseudo perilune longitude is λprl=181.96∘ and the pseudo argument of latitude is uL=−0.63∘, which verifies the assumption of the initial values in Equations (6) and (7). The flight time Δtpri‐peri is determined to be 3.2 days, reaching the upper boundary of the time constraint because the velocity cost at TEI decreases with flight time for a Moon-returning orbit [11].


### 5.1.2. Aerocapture Trajectory near the Earth

The initial states of the aerocapture phase at the atmospheric entry injection point are listed in Table4. Taking the Apollo spacecraft as the reference [29], the drag coefficient is set to 1.2891 and the lift coefficient to 0.3877. The radius of the reference area is 3.0 m. The other spacecraft parameters are listed in Table 4. The initial states of the aerocapture process are shown in Table 5. The time history of the altitude, speed, longitude, latitude, heading angle, and flight path angle is shown in Figures 8–10.Table 4
Aerocapture parameters for the PM.
ParameterOptimal valueMassm (kg)6000Reference area radiusS (m)3.0Lift coefficientCL0.3877Drag coefficientCD1.2891Reentry altitudeh (km)122Nose radiusRd (m)2.5Table 5
Initial states of the aerocapture phase.
h2,imθ2,i∘φ2,i∘v2,im/sγ2,i∘ψ2,i∘122241.8236-52.310518.958010639.7154-5.337949.9770Figure 8
Time history of (a) geocentric altitudeh and (b) velocity v.
(a)(b)Figure 9
Time history of (a) geocentric longitudeθ and (b) latitude φ.
(a)(b)Figure 10
Time history of (a) heading angleψ and (b) flight-path angle γ.
(a)(b)It can be observed in Figure8 that the altitude of the PM decreases sharply to h=67.48 km at t=100.2 s, satisfying the altitude constraint in Equation (17), and then increases modestly to h=122 km. The velocity v increases by approximately 50 m/s at the beginning of the airbrake but then consistently decreases to a minimum of 7.888 km/s before the PM exits the atmosphere. The total deceleration is approximately 2.7517 km/s. The entire flight time of the aerocapture phase is about 782.8 s.The geocentric longitudeλ and latitude φ are shown in Figure 9 to both increase with time. The longitude increases from 18.96∘ to 43.52∘ following a nearly logarithmic trend whereas the latitude increases almost linearly from −52.31∘ to 12.51∘.In Figure10, the heading angle ψ increases steadily from 50.00∘ at t=0 s to 87.22∘ at t=100.2 s. The flight-path angle γ first increases sharply from −5.33∘ to 0.57∘ prior to t=107.3 s, preventing the PM from hitting the Earth, then varies modestly from 0.57∘ to 0.85∘ at the end of the flight to achieve the postatmosphere apoapsis.The optimization result of the bank angle profileσ is shown in Figure 11. The PM moves with a bank angle of σ=0∘ during the early part of its flight. At t=108.9 s, the angle switches from 0∘ to 90∘ before immediately switching from 90∘ to −90∘ following a “bang-bang” profile. Note that it was theoretically proven in [15] that an optimal solution with such a “bang-bang” profile is reflective of the aerocapture problem.Figure 11
Bank angleσ with time t.The time histories of the constraint parameters including the load factorn, heat rate Q·, and dynamic pressure q are shown in Figure 11, in which it can be observed that all three constraints initially increase rapidly but then slowly decrease with time. The maximum values of all three constraints are within the allowable range, with a load factor of 2.5024 g, heat rate of 708.2966 kw/m2, and dynamic pressure of 5566.9255 kPa as demonstrated in Figure 12. Furthermore, the integrated heat minimized in Equation (18) by the optimization process is 123.4263 MJ/m2.Figure 12
(a) Load, (b) heat rate, and (c) dynamic pressure vs. timet.
(a)(b)(c)


### 5.1.3. Postatmospheric Trajectory

At its atmospheric exit point, the spacecraft reaches a postatmospheric orbit with an apoapsis altitudeha=343.00km and a periapsis altitude hp=78.95km. After a flight time of 1966.88 s, the spacecraft transfers from its atmospheric exit point to its apoapsis. Using Equation (29), the velocity increment that inserts the spacecraft into the station LEO is about vpost‐atm=77.66m/s, which is less than the result vpost‐atm=95.9623m/s obtained in Ref. [21]. The resulting three-dimensional trajectory of the postatmospheric orbit is shown in Figure 13.Figure 13
Three-dimensional demonstration of the three-phase trajectory near the Earth.After the transfer of the trans-Earth phase, aerocapture phase, and postatmospheric phase, the mission constraints are recalculated to determine whether the spacecraft could return to the target LEO space station. The errors of the coplanar constraint and orbital height areΔiS=0.0020∘, ΔΩS=0.0022∘, and ΔhS=0.51m, respectively, which are acceptable for the engineering requirements. It can be concluded that the proposed algorithm is effective and satisfies the nonlinear constraints in the whole mission.


## 5.1.1. Trans-Earth Trajectory in Cislunar Space

The TEI departure timetpri is set between 15 April 2028 08:00:00 and 15 May 2028 08:00:00. The orbital elements of the LEO space station are assumed as given in Table 1. Using the pseudo perilune longitude rule, the initial value of tpri is determined to be 19 Apr 2028 15:29:1252. The altitude and inclination constraints of LLO are set to hLLO∗=200km and iLLO∗=90∘. A desktop PC with a 3.6 GHz Intel Core i7-4790 and 4 GB of RAM was used to run this simulation.Table 1
Orbital elements of an LEO space station.
hSkmeiSoΩS∘343.0042.430.0The four design parameters are constrained by the following inequality:(36)10500m/s<vvcp<11500m/s,−10∘<Δuvcp<10∘,−10h<Δtvcp<10h,2.8day<Δtperi‐vcp<3.2day.For the top level in the optimization scheme, the optimal VCP height is determined to behvcp=76.08km with a total velocity increment of 914.8666 m/s. It consists of 837.20 m/s at the perilune and 77.66 m/s in the postatmospheric insertion. The total computation time is 261.68 s. The other numerical results of the trans-Earth subproblem are listed in Tables 2 and 3. A three-dimensional demonstration of the trans-Earth trajectory in cislunar space is given in Figure 7.Table 2
Solutions for ETOs.
Design parameterOptimal valuevvcpm/s11030.9310uvcp∘39.2871uL∘-0.6291tpri19 April 2028 19:8:31.3362λprl∘181.96Δtperi‐vcpday3.20Table 3
Orbital parameters at TEI.
Design parameterOptimal valueiLLO∘89.96ΩLLO∘354.10hLLOkm199585.42uLLO∘50.31vBm/s837.20Figure 7
Three-dimensional demonstration of ETO in cislunar space.It can be observed in the optimization results in Table2 that the pseudo perilune longitude is λprl=181.96∘ and the pseudo argument of latitude is uL=−0.63∘, which verifies the assumption of the initial values in Equations (6) and (7). The flight time Δtpri‐peri is determined to be 3.2 days, reaching the upper boundary of the time constraint because the velocity cost at TEI decreases with flight time for a Moon-returning orbit [11].


## 5.1.2. Aerocapture Trajectory near the Earth

The initial states of the aerocapture phase at the atmospheric entry injection point are listed in Table4. Taking the Apollo spacecraft as the reference [29], the drag coefficient is set to 1.2891 and the lift coefficient to 0.3877. The radius of the reference area is 3.0 m. The other spacecraft parameters are listed in Table 4. The initial states of the aerocapture process are shown in Table 5. The time history of the altitude, speed, longitude, latitude, heading angle, and flight path angle is shown in Figures 8–10.Table 4
Aerocapture parameters for the PM.
ParameterOptimal valueMassm (kg)6000Reference area radiusS (m)3.0Lift coefficientCL0.3877Drag coefficientCD1.2891Reentry altitudeh (km)122Nose radiusRd (m)2.5Table 5
Initial states of the aerocapture phase.
h2,imθ2,i∘φ2,i∘v2,im/sγ2,i∘ψ2,i∘122241.8236-52.310518.958010639.7154-5.337949.9770Figure 8
Time history of (a) geocentric altitudeh and (b) velocity v.
(a)(b)Figure 9
Time history of (a) geocentric longitudeθ and (b) latitude φ.
(a)(b)Figure 10
Time history of (a) heading angleψ and (b) flight-path angle γ.
(a)(b)It can be observed in Figure8 that the altitude of the PM decreases sharply to h=67.48 km at t=100.2 s, satisfying the altitude constraint in Equation (17), and then increases modestly to h=122 km. The velocity v increases by approximately 50 m/s at the beginning of the airbrake but then consistently decreases to a minimum of 7.888 km/s before the PM exits the atmosphere. The total deceleration is approximately 2.7517 km/s. The entire flight time of the aerocapture phase is about 782.8 s.The geocentric longitudeλ and latitude φ are shown in Figure 9 to both increase with time. The longitude increases from 18.96∘ to 43.52∘ following a nearly logarithmic trend whereas the latitude increases almost linearly from −52.31∘ to 12.51∘.In Figure10, the heading angle ψ increases steadily from 50.00∘ at t=0 s to 87.22∘ at t=100.2 s. The flight-path angle γ first increases sharply from −5.33∘ to 0.57∘ prior to t=107.3 s, preventing the PM from hitting the Earth, then varies modestly from 0.57∘ to 0.85∘ at the end of the flight to achieve the postatmosphere apoapsis.The optimization result of the bank angle profileσ is shown in Figure 11. The PM moves with a bank angle of σ=0∘ during the early part of its flight. At t=108.9 s, the angle switches from 0∘ to 90∘ before immediately switching from 90∘ to −90∘ following a “bang-bang” profile. Note that it was theoretically proven in [15] that an optimal solution with such a “bang-bang” profile is reflective of the aerocapture problem.Figure 11
Bank angleσ with time t.The time histories of the constraint parameters including the load factorn, heat rate Q·, and dynamic pressure q are shown in Figure 11, in which it can be observed that all three constraints initially increase rapidly but then slowly decrease with time. The maximum values of all three constraints are within the allowable range, with a load factor of 2.5024 g, heat rate of 708.2966 kw/m2, and dynamic pressure of 5566.9255 kPa as demonstrated in Figure 12. Furthermore, the integrated heat minimized in Equation (18) by the optimization process is 123.4263 MJ/m2.Figure 12
(a) Load, (b) heat rate, and (c) dynamic pressure vs. timet.
(a)(b)(c)


## 5.1.3. Postatmospheric Trajectory

At its atmospheric exit point, the spacecraft reaches a postatmospheric orbit with an apoapsis altitudeha=343.00km and a periapsis altitude hp=78.95km. After a flight time of 1966.88 s, the spacecraft transfers from its atmospheric exit point to its apoapsis. Using Equation (29), the velocity increment that inserts the spacecraft into the station LEO is about vpost‐atm=77.66m/s, which is less than the result vpost‐atm=95.9623m/s obtained in Ref. [21]. The resulting three-dimensional trajectory of the postatmospheric orbit is shown in Figure 13.Figure 13
Three-dimensional demonstration of the three-phase trajectory near the Earth.After the transfer of the trans-Earth phase, aerocapture phase, and postatmospheric phase, the mission constraints are recalculated to determine whether the spacecraft could return to the target LEO space station. The errors of the coplanar constraint and orbital height areΔiS=0.0020∘, ΔΩS=0.0022∘, and ΔhS=0.51m, respectively, which are acceptable for the engineering requirements. It can be concluded that the proposed algorithm is effective and satisfies the nonlinear constraints in the whole mission.


## 5.2. Performance of the Hybrid Optimization Algorithms

To validate the trajectory design method based on the “perilune” and “argument” pseudo estimation rules shown in Equations (6) and (7), respectively, 1000 space station LEOs with random inclinations and RAANs were generated and optimized. The results of the subsequent convergence analysis are shown in Table 6. It can be observed that all 1000 generated orbits optimized using the proposed orbit design method based on the combined perilune and argument pseudo estimation rules converge to corresponding optimal solutions, a convergence rate of 100.0%. However, note that if no good initial value is provided to the algorithm, none of these orbits will converge. It can also be observed that the pseudo perilune rule plays a more important role than the pseudo argument of latitude rule, as they exhibit convergence rates of 87.5% and 21.9%, respectively. The results in Table 6 thus show that the proposed orbital design method provides good convergence and can effectively solve the ETO to a space station in LEO with an arbitrary orbital plane.Table 6
Convergence analysis.
Initial estimated valueOrbitConvergenceAgreementPerilune+argument10001000100.0%Perilune100087587.5%Argument100021921.9%None100000.00%The variation in the velocity incrementvpost‐atm with the reentry angle γreen can be observed in Figure 14 to be uniformly dispersed between −6.5∘ and −4.3∘. It is interesting that vpost‐atm decreases dramatically from 900 m/s to 100 m/s for γreen<−6.0∘ and that vpost‐atm is greater than 100 m/s when γreen>−4.5∘. According to the relationship between the VCP height and reentry angle shown in Figure 15, it can be determined that the preliminary VCP height corresponding to a velocity increment of less than 100 m/s is 50–80 km, which is acceptable from the perspective of the propellant consumption of the reusable PM.Figure 14
Variation in velocity increment with reentry angle.Figure 15
Relationship between VCP height and reentry angle.


## 6. Conclusions

This paper presented an entire flight trajectory design method for a reusable spacecraft flying from the Moon to an LEO space station. A hybrid optimization scheme was proposed to solve the problem by minimizing the total velocity increment from a propellant consumption perspective. It was determined that the trans-Earth phase design method is effective in satisfying the LEO coplanar constraints of the space station and that the two pseudo estimation rules evaluated in this study produced good initial guesses of the design variables. In the aerocapture phase, the maximum value of the heat rate and dynamic pressure are all within the given mission constraints, and in the postatmospheric phase, the maneuvering impulse required for the spacecraft to enter the station LEO is approximately 77.66 m/s, which is less than that of the previous study. Additionally, a robustness analysis indicated that the proposed hybrid optimization method is appropriate for an LEO space station with arbitrary orbital planes. Finally, it was found that the selection of VCP height should be limited to within 50–80 km to save propellant consumption. The conclusions of this study are intended to be used as a reference for the design of the future Chinese manned lunar mission.



---

*Source: 8876489-2021-01-09.xml*