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.

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 [310], indirect methods [1114], and hybrid methods [1517] 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 [2022]. 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 (Figure 1). 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.

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]where 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), is the dimensionless geocentric distance, is the longitude, is the latitude, is the engine thrust, is the gas jet velocity, is the dimensionless spacecraft quality, is the dimensionless time, is the dimensionless resistance, is the dimensionless lift, is the dimensionless lateral force, is the angle of attack, is the sideslip angle, is the dimensionless earth rotation angular rate, and and 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 , , and are zero. In the reentry trajectory segment, the value of is zero and the atmospheric model is US standard atmosphere (1976).

The equations for the dimensionless parameters arewhere is the coefficient of second-order principal spherical harmonic function (the value of is 1.08263e-3.); and , , , , and arewhere (the value of is 6.37100e+6 m.) is the reference radius of the earth, is the gravitation constant (the value of is 3.986005e+14), and 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) [2628]. At last, we obtain the minimum duration transfer trajectory.

The minimum duration optimization index is written aswhere , , and 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

From (5), the covariate variables satisfy the differential equations (only the transition trajectory with thrust control is considered):where the expressions of , , , and are given by

From the sufficient conditions of optimality and formula (8), the optimal thrust directions are obtained bywhere is the sign function.

The minimum duration transfer trajectory consists of the transition trajectory segment and the reentry trajectory segment. The dimensionless geocentric distance at the reentry point is a constant. The dimensionless velocity and the velocity inclination angle 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 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:where , , , , , , and 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:

Therefore, the cross-section conditions are given by

The reentry trajectory satisfies the constraints

Because the index includes time, we havewhere , , and are the dimensionless speed, velocity inclination angle, and geocentric distance at the terminal of transition trajectory, respectively; and are the longitude and latitude of the landing point, respectively; and 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 (, and ) and the values of four state variables (, , , and ). Since the initial value of an adjoint variable can be obtained from (5) and (13), and (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. , , , , , , , , and are encoded. Twenty chromosomes are generated randomly as the initial population. The crossover probability , the mutation probability , 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. Steps 2 and 3 are repeatedly executed until the GA converges 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.

Figures 2, 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 , and the grid points are the intersections of the set longitude and tracks of subsatellite points. When , , , and are set in formula (1) and an arbitrary longitude value is given, the latitude value at the grid point can be obtained.

After the region (shown as Figure 5) composed of maneuverable points is gridded, gridded points are obtained. The function value at the gridded point ( is the set longitude and is latitude value of the track of subsatellite) is the vector :

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 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 are inserted in formula (14) to obtain the ten constant coefficients . The interpolated values (the predicted values) of the initial values of the adjoint variables and the values of partial state variables at the point can be computed by using :where are the longitude and the latitude of the maneuver point, respectively, and is the value, which is obtained by interpolation, in the vector .

The neighborhoods 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 Figure 6) 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 is [7.5km/s, 8.8 km/s]. The range of is [-9°, -6°]. All maneuver points are in the region S. The position of one hundred maneuver points is shown in Table 1.

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 , , , , , , , , and at grid points are shown in Table 2. Landing speed and final mass of vehicles in one hundred experiments are shown in Table 3. Figures 715 are the spatial distribution of , , , , , , , , and 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.

Figures 1927 are the fitting errors of 9 parameters in 100 sets of comparative experiments. We can see from Figures 1927 that the fitting error of is less than 5e-4, the fitting error of is less than 1.5e-4, the fitting error of is less than 1e-3, the fitting error of is less than 8e-3, the fitting error of is less than 1.5e-5, the fitting error of is less than 8e-5, the fitting error of is less than 5e-8 km/s, the fitting error of is less than 4e-5, and the fitting error of 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 1930 that the proposed method satisfies the standard accuracy requirements and has computational efficiency.

Figure 31 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.

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.

Data Availability

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

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.