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

Lunar exploration architecture can be made more ﬂ exible 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 satis ﬁ es 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 e ﬀ ective 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 satis ﬁ ed. 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.


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 Moonto-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-layerhybrid 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][22][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 characteristics The remainder is organized as follows. Section 2 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.

Problem Statement
The manned lunar mission architecture employing an LEO space station is shown in Figure 1 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. 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 Figure 2, the entire trajectory of this mission is separated into three phases: trans-Earth phase, aerocapture phase, and postatmospheric phase.
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 time t prl , an impulse v prl 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 height h, 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.

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  3 International Journal of Aerospace Engineering 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.

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.

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.

Moon-to-Earth Trajectory Optimization
Method. The station LEO should be near-circular with a geocentric altitude of h S , orbital inclination of i S , and RAAN of Ω S . The vacuum perigee (VCP) is defined as the perigee without considering atmospheric perturbations, as shown in Figure 3. Supposing that i vcp = i S and Ω vcp = Ω S , the coplanar constraint can be satisfied. Therefore, the design variables of the ETO are selected as where v vcp is the velocity, u vcp is the argument of latitude, Δ t peri-vcp is the transfer time, and t vcp is the time at VCP. Given these four parameters and the VCP height h vcp , the position vector r vcp and the velocity vector v vcp at VCP can be obtained by where R E is the radius of the Earth. The r vcp and v vcp 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: where r * LLO and i * LLO are the selenocentric radius and orbital inclination of LLO, respectively, r prl and i prl are the perilune altitude and inclination of ETO, respectively, and R M is the radius of the Moon.
In the trans-Earth phase, the objective function is to minimize TEI velocity increment as follows: 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 Earthcentered 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.

Trans-Earth Window Estimation by Analytical
Methods. The geocentric relationship between the Earth, the Moon, ETO, and LEO is illustrated using the ECMEPC in  International Journal of Aerospace Engineering Moon-to-Earth plane at point E and the Earth equatorial plane at point D. At TEI time t prl , 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 Similarly, for the preliminary analysis, the position of the VCP on the Moon-to-Earth plane can be estimated by 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 i L , respectively, while in the Earth-centered J2000 inertial coordinate system, the RAAN and inclination of the Moon orbit are Ω M and i M , respectively.
In the spherical triangle ΔCDE, CE = Ω M − Ω S , ∠CDE = π − i S , and ∠CED = i L . The inclination i L can then be obtained using the cosine formula as follows: Applying the sine formula to the spherical triangle ΔCDE yields By substituting i L into Equation (9), the Ω L can be obtained by In the ECMEPC, the longitude of the Moon at TEI is Therefore, the perilune time t prl can be computed by interpolation using Ω prl and the VCP point time t vcp as follows: In Figure 4, u A is expressed by DA and u L is expressed by EA, which satisfiesDA + AE = DE. Combining Equations (9) and (10), the argument of latitude at VCP is The velocity v vcp is mainly decided by the Moon-to-Earth distance d M and can be estimated by In the above analysis, the initial values of t vcp and u vcp 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.
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. The L and D can be obtained by where C D is the drag coefficient, C L 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 pressure q, 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 h min . All these constraints are formulized as where C 1 is a constant coefficient based on the PM characteristics, R d is the nose radius of the PM, ρ 0 is the density of the atmosphere at sea level, and v c is the reference velocity of a circular orbit with an altitude equal to the Earth's radius. In a lunar exploration mission, q max , n max , _ Q max , and h min are usually set to 30 kPa, 4.5 g, 2.6 MW/m 2 , and 40 km, respectively [17].
In order to prevent the spacecraft's heating shield system from being burned, the integrated heat Q S should be minimized as follows: where t 2,0 and t 2,f are the starting and ending times, respectively, of the aerocapture phase.  International Journal of Aerospace Engineering where the subscript "reen" denotes the reentry orbital states of the trans-Earth trajectory. The initial states are specified as follows.

Boundary
In the trans-Earth trajectory, the true anomaly f reen at the atmospheric reentry point is obtained by The corresponding eccentric anomaly E reen is E reen = 2 arctan The flight time between VCP and the reentry injection, Δt reen-vcp , can be obtained by where a reen is the semimajor axis of ETO and e the eccentricity. The reentry time t reen is By backward integrating ðr vcp , v vcp Þ in Equations (2) and (3), ðr reen , v reen Þ 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 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.

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 At the atmospheric exit point, let r 2,f and v 2,f be the position and velocity vectors, respectively, of the aerocapture tra-jectory. It should be noted that these two vectors are relative to the rotating Earth, and given t 2,f , they can be transformed into the inertial vectors r exit and v exit in the Earth-centered inertial coordinate system.
The semimajor axis of the postatmospheric orbit is and the apoapsis and periapsis radius is To reach the space station LEO, two engine burns are planned that are similar to a Hohmann transfer. The first tangential impulse is applied at r a to raise r p to r * p = ðr a + r S Þ/2. The second impulse is applied at r * p to circularize the orbit and bring the apoapsis radius to r S .
The total two-velocity increment is From Equation (28), it can be found that v post-atm is minimized when r a = r S , 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 apoapsistargeting problem. Therefore, the velocity increment in Equation (28) can be rewritten as Another set of boundary conditions in atmospheric flight is given by Substituting Equation (29) into Equation (25), the minimum velocity increment of the objective function can be specified.

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]. Let x i denote the optimization variable in phase i. Each subproblem is then solved in sequence as follows: where J i is the goal function, g i is the dynamic equation, eq i is the equality constraint, and ieq i 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: 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 variable X is defined to minimize the total mission objective function J as follows: 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 Section 3. 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 Δv post-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 where N is the number of maneuver times and Δv i is the velocity increment in each phase. The objective function J includes the impulse Δv peri at the perilune in the first phase and Δv post-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 height h vcp 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 h vcp 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 h vcp and γ can be estimated as 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 Figure 6. The detailed implementation process is described as follows.
Step 1. According to the ði S , Ω 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 of h vcp randomly between 0 and 100 km. 8 International Journal of Aerospace Engineering Step 3. Estimate the initial values of the design variables fv vcp , f vcp , t vcp , Δt peri-vcp g, employing the SQP algorithm to optimize the trans-Earth phase in cislunar space.
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 states ðr, θ, φ, v, γ, ψÞ 3,i and optimize the insertion velocity increment v post-atm to the target space station LEO which is similar to a Hohmann transfer.
Step 7. Apply the golden section method to optimize h vcp according to the total velocity increment J = min ðΔv prl + v post-atm Þ and repeat Steps 3-7 until the mission constraints are all satisfied and the algorithm at the top-level is converged.

Simulation
For the top level in the optimization scheme, the optimal VCP height is determined to be h  Tables 2 and 3. A three-dimensional demonstration of the trans-Earth trajectory in cislunar space is given in Figure 7.
It can be observed in the optimization results in Table 2 that the pseudo perilune longitude is λ prl = 181:96 ∘ and the pseudo argument of latitude is u L = −0:63 ∘ , which verifies the assumption of the initial values in Equations (6) and (7). The flight time Δt pri-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].  Table 4. 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.
It can be observed in Figure 8 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 Figure 10, 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.
The time histories of the constraint parameters including the load factor n, 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/m 2 , 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/m 2 .

Postatmospheric Trajectory.
At its atmospheric exit point, the spacecraft reaches a postatmospheric orbit with an apoapsis altitude h a = 343:00 km and a periapsis altitude h p = 78:95 km. 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 v post-atm = 77:66 m/s , which is less than the result v post-atm = 95:9623 m/s obtained in Ref. [21]. The resulting three-dimensional trajectory of the postatmospheric orbit is shown in Figure 13.
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 Δi S = 0:0020 ∘ , ΔΩ S = 0:0022 ∘ , and Δh S = 0:51 m, 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.

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.
The variation in the velocity increment v post-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 v post-atm decreases dramatically from 900 m/s to 100 m/s for γ reen < −6:0 ∘ and that v post-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.

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   International Journal of Aerospace Engineering 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.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.   14 International Journal of Aerospace Engineering