1st ACT global trajectory optimization competition: Tsinghua University results

Acta Astronautica 61 (2007) 735 – 741 www.elsevier.com/locate/actaastro

1st ACT global trajectory optimization competition: Tsinghua University results Yanheng Li∗ , Hexi Baoyin, Shengping Gong, Yunfeng Gao, Junfeng Li School of Aerospace, Tsinghua University, Beijing 100084, China Available online 18 May 2007

Abstract We present the results obtained for the 1st ACT Global Trajectory Optimization Competition. To simplify this complicated interplanetary trajectory design problem, the escape velocity, the thrust history and the impact point were imposed, and gravity assists were not considered. A genetic algorithm was employed to implement the optimization. The initialization of the algorithm and our result are given in detail. Though the method proved practicable for the imposed problem, final rank of the competition showed that our result was far from global optimum due to these superfluous simplifications. Our further work after the competition is also mentioned in brief. © 2007 Elsevier Ltd. All rights reserved.

1. Introduction The possibility of an asteroid impacting with our planet in the near future has been widely discussed. For instance the asteroid 2004VD17 may collide with the Earth within the next century [1] with catastrophic consequences. Compared with explosion or vaporization, an orbital deflection seems to be a more promising approach to eliminate this threat [2]. All the deflection strategies act on the asteroid by means of some force interactions (gravitational [3], propulsion [4] and impact [5]) to increase the miss distance between the Earth and the asteroid to some predefined values. The concept of asteroid deflection was embodied in the 1st ACT Competition on Global Trajectory Optimization [6]. The goal was to intercept a target asteroid, TW229, with a low-thrust impactor while maximizing the impact energy. The complexity of this global ∗ Corresponding author.

E-mail address: [email protected] (Y. Li). 0094-5765/$ - see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.actaastro.2007.03.008

optimization problem is considerable given the length of the flight duration, the width of the launch window and the possibility of multiple planet swing-bys. All these facts contribute to an increasing number of local optima. We decided to participate in this highly challenging competition giving the good learning and exchange opportunity. Since our team lacks experience and commercial software tools for this problem, a direct Earth to asteroid strategy was employed. Several significant assumptions were made to simplify the problem, which also narrowed the solution space. Our methods and results during the competition are presented in detail in this paper. After the competition and the subsequent workshop, we did more work based on a further understanding of the problem. The results of our improved method are also given to illustrate the potential of our idea. 2. General description of the method Genetic algorithms, based on principles of natural selection and genetics, are a class of adaptive stochastic

736

Y. Li et al. / Acta Astronautica 61 (2007) 735 – 741

global optimization algorithms involving search and optimization. They have some fascinating features [7,8]. Firstly, derivative information is not needed; therefore the objective functions are not required to be continuous over the search domain. Secondly, Genetic algorithms do not need to know the mechanism of the problem or the solution structure; they just need to recognize a good solution when they see it. And lastly, the parallel and stochastic nature of Genetic algorithms makes them more likely to locate a global optimum than traditional techniques, and the performance is much less sensitive to initial conditions. However, Genetic algorithms also have some apparent disadvantages, for example, the convergence of Genetic algorithms is usually slower and the solutions are less accurate than other numerical methods. Genetic algorithms had been used in low-thrust interplanetary trajectory optimization before the competition [9–11]. They were chosen as our algorithm mainly because their features fit this problem well. It is a global optimization problem which is discontinuous about the control variables and we had no idea about the solution characteristics. All these difficulties can be handled in principle by Genetic algorithms. Besides, the only professional optimization software at our hand is Matlab姠7.1 optimization toolbox. It provides a genetic algorithm toolbox called “gatool” which had been used for optimal control problems in our research not long ago and we were familiar with it. Therefore it was employed to tackle the problem. The approach of incorporating Genetic algorithms into the problem was inspired by the work of Wuerl et al. [10] and Wall et al. [11] on the Earth to Mars trajectory optimization. They divided the trajectory into several coasting and thrusting arcs in which the structure of the thrust switch is pre-defined. We decided to impose this on the problem of the trajectory. In addition, the thrust vector was fixed along the tangent of the trajectory. And last, the problem was reduced to be a planar one. Due to these simplifying assumptions the final solution will be, inevitably, relatively far from the global optimum. 3. Application of the method to the problem 3.1. Analysis of the problem The objective function of the problem is J = mf |Uf · Vaf |,

(1)

where mf is the final mass of the spacecraft, Vaf is the velocity of the asteroid at the impact time, Vsf is the

Vsf

Uf

Vaf Fig. 1. Geometrical relationship of Uf ,Vsf and Vaf .

velocity of the spacecraft at the impact time, and Uf equals Vsf − Vaf . The geometrical relationship between Uf , Vsf and Vaf is shown in Fig. 1 together with some important angle definition. According to the definition of J, we should (a) reduce the fuel consumption (increase mf ), (b) impact near the perihelion of the asteroid (increase Vaf ), (c) increase the impact angle (the angle  between Vsf and Vaf ). The first two tasks are easy to understand. The following formula can better explain the third one:  Jv = |Uf · Vaf | = Vaf2 + Vsf2 − 2Vaf Vsf cos | cos |. (2) Fig. 2 gives the –Jv relationship, where Uf , Vf and Jv are dimensionless. From Fig. 2 we see that whether Vaf is greater than Vsf or not, a retrograde orbit ( >90◦ ) is needed to get a high J. Note that the use of gravity assist would be extremely beneficial in achieving large values of  and mf at the same time, so at the end we can expect that the result of direct Earth to asteroid strategy will be remarkably lower than that of a swing-by strategy. 3.2. Description of our method At a preliminary stage, a pre-imposed flying-andcontrol strategy without any swing-bys is adopted. It is a simple four-phase procedure composed of a pre-launch phase, an escape-orbit coasting phase, a thrust phase and a final coasting phase. All the thrust and the escape velocity are assumed to lie in the ecliptic plane and be fixed along the tangent of the trajectory. The problem is therefore reduced to a two-dimensional one. With the above impositions, the impact point is confined to

Y. Li et al. / Acta Astronautica 61 (2007) 735 – 741

2.5

737

2.5 2

2 Jv

Jv

1.5 1.5

1 1

0.5 0

0.5 0

50 100 150 Impact angle (deg)

0

50 100 150 Impact angle (deg)

Fig. 2. Relationship between impact angle  and objective function J: (a) Vaf > Vsf (Vaf = 1.5, Vsf = 1); (b) Vaf < Vsf (Vaf = 1, Vsf = 1.5).

impact condition [11]. The compound fitness function for Genetic algorithm is J¯ = −J + kd f ,

(4)

where k is the penalty coefficient. 3.3. Initialization of the problem Although Genetic algorithms have good performances in global optimum location, sensible initial guess for the optimization variables should be made to narrow the search scope and reduce the cpu time. Fig. 3. Illustration of the impact points.

be the intersection of the ecliptic plane and the asteroid orbit. Two possible impact points can be expected (Fig. 3). According to the Vaf –J relationship, we choose the point that is nearer to the sun (noted as P ), whose heliocentric ecliptic position vector is XP = [1.4505, −1.8336, 0] AU. T

(3)

There are five important time nodes in our strategy, labeled as t0 , t1 , t2 , t3 and t4 . Events on the nodes are listed in Table 1. The whole mission duration is divided into four phases with these five time nodes. Actions which the spacecraft should take in each phase are listed in Table 2. The problem is now reduced to a simple parameter optimization problem with only four parameters t01 , t12 , t23 and t34 to optimize. Genetic algorithm can be employed to carry out the optimization. The final distance df between the spacecraft and asteroid at t4 is introduced as a penalty term to meet the terminal

3.3.1. Initial guess for total mission duration (t04 ) Viewed from t0 , the asteroid will be on the impact point P in 1458.752 days, after that it will pass P every 1522.226 days (the period of the asteroid), therefore t04 = 1458.752 + 1522.226N days (N = 1, 2, 3, . . . , 12).

(5)

Due to the restriction of the problem (launch between 2010 and 2030, impact before 2060), N ranges from 1 to 12. 3.3.2. Initial guess for t01 According to the launch time restriction, the range of t01 can be given by t01 ∈ [1597, 9266] days.

(6)

The distance between the sun and the impact point P is 2.338 AU, which is outbound for the Earth, and for maneuver without gravity assist, the most effective way to increase the aphelion distance of an orbiter is to give it a delta V on perihelion. Therefore the position of the Earth at launch should be around the opposite of P about the Sun, which is noted as point L (Fig. 4). Through

738

Y. Li et al. / Acta Astronautica 61 (2007) 735 – 741

Table 1 Events on time nodes t0

t1

t2

t3

t4

MJD53600 (08/18/2005)

Launch

Thrust beginning

Thrust end

Impact

Table 2 Action in each phase of the mission Phase 1 (duration t01 )

Phase 2 (duration t12 )

Phase 3 (duration t23 )

Phase 4 (duration t34 )

Wait for launch

Coast along the escape orbit

Start up the propulsion system

Shut down the propulsion and coast until impact

and the initial mass of the spacecraft m0 = 1500 kg, we have t23 <

m0 = 10633.681 days = 29.1126 years. m ˙ max

(9)

The impact point P is even outbound for the spacecraft with the maximum escape velocity (2.5 km/s). Therefore the thrust duration t23 must have a lower boundary to make the spacecraft reach this heliocentric distance. This boundary can be got through a simple line search, with t23 as the search variable and t12 fixed to be some small constant (within the period of escape orbit). The condition that the aphelion distance of the final orbit equals 2.338 AU should be met in the search. The lower boundary is t23 > 5.97 years. Fig. 4. Illustration of launch window.

preliminary implementations of our method, we found that the launch position of the Earth in the best solutions laid in an interval whose length is about 90 days and was approximately symmetrical about point L. In later implementations the width of the window was extended to be 120 days lest some solutions should be missed. It was supposed to be symmetrical about L, that is, t01 ∈ [1597, 1685] ∪ [1990 + 365.25M − 60, 1990 + 365.25M + 60] ∪ [9235, 9266] days (M = 0, 1, . . . , 19). (7) 3.3.3. Initial guess for t23 From the maximum mass flow rate m ˙ max =

Tmax = 1.632 × 10−6 kg/s, Isp · g

In practice, this boundary is approximately set to be 6 years, which equals 2191 days. With the above initializations N can be re-ranged from 2 to 12. 3.3.4. Initial guess for t12 Apparently, the lower boundary of t12 is 0, which means that thrust will begin right after launch. The upper boundary of t12 can be determined according to the lower boundary of t01 and t23 , i.e. t12 < t04 − 1597 − 2191 days.

(11)

3.3.5. Other initialization Duration t34 is not an independent variable which can be determined by N, i.e. t34 = t04 − t01 − t12 − t23 .

(8)

(10)

(12)

Because N is an integer within a small interval, we use enumeration rather than optimization to find the best N.

Y. Li et al. / Acta Astronautica 61 (2007) 735 – 741

739

Table 3 Gatool setup Population type

Population size

Selection operator

Mutation operator

Crossover operator:

Stop criteria

Double vector

20

Stochastic uniform

Gaussian

Heuristic

1500 generations/df < 0.001 au

Table 4 Best solutions for each N N Best J

(109 J)

2

3

4

5

6

7

8

9

10

11

12

80.8

70.1

72.5

89.0

86.5

83.6

75.2

75.0

81.4

86.1

74.6

This reduces the Genetic algorithm optimization variables to t01 , t12 and t23 , which are all real numbers. This treatment can facilitate the coding and the implementation of the algorithm. 3.4. Gatool parameter setup and other techniques

Table 5 Best solution of the problem found during the competition N

t01 (day)

t12 (day)

t23 (day)

t34 (day)

5

4232.133

2263.369

2424.547

149.808

(109 J)

mf (kg)

df (km)

J

1157.99

2026

89.0

 (deg) 24.2

The detailed setup of gatool is listed in Table 3. These choices are all provided by gatool and the setup was made according to our experience in the optimal control problems. Since the simplified problem is a planar trajectory optimization problem, polar coordinates are used as the state variables of the equation of motion. The penalty coefficient k in the fitness function (4) is initialized with the value 0.0025. It was found that if k <0.0024, the method will fail to converge with a desirable df , which can be up to the order of 108 km and far beyond tolerance (0.01 au). After each implementation, k is multiplied by 10 to decrease df . Usually two to three implementations of the algorithm are needed to make df less than the order of 104 km. 4. Competition results After a week of debugging and computation, we acquired the best solution of the simplified problem. The best solutions for each N are listed in Table 4. The particular data of the best solution for all N are given in Table 5. The spacecraft will be launched in March, 2017, and impact in June, 2030. The total flight duration is 13.25 years. The trajectory of the optimal solution is showed in Fig. 5. The results in Table 4 cannot be insured to be global optima. Since our method introduces penalty term, the solution is sensitive to the penalty coefficient k and tends to prematurely converge to local optima if k is not

Fig. 5. Trajectory of the best solution.

optimal [12]. In addition, different implementations can lead to different results with the same initialization. Therefore optimum could possibly be left out. For instance, after the competition, we found out that for N = 6, there is a solution in which J = 92.0 × 109 J and for N = 11, there is a solution in which J = 92.7 × 109 J. Though these results are slightly better than the submitted result, they were missed by our method. It

740

Y. Li et al. / Acta Astronautica 61 (2007) 735 – 741

Table 6 Parameters of the best result t01 (day)

t12 (day)

t23 (day)

N

esc (deg)

1 (deg)

2 (deg)

3 (deg)

4 (deg)

3485.647

789.261

8962.530

8

−34.57

48.64

207.14

289.580

334.735

mf (kg)

df (km)

J

(109 J)

Launch date

Flight time

Impact date

847.27

1184968

244

03/04/2015

27.79 years

18/12/2042

indicates that to get a solution closer to the global optimum, more restarts should be taken and more precise initialization is needed. For this problem, five restarts for each N were taken. Based on the analysis above, it can be expected that this result is far from globally optimal according to the low impact angle. The one-thrust-arc strategy causes the low eccentricity of the osculating orbit of the spacecraft and therefore its flying path angle is small. The impact angle is determined by the inclination between the spacecraft and asteroid’s velocity on impact. With the asteroid’s flying path angle at P being only −14◦ , it is unrealistic to expect a large impact angle. We had some plans to improve our strategy in the last couple of days during the competition, for instance, to adjust the escape velocity angles and to increase the number of the thrust arcs, but they cost a longer time to implement on a PC and were not finished before the deadline. 5. Post-competition results After the release of the competition results and the workshop, we knew about other teams’ work, especially from teams who employed direct Earth to asteroid strategy [13–15]. The best result for direct transfer strategies1 was 350 × 109 J achieved by Lomonosov’s Moscow State University [14]. Since we lack the software tools and solvers they own, their results could not be repeated. However, they inspired us to exploit the potential of our own method. Take Lomonosov’s Moscow State University results for example, they acquire a highly eccentric final orbit through acceleration and deceleration near the apsides. The control law is acquired through optimal control theory [14]. Enlightened by this type of trajectory, we got a better result using direct method. Some simplifications, such as planar escape velocity 1 The best direct transfer was actually found by Moscow Aviation Institute, but not listed in the final rankings as they provided a better trajectory using fly-bys.

and planar tangential thrust, are kept. But a multiple acceleration–deceleration arc is employed instead of the one thrust arc strategy used in the competition. The start and end angles of the thrust in all revolutions are unified and they should be optimized. There are four angles in total: (1) (2) (3) (4)

acceleration start angle 1 , acceleration end angle 2 , deceleration start angle 3 , deceleration end angle 4 .

All of these angles are described in the heliocentric ecliptic frame. Besides, the direction of the escape velocity esc is also incorporated as an optimization variable. Now there are nine variables to optimize which are t01 , t12 , t23 , N, 1 , 2 , 3 , 4 and esc . Here t23 starts from the beginning of the first thrust arc and ends at the end of the last thrust arc. In the process of implementation, minimum heliocentric distance constraint (0.02 AU) is not considered. After the best solution is found, this restriction is incorporated into the fitness function of Genetic algorithm as a penalty term to meet the problem’s requirement. Necessary initial guesses are also made, and the best result can be acquired with J = 244 × 109 J. Solution parameters of the improved strategy are listed in Table 6. The trajectory of the result is shown in Fig. 6. The impact angle is 45.9◦ and the minimum heliocentric distance is 0.2003 AU. It can be seen that for this improved strategy, though more fuel is consumed, the larger impact angle compared with that of the original method leads to a much better solution. 6. Conclusions For a direct Earth to asteroid strategy, it is hard and time consuming to achieve a retrograde orbit without swing-bys. From the analysis and results of the competition participants who employ non-swing-by strategy, it can be concluded that the key point to achieve a large J is to increase the impact angle, which

Y. Li et al. / Acta Astronautica 61 (2007) 735 – 741

References

x 108 3 2 1 Y(km)

741

0 -1 -2 -3 -4 -4

-2

0 X (km)

2 x108

Fig. 6. Trajectory of the best result of the improved method.

can be achieved through increasing the eccentricity of the spacecraft’s final orbit. For our original method, the one thrust arc strategy leads to low eccentricity of the final orbit and thus causes the low value of objective function. An effective way to get a more eccentric orbit is to accelerate the spacecraft near the perihelion and decelerate near the aphelion of its osculating orbit. This can be easily realized in a direct method and was employed in our post-competition work. The improved method, though consumes more fuel, can get a much larger impact angle and the result is significantly improved. The Genetic Algorithm, which is employed by our method, basically shows its competence in solution location. However, even in our simplified strategy it cannot converge to the global optimum every time according to our results. Good initialization and several restarts are needed to avoid convergence to local optima, which also increases the workload of the method and reduces its reliability.

[1] http://impact.arc.nasa.gov/news_detail.cfm?ID = 167, 24/04/ 2006. [2] T.J. Ahrens, A.W. Harris, Deflection and fragmentation of nearEarth asteroids, Nature 360 (1992) 429–433. [3] E.T. Lu, S.G. Love, Gravitational tractor for towing asteroids, Nature 438 (7065) (2005) 177–178. [4] D.J. Scheeres, R.L. Schweickart, The mechanics of moving asteroids, AIAA paper 2004-1446 (2004). [5] C. McInnes, Deflection of near-Earth asteroids by kinetic energy impacts from retrograde orbits, Planetary Space Science 52 (2004). [6] D. Izzo, 1st ACT global trajectory optimisation competition: problem description and summary of the results, Acta Astronautica (2007), this issue, doi: 10.1016/j.actaastro.2007. 03.003. [7] D.E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, Reading, MA, 1989. [8] M. Mitchell, An Introduction to Genetic Algorithms, The MIT Press, Cambridge, MA, 1996. [9] G.A. Rauwolf, V.L. Coverstone-Carroll, Near-optimal lowthrust orbit transfers generated by a genetic algorithm, Journal of Spacecraft and Rockets 33 (6) (1996). [10] A. Wuerl, T. Crain, E.B. Braden, Genetic algorithm and calculus of variations-based trajectory optimization technique, Journal of Spacecraft and Rockets 40 (6) (2003). [11] B. Wall, B. Conway, Near-optimal low-thrust Earth–Mars trajectories via a genetic algorithm, Journal of Guidance, Control and Dynamics 28 (5) (2005). [12] N. Yokoyama, S. Suzuki, Modified genetic algorithm for constrained trajectory optimization, Journal of Guidance, Control and Dynamics 28 (1) (2005). [13] B. Dachwald, A. Ohndorf, 1st ACT global trajectory optimisation competition: results found at DLR, Acta Astronautica (2007), this issue, doi: 10.1016/j.actaastro.2007.03.011. [14] I.S. Grigoriev, M.P. Zapletin, E.V. Zapletina, 1st ACT global trajectory optimisation competition: results found at Moscow State University, Acta Astronautica (2007), this issue, doi: 10.1016/j.actaastro.2007.03.002. [15] T. Dargent, V. Martinot, 1st ACT global trajectory optimisation competition: results found at Alcatel/Alenia, Acta Astronautica (2007), this issue, doi: 10.1016/j.actaastro.2007.03.010.