# A Fast Computational Method for the Minimum Duration Transfer Trajectories of Space-to-Ground Vehicles

Authors: Qingguo Liu; Xinxue Liu; Jian Wu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

On the conditions that the spacecraft engine is in finite thrust mode and the maneuver time is given, it takes a long time to compute the minimum duration transfer trajectories of space-to-ground vehicles, which is mainly because the initial values of the adjoint variables involved in the optimization model have no definite physical meanings and the model is sensitive to them. In order to develop space-to-ground transfer trajectory programmes in real time in an uncertain environment for the decision makers, we propose a fast method for computing the minimum duration transfer trajectories of space-to-ground vehicles with the given position of the landing point and the arbitrary maneuver point. First, the optimization model based on the hybrid method is established to compute the minimum duration transfer trajectory. Then, the region composed of maneuverable points is gridded and the initial values of the adjoint variables and the values of partial state variables of the minimum duration transfer trajectories at all gridded points are computed and saved to a database. Finally, the predicted values of the initial values of the adjoint variables and the values of partial state variables at any maneuver point within the region composed of maneuverable points are computed by using a binary cubic interpolation method. Finally, the minimum duration transfer trajectory is obtained by the hybrid method which takes the neighborhood of the predicted values as the search ranges of the initial values of the adjoint variables and the values of partial state variables. Simulation results demonstrate that the proposed method, which requires only 2.93% of the computational time of the hybrid method, can improve substantially the computational time of the minimum duration transfer trajectory of a space-to-ground vehicle under the guarantee of ensuring accuracy. The methodology of converting the time domain into the space domain is well applied in this paper.


---

## Body

## 1. Introduction

The transfer trajectories of spacecraft, such as return satellites, manned spacecraft, space shuttles, space-to-ground kinetic weapons, and other kinds of space-to-ground vehicles, from their orbit to the landing point is uniformly referred to as space-to-ground transfer trajectories. A typical space-to-ground transfer trajectory consists of a transition trajectory segment and a reentry trajectory segment. Depending on the resistance and lift of the spacecraft in the reentry segment, the re-entry trajectories are classified into ballistic, semiballistic, and gliding [1, 2].In the existing literature, direct methods [3–10], indirect methods [11–14], and hybrid methods [15–17] are used to study the optimization of transfer trajectories of different types of spacecraft. Space-to-ground transfer trajectory is one kind of transfer trajectories. The direct methods are computationally efficient but suffer from the lack of high precision, and the indirect methods, which are based on Pontryagin’s minimum principle [18, 19], benefit from a robust theoretical foundation but require initial values for the adjoint variables in the conjugate equations that are not easy to guess and on which the system is sensitive. The hybrid methods combine the advantages of the two previous methods and consist of converting the transfer trajectory optimization problem into a parametric optimization problem and obtaining the transfer trajectory with higher optimization index value by adjusting the initial values of the adjoint variables and the values of partial state variables. The space-to-ground transfer trajectories are studied by hybrid methods [20–22]. The advantages and disadvantages of direct, indirect, and hybrid methods are discussed in solving the transfer trajectory of the space-to-ground kinetic weapon [23]. Computing a minimum duration transfer trajectory based on the hybrid method can obtain an accurate solution, but still takes a relatively long time. In general, the region composed of maneuverable points is continuous, so the minimum duration transfer trajectories of all maneuverable points cannot be computed and saved to a database in advance. According to the research status, we propose a fast method based on the interpolation scheme and the hybrid method for computing the minimum duration transfer trajectories. The methodology of this paper is converting the time domain into the space domain. Rapid and efficient computation of space-to-ground transfer trajectories can lay the foundation for decision makers to develop space-to-ground transfer trajectory programmes in real time in an uncertain environment.We study a typical space-to-ground transfer trajectory which only takes into consideration of the energy but not the heat rate, normal load, and dynamic pressure. The typical space-to-ground transfer trajectory adopts the zero angle of attack and ballistic reentry mode (Figure1). In this paper, the main content is as follows: First, we establish the optimization model for the minimum duration transfer trajectory by using hybrid method. Then, we propose the fast method for computing the minimum duration transfer trajectory. Finally, we carry on the simulations and the result analysis to demonstrate the high computational efficiency of the proposed method.Figure 1
A sketch diagram of a space-to-ground transfer trajectory.


## 2. Optimization for the Minimum Duration Transfer Trajectory

In this section, we provide the differential equations of motion in the earth-fixed coordinate system and then establish the optimization model for the minimum duration transfer trajectory.


### 2.1. Motion Differential Equations in the Earth-Fixed Coordinate System

The differential equations of motion are [24, 25](1)dV-dt-=T-cos⁡αcos⁡β-X-m--g-rsin⁡γ-g-Φcos⁡γsin⁡ψ+ω-e2r-cos⁡Φsin⁡γcos⁡Φ-cos⁡γsin⁡Φcos⁡ψdγdt-=T-sin⁡α+Y-m-V-+V-cos⁡γr--g-rcos⁡γV-+g-Φsin⁡γsin⁡ψV-+2ω-ecos⁡Φsin⁡ψ+ω-e2r-cos⁡ΦV-cos⁡γcos⁡Φ-sin⁡γsin⁡Φcos⁡ψdψdt-=T-cos⁡αsin⁡β-Z-m-V-cos⁡γ-V-cos⁡γcos⁡ψtan⁡Φr--g-θcos⁡ψV-cos⁡γ-2ω-ecos⁡Φcos⁡ψtan⁡γ-sin⁡Φ+ω-e2r-V-cos⁡γsin⁡ψsin⁡Φcos⁡Φdr-dt-=V-sin⁡γdΘdt-=V-cos⁡γcos⁡ψr-cos⁡ΦdΦdt-=V-cos⁡γsin⁡ψr-dm-dt-=-T-V-ewhere V- is the dimensionless velocity, γ is the velocity inclination angle, ψ is the course angle (the angle between the projection of the velocity vector on the local horizontal plane and the latitude tangent), r- is the dimensionless geocentric distance, Θ is the longitude, Φ is the latitude, T- is the engine thrust, V-e is the gas jet velocity, m- is the dimensionless spacecraft quality, t- is the dimensionless time, X- is the dimensionless resistance, Y- is the dimensionless lift, Z- is the dimensionless lateral force, α is the angle of attack, β is the sideslip angle, ω-e is the dimensionless earth rotation angular rate, and g-r and g-Φ are the dimensionless gravitational components when only the first three terms of the spherical harmonic expansion are considered. In the transition trajectory segment, the values of X-, Y-, and Z- are zero. In the reentry trajectory segment, the value of T- is zero and the atmospheric model is US standard atmosphere (1976).The equations for the dimensionless parameters are(2)V-=VVref,r-=rrref,T-=Tmrefgref,V-e=VeVref,m-=mmreft-=ttref,X-=Xmrefgref,Y-=Ymrefgref,Z-=Zmrefgrefω-e=ωetref,g-r=1r-21-3J22r-23sin2⁡Φ-1,g-Φ=3J2sin⁡2Φ2r-4where J2 is the coefficient of second-order principal spherical harmonic function (the value of J2 is 1.08263e-3.); and rref, mref, Vref, tref, and gref are(3)rref=RE,mref=m0,Vref=μrreftref=rrefVref,gref=μr2refwhere RE (the value of RE is 6.37100e+6 m.) is the reference radius of the earth, μ is the gravitation constant (the value of μ is 3.986005e+14m3/s2), and m0 is the initial mass of the spacecraft.


### 2.2. Minimum Duration Transfer Trajectory Optimization Model Based on the Hybrid Method

The minimum duration transfer trajectory optimization problem based on the hybrid method is a high-dimensional nonlinear optimization problem, which is hard to obtain global optimal solutions. GA has strong global search ability and simple implementation steps, which is suitable for solving the minimum duration transfer trajectory optimization problem. The main idea of the hybrid method for computing the minimum duration transfer trajectory can be summarized as follows. The optimization problem of the transition trajectory segment is converted into a two-point boundary value problem according to Pontryagin’s minimum principle and the reentry trajectory segment needs to satisfy constraints. Then, the initial values of the adjoint variables and the values of partial state variables are adjusted by means of the genetic algorithm (GA) [26–28]. At last, we obtain the minimum duration transfer trajectory.The minimum duration optimization index is written as(4)minJ=t-=t-a+t-bwhere t-, t-a, and t-b are the durations to complete the transfer trajectory, the transition trajectory segment, and the reentry trajectory segment of the vehicle.According to Pontryagin’s minimum principle, the Hamiltonian function is given by(5)H=λV-V-˙+λγγ˙+λψψ˙+λr-r-˙+λΘΘ˙+λΦΦ˙+λm-m-˙From (5), the covariate variables satisfy the differential equations (only the transition trajectory with thrust control is considered):(6)
d
λ
V
-
d
t
-
=
λ
γ
g
-
Φ
sin
⁡
γ
sin
⁡
ψ
V
-
2
-
1
r
-
+
g
-
r
V
-
2
cos
⁡
γ
+
ω
-
e
2
r
-
cos
⁡
Φ
V
-
2
cos
⁡
γ
cos
⁡
Φ
-
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
T
-
sin
⁡
α
m
-
V
-
2
-
λ
Θ
cos
⁡
γ
cos
⁡
ψ
r
-
cos
⁡
Φ
-
λ
r
-
sin
⁡
γ
-
λ
Φ
cos
⁡
γ
sin
⁡
ψ
r
-
+
λ
ψ
T
-
cos
⁡
α
sin
⁡
β
-
Z
-
m
-
V
-
2
cos
⁡
γ
+
cos
⁡
γ
cos
⁡
ψ
tan
⁡
Φ
r
-
-
g
-
Φ
cos
⁡
ψ
V
-
2
cos
⁡
γ
+
ω
-
e
2
r
-
V
-
2
cos
⁡
γ
sin
⁡
ψ
sin
⁡
Φ
cos
⁡
Φ
d
λ
γ
d
t
-
=
λ
V
-
g
-
r
cos
⁡
γ
-
g
-
Φ
sin
⁡
γ
sin
⁡
ψ
-
ω
-
e
2
r
-
cos
⁡
Φ
cos
⁡
γ
cos
⁡
Φ
+
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
λ
γ
V
-
sin
⁡
γ
r
-
-
g
-
r
sin
⁡
γ
+
g
-
Φ
cos
⁡
γ
sin
⁡
ψ
V
-
+
ω
-
e
2
r
-
cos
⁡
Φ
V
-
sin
⁡
γ
cos
⁡
Φ
+
cos
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
λ
Θ
V
-
sin
⁡
γ
cos
⁡
ψ
r
-
cos
⁡
Φ
+
λ
Φ
V
-
sin
⁡
γ
sin
⁡
ψ
r
-
-
λ
ψ
T
-
cos
⁡
α
sin
⁡
β
tan
⁡
γ
m
-
V
-
cos
⁡
γ
+
V
-
sin
⁡
γ
cos
⁡
ψ
tan
⁡
Φ
r
-
-
g
-
Φ
tan
⁡
γ
cos
⁡
ψ
V
-
cos
⁡
γ
-
2
ω
-
e
cos
⁡
Φ
cos
⁡
ψ
s
e
c
2
γ
+
ω
-
e
2
r
-
tan
⁡
γ
V
-
cos
⁡
γ
sin
⁡
ψ
sin
⁡
Φ
cos
⁡
Φ
-
λ
r
-
V
-
cos
⁡
γ
d
λ
ψ
d
t
-
=
λ
V
-
g
-
Φ
cos
⁡
γ
cos
⁡
ψ
-
ω
-
e
2
r
-
cos
⁡
Φ
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
λ
Θ
V
-
cos
⁡
γ
sin
⁡
ψ
r
-
cos
⁡
Φ
-
λ
γ
g
-
Φ
sin
⁡
γ
cos
⁡
ψ
V
-
+
2
ω
-
e
cos
⁡
Φ
cos
⁡
ψ
+
ω
-
e
2
r
-
cos
⁡
Φ
V
-
sin
⁡
γ
sin
⁡
Φ
sin
⁡
ψ
-
λ
Φ
V
-
cos
⁡
γ
cos
⁡
ψ
r
-
-
λ
ψ
V
-
cos
⁡
γ
sin
⁡
ψ
tan
⁡
Φ
r
-
+
g
-
Φ
sin
⁡
ψ
V
-
cos
⁡
γ
+
2
ω
-
e
cos
⁡
Φ
sin
⁡
ψ
tan
⁡
γ
+
ω
-
e
2
r
-
V
-
cos
⁡
γ
cos
⁡
ψ
sin
⁡
Φ
cos
⁡
Φ
d
λ
r
-
d
t
-
=
-
λ
V
-
g
-
r
r
sin
⁡
γ
+
g
-
Φ
r
cos
⁡
γ
sin
⁡
ψ
+
ω
-
e
2
cos
⁡
Φ
sin
⁡
γ
cos
⁡
Φ
-
cos
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
-
λ
γ
g
-
r
r
cos
⁡
γ
V
-
-
g
-
Φ
r
sin
⁡
γ
sin
⁡
ψ
V
-
+
ω
-
e
2
cos
⁡
Φ
V
-
cos
⁡
γ
cos
⁡
Φ
-
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
λ
Θ
V
-
cos
⁡
γ
sin
⁡
ψ
r
-
2
cos
⁡
Φ
+
λ
Φ
V
-
cos
⁡
γ
cos
⁡
ψ
r
-
2
-
λ
ψ
V
-
sin
⁡
γ
cos
⁡
ψ
tan
⁡
Φ
r
-
2
+
g
-
Φ
r
cos
⁡
ψ
V
-
cos
⁡
γ
+
ω
-
e
2
V
-
cos
⁡
γ
sin
⁡
ψ
sin
⁡
Φ
cos
⁡
Φ
d
λ
Θ
d
t
-
=
0
d
λ
Φ
d
t
-
=
-
λ
V
-
g
-
r
Φ
sin
⁡
γ
+
g
-
Φ
Φ
cos
⁡
γ
sin
⁡
ψ
-
ω
-
e
2
r
-
sin
⁡
Φ
sin
⁡
γ
cos
⁡
Φ
-
cos
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
-
ω
-
e
2
r
-
cos
⁡
Φ
sin
⁡
γ
sin
⁡
Φ
+
cos
⁡
γ
cos
⁡
Φ
cos
⁡
ψ
-
λ
γ
g
-
r
Φ
cos
⁡
γ
-
g
-
Φ
Φ
sin
⁡
γ
sin
⁡
ψ
V
-
+
2
ω
-
e
sin
⁡
Φ
sin
⁡
ψ
-
ω
-
e
2
r
-
sin
⁡
Φ
V
-
cos
⁡
γ
cos
⁡
Φ
-
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
-
ω
-
e
2
r
-
cos
⁡
Φ
V
-
cos
⁡
γ
sin
⁡
Φ
+
sin
⁡
γ
cos
⁡
Φ
cos
⁡
ψ
+
λ
ψ
V
-
cos
⁡
γ
cos
⁡
ψ
r
-
c
o
s
2
Φ
-
g
-
Φ
Φ
cos
⁡
ψ
V
-
cos
⁡
γ
-
2
ω
-
e
sin
⁡
Φ
cos
⁡
ψ
tan
⁡
γ
-
cos
⁡
Φ
+
ω
-
e
2
r
-
V
-
cos
⁡
γ
sin
⁡
ψ
s
i
n
2
Φ
-
c
o
s
2
Φ
-
λ
Θ
V
-
cos
⁡
γ
cos
⁡
ψ
tan
⁡
Φ
r
-
cos
⁡
Φ
d
λ
m
-
d
t
-
=
T
-
m
-
2
λ
V
-
cos
⁡
α
cos
⁡
β
+
λ
γ
sin
⁡
α
V
-
+
λ
ψ
cos
⁡
α
sin
⁡
β
V
-
cos
⁡
γwhere the expressions of g-rr, g-Φr, g-rΦ, and g-ΦΦ are given by(7)g-rr=-2g-rr-+3J2r-53sin2⁡Φ-1,g-Φr=-6J2sin⁡2Φr-5g-rΦ=6J2sin⁡2Φ2r-4,g-ΦΦ=3J2cos⁡2Φr-4From the sufficient conditions of optimality and formula (8), the optimal thrust directions are obtained by(8)∂H∂α=0=-λV-T-sin⁡αcos⁡βm-+λγT-cos⁡αm-V--λψT-sin⁡αsin⁡βm-V-cos⁡γ∂H∂β=0=λV-T-cos⁡αsin⁡βm-+λψT-cos⁡αcos⁡βm-V-cos⁡γ(9)α=arctan⁡sgnλV-·λγcos⁡γλV-V-cos⁡γ2+λψ2β=arctan⁡λψλV-V-cos⁡γwhere sgn() is the sign function.The minimum duration transfer trajectory consists of the transition trajectory segment and the reentry trajectory segment. The dimensionless geocentric distancer-a at the reentry point is a constant. The dimensionless velocity V-a and the velocity inclination angle γa at the reentry point are two key parameters that affect the reentry heat flow, dynamic pressure, fuel consumption, and the duration. Thus, we restrict the search ranges of parameters V- and γ, which also can improve the computational efficiency. Therefore, the terminal constraints and the cross-section constraints are shown in (10), (11), and (12).The transition trajectory satisfies initial boundary conditions:(10)V-t-0-V-0=0,γt-0-γ0=0,r-t-0-r-0=0,ψt-0-ψ0=0Θt-0-Θ0=0,Φt-0-Φ0=0,m-t-0-m-0=0where V-0, γ0, ψ0, r-0, Θ0, Φ0, and m-0 are the dimensionless velocity, the velocity inclination angle, the course angle, the dimensionless geocentric distance, the longitude, the latitude, and the dimensionless spacecraft quality at the maneuver point of transition trajectory, respectively.The transition trajectory satisfies the terminal boundary constraints:(11)V-t-a-V-a=0γt-a-γa=0r-t-a-r-a=0Therefore, the cross-section conditions are given by(12)λψt-a=0,λΘt-a=0λΦt-a=0,λm-t-a=0The reentry trajectory satisfies the constraints(13)λΘt-b=Θb,λΦt-b=Φbr-t-b=RErref=1Because the index includes time, we have(14)Ht-0=Ht-=Ht-a=-1where V-a, γa, and r-a are the dimensionless speed, velocity inclination angle, and geocentric distance at the terminal of transition trajectory, respectively; Θb and Φb are the longitude and latitude of the landing point, respectively; and t-0 is the dimensionless time at the maneuver point.In order to achieve the minimum index under the given constraints, eleven parameters are taken as optimization variables including the initial values of seven adjoint variables (λV-(t-0),λγ(t-0),λψ(t-0),λr-(t-0),λΘ(t-0),λΦ(t-0), and λm-(t-0)) and the values of four state variables (V-a, γa, r-a, and t-a). Since the initial value of an adjoint variable can be obtained from (5) and (13), and r-a=RE+120000/RE (the height of reentry point is 120km), only nine variables need to be optimized. A fourth-order Runge-Kutta method is used to compute the starts of formula (1) and (6) integral equations and the Adams predictor-corrector method is adopted to compute the remaining integral equations. The implementation steps of the GA are as follows.Step 1.
λ V - ( t - 0 ), λγ(t-0), λψ(t-0), λr-(t-0), λΘ(t-0), λm-(t-0), V-a, γa, and t-a are encoded. Twenty chromosomes are generated randomly as the initial population. The crossover probability Pc, the mutation probability Pm, and the maximum number of iterations N are set.Step 2.
The fitness function is shown in formula (4). When a chromosome is determined, formulas (6) and (8) are integrated to get the optimal thrust directions under the constraints (11), (12), and (13). And then formula (1) is integrated to obtain the transfer trajectory. The fitness values of all the chromosomes in the current generation are computed.Step 3.
The chromosomes are selected using roulette wheel selection. Some genes on two different chromosomes reciprocally cross according to the crossover probability and others mutate according to the mutation probability. The execution of selection, crossover, and mutation leads to the next generation population.Step 4.
Steps2 and 3 are repeatedly executed until the GA converges ε=10-6 or the maximum number of iterations N is reached.The minimum duration transfer trajectory at a given maneuver point can be computed through these steps.


## 2.1. Motion Differential Equations in the Earth-Fixed Coordinate System

The differential equations of motion are [24, 25](1)dV-dt-=T-cos⁡αcos⁡β-X-m--g-rsin⁡γ-g-Φcos⁡γsin⁡ψ+ω-e2r-cos⁡Φsin⁡γcos⁡Φ-cos⁡γsin⁡Φcos⁡ψdγdt-=T-sin⁡α+Y-m-V-+V-cos⁡γr--g-rcos⁡γV-+g-Φsin⁡γsin⁡ψV-+2ω-ecos⁡Φsin⁡ψ+ω-e2r-cos⁡ΦV-cos⁡γcos⁡Φ-sin⁡γsin⁡Φcos⁡ψdψdt-=T-cos⁡αsin⁡β-Z-m-V-cos⁡γ-V-cos⁡γcos⁡ψtan⁡Φr--g-θcos⁡ψV-cos⁡γ-2ω-ecos⁡Φcos⁡ψtan⁡γ-sin⁡Φ+ω-e2r-V-cos⁡γsin⁡ψsin⁡Φcos⁡Φdr-dt-=V-sin⁡γdΘdt-=V-cos⁡γcos⁡ψr-cos⁡ΦdΦdt-=V-cos⁡γsin⁡ψr-dm-dt-=-T-V-ewhere V- is the dimensionless velocity, γ is the velocity inclination angle, ψ is the course angle (the angle between the projection of the velocity vector on the local horizontal plane and the latitude tangent), r- is the dimensionless geocentric distance, Θ is the longitude, Φ is the latitude, T- is the engine thrust, V-e is the gas jet velocity, m- is the dimensionless spacecraft quality, t- is the dimensionless time, X- is the dimensionless resistance, Y- is the dimensionless lift, Z- is the dimensionless lateral force, α is the angle of attack, β is the sideslip angle, ω-e is the dimensionless earth rotation angular rate, and g-r and g-Φ are the dimensionless gravitational components when only the first three terms of the spherical harmonic expansion are considered. In the transition trajectory segment, the values of X-, Y-, and Z- are zero. In the reentry trajectory segment, the value of T- is zero and the atmospheric model is US standard atmosphere (1976).The equations for the dimensionless parameters are(2)V-=VVref,r-=rrref,T-=Tmrefgref,V-e=VeVref,m-=mmreft-=ttref,X-=Xmrefgref,Y-=Ymrefgref,Z-=Zmrefgrefω-e=ωetref,g-r=1r-21-3J22r-23sin2⁡Φ-1,g-Φ=3J2sin⁡2Φ2r-4where J2 is the coefficient of second-order principal spherical harmonic function (the value of J2 is 1.08263e-3.); and rref, mref, Vref, tref, and gref are(3)rref=RE,mref=m0,Vref=μrreftref=rrefVref,gref=μr2refwhere RE (the value of RE is 6.37100e+6 m.) is the reference radius of the earth, μ is the gravitation constant (the value of μ is 3.986005e+14m3/s2), and m0 is the initial mass of the spacecraft.


## 2.2. Minimum Duration Transfer Trajectory Optimization Model Based on the Hybrid Method

The minimum duration transfer trajectory optimization problem based on the hybrid method is a high-dimensional nonlinear optimization problem, which is hard to obtain global optimal solutions. GA has strong global search ability and simple implementation steps, which is suitable for solving the minimum duration transfer trajectory optimization problem. The main idea of the hybrid method for computing the minimum duration transfer trajectory can be summarized as follows. The optimization problem of the transition trajectory segment is converted into a two-point boundary value problem according to Pontryagin’s minimum principle and the reentry trajectory segment needs to satisfy constraints. Then, the initial values of the adjoint variables and the values of partial state variables are adjusted by means of the genetic algorithm (GA) [26–28]. At last, we obtain the minimum duration transfer trajectory.The minimum duration optimization index is written as(4)minJ=t-=t-a+t-bwhere t-, t-a, and t-b are the durations to complete the transfer trajectory, the transition trajectory segment, and the reentry trajectory segment of the vehicle.According to Pontryagin’s minimum principle, the Hamiltonian function is given by(5)H=λV-V-˙+λγγ˙+λψψ˙+λr-r-˙+λΘΘ˙+λΦΦ˙+λm-m-˙From (5), the covariate variables satisfy the differential equations (only the transition trajectory with thrust control is considered):(6)
d
λ
V
-
d
t
-
=
λ
γ
g
-
Φ
sin
⁡
γ
sin
⁡
ψ
V
-
2
-
1
r
-
+
g
-
r
V
-
2
cos
⁡
γ
+
ω
-
e
2
r
-
cos
⁡
Φ
V
-
2
cos
⁡
γ
cos
⁡
Φ
-
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
T
-
sin
⁡
α
m
-
V
-
2
-
λ
Θ
cos
⁡
γ
cos
⁡
ψ
r
-
cos
⁡
Φ
-
λ
r
-
sin
⁡
γ
-
λ
Φ
cos
⁡
γ
sin
⁡
ψ
r
-
+
λ
ψ
T
-
cos
⁡
α
sin
⁡
β
-
Z
-
m
-
V
-
2
cos
⁡
γ
+
cos
⁡
γ
cos
⁡
ψ
tan
⁡
Φ
r
-
-
g
-
Φ
cos
⁡
ψ
V
-
2
cos
⁡
γ
+
ω
-
e
2
r
-
V
-
2
cos
⁡
γ
sin
⁡
ψ
sin
⁡
Φ
cos
⁡
Φ
d
λ
γ
d
t
-
=
λ
V
-
g
-
r
cos
⁡
γ
-
g
-
Φ
sin
⁡
γ
sin
⁡
ψ
-
ω
-
e
2
r
-
cos
⁡
Φ
cos
⁡
γ
cos
⁡
Φ
+
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
λ
γ
V
-
sin
⁡
γ
r
-
-
g
-
r
sin
⁡
γ
+
g
-
Φ
cos
⁡
γ
sin
⁡
ψ
V
-
+
ω
-
e
2
r
-
cos
⁡
Φ
V
-
sin
⁡
γ
cos
⁡
Φ
+
cos
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
λ
Θ
V
-
sin
⁡
γ
cos
⁡
ψ
r
-
cos
⁡
Φ
+
λ
Φ
V
-
sin
⁡
γ
sin
⁡
ψ
r
-
-
λ
ψ
T
-
cos
⁡
α
sin
⁡
β
tan
⁡
γ
m
-
V
-
cos
⁡
γ
+
V
-
sin
⁡
γ
cos
⁡
ψ
tan
⁡
Φ
r
-
-
g
-
Φ
tan
⁡
γ
cos
⁡
ψ
V
-
cos
⁡
γ
-
2
ω
-
e
cos
⁡
Φ
cos
⁡
ψ
s
e
c
2
γ
+
ω
-
e
2
r
-
tan
⁡
γ
V
-
cos
⁡
γ
sin
⁡
ψ
sin
⁡
Φ
cos
⁡
Φ
-
λ
r
-
V
-
cos
⁡
γ
d
λ
ψ
d
t
-
=
λ
V
-
g
-
Φ
cos
⁡
γ
cos
⁡
ψ
-
ω
-
e
2
r
-
cos
⁡
Φ
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
λ
Θ
V
-
cos
⁡
γ
sin
⁡
ψ
r
-
cos
⁡
Φ
-
λ
γ
g
-
Φ
sin
⁡
γ
cos
⁡
ψ
V
-
+
2
ω
-
e
cos
⁡
Φ
cos
⁡
ψ
+
ω
-
e
2
r
-
cos
⁡
Φ
V
-
sin
⁡
γ
sin
⁡
Φ
sin
⁡
ψ
-
λ
Φ
V
-
cos
⁡
γ
cos
⁡
ψ
r
-
-
λ
ψ
V
-
cos
⁡
γ
sin
⁡
ψ
tan
⁡
Φ
r
-
+
g
-
Φ
sin
⁡
ψ
V
-
cos
⁡
γ
+
2
ω
-
e
cos
⁡
Φ
sin
⁡
ψ
tan
⁡
γ
+
ω
-
e
2
r
-
V
-
cos
⁡
γ
cos
⁡
ψ
sin
⁡
Φ
cos
⁡
Φ
d
λ
r
-
d
t
-
=
-
λ
V
-
g
-
r
r
sin
⁡
γ
+
g
-
Φ
r
cos
⁡
γ
sin
⁡
ψ
+
ω
-
e
2
cos
⁡
Φ
sin
⁡
γ
cos
⁡
Φ
-
cos
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
-
λ
γ
g
-
r
r
cos
⁡
γ
V
-
-
g
-
Φ
r
sin
⁡
γ
sin
⁡
ψ
V
-
+
ω
-
e
2
cos
⁡
Φ
V
-
cos
⁡
γ
cos
⁡
Φ
-
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
+
λ
Θ
V
-
cos
⁡
γ
sin
⁡
ψ
r
-
2
cos
⁡
Φ
+
λ
Φ
V
-
cos
⁡
γ
cos
⁡
ψ
r
-
2
-
λ
ψ
V
-
sin
⁡
γ
cos
⁡
ψ
tan
⁡
Φ
r
-
2
+
g
-
Φ
r
cos
⁡
ψ
V
-
cos
⁡
γ
+
ω
-
e
2
V
-
cos
⁡
γ
sin
⁡
ψ
sin
⁡
Φ
cos
⁡
Φ
d
λ
Θ
d
t
-
=
0
d
λ
Φ
d
t
-
=
-
λ
V
-
g
-
r
Φ
sin
⁡
γ
+
g
-
Φ
Φ
cos
⁡
γ
sin
⁡
ψ
-
ω
-
e
2
r
-
sin
⁡
Φ
sin
⁡
γ
cos
⁡
Φ
-
cos
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
-
ω
-
e
2
r
-
cos
⁡
Φ
sin
⁡
γ
sin
⁡
Φ
+
cos
⁡
γ
cos
⁡
Φ
cos
⁡
ψ
-
λ
γ
g
-
r
Φ
cos
⁡
γ
-
g
-
Φ
Φ
sin
⁡
γ
sin
⁡
ψ
V
-
+
2
ω
-
e
sin
⁡
Φ
sin
⁡
ψ
-
ω
-
e
2
r
-
sin
⁡
Φ
V
-
cos
⁡
γ
cos
⁡
Φ
-
sin
⁡
γ
sin
⁡
Φ
cos
⁡
ψ
-
ω
-
e
2
r
-
cos
⁡
Φ
V
-
cos
⁡
γ
sin
⁡
Φ
+
sin
⁡
γ
cos
⁡
Φ
cos
⁡
ψ
+
λ
ψ
V
-
cos
⁡
γ
cos
⁡
ψ
r
-
c
o
s
2
Φ
-
g
-
Φ
Φ
cos
⁡
ψ
V
-
cos
⁡
γ
-
2
ω
-
e
sin
⁡
Φ
cos
⁡
ψ
tan
⁡
γ
-
cos
⁡
Φ
+
ω
-
e
2
r
-
V
-
cos
⁡
γ
sin
⁡
ψ
s
i
n
2
Φ
-
c
o
s
2
Φ
-
λ
Θ
V
-
cos
⁡
γ
cos
⁡
ψ
tan
⁡
Φ
r
-
cos
⁡
Φ
d
λ
m
-
d
t
-
=
T
-
m
-
2
λ
V
-
cos
⁡
α
cos
⁡
β
+
λ
γ
sin
⁡
α
V
-
+
λ
ψ
cos
⁡
α
sin
⁡
β
V
-
cos
⁡
γwhere the expressions of g-rr, g-Φr, g-rΦ, and g-ΦΦ are given by(7)g-rr=-2g-rr-+3J2r-53sin2⁡Φ-1,g-Φr=-6J2sin⁡2Φr-5g-rΦ=6J2sin⁡2Φ2r-4,g-ΦΦ=3J2cos⁡2Φr-4From the sufficient conditions of optimality and formula (8), the optimal thrust directions are obtained by(8)∂H∂α=0=-λV-T-sin⁡αcos⁡βm-+λγT-cos⁡αm-V--λψT-sin⁡αsin⁡βm-V-cos⁡γ∂H∂β=0=λV-T-cos⁡αsin⁡βm-+λψT-cos⁡αcos⁡βm-V-cos⁡γ(9)α=arctan⁡sgnλV-·λγcos⁡γλV-V-cos⁡γ2+λψ2β=arctan⁡λψλV-V-cos⁡γwhere sgn() is the sign function.The minimum duration transfer trajectory consists of the transition trajectory segment and the reentry trajectory segment. The dimensionless geocentric distancer-a at the reentry point is a constant. The dimensionless velocity V-a and the velocity inclination angle γa at the reentry point are two key parameters that affect the reentry heat flow, dynamic pressure, fuel consumption, and the duration. Thus, we restrict the search ranges of parameters V- and γ, which also can improve the computational efficiency. Therefore, the terminal constraints and the cross-section constraints are shown in (10), (11), and (12).The transition trajectory satisfies initial boundary conditions:(10)V-t-0-V-0=0,γt-0-γ0=0,r-t-0-r-0=0,ψt-0-ψ0=0Θt-0-Θ0=0,Φt-0-Φ0=0,m-t-0-m-0=0where V-0, γ0, ψ0, r-0, Θ0, Φ0, and m-0 are the dimensionless velocity, the velocity inclination angle, the course angle, the dimensionless geocentric distance, the longitude, the latitude, and the dimensionless spacecraft quality at the maneuver point of transition trajectory, respectively.The transition trajectory satisfies the terminal boundary constraints:(11)V-t-a-V-a=0γt-a-γa=0r-t-a-r-a=0Therefore, the cross-section conditions are given by(12)λψt-a=0,λΘt-a=0λΦt-a=0,λm-t-a=0The reentry trajectory satisfies the constraints(13)λΘt-b=Θb,λΦt-b=Φbr-t-b=RErref=1Because the index includes time, we have(14)Ht-0=Ht-=Ht-a=-1where V-a, γa, and r-a are the dimensionless speed, velocity inclination angle, and geocentric distance at the terminal of transition trajectory, respectively; Θb and Φb are the longitude and latitude of the landing point, respectively; and t-0 is the dimensionless time at the maneuver point.In order to achieve the minimum index under the given constraints, eleven parameters are taken as optimization variables including the initial values of seven adjoint variables (λV-(t-0),λγ(t-0),λψ(t-0),λr-(t-0),λΘ(t-0),λΦ(t-0), and λm-(t-0)) and the values of four state variables (V-a, γa, r-a, and t-a). Since the initial value of an adjoint variable can be obtained from (5) and (13), and r-a=RE+120000/RE (the height of reentry point is 120km), only nine variables need to be optimized. A fourth-order Runge-Kutta method is used to compute the starts of formula (1) and (6) integral equations and the Adams predictor-corrector method is adopted to compute the remaining integral equations. The implementation steps of the GA are as follows.Step 1.
λ V - ( t - 0 ), λγ(t-0), λψ(t-0), λr-(t-0), λΘ(t-0), λm-(t-0), V-a, γa, and t-a are encoded. Twenty chromosomes are generated randomly as the initial population. The crossover probability Pc, the mutation probability Pm, and the maximum number of iterations N are set.Step 2.
The fitness function is shown in formula (4). When a chromosome is determined, formulas (6) and (8) are integrated to get the optimal thrust directions under the constraints (11), (12), and (13). And then formula (1) is integrated to obtain the transfer trajectory. The fitness values of all the chromosomes in the current generation are computed.Step 3.
The chromosomes are selected using roulette wheel selection. Some genes on two different chromosomes reciprocally cross according to the crossover probability and others mutate according to the mutation probability. The execution of selection, crossover, and mutation leads to the next generation population.Step 4.
Steps2 and 3 are repeatedly executed until the GA converges ε=10-6 or the maximum number of iterations N is reached.The minimum duration transfer trajectory at a given maneuver point can be computed through these steps.


## 3. The Fast Computational Method

Given the landing point position, the minimum duration transfer trajectory at any one maneuver point can be obtained by the hybrid method, but the computational time is still too long for real-time applications. On the basis of the characteristics of the problem and the hybrid method, the region composed of maneuverable points is gridded. Then, the initial values of the adjoint variables and the values of partial state variables of the minimum duration transfer trajectory at all gridded points are computed and saved to a database. The predicted values of the initial values of the adjoint variables and the values of partial state variable at any one maneuver point in the region composed of maneuverable points are computed by using a binary cubic interpolation method [29, 30], and the minimum duration transfer trajectory is computed by the hybrid method which takes neighborhood of predicted values as the search ranges of the initial values of the adjoint variables and the values of partial state variables. The proposed method can be summarized as combining an interpolation scheme with hybrid method.


### 3.1. Discrete Grids

After analysing a large number of computational results, we find that the initial values of the adjoint variables and the values of partial state variables of the minimum duration transfer trajectories have characteristics of regular change. Therefore, we come up with the idea of gridding the region composed of maneuverable points and the sizes of discrete grids are uniform. The region composed of maneuverable points is gridded according to the rules described below. The initial values of the adjoint variables and the values of partial state variables of the minimum duration transfer trajectory of at all gridded points are computed and saved to a database.Figures2, 3, and 4 show the track of subsatellite points after the spacecraft runs during one day, five days, and twenty days. It can be seen that as the spacecraft running time gets longer, the track of subsatellite points gets denser. In general, the track of subsatellite points does not overlap during the least common multiple of the spacecraft period of revolution and the earth rotation period. The longitude Θ and latitude Φ of the track of subsatellite points are taken as the grid discrete parameters. The discrete granularity of the longitude is set as Θd, and the grid points are the intersections of the set longitude and n tracks of subsatellite points. When T-=0, X-=0, Y-=0, and Z-=0 are set in formula (1) and an arbitrary longitude value is given, the latitude value at the grid point can be obtained.Figure 2
Track of sub-satellite points after the spacecraft runs one day.Figure 3
Track of sub-satellite points after the spacecraft runs five days.Figure 4
Track of sub-satellite points after the spacecraft runs twenty days.After the region (shown as Figure5) composed of maneuverable points is gridded, Dsum gridded points are obtained. The function value at the gridded point (Θi,Φi) (Θi is the ith set longitude and Φi is latitude value of the jth track of subsatellite) is the vector Z(Θi,Φi):(15)ZΘi,Φi=λV-t-0,λγt-0,λψt-0,λr-t-0,λΘt-0,λm-t-0,V-a,γa,t-aΘi,ΦiFigure 5
Sketch diagram of the maneuverable region of the spacecraft.


### 3.2. The Fast Computational Method Based on Binary Cubic Interpolation

The predicted values of the initial values of the adjoint variables and the values of partial state variables at any maneuver point in the region composed of maneuverable points are computed using a binary cubic interpolation method. A neighborhoodδ(k)(k=1,2,⋯,9) of the predicted values is taken as search ranges of the initial values of the adjoint variables and the values of partial state variables, and the minimum duration transfer trajectory is computed by the hybrid method.Given the time at any maneuver point, the initial values of the adjoint variables and the values of partial state variables at the ten points which are the nearest to the point(u,v) are inserted in formula (14) to obtain the ten constant coefficients ak0,ak1,⋯,ak10(k=1,2,⋯,9). The interpolated values (the predicted values) of the initial values of the adjoint variables and the values of partial state variables at the point (u,v) can be computed by using I(k)(u,v):(16)Iku,v=ak0+ak1u+ak2v+ak3u2+ak4v2+ak5u3+ak6v3+ak7uv+ak8uv2+ak9u2vk=1,2,⋯,9where u,v are the longitude and the latitude of the maneuver point, respectively, and I(k)(u,v) is the kth value, which is obtained by interpolation, in the vector Z(Θi,Φi).The neighborhoodsδ(k)(k=1,2,⋯,9) of the initial values of the adjoint variables and the values of partial state variables are taken as search ranges of optimization variables, and the minimum duration transfer trajectory is computed by the hybrid method.


## 3.1. Discrete Grids

After analysing a large number of computational results, we find that the initial values of the adjoint variables and the values of partial state variables of the minimum duration transfer trajectories have characteristics of regular change. Therefore, we come up with the idea of gridding the region composed of maneuverable points and the sizes of discrete grids are uniform. The region composed of maneuverable points is gridded according to the rules described below. The initial values of the adjoint variables and the values of partial state variables of the minimum duration transfer trajectory of at all gridded points are computed and saved to a database.Figures2, 3, and 4 show the track of subsatellite points after the spacecraft runs during one day, five days, and twenty days. It can be seen that as the spacecraft running time gets longer, the track of subsatellite points gets denser. In general, the track of subsatellite points does not overlap during the least common multiple of the spacecraft period of revolution and the earth rotation period. The longitude Θ and latitude Φ of the track of subsatellite points are taken as the grid discrete parameters. The discrete granularity of the longitude is set as Θd, and the grid points are the intersections of the set longitude and n tracks of subsatellite points. When T-=0, X-=0, Y-=0, and Z-=0 are set in formula (1) and an arbitrary longitude value is given, the latitude value at the grid point can be obtained.Figure 2
Track of sub-satellite points after the spacecraft runs one day.Figure 3
Track of sub-satellite points after the spacecraft runs five days.Figure 4
Track of sub-satellite points after the spacecraft runs twenty days.After the region (shown as Figure5) composed of maneuverable points is gridded, Dsum gridded points are obtained. The function value at the gridded point (Θi,Φi) (Θi is the ith set longitude and Φi is latitude value of the jth track of subsatellite) is the vector Z(Θi,Φi):(15)ZΘi,Φi=λV-t-0,λγt-0,λψt-0,λr-t-0,λΘt-0,λm-t-0,V-a,γa,t-aΘi,ΦiFigure 5
Sketch diagram of the maneuverable region of the spacecraft.


## 3.2. The Fast Computational Method Based on Binary Cubic Interpolation

The predicted values of the initial values of the adjoint variables and the values of partial state variables at any maneuver point in the region composed of maneuverable points are computed using a binary cubic interpolation method. A neighborhoodδ(k)(k=1,2,⋯,9) of the predicted values is taken as search ranges of the initial values of the adjoint variables and the values of partial state variables, and the minimum duration transfer trajectory is computed by the hybrid method.Given the time at any maneuver point, the initial values of the adjoint variables and the values of partial state variables at the ten points which are the nearest to the point(u,v) are inserted in formula (14) to obtain the ten constant coefficients ak0,ak1,⋯,ak10(k=1,2,⋯,9). The interpolated values (the predicted values) of the initial values of the adjoint variables and the values of partial state variables at the point (u,v) can be computed by using I(k)(u,v):(16)Iku,v=ak0+ak1u+ak2v+ak3u2+ak4v2+ak5u3+ak6v3+ak7uv+ak8uv2+ak9u2vk=1,2,⋯,9where u,v are the longitude and the latitude of the maneuver point, respectively, and I(k)(u,v) is the kth value, which is obtained by interpolation, in the vector Z(Θi,Φi).The neighborhoodsδ(k)(k=1,2,⋯,9) of the initial values of the adjoint variables and the values of partial state variables are taken as search ranges of optimization variables, and the minimum duration transfer trajectory is computed by the hybrid method.


## 4. Simulations and Analyses of Results

### 4.1. Simulations

The parameters of the vehicle are as follows: the orbital radius is 6678.14km, the flattening is 0, the value of the orbital inclination is 50°, the right ascension of ascending node is 0°, the argument of perigee is 105°, the true anomaly is 0°, the mass of the spacecraft is 500kg, the engine thrust is 100N, the gas jet velocity is 3000m/s, the shape of the vehicle is an axisymmetric cone, and the characteristic area of vehicle is 0.02m2.The position of the landing point is (-158°, 38°) (the first value is the longitude and the second is the latitude). The region S composed of maneuverable points (shown as Figure6) is bounded by the four points A(13°, 40°), B(15°, 38°), C(19°, 39°), and D(16°, 41°) in turn on the surface of the earth. The discrete granularity of the longitude is 0.5° and the longitude range is between 13° and 19°. The grid points are the intersections of the set longitudes and the tracks of subsatellite points with a spacing distance between 29 km and 31 km in the direction of AC. The neighborhood ranges are the predicted values±2%. The range of V-a is [7.5km/s, 8.8 km/s]. The range of γa is [-9°, -6°]. All maneuver points are in the region S. The position of one hundred maneuver points is shown in Table 1.Table 1
Position of one hundred maneuver points.
No. Θ Φ 1 13.1 40.1 2 13.2 40.1 3 13.3 40.1 4 13.4 40.1 5 13.6 40.1 6 13.7 40.1 7 13.7 40.1 8 13.8 40.1 9 13.9 40.1 10 14.1 40.2 11 14.2 40.2 12 14.3 40.2 13 14.4 40.2 14 14.6 40.2 15 14.7 40.2 16 14.8 40.2 17 14.9 40.2 18 15.1 40.2 19 15.2 40.2 20 15.3 40.2 21 15.4 40.3 22 15.6 40.3 23 15.7 40.3 24 15.8 40.3 25 15.9 40.3 26 16.1 40.3 27 16.2 39.6 28 16.3 39.6 29 16.4 39.6 30 16.6 39.6 31 16.7 39.6 32 16.8 39.6 33 16.9 39.6 34 17.1 39.3 35 17.2 39.3 36 17.3 39.3 37 17.4 39.3 38 17.6 39.3 39 17.7 39.3 40 17.8 39.3 41 17.9 39.3 42 18.1 39.1 43 18.2 39.1 44 18.3 39.1 45 18.4 39.1 46 18.6 39.1 47 18.7 39.1 48 18.8 39.1 49 18.9 39.1 50 14.1 40.1 51 14.2 40.1 52 14.3 40.1 53 14.4 40.1 54 14.6 40.1 55 14.7 40.1 56 14.8 40.1 57 14.9 40.1 58 15.1 40.1 59 15.2 40.1 60 15.3 40.1 61 15.4 39.9 62 15.6 39.9 63 15.7 39.9 64 15.8 39.9 65 15.9 39.9 66 16.1 39.9 67 16.2 39.9 68 16.3 39.9 69 16.4 39.9 70 16.6 39.9 71 16.7 39.9 72 16.8 39.9 73 16.9 39.9 74 17.1 39.4 75 17.2 39.4 76 17.3 39.4 77 17.4 39.4 78 17.6 39.4 79 17.7 39.4 80 17.8 39.4 81 15.4 38.9 82 15.6 38.9 83 15.7 38.9 84 15.8 38.9 85 15.9 38.9 86 16.1 38.9 87 16.2 38.9 88 16.3 38.9 89 16.4 38.9 90 16.6 38.9 91 16.7 38.9 92 16.8 38.9 93 16.9 38.9 94 17.1 39.2 95 17.2 39.2 96 17.3 39.2 97 17.4 39.2 98 17.6 39.2 79 17.7 39.2 100 17.8 39.2 In Table1, the units of Θ and Φ are ° and ° respectivelyFigure 6
The area S where the spacecraft can maneuver.The idea of the proposed method is as follows. First, the minimum duration transfer trajectories at all grid points are obtained by using the traditional hybrid method. Second, the predicted values of minimum duration transfer trajectories at one hundred maneuver points, excluding the grid points, are computed by using the proposed method. Finally, by taking the neighborhood ranges as the search range, the minimum duration transfer trajectories at one hundred maneuver points are obtained by using the traditional hybrid method. The results obtained by the proposed method are used for comparison with the traditional hybrid method.


### 4.2. Analyses of Results

TheλV-(t-0), λγ(t-0), λψ(t-0), λr-(t-0), λΘ(t-0), λm-(t-0), V-a, γa, and t-a at grid points are shown in Table 2. Landing speed Vb and final mass mb of vehicles in one hundred experiments are shown in Table 3. Figures 7–15 are the spatial distribution of λV-(t-0), λγ(t-0), λψ(t-0), λr-(t-0), λΘ(t-0), λm-(t-0), V-a, γa, and t-a at all grid points. We can see that values of each parameter at all grid points are on a smooth surface, satisfying the interpolation conditions. In order to further demonstrate the effectiveness of our proposed method, we conducted 100 sets of comparative experiments between our method and the traditional hybrid method. Figure 16 shows the transfer trajectory of No. 5 experiment. Figure 17 shows the transfer trajectory of No. 45 experiment. Figure 18 shows the transfer trajectory of No. 85 experiment.Table 2
Initial values of the adjoint variables (partial) and the values of partial state variables at grid points.
Θ Φ λ V - t - 0 λ γ t - 0 λ ψ t - 0 λ r - t - 0 λ Θ t - 0 λ m - t - 0 V - a γ a t - a 13 40.000 11.308 3.217 -19.323 14.287 0.002 1.214 8.383 -7.216 4505.6  13.5 40.159 11.209 3.289 -19.306 14.293 0.154 1.361 8.352 -7.142 4486.5 39.732 11.176 3.295 -19.342 14.248 -0.283 1.792 8.439 -7.357 4509.3  14 40.280 11.159 3.348 -19.289 14.204 0.259 1.472 8.330 -7.092 4466.1 39.904 11.111 3.357 -19.308 14.165 -0.123 1.848 8.407 -7.282 4486.2 39.137 11.042 3.374 -19.383 14.086 -0.887 2.618 8.559 -7.667 4524.4  14.5 40.603 11.103 3.412 -19.167 14.136 0.571 1.139 8.268 -6.943 4438.9 40.216 11.057 3.419 -19.216 14.095 0.181 1.527 8.347 -7.137 4459.7 39.757 11.001 3.432 -19.274 14.043 -0.283 1.990 8.439 -7.366 4483.0 39.059 10.931 3.447 -19.351 13.966 -0.991 2.688 8.580 -7.718 4518.3  15 40.643 10.906 3.482 -19.418 14.044 0.602 1.089 8.261 -6.933 4425.6 40.311 10.875 3.488 -19.456 14.004 0.269 1.424 8.328 -7.101 4443.3 39.878 10.833 3.451 -19.502 13.955 -0.174 1.861 8.415 -7.319 4466.2 39.125 10.752 3.469 -19.587 13.870 -0.928 2.615 8.568 -7.697 4505.3 38.245 10.659 3.491 -19.685 13.774 -1.818 3.494 8.746 -8.119 4550.3  15.5 40.727 10.518 3.540 -19.775 13.948 0.678 0.996 8.246 -6.902 4410.2 40.413 10.471 3.579 -19.824 13.914 0.362 1.313 8.309 -7.061 4427.3 39.758 10.401 3.604 -19.896 13.847 -0.301 1.969 8.442 -7.391 4460.4 39.218 10.342 3.619 -19.961 13.780 -0.845 2.511 8.553 -7.664 4488.8 38.365 10.248 3.642 -20.060 13.686 -1.699 3.364 8.724 -8.089 4531.4  16 40.830 10.414 3.651 -19.836 13.859 0.769 0.880 8.227 -6.862 4393.6 40.487 10.377 3.659 -19.871 13.823 0.423 1.225 8.297 -7.035 4411.5 39.873 10.313 3.674 -19.941 13.756 -0.188 1.843 8.422 -7.345 4442.5 39.258 10.246 3.689 -20.012 13.685 -0.796 2.461 8.546 -7.654 4474.4 38.607 10.174 3.706 -20.091 13.613 -1.452 3.117 8.677 -7.983 4506.7  16.5 40.578 10.287 3.730 -19.917 13.734 0.504 1.120 8.280 -7.000 4394.8 40.192 10.247 3.739 -19.964 13.690 0.115 1.508 8.358 -7.196 4414.6 39.558 10.181 3.754 -20.037 13.618 -0.517 2.146 8.486 -7.515 4446.3 38.879 10.109 3.770 -20.117 13.543 -1.199 2.832 8.623 -7.856 4480.8  17 40.260 10.157 3.804 -20.013 13.596 0.175 1.431 8.346 -7.171 4379.4 39.821 10.108 3.810 -20.060 13.553 -0.269 1.873 8.434 -7.390 4402.5 39.372 10.059 3.822 -20.114 13.508 -0.721 2.324 8.525 -7.615 4425.7 38.877 10.007 3.832 -20.171 13.454 -1.220 2.824 8.624 -7.864 4451.6  17.5 39.827 10.010 3.881 -20.119 13.459 -0.278 1.860 8.435 -7.399 4329.4 39.546 9.985 3.888 -20.150 13.427 -0.565 2.143 8.492 -7.540 4343.6 38.931 9.913 3.903 -20.223 13.364 -1.012 2.763 8.616 -7.849 4374.5  18 39.588 9.882 3.957 -20.196 13.333 -0.513 2.090 8.485 -7.528 4310.4 38.896 9.814 3.973 -20.272 13.262 -1.210 2.785 8.624 -7.876 4345.7  18.5 39.274 9.753 4.038 -20.459 13.201 -0.838 2.396 8.551 -7.697 4315.3  19 39.072 9.726 4.096 -20.613 13.076 -1.051 2.589 8.594 -7.808 4315.7 In Table2, the units of Θ, Φ, V-a, γa and t-a are °, °, km/s, ° and s respectively.Table 3
Landing speed and final mass of vehicles in one hundred experiments.
No. V b m b 1 7.641 404.2 2 7.610 404.6 3 7.554 404.9 4 7.502 405.1 5 7.445 405.3 6 7.410 405.4 7 7.399 405.6 8 7.373 405.8 9 7.300 405.9 10 7.276 406.1 11 7.223 406.2 12 7.189 406.4 13 7.067 406.6 14 7.001 406.7 15 6.948 406.7 16 6.902 406.9 17 6.852 407.1 18 6.821 407.2 19 6.794 407.3 20 6.756 407.5 21 6.713 407.7 22 6.682 407.7 23 6.653 407.8 24 6.629 408.1 25 6.206 408.2 26 6.197 408.3 27 7.182 406.5 28 7.154 406.7 29 7.001 406.9 30 6.922 407.0 31 6.879 407.0 32 6.798 407.1 33 6.740 406.9 34 7.165 406.6 35 7.133 406.8 36 7.095 406.9 37 7.010 407.4 38 6.921 407.7 39 6.833 408.1 40 6.795 408.2 41 6.721 408.5 42 7.109 407.1 43 7.056 407.3 44 7.001 407.4 45 6.950 407.6 46 6.923 407.6 47 6.887 407.7 48 6.832 407.9 49 6.821 408.0 50 7.229 406.4 51 7.155 406.5 52 7.102 406.7 53 7.079 406.8 54 7.001 407.0 55 6.938 407.3 56 6.899 407.5 57 6.723 407.8 58 6.686 407.9 59 6.569 408.1 60 6.430 408.3 61 6.442 408.1 62 6.402 408.3 63 6.368 408.4 64 6.321 408.5 65 6.300 408.5 66 6.279 408.6 67 6.225 408.8 68 6.198 408.9 69 6.158 409.1 70 6.113 409.2 71 6.083 409.2 72 6.053 409.3 73 6.002 409.5 74 7.102 406.5 75 7.086 406.6 76 7.026 406.8 77 6.956 407.1 78 6.823 407.4 79 6.791 407.5 80 6.722 407.7 81 6.865 407.0 82 6.823 407.2 83 6.784 407.4 84 6.736 407.5 85 6.695 407.6 86 6.652 407.8 87 6.601 407.9 88 6.539 408.2 89 6.493 408.3 90 6.457 408.5 91 6.401 408.7 92 6.380 408.7 93 6.333 408.9 94 6.555 408.3 95 6.512 408.4 96 6.478 408.6 97 6.402 408.8 98 6.373 408.9 79 6.345 409.0 100 6.291 409.2 In Table3,  Vb is landing speed and mb is final mass. The units of Vb and mb are km/s and kg respectively.Figure 7
Spatial distribution ofλV-(t-0) at all mesh points.Figure 8
Spatial distribution ofλγ(t-0) at all mesh points.Figure 9
Spatial distribution ofλψ(t-0) at all mesh points.Figure 10
Spatial distribution ofλr-(t-0) at all mesh points.Figure 11
Spatial distribution ofλΘ(t-0) at all mesh points.Figure 12
Spatial distribution ofλm-(t-0) at all mesh points.Figure 13
Spatial distribution ofV-a at all mesh points.Figure 14
Spatial distribution ofγa at all mesh points.Figure 15
Spatial distribution oft-a at all mesh points.Figure 16
Transfer trajectory of No.5 experiment.Figure 17
Transfer trajectory of No.45 experiment.Figure 18
Transfer trajectory of No.85 experiment.Figures19–27 are the fitting errors of 9 parameters in 100 sets of comparative experiments. We can see from Figures 19–27 that the fitting error of λV-(t-0) is less than 5e-4, the fitting error of λγ(t-0) is less than 1.5e-4, the fitting error of λψ(t-0) is less than 1e-3, the fitting error of λr-(t-0) is less than 8e-3, the fitting error of λΘ(t-0) is less than 1.5e-5, the fitting error of λm-(t-0) is less than 8e-5, the fitting error of V-a is less than 5e-8 km/s, the fitting error of γa is less than 4e-5, and the fitting error of t-a is less than 5e-4s. Figure 28 shows the landing point error, from which we can see that the landing point error is less than 2e-3m. Figure 29 shows fitting error of final mass, from which we can see that the fitting error of final mass is less than 3e-3m. Figure 30 shows the fitting error of total time, from which we can see that the fitting error of total time is less than 1.8e-3m. We can see from Figures 19–30 that the proposed method satisfies the standard accuracy requirements and has computational efficiency.Figure 19
Fitting error ofλV-(t-0).Figure 20
Fitting error ofλγ(t-0).Figure 21
Fitting error ofλψt-0.Figure 22
Fitting error ofλr-(t-0).Figure 23
Fitting error ofλΘ(t-0).Figure 24
Fitting error ofλm-(t-0).Figure 25
Fitting error ofV-a.Figure 26
Fitting error ofγa.Figure 27
Fitting error oft-a.Figure 28
Landing point error.Figure 29
Fitting error of final mass.Figure 30
Fitting error of total time.Figure31 is the comparison between the proposed method (interpolation scheme + hybrid method) and the traditional hybrid method, from which we can see that the computational time of our method is significantly shorter than that of the traditional hybrid method. The computational time in each of the 100 experiments, when using the traditional hybrid method, is between 1677.7s and 1682.3s, and the average computational time is 1680.4s. In contrast, the computational time with our method is between 49.13s and 49.47s, and the average computational time is 49.34s. Thus, the proposed method requires only 2.93% of the computational time of the traditional hybrid method.Figure 31
Comparison of the computational time of two methods.The simulation results demonstrate that the proposed method is feasible. Although the proposed method still needs some computational time to obtain the accurate transfer trajectories, yet the proposed method improves substantially the computational time. Fast computational of the space-to-ground transfer trajectories lays the foundation for the decision makers to develop space-to-ground transfer trajectory programmes in real time in an uncertain environment.


## 4.1. Simulations

The parameters of the vehicle are as follows: the orbital radius is 6678.14km, the flattening is 0, the value of the orbital inclination is 50°, the right ascension of ascending node is 0°, the argument of perigee is 105°, the true anomaly is 0°, the mass of the spacecraft is 500kg, the engine thrust is 100N, the gas jet velocity is 3000m/s, the shape of the vehicle is an axisymmetric cone, and the characteristic area of vehicle is 0.02m2.The position of the landing point is (-158°, 38°) (the first value is the longitude and the second is the latitude). The region S composed of maneuverable points (shown as Figure6) is bounded by the four points A(13°, 40°), B(15°, 38°), C(19°, 39°), and D(16°, 41°) in turn on the surface of the earth. The discrete granularity of the longitude is 0.5° and the longitude range is between 13° and 19°. The grid points are the intersections of the set longitudes and the tracks of subsatellite points with a spacing distance between 29 km and 31 km in the direction of AC. The neighborhood ranges are the predicted values±2%. The range of V-a is [7.5km/s, 8.8 km/s]. The range of γa is [-9°, -6°]. All maneuver points are in the region S. The position of one hundred maneuver points is shown in Table 1.Table 1
Position of one hundred maneuver points.
No. Θ Φ 1 13.1 40.1 2 13.2 40.1 3 13.3 40.1 4 13.4 40.1 5 13.6 40.1 6 13.7 40.1 7 13.7 40.1 8 13.8 40.1 9 13.9 40.1 10 14.1 40.2 11 14.2 40.2 12 14.3 40.2 13 14.4 40.2 14 14.6 40.2 15 14.7 40.2 16 14.8 40.2 17 14.9 40.2 18 15.1 40.2 19 15.2 40.2 20 15.3 40.2 21 15.4 40.3 22 15.6 40.3 23 15.7 40.3 24 15.8 40.3 25 15.9 40.3 26 16.1 40.3 27 16.2 39.6 28 16.3 39.6 29 16.4 39.6 30 16.6 39.6 31 16.7 39.6 32 16.8 39.6 33 16.9 39.6 34 17.1 39.3 35 17.2 39.3 36 17.3 39.3 37 17.4 39.3 38 17.6 39.3 39 17.7 39.3 40 17.8 39.3 41 17.9 39.3 42 18.1 39.1 43 18.2 39.1 44 18.3 39.1 45 18.4 39.1 46 18.6 39.1 47 18.7 39.1 48 18.8 39.1 49 18.9 39.1 50 14.1 40.1 51 14.2 40.1 52 14.3 40.1 53 14.4 40.1 54 14.6 40.1 55 14.7 40.1 56 14.8 40.1 57 14.9 40.1 58 15.1 40.1 59 15.2 40.1 60 15.3 40.1 61 15.4 39.9 62 15.6 39.9 63 15.7 39.9 64 15.8 39.9 65 15.9 39.9 66 16.1 39.9 67 16.2 39.9 68 16.3 39.9 69 16.4 39.9 70 16.6 39.9 71 16.7 39.9 72 16.8 39.9 73 16.9 39.9 74 17.1 39.4 75 17.2 39.4 76 17.3 39.4 77 17.4 39.4 78 17.6 39.4 79 17.7 39.4 80 17.8 39.4 81 15.4 38.9 82 15.6 38.9 83 15.7 38.9 84 15.8 38.9 85 15.9 38.9 86 16.1 38.9 87 16.2 38.9 88 16.3 38.9 89 16.4 38.9 90 16.6 38.9 91 16.7 38.9 92 16.8 38.9 93 16.9 38.9 94 17.1 39.2 95 17.2 39.2 96 17.3 39.2 97 17.4 39.2 98 17.6 39.2 79 17.7 39.2 100 17.8 39.2 In Table1, the units of Θ and Φ are ° and ° respectivelyFigure 6
The area S where the spacecraft can maneuver.The idea of the proposed method is as follows. First, the minimum duration transfer trajectories at all grid points are obtained by using the traditional hybrid method. Second, the predicted values of minimum duration transfer trajectories at one hundred maneuver points, excluding the grid points, are computed by using the proposed method. Finally, by taking the neighborhood ranges as the search range, the minimum duration transfer trajectories at one hundred maneuver points are obtained by using the traditional hybrid method. The results obtained by the proposed method are used for comparison with the traditional hybrid method.


## 4.2. Analyses of Results

TheλV-(t-0), λγ(t-0), λψ(t-0), λr-(t-0), λΘ(t-0), λm-(t-0), V-a, γa, and t-a at grid points are shown in Table 2. Landing speed Vb and final mass mb of vehicles in one hundred experiments are shown in Table 3. Figures 7–15 are the spatial distribution of λV-(t-0), λγ(t-0), λψ(t-0), λr-(t-0), λΘ(t-0), λm-(t-0), V-a, γa, and t-a at all grid points. We can see that values of each parameter at all grid points are on a smooth surface, satisfying the interpolation conditions. In order to further demonstrate the effectiveness of our proposed method, we conducted 100 sets of comparative experiments between our method and the traditional hybrid method. Figure 16 shows the transfer trajectory of No. 5 experiment. Figure 17 shows the transfer trajectory of No. 45 experiment. Figure 18 shows the transfer trajectory of No. 85 experiment.Table 2
Initial values of the adjoint variables (partial) and the values of partial state variables at grid points.
Θ Φ λ V - t - 0 λ γ t - 0 λ ψ t - 0 λ r - t - 0 λ Θ t - 0 λ m - t - 0 V - a γ a t - a 13 40.000 11.308 3.217 -19.323 14.287 0.002 1.214 8.383 -7.216 4505.6  13.5 40.159 11.209 3.289 -19.306 14.293 0.154 1.361 8.352 -7.142 4486.5 39.732 11.176 3.295 -19.342 14.248 -0.283 1.792 8.439 -7.357 4509.3  14 40.280 11.159 3.348 -19.289 14.204 0.259 1.472 8.330 -7.092 4466.1 39.904 11.111 3.357 -19.308 14.165 -0.123 1.848 8.407 -7.282 4486.2 39.137 11.042 3.374 -19.383 14.086 -0.887 2.618 8.559 -7.667 4524.4  14.5 40.603 11.103 3.412 -19.167 14.136 0.571 1.139 8.268 -6.943 4438.9 40.216 11.057 3.419 -19.216 14.095 0.181 1.527 8.347 -7.137 4459.7 39.757 11.001 3.432 -19.274 14.043 -0.283 1.990 8.439 -7.366 4483.0 39.059 10.931 3.447 -19.351 13.966 -0.991 2.688 8.580 -7.718 4518.3  15 40.643 10.906 3.482 -19.418 14.044 0.602 1.089 8.261 -6.933 4425.6 40.311 10.875 3.488 -19.456 14.004 0.269 1.424 8.328 -7.101 4443.3 39.878 10.833 3.451 -19.502 13.955 -0.174 1.861 8.415 -7.319 4466.2 39.125 10.752 3.469 -19.587 13.870 -0.928 2.615 8.568 -7.697 4505.3 38.245 10.659 3.491 -19.685 13.774 -1.818 3.494 8.746 -8.119 4550.3  15.5 40.727 10.518 3.540 -19.775 13.948 0.678 0.996 8.246 -6.902 4410.2 40.413 10.471 3.579 -19.824 13.914 0.362 1.313 8.309 -7.061 4427.3 39.758 10.401 3.604 -19.896 13.847 -0.301 1.969 8.442 -7.391 4460.4 39.218 10.342 3.619 -19.961 13.780 -0.845 2.511 8.553 -7.664 4488.8 38.365 10.248 3.642 -20.060 13.686 -1.699 3.364 8.724 -8.089 4531.4  16 40.830 10.414 3.651 -19.836 13.859 0.769 0.880 8.227 -6.862 4393.6 40.487 10.377 3.659 -19.871 13.823 0.423 1.225 8.297 -7.035 4411.5 39.873 10.313 3.674 -19.941 13.756 -0.188 1.843 8.422 -7.345 4442.5 39.258 10.246 3.689 -20.012 13.685 -0.796 2.461 8.546 -7.654 4474.4 38.607 10.174 3.706 -20.091 13.613 -1.452 3.117 8.677 -7.983 4506.7  16.5 40.578 10.287 3.730 -19.917 13.734 0.504 1.120 8.280 -7.000 4394.8 40.192 10.247 3.739 -19.964 13.690 0.115 1.508 8.358 -7.196 4414.6 39.558 10.181 3.754 -20.037 13.618 -0.517 2.146 8.486 -7.515 4446.3 38.879 10.109 3.770 -20.117 13.543 -1.199 2.832 8.623 -7.856 4480.8  17 40.260 10.157 3.804 -20.013 13.596 0.175 1.431 8.346 -7.171 4379.4 39.821 10.108 3.810 -20.060 13.553 -0.269 1.873 8.434 -7.390 4402.5 39.372 10.059 3.822 -20.114 13.508 -0.721 2.324 8.525 -7.615 4425.7 38.877 10.007 3.832 -20.171 13.454 -1.220 2.824 8.624 -7.864 4451.6  17.5 39.827 10.010 3.881 -20.119 13.459 -0.278 1.860 8.435 -7.399 4329.4 39.546 9.985 3.888 -20.150 13.427 -0.565 2.143 8.492 -7.540 4343.6 38.931 9.913 3.903 -20.223 13.364 -1.012 2.763 8.616 -7.849 4374.5  18 39.588 9.882 3.957 -20.196 13.333 -0.513 2.090 8.485 -7.528 4310.4 38.896 9.814 3.973 -20.272 13.262 -1.210 2.785 8.624 -7.876 4345.7  18.5 39.274 9.753 4.038 -20.459 13.201 -0.838 2.396 8.551 -7.697 4315.3  19 39.072 9.726 4.096 -20.613 13.076 -1.051 2.589 8.594 -7.808 4315.7 In Table2, the units of Θ, Φ, V-a, γa and t-a are °, °, km/s, ° and s respectively.Table 3
Landing speed and final mass of vehicles in one hundred experiments.
No. V b m b 1 7.641 404.2 2 7.610 404.6 3 7.554 404.9 4 7.502 405.1 5 7.445 405.3 6 7.410 405.4 7 7.399 405.6 8 7.373 405.8 9 7.300 405.9 10 7.276 406.1 11 7.223 406.2 12 7.189 406.4 13 7.067 406.6 14 7.001 406.7 15 6.948 406.7 16 6.902 406.9 17 6.852 407.1 18 6.821 407.2 19 6.794 407.3 20 6.756 407.5 21 6.713 407.7 22 6.682 407.7 23 6.653 407.8 24 6.629 408.1 25 6.206 408.2 26 6.197 408.3 27 7.182 406.5 28 7.154 406.7 29 7.001 406.9 30 6.922 407.0 31 6.879 407.0 32 6.798 407.1 33 6.740 406.9 34 7.165 406.6 35 7.133 406.8 36 7.095 406.9 37 7.010 407.4 38 6.921 407.7 39 6.833 408.1 40 6.795 408.2 41 6.721 408.5 42 7.109 407.1 43 7.056 407.3 44 7.001 407.4 45 6.950 407.6 46 6.923 407.6 47 6.887 407.7 48 6.832 407.9 49 6.821 408.0 50 7.229 406.4 51 7.155 406.5 52 7.102 406.7 53 7.079 406.8 54 7.001 407.0 55 6.938 407.3 56 6.899 407.5 57 6.723 407.8 58 6.686 407.9 59 6.569 408.1 60 6.430 408.3 61 6.442 408.1 62 6.402 408.3 63 6.368 408.4 64 6.321 408.5 65 6.300 408.5 66 6.279 408.6 67 6.225 408.8 68 6.198 408.9 69 6.158 409.1 70 6.113 409.2 71 6.083 409.2 72 6.053 409.3 73 6.002 409.5 74 7.102 406.5 75 7.086 406.6 76 7.026 406.8 77 6.956 407.1 78 6.823 407.4 79 6.791 407.5 80 6.722 407.7 81 6.865 407.0 82 6.823 407.2 83 6.784 407.4 84 6.736 407.5 85 6.695 407.6 86 6.652 407.8 87 6.601 407.9 88 6.539 408.2 89 6.493 408.3 90 6.457 408.5 91 6.401 408.7 92 6.380 408.7 93 6.333 408.9 94 6.555 408.3 95 6.512 408.4 96 6.478 408.6 97 6.402 408.8 98 6.373 408.9 79 6.345 409.0 100 6.291 409.2 In Table3,  Vb is landing speed and mb is final mass. The units of Vb and mb are km/s and kg respectively.Figure 7
Spatial distribution ofλV-(t-0) at all mesh points.Figure 8
Spatial distribution ofλγ(t-0) at all mesh points.Figure 9
Spatial distribution ofλψ(t-0) at all mesh points.Figure 10
Spatial distribution ofλr-(t-0) at all mesh points.Figure 11
Spatial distribution ofλΘ(t-0) at all mesh points.Figure 12
Spatial distribution ofλm-(t-0) at all mesh points.Figure 13
Spatial distribution ofV-a at all mesh points.Figure 14
Spatial distribution ofγa at all mesh points.Figure 15
Spatial distribution oft-a at all mesh points.Figure 16
Transfer trajectory of No.5 experiment.Figure 17
Transfer trajectory of No.45 experiment.Figure 18
Transfer trajectory of No.85 experiment.Figures19–27 are the fitting errors of 9 parameters in 100 sets of comparative experiments. We can see from Figures 19–27 that the fitting error of λV-(t-0) is less than 5e-4, the fitting error of λγ(t-0) is less than 1.5e-4, the fitting error of λψ(t-0) is less than 1e-3, the fitting error of λr-(t-0) is less than 8e-3, the fitting error of λΘ(t-0) is less than 1.5e-5, the fitting error of λm-(t-0) is less than 8e-5, the fitting error of V-a is less than 5e-8 km/s, the fitting error of γa is less than 4e-5, and the fitting error of t-a is less than 5e-4s. Figure 28 shows the landing point error, from which we can see that the landing point error is less than 2e-3m. Figure 29 shows fitting error of final mass, from which we can see that the fitting error of final mass is less than 3e-3m. Figure 30 shows the fitting error of total time, from which we can see that the fitting error of total time is less than 1.8e-3m. We can see from Figures 19–30 that the proposed method satisfies the standard accuracy requirements and has computational efficiency.Figure 19
Fitting error ofλV-(t-0).Figure 20
Fitting error ofλγ(t-0).Figure 21
Fitting error ofλψt-0.Figure 22
Fitting error ofλr-(t-0).Figure 23
Fitting error ofλΘ(t-0).Figure 24
Fitting error ofλm-(t-0).Figure 25
Fitting error ofV-a.Figure 26
Fitting error ofγa.Figure 27
Fitting error oft-a.Figure 28
Landing point error.Figure 29
Fitting error of final mass.Figure 30
Fitting error of total time.Figure31 is the comparison between the proposed method (interpolation scheme + hybrid method) and the traditional hybrid method, from which we can see that the computational time of our method is significantly shorter than that of the traditional hybrid method. The computational time in each of the 100 experiments, when using the traditional hybrid method, is between 1677.7s and 1682.3s, and the average computational time is 1680.4s. In contrast, the computational time with our method is between 49.13s and 49.47s, and the average computational time is 49.34s. Thus, the proposed method requires only 2.93% of the computational time of the traditional hybrid method.Figure 31
Comparison of the computational time of two methods.The simulation results demonstrate that the proposed method is feasible. Although the proposed method still needs some computational time to obtain the accurate transfer trajectories, yet the proposed method improves substantially the computational time. Fast computational of the space-to-ground transfer trajectories lays the foundation for the decision makers to develop space-to-ground transfer trajectory programmes in real time in an uncertain environment.


## 5. Conclusion and Future Research

Fast computational method of the space-to-ground transfer trajectories can lay the foundation for decision makers to develop space-to-ground transfer trajectory programmes in real time in an uncertain environment. In order to overcome the computational cost encountered with the hybrid method, we solve the transfer trajectory by using interpolation scheme and hybrid method. The main idea of this paper is to convert the problem from the time domain into the spatial domain. The simulation results demonstrate that the proposed method has high computational efficiency.Although the proposed method has high computational efficiency, it still has some limitations that some time is still needed to obtain the accurate transfer trajectory. To further accelerate the computational time of the transfer trajectories, investigating more interpolation methods as well as optimization algorithms in future research would be better. The methods that combine interpolation schemes and hybrid methods can be applied to other problems related to fixed transfer trajectories and finite thrust maneuver mode. Furthermore, the methodology of converting the time domain into the space domain can also be applied to other research areas.



---

*Source: 5319530-2019-05-20.xml*