Abstract
The paper deals with the multiple gravity assist trajectories design. In order to improve the performance of the heuristic algorithms, such as differential evolution algorithm, in multiple gravity assist trajectories design optimization, a method combining BFS (breadth-first search) and EP_DE (differential evolution algorithm based on search space exploring and principal component analysis) is proposed. In this method, firstly find the possible multiple gravity assist planet sequences with pruning based BFS and use standard differential evolution algorithm to judge the possibility of all the possible trajectories. Then select the better ones from all the possible solutions. Finally, use EP_DE which will be introduced in this paper to find an optimal decision vector of spacecraft transfer time schedule (launch window and transfer duration) for each selected planet sequence. In this paper, several cases are presented to prove the efficiency of the method proposed.
1. Introduction
Multiple gravity assist in deep space exploration mission design can greatly reduce the fuel expense. In general, all the planets in solar system can be used to design the gravity assist trajectories. At the same time, each planet can be used more than once. So, there are so many possible planet sequences for the gravity assist trajectories design. For each possible gravity assist planet sequence, there also must be an outstanding algorithm which can find an optimal decision vector of spacecraft transfer time schedule. A global optimization method, such as differential evolution, is suggested for seeking solutions which are likely to be close to the global optimum. However, the optimal solution may not be found only with standard differential evolution, so there must be some improvements on it. How to design the gravity assist trajectories has attracted a lot of researchers to work on it, and lots of research achievements have been published [1–5]. Many multiple gravity assist missions such as Cassini, Rosetta, and Messenger have been the benchmarks announced by ACT (Advanced Concepts Team of European Space Agency) [1]. In addition, a combination of stochastic search algorithms and search algorithm is proposed to solve the problem of multiple gravity assist trajectory design by Vasile [2]. In his recent research, he also applies bionics to solving the problem [3]. Izzo et al. put forward the method of pruning search space based on reducing the search space with a series of constraints [4]. Hennes and Izzo also use Monte Carlo Tree to find a solution [5]; at the same time, Monte Carlo Tree search method has been the hottest technology for optimization, especially after the excellent performance of AlphaGo [6]. Hartmann et al. use stochastic search algorithms to optimize the space trajectories design [7]. However, search space pruning proposed by Izzo et al. also might delete a better solution, though it is an effective algorithm. In paper [8], clustering was used to decide where to branch. In the study of the problem, it is difficult to get a better result if we only use DE and other search algorithms to search the optimal trajectory. In Cassini, it seems to be hard when we try to find a better solution than 5.3 km/s if only using DE algorithm [9]. So, we also use combinatorial method to deal with multiple gravity assist trajectories design.
In addition, in this paper, breadth-first search was also mentioned. It is true that BFS can lead to an optimal solution if there is enough calculation expense. Frequently, some pruning methods must be used for BFS. However, usually there is no suitable pruning method that can be found easily. So some computer calculating methods are used here. Korf and Schultze propose the method of how to design large-scale parallel breadth-first search which can get a result in less calculation time [10]. Similarly, Zhu and Cheung and Awerbuch and Gallager are devoted to a new distributed breadth-first search algorithm [11, 12]. Bisson et al. present the results obtained by using an evolution of CUDA-based solution, via a breadth-first search, for the exploration [13]. Katsuta et al. propose a tightly coupled accelerator which can accelerate breadth-first search [14]. Of course, in order to improve the efficiency of BFS, Liu et al. finish some improvements to it [15].
When designing the multiple gravity assist trajectories, we must find the most possible gravity assist planets sequence and furthermore find the optimal solution of gravity assist time sequence [16]. In this progress, one possible planets sequence is usually determined by deterministic algorithms, such as BFS algorithm, and the optimal gravity assist time solution is frequently obtained by using heuristic algorithms [2]. But for heuristic algorithms including DE, GA, and PSO, because of the complexity of search space, there is no universally valid method which can deal with the multiple gravity assist trajectories optimization. Therefore, in this paper, the EP_DE algorithm (we also can call it EP_DE I algorithm, because there may be some improved variants in the future research) is proposed based on DE algorithm, whose purpose is to reduce the global search space and improve the performance of DE algorithm. Also, this method of reducing search space is also available for other heuristic algorithms, which we will prove from the experiments presented in this paper.
We are devoted to finding a heuristic search algorithm based on machine learning, and EP_DE algorithm proposed in this paper can be seen as a basic version of this idea. And the feasibility of the algorithm has been proved by the benchmarks Cassini and GTOC1. With the development of hardware technology, heuristic algorithms based on machine learning will become a new research direction. So, EP_DE algorithm has important reference meaning for the development of heuristic search algorithm based on machine learning.
Obviously, compared to previous studies, the advantages of the method proposed in this paper are the guidance strategy based on data analysis. In prior studies, basic aerospace dynamics knowledge is usually used to reduce the global search space [4], while, in this paper, statistical features of a large number of data pieces are used to guide seeking the local space of optimal value. Therefore, EP_DE algorithm makes full use of the historical data generated in the process of optimization, rather than simply depending on some basic pruning method, such as angle constraint.
In this paper, a method combining BFS (breadth-first search) and EP_DE (differential evolution algorithm based on search space exploring and principal component analysis) is proposed. Firstly, find the possible multiple gravity assist planet sequences with pruning BFS; here, some thresholds have been set to prune the search space. In order to judge the efficiency of the planet sequence, standard differential evolution algorithm is applied to find the initial decision of the planet sequence. And it should be noted that the reason why we use standard differential evolution algorithm here is that it can cost less time. Then, select out the better planet sequence from all the possible sequences. Finally, find the optimal decision vector of spacecraft transfer time schedule for the planet sequence selected with EP_DE algorithm.
The rest of the paper is organized as follows. Section 2 briefly describes multiple gravity assist calculation model. In Section 3, the proposed method is presented in detail. Two experiments based on Cassini and GTOC1 are introduced in Section 4 which can prove the efficiency of EP_DE algorithm. Experiment result analysis and trajectory simulation are in Section 5. Finally, Section 6 is devoted to conclusion and future work.
2. Multiple Gravity Assist Calculation Model
In general, there are two kinds of MGA (multiple gravity assist) trajectory design: PMGA (Powered Multiple Gravity Assist) and Unpowered Multiple Gravity Assist. In this paper, assume that the final planet is a fixed destination, and dynamical model in the numerical test cases is MGA with powered flyby. PMGA is the method of applying a pulse to the spacecraft when it flies by a planet. And Unpowered Multiple Gravity Assist is the method of applying a pulse to spacecraft in heliocentric coordinate system.
In the trajectory design of PMGA, the trajectory can be seen as a combination of several partial trajectories which are solved as Lambert problems. It can be seen from Figure 1. In order to introduce the calculation model, we now define some mathematical symbols.
The position vector of spacecraft in the heliocentric coordinate system is . The velocity vector of spacecraft in the heliocentric coordinate system is . The position vector of gravity assist planet in the heliocentric coordinate system is . The velocity vector of gravity assist planet in the heliocentric coordinate system is . The time before spacecraft flies by the gravity assist planet is , and the time after spacecraft flies by the gravity assist planet is .
Before the spacecraft flies by the gravity assist planet, the velocity of the spacecraft relative to the velocity of planet is
After the spacecraft flies by the gravity assist planet, the relative velocity is
Of course, the value of the two relative velocities is the same. The value can be represented to be
If the angle between and is represented by , the relationship between them is
Assume that the planetary gravitational constant is , and the distance between spacecraft and planet center is . can be obtained by
The value range of is [0°, 180°].
Multiple gravity assist spacecraft trajectory design relies on two partial decision vectors. One is the sequence of planet; the other is decision vector of spacecraft transfer time schedule. Determine the sequence of planets and then solve the sequence of launch window and transfer duration. The calculation model is described as follows.
One partial decision vector is
is Earth; is each flyby planet. is the objective planet. And another partial decision vector is
is the time when spacecraft is launched from Earth, is the time when it flies by a planet, and is the time when arrives at the objective planet. is the transfer time interval from to .
Decision vector is
In each partial trajectory solved as a Lambert problem, there is a delta-. Based on the decision vector, the objective function can be expressed as follows:where is the delta- which can help spacecraft escape from Earth. is the max velocity constraint.
Meanwhile the delta applied here is
So we can know that is decided by , also can find the answer to , and go on to obtain the minimal .
The partial trajectory between two planets is solved as a Lambert problem. It is a two-point boundary value problem,
When time is , the position of spacecraft is , and the position of gravity assist planet is .
If there are partial trajectories, the problem can be solved by using Lambert calculations.
Finally, when combining two partial trajectories, we should find the suitable and subject to
If , there will be
3. Trajectory Design Based on BFS and EP_DE Algorithm
In this section, we will introduce the method proposed in detail. The introduction of the method is divided into two parts. One is how to use pruning based BFS to seek the gravity assist planet sequence, and the others are the process and efficiency of EP_DE algorithm.
There will be many possible planet sequences when optimizing the gravity assist trajectory if we do not prune the search space. So, some algorithm improvements about BFS must be proposed.
Assume that m planets are available and can be used for the gravity assist trajectory design, and the maximum times of gravity assist are . So the search space will be a search tree with layers and branches in each layer. It can be presented in Figure 2. In this paper, the experiment is based on the trajectory from Earth to Saturn. So, in each possible planet sequence, the objective planet is always Saturn. In general, the available planets are Venus, Earth, Mars, and Jupiter. That is to say, is 5 when counting in Saturn.
Pruning methods appear to be especially important while there are so many branches. In addition, when the pruning method is not suitable, the optimal solution obtained will be not ideal. The reason why we only obtain some suboptimal solutions will be introduced in this paper.
As for how to prune the search space, the most direct method is to set a threshold value, although it is difficult to determine the threshold. In the experiment, the threshold is determined by experience. The value of each solution is obtained by standard DE algorithm which can estimate the efficiency of the gravity assist planet sequence. Then, pick out the most possible planet sequences, and go on to calculate the optimal value of each planet sequence picked out with EP_DE algorithm. Next, there is the introduction of the EP_DE proposed in this paper. When calculating the optimal value of the mission Cassini at which gravity assist planet sequence is EVVEJS (Earth, Venus, Venus, Earth, Jupiter, and Saturn), we usually get 5.3 km/s.
When using standard differential evolution algorithm to optimize trajectory design problem, it can be found that usually the suboptimal value valley is much larger than a better value valley in spatial size, so it is hard to find a better solution. In order to solve this special problem, the algorithm proposed in this paper is showed in Algorithm 1 EP_DE. An earlier version of this algorithm was presented at the International Congress on Evolutionary Computation. The algorithm uses differential evolution algorithm based on search space exploring and PCA (principal component analysis), so it is called “EP_DE.”
Algorithm 1 (EP_DE).
Inputs Number of initialization samples: Pop0. Number of samples retained: Pop1. Number of clusters: Pop2. Number of dimension division: Pop3.Outputs Solution vector: .Step 1. Initialize Pop0 population samples and calculate the function values of all samples. Rank all the values of the function ascending order. Reserve the top Pop1 samples.
Step 2. Use -mean clustering algorithm to cluster Pop1 samples into Pop2 communities. In general, we need C communities, and the value range of it is usually from 3 to 5.
Step 3. Choose the community where we can find a better value. And continue to find the densest direction of samples distribution by PCA.
Step 4. Divide the community reserved into Pop3 parts in the densest direction of samples distribution. Use differential evolution algorithm to search optimal value of the small space in every partition obtained.
The main idea of algorithm is to explore the relationship between function value and variables. For one thing, estimate the general position of function value valley by initializing P0 samples in global search space; for another thing, learn the densest dimension of samples distribution in local search space.
In order to show the validity of method proposed in this paper, we use it to find the optimal solution in benchmark Cassini1 and GTOC1.
4. Efficiency of EP_DE Algorithm
4.1. Search Space Exploring
The meaning of search space exploring is to know the change of the function value in the search space by mathematical statistics.
There are two methods of search space exploring: method one is grid-initializing by using a large number of samples to know the relationship between search space function value and the variables change, which we call GI method. Method two randomly initializes a few parts of samples to know about the search space well, which is called method RI. For example, we initialize 25682962 samples in the experiment when using method RI and select the 500 best samples eventually. Apparently, the feature of method one in experiment results is low speed and high repeatability. When using method two, randomly initialize 2000 samples and pick out 500 excellent samples whose function values are outstanding. Apparently, the feature of method two in experiment results is high speed and low repeatability. In addition, the time system is JMD2000.
In the experiment of Cassini1, we use method one to explore the search space, and value range for each variable is in Table 1. Grid-initialization number on each variable is in Table 2. The grid number on each variable is determined according to the length of the variable. In the experiment of GTOC1, we use method two to explore the search space, and value range for each variable is in Table 3.
4.2. Clustering
In the experiment of Cassini, cluster selected 500 samples into five communities. Principal component analysis was performed on the samples in each community we selected. Also there is a mechanism to ensure that there are enough samples in each community for the PCA to be statistically meaningful. If, in a community, there are no enough meaningful samples, we must produce some new samples around the samples reserved by normal distribution. Because the best solution was found in the fourth community, experiment result of the fourth community is showed here. Its samples number is 67, and the value range of community is showed in Table 4.
In the experiment of GTOC1, the number of clusters is 3. Table 5 is spatial range of each community after clustering. Because we use method one to solve Cassini and method two to deal with GTOC1, therefore, in each experiment of GTOC1, the results will be different due to the instability of method two. But for it, if there are a suitable number of random initialization samples, the results also tend to be stable. So, in Table 5, we show two results in (a) and (b) for different experiments when the samples number is suitable. In order to know the location of the selected community clearly, we show the community range in Figures 3 and 4. In addition, the range in Figure 4 is for the community in Table 5(a).
4.3. Principal Component Analysis
Gain the covariance matrix by the fourth community samples, and eigenvalues and eigenvectors of the covariance matrix are calculated. From left to right, the characteristic value is arranged in the successive increasing order. The feature vector is represented as a column vector, and the feature vector corresponding to the characteristic value is also arranged from left to right in Table 6. Here, the sixth component vector with lowest variance and highest “density” is used to divide the fourth community. In fact, the direction of the sixth component is similar to . Because of the calculation convenience, we can directly divide the community space into parts in the direction of .
Optimal value found is in the second community, and two PCA results which mean there are two experiments related to Table 5 are showed in Table 7. For each range of the community we obtain that in every experiment there must be one only PCA result for it. So, Tables 7(a) and 7(b) are the results for Tables 5(a) and 5(b), respectively. We still use feature vector matrix to show the results. The eighth component is the densest dimension of samples distribution. From this direction we might find several function value channels.
4.4. Computational Results
In order to prove the efficiency of the EP_DE algorithm, the different calculation results are compared here by using DE and EP_DE. Based on the idea of improvement for DE, we also try some other improvements for PSO algorithm and GA algorithm, and they are called EP_PSO and EP_GA. 20 independent runs are performed by each algorithm. In Table 8, there are the best results obtained by standard algorithms including DE, PSO, and GA, while in Table 9, there are the best results gained by using improved algorithms including EP_DE, EP_PSO, and EP_GA.
4.5. The Relation between the Optimal Value and the Step Size
When dividing the community into several partitions on the selected dimension, the step size of the partition determines the optimal value we can find. Set different width of the partition and study the relationship between width of partition and the optimal value. The experiment results are in Table 10. From it, we can find that the calculation result of differential evolution algorithm is more stable than particle swarm optimization algorithm and genetic algorithm by the change of partition width.
Dividing the space in the selected dimension, we can find that step size of partition determines the optimal value. Set different partition width and study the relationship between partition width and optimal value. Experiment result is in Table 11. First of all, divide each community into several partitions with grid width of 1000 days and use differential evolution algorithm to find optimal value in each partition of community. Then because the best solution is found in the second community, we go on dividing the second community with grid width of 500 days, which can contribute to a better solution.
4.6. Verification
In the experiment of Cassini1, we ever tried to use second screening to find the optimal value, but the optimal value was lost. So there might be rapid growing of function value near the optimal value. When exploring search space, the optimal value is also removed. We fix four of all variables and draw the relationship between function value and the remaining two variables and in Figure 5.
As can be seen from Figure 5, there is rapid growing of function value. When value range of is set from 449.3857 to 449.3860, value range of is set from 54.7115 to 54.7119, function value is mainly affected by . Still use this method to view the relationship between the function value and the other variables in Figures 6, 7, 8, and 9. These graphs can also verify the analytical results in PCA.
Differential evolution algorithm factors have a great influence on results. The results of experiments are summarized and analyzed here. When there are 200 samples and 800 generations, best result we found is −473311.67197 km. When there are 2000 samples and 2000 generations, best result we found is −885276.164762 km. After clustering, use differential evolution algorithm to find the optimal value in each community. Each optimal value is showed in Table 12. The optimal value we found with the method proposed is in Table 13. Use ten experiments to count the repeatability of the algorithm; result is in Table 14. If there are more variables, the repeatability of the algorithm needs to be improved.
The global search space is estimated by the function value of samples in the global space. The community is estimated by PCA in each community. EP_DE algorithm can obtain the distribution of search space through this statistical method. Also, its efficiency in different problems is also different. In the calculation of the Cassini1, the results are relatively stable, while, in the calculation of the GTOC1, experiment result is still far from the optimal value that had already been announced. We will test the efficiency of EP_DE by other benchmarks in the next work.
4.7. Trajectory Simulation
Show the trajectory simulation based on the solution we have found. The trajectories simulations of Cassini and GTOC1 are in Figure 10.
(a)
(b)
5. Seek the Optimal Gravity Assist Planet Sequence for Cassini
5.1. Set Transfer Durations
The gravity assist planet sequence of above Cassini experiment mentioned is EVVEJS announced by ACT in ESA. And the best value found in Cassini with EVVEJS is 4.93 km/s. Of course, there may be a better gravity assist planet sequence. Next, we will try to seek a better solution with the method proposed in this paper.
In general, the value of the solution could be affected by transfer duration between two planets. In Tables 15 and 16, there is transfer duration between random two planets.
5.2. Set Pruning Threshold
The search computation expense will be high if there is no pruning threshold. However how to find the suitable threshold is difficult. But there is some prior knowledge that can be used by us. Firstly, while seeking the optimal value in Cassini with DE algorithm, the result found usually comes to 5.3 km/s. Then, the maximum energy expense of direct transfer trajectory is no more than 15 km/s. So, in the former search process, in order to prevent the pruning of possible solution, the threshold we set is 20 km/s while in the later search process, in order to reduce the computation expense, the threshold we set is 6 km/s. In addition, sometimes, it cannot hold true that more gravity assist planets can reduce the total fuel expense. For example, in the experiment, we find that the fuel expense of planet sequence EVVS is higher than that EJS, while the fuel expense of planet sequence EVVEJS is lower than that EJS.
5.3. Experiment Result
The planet sequence we gain with method proposed is EVVEEJS. Every time the experiment result obtained by EP_DE algorithm seems to be not stable. In several experiments we have finished, there are three different results shown in Table 17.
5.4. Analysis
The feature of the optimal solution is obvious when seeking the optical solution of multiple gravity assist trajectory. It looks like looking for a needle in a sea; that is to say, it is so hard to find the optimal solution. In this paper, when fixing the gravity assist planet sequence and looking for the optimal solution with EP_DE algorithm, the method proposed was proved to be effective, although sometimes there may not be perfect performance. The value of Cassini obtained in this paper is better than that announced by ACT in ESA. Likewise, this method can be used for any problems of multiple gravity assist trajectory design.
6. Conclusions
All of the above experiments can illustrate the efficiency of the method proposed in this paper, especially the availability of EP_DE algorithm. Experimental result of Cassini has the same function value as that announced in the database of ACT in ESA, though the specific solution is different from that. Also, the experimental result of GTOC1 also performs well, even if the accuracy of the calculation depends on the step size. For Cassini, if minimum energy is used as the optimization objective, the planets sequence EVVEJS in the benchmark may not be the best selection, while the function value 4.425189 km/s based on the planets sequence EVVEEJS found in this paper is better than 4.93 km/s which is announced by ACT.
It is a new research trend for heuristic algorithms based on machine learning. EP_DE algorithm combining search space exploring and PCA proposed in this paper is the basic research work for it. But, in general, there are also some limitations of EP_DE algorithm. For example, the processes of data statistics and analysis will cost lots of time and computing resources. However, this problem can also be solved by means of the distributed computing and development of hardware technology.
The purpose of the machine learning here is to obtain more information from less data. In the future research work, in order to improve the calculation accuracy of EP_DE algorithm, deep learning, decision tree, neural network, and reinforcement learning can also be added into the EP_DE algorithm.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by Natural Science Foundation of China under Grant nos. 41571403, 61472375, and 61103144.