Jules Verne

In the last blog entry we introduced a new optimization algorithm - GCMA-ES  - and a new method to apply Lambert transfers for trajectory optimization - the Stochastic Lambert Solver. But we applied both methods only to variants of existing problems from the ESA GTOP optimization benchmark database GTOP Benchmarks . 

As I showed Ingo Althöfer these results, he had the idea to apply the methods to something new. Ingo is Team-Manager of Team Jena which participated in all GTOC competitions since GTOC 5. This team has two special properties: it is both the smallest and the least successful team in the top ten. Probably these properties correlate. Beside its manager Team Jena has only one worker, the author of this blog.

GTOC 10 is near, so if my manager has a training task for me, I should start working on it. And it is getting cold in German November 2018, but not in my flat as long as my AMD TR 2990WX is running full thrust utilizing all 64 threads.

Around the World in Eighty Years 

If you have a car, do you know how far you can drive using 550 kg fuel? May be about 6000-8000 km. That is not exactly the distance Ingo had in mind. In reminiscence to Around the World in Eighty Days (1872) from Jules Verne he thought about the following task:

Start with a small spaceship having a mass of 1500 kg from earth with a low start delta velocity of maximal 2.5 km/s, visit all nine planets of the solar system using only about 550 kg propellant in at most 80 years. Sounds crazy? It is, but this is not all: Unfortunately Ingo has seen Pluto Orbiter Mission .



You can do all kinds of exciting stuff as soon as you arrive at Pluto. Only problem: Delta velocity relative to Pluto at arrival time should be quite low so that you can start gravity assist maneuvers around Pluto's moon Charon. 80 years sounds like a long time, but look at Solar System to realize how far we have to travel. Ingo also thought about using solar sails but fortunately he has other friends he will ask about this. I have no clue how to slow down for the Pluto orbit insertion if the “wind” blows almost directly from behind. Maybe you have?
Remark: In the discussion whether Pluto is a planet we agree with Pluto Planet Definition .

Trajectory found

Finally we found a trajectory which can be viewed at Jules Verne Tour .

Start from Earth, head for Venus and Mercury

Jupiter, then Mars and Saturn

Neptune, then Uranus and heading for Pluto

Smooth Arrival at Pluto 

Trajectory Planning and Optimization 

How to start? There are two main options:

1. Implement a (beam-)search algorithm to search for planet sequences / timings based on Lambert arc / gravity assist sequences. We tried this successfully in the past to solve GTOC1 - the real one, not its GTOP variant which uses a fixed planet sequence. This approach relied on precomputed Lambert arc tables. For all 9 planets these tables wouldn’t fit in our memory. Therefore we opted for the second approach: 
2. Apply a global optimization algorithm for different fixed planet sequences to determine only the timings. Planet sequences have to be chosen using experience, for instance gained be applying the search approach to GTOC1. We can use ESAs AstroToolbox the GTOP benchmarks are based on to implement the function to be optimized. 

We already extended this code by a Stochastic Lambert Solver but we need some more changes:

• The GTOP code needs to be extended to cover GA maneuvers at Uranus and Neptune. 
• We have to add the 80 years overall time of flight constraint. 
• We have to check out and evaluate different planet sequences and transfer time limits. 
• Finally we have to convert the final mission into one using continuous low thrust. An impulse engine is by far not efficient enough for this tour. 

Global Optimization

To our surprise we found out, that this optimization task is significantly easier than both messenger full and the modified GTOP GTOC1 task. No impact as second objective as for GTOC1, and no DSM-1 maneuvers which complicate messenger full.

As usual (see How to utilize the TR 2990WX) we configure the AMD TR-2990WX to run 8 GCMA-ES optimization tasks in parallel each using 8 parallel threads. This way about 40x10^9 Lambert transfers are computed each hour. After 1-2 hours we collect the GCMA-ES results and use this as basis for evaluation of the planet sequences tested. We use the linux numactl command to minimize the negative impact imposed by the two NUMA nodes without direct memory access. Without numactl all cores will still show 100% usage, but both performance and power consumption will be reduced. More importantly also the room temperature will drop, something you want to avoid during a cold November.

Remark: Overclocking the TR2990WX is a very good method to further increase the room temperature since it reduces the power efficiency. In November an all cores overclock to 3.85 Ghz proved sufficient, something to be reconsidered during winter or if you plan to heat up a larger room.

We tried sequentially the the following planet sequences (m=Mercury, M=Mars):

• E m V E J S U N P 
• E V m V E M J S U N P 
• E V V m V V E M J S U N P 
• E V V m V V E M S J U N P 
• E V V m V V E M S J N U P 
• E V V m V V E J M S N U P 

in overall less than 12 h CPU time until we were satisfied with the result:

Less than 4 km/s overall delta velocity budget plus about 8 km/s for the Pluto orbit insertion. Visiting more planets makes no sense because of the 80 years limit although

E V V V m V V E J M S N U P also delivered reasonable results.

As limits for the global optimization we used:

• lower limit = [5000, 100, 100, 100, 100, 100, 100, 300, 300, 300, 300, 300, 300] 
• upper limit = [15000, 1000, 1000, 1000, 1200, 1000, 1000, 3000, 3000, 9000, 9000, 9000, 9000] 

Multiple Venus encounters are required to first slow down for Mercury and then accelerate for the outer planets. At Jupiter we have to keep distance because of its radiation. Nearly all of the needed delta velocity is provided using gravity assists. 

Lambert Arcs Resulting from Global Optimization

For the resulting Lambert arcs we list the dates of the encounters and the delta velocity values between these arcs.

• Time of flight = 80 years 
• Start at 2035-10-03T11:36:18.524 from Earth with dv = 3.354 
• Arrival at 2115-10-04T11:35:55.196 at Pluto with dv = 9.824 
• Orbit insertion dv at Pluto = 8.889 

Planet encounters:

• dv: absolute delta velocity at the planet 
• dvE: minimal impulse delta velocity required from the engine 
• rp: distance as factor of the planet radius. 

1. Venus: dv = 6.164 dvE = 0.455 rp = 1.437 at 2036-12-22T15:12:40.013 
2. Venus: dv = 5.278 dvE = 0.0 rp = 1.876 at 2038-06-12T03:41:05.304 
3. Mercury: dv = 1.738 dvE = 0.0 rp = 1.05 at 2039-12-07T09:23:54.149 
4. Venus: dv = 6.838 dvE = 0.0 rp = 1.058 at 2041-06-06T04:45:25.432 
5. Venus: dv = 6.437 dvE = 0.0 rp = 1.161 at 2043-08-23T21:22:10.182 
6. Earth: dv = 6.275 dvE = 0.374 rp = 1.063 at 2045-06-24T01:39:10.862 
7. Jupiter: dv = 2.567 dvE = 0.0 rp = 133.43 at 2047-01-17T04:12:20.721 
8. Mars: dv = 1.4 dvE = 0.293 rp = 1.766 at 2053-10-08T18:57:14.376 
9. Saturn: dv = 6.001 dvE = 0.0 rp = 17.41 at 2058-05-29T23:06:06.976 
10. Neptune: dv = 11.664 dvE = 1.81 rp = 1.1 at 2081-07-23T13:50:49.293 
11. Uranus: dv = 4.904 dvE = 0.0 rp = 6.583 at 2100-10-23T21:41:21.983 

We calculate a year as 365.25 days. Since 2100 there is no leap year, we arrive one day later as expected. Likely some better tours are possible. Who finds some?

Continuous Thrust Conversion

If we take the efficient low thrust engine from the GTOC 8 task - Isp = 5000 and maximal thrust = 0.1 N - these are very reasonable delta velocity values. We performed the conversion into a low thrust trajectory on a leg by leg basis. A global continuous thrust optimization would not result in a significant improvement, since the very high time of flight values simplifies the conversion.

1. In a first step we adjust the incoming and outgoing velocity vectors at each planet encounter obtained from the Lambert transfers so that the resulting delta velocity is fully covered by the gravity assist maneuver. This adjustment results in an incoming and outgoing delta velocity between the new (corrected) and the old (Lambert) velocities which have to be covered by continous thrust. A factor at each encounter determines how "much" of this correction is covered by the incoming velocity change, and how much by the outgoing one. A global CMA-ES based optimization is applied to determine the globally optimal set of factors. Value to optimize is the sum of the (incoming dv + outgoing dv) / (time of flight) for each leg. The optimal factors are then applied, the resulting velocities - and the timings of the lambert arcs-  are the basis of the following continous thrust optimization.
2. Then we use CMA-ES (GCMA-ES would be overkill here) to compute 25 optimal thrust vectors associated to 25 equally timed leg segments for each leg so that we arrive at correct time, position and velocity thereby minimizing fuel consumption. Using more than 25 segments only provides minimal fuel consumption gains not worth the effort for the given task. 

Processor Configuration for Continuous Thrust Conversion

Our processor is utilized as usual performing 8 low thrust trajectory conversions in parallel each using 8 threads allocated to a parallelized CMA-ES implementation. Then we chose the one with lowest fuel consumption. The whole process takes less than 30 minutes.

Note: Conversion of multiple legs of the same trajectory cannot be performed in parallel since we don't know the start mass of each leg in advance.

Continuous Thrust Conversion Results

The results scale with start mass and maximal thrust, for example you get equivalent results with start mass 3000 kg and max thrust 0.2 N.

But what happens if you try to save engine mass by equipping the 3000 kg ship with a weak 0.1 N engine? May be a bad idea, since now the slow down at Pluto is significantly harder and the trajectory needs 796 kg more propellant.

Additionally you need a more expensive rocket at start since now a start delta velocity of 2.8 km/s is required.

Result for start mass = 1500 kg and maximal thrust = 0.1 N

Start delta velocity from Earth = {0,19650; -2,14745; -1,26485}

• m: mass at end of the leg 
• f: fuel/propellant consumption for the leg 
• tf: time of flight of the leg 
• th: relative thrust level at each segment of the leg 

1. m = 1461.682 f = 38.318 tf = 446.15 th = 1.0 1.0 1.0 1.0 0.438 0.0 0.0 0.868 0.0 0.0 0.0 0.0 0.0 0.999 1.0 1.0 1.0 0.981 0.0 0.0 0.0 0.0 0.0 0.996 0.903  
2. m = 1449.236 f = 12.446 tf = 536.52 th = 0.0 0.0 0.409 0.0 0.0 0.0 0.0 0.873 0.0 0.0 0.0 0.0 0.657 0.0 0.0 0.0 0.0 0.468 0.0 0.0 0.0 0.0 0.885 0.0 0.0 
3. m = 1449.235 f = 0.0 tf = 543.238 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
4. m = 1449.235 f = 0.0 tf = 546.807 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
5. m = 1449.235 f = 0.0 tf = 808.692 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
6. m = 1440.855 f = 8.38 tf = 670.178 th = 0.095 0.0 0.0 0.0 0.0 0.0 0.014 0.235 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.023 0.23 0.0 0.0 0.0 0.0 0.0 0.532 0.64 0.0 
7. m = 1440.855 f = 0.0 tf = 572.106 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
8. m = 1428.814 f = 12.04 tf = 2456.615 th = 0.012 0.0 0.0 0.012 0.007 0.007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.659 
9. m = 1428.814 f = 0.0 tf = 1694.173 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
10. m = 1329.229 f = 99.585 tf = 8455.614 th = 0.033 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.637 1.0 
11. m = 1234.09 f = 95.139 tf = 7031.327 th = 1.0 0.832 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.088 
12. m = 950.533 f = 283.556 tf = 5458.58 th = 0.827 0.04 0.01 0.009 0.008 0.0 0.002 0.003 0.002 0.001 0.006 0.003 0.004 0.006 0.006 0.01 0.02 0.009 0.418 0.993 0.998 0.999 0.999 0.999 0.999 

Result for start mass = 3000 kg and maximal thrust = 0.1 N 

Start delta velocity from Earth = {0,270656; -2,35800; -1,485455} 

1. m = 2946.222 f = 53.778 tf = 446.15 th = 1.0 1.0 1.0 1.0 1.0 0.229 0.548 1.0 1.0 0.0 0.0 0.0 0.324 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0 1.0  
2. m = 2920.298 f = 25.925 tf = 536.52 th = 0.0 0.0 1.0 0.0 0.0 0.0 0.823 1.0 0.0 0.0 0.0 0.095 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.155 1.0 0.783 0.0 
3. m = 2920.296 f = 0.001 tf = 543.238 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
4. m = 2920.296 f = 0.0 tf = 546.807 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
5. m = 2920.296 f = 0.0 tf = 808.692 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
6. m = 2903.242 f = 17.054 tf = 670.178 th = 0.198 0.0 0.0 0.0 0.0 0.0 0.0 0.246 0.23 0.0 0.0 0.0 0.0 0.0 0.0 0.444 0.0 0.0 0.0 0.0 0.0 0.138 1.0 1.0 0.355 
7. m = 2903.242 f = 0.0 tf = 572.106 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
8. m = 2876.396 f = 26.846 tf = 2456.615 th = 0.128 0.001 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.42 1.0 
9. m = 2876.396 f = 0.0 tf = 1694.173 th = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
10. m = 2662.053 f = 214.343 tf = 8455.614 th = 0.138 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.459 1.0 1.0 1.0 
11. m = 2451.19 f = 210.863 tf = 7031.327 th = 1.0 1.0 1.0 0.88 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.375 
12. m = 1654.509 f = 796.681 tf = 5458.58 th = 1.0 1.0 1.0 1.0 1.0 0.68 0.0 0.0 0.0 0.027 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0