How to solve GTOC1 today

In this blog entry we review methods described in (1) "1st ACT global trajectory optimization competition: Results found at the Jet Propulsion Laboratory" (Acta Astronautica 61 (2007) 806–815)" to solve the GTOC 1 problem back in 2006. Different aspects of the methods used back then are compared to what we regard more suitable to the set of tools available today. 

Analysis of the JPL method

How did JPL determine the planet sequence ?

JPL used forward search using a sims flanagan approximation of the required continuous thrust transfers between the GA maneuvers. Forward search is motivated by the fact that there are only very few valid transfers with v_inf <= 2.5 km/s because of the very low max_thrust constraint of 0.04N. Problem is to find a search heuristics to be used for pruning targeting at a high impact velocity at the asteroid. Searching backwards from the impact can be motivated in a similar way: Very few valid transfers to the asteroid with impact velocity >= 52.0 km/s, but we have to find a heuristics targeting a low v_inf at earth. Advantage of the backwards search is that if we tighten the impact requirement even more (>= 53 km/s), we further reduce the size of the search tree. Below we will compare sequences obtained by forward and backward search and compare the resulting J-values by applying identical methods after the determination of the planet sequence. 

Lambert coast arcs or multiple impulses per leg during search/optimization?

Next question is how accurate our continuous thrust approximation needs to be. Sims Flanagan (multiple impulses per leg) requires an expensive local optimization step to determine the optimal impulses. Using a single coast arc between the GAs and impulses only before and after the GA maneuvers is much cheaper, requires a single application of the Lambert algorithm. The final optimized results presented in (1) use mostly coast arcs between the GAs caused by the very low max_thrust requirement. This could be viewed as a hint coast arc approximations are sufficient for the GTOC1 planet sequence search. 

Which Lambert algorithm should be applied for the search/optimization

As already mentioned in a previos blog entry GTOC1 Lambert the method used in ESAs GTOC1-GTOP benchmark code AstroToolbox is oversimplifying, it is better to use a lambert solver providing all solutions like PYKEP Lambert Solver which would convert the GTOC1-GTOP benchmark into a mixed integer optimization problem. When used for search, each solution of the Lambert solver opens a separate sub branch. 

Beam search or global optimization?

Application of a search algorithm for the determination of suitable planet sequences seems undisputed, but what about the phasing / timing of the transfers? JPL used search also in this area, where the GTOP GTOC1 benchmark uses a global optimization algorithm instead. JPLs decision back in 2006 was justified by the limited availability of good optimization algorithms at that time. Another advantage is that the search for planet sequences and for good transfer timings can be performed in a single step. 

But replacing search by optimization has important advantages:

• Memory consumption: Depth first search is not suitable to trajectory optimization and beam search requires a lot of memory as soon as you increase the number of branches. This is not the case with optimization algorithms. 
• Ease of use: Beam search requires sophisticated heuristics required for branch pruning where optimization algorithms only require a target function as long as this target function is deterministic. But even with stochastic - heuristics driven - target functions optimization is easier to apply. 
• Parallelization: Both methods (search and optimization) can distribute workload on separate threads, but if the load needs to be distributed on multiple processors, synchronization of the results is easier with optimization. We will not present details here, because current processors support up to 64 parallel threads so it is much less attractive to use multiple processors these days.  
• Separation of the phasing/timing optimization from the search for good planet sequences simplifies the latter. We used fixed planet orbits for this search which means we have to compute lambert transfer tables only for a single planet revolution. The timings we obtain this way are not accurate, but since our aim was only to identify good planet sequences this doesn't hurt.   

Are old J-values comparable to actual ones?

Simple answer: Yes, if the same method to compute the planet orbits is used. The GTOP benchmarks uses ESAs Astrotoolbox which uses the following approximation AstroToolbox . JPL doesn't elaborate on the method they used, today it probably would be data from NASA Horizons  . We tried to apply horizons data to the Earth and Venus encounters of the first trajectory leg of JPLs winning solution, but couldn't find a valid low thrust transfer with v_inf <= 4.5 km/s independent from the velocity at Venus. Same result with ESAs Astrotoolbox planet orbit approximation. Our optimization led to inferior results for both JPL planet sequences, but we could reproduce JPLs J-value for the DEIMOS sequence - which is better than the value from the solution submitted by DEIMOS. Maybe someone could comment on the method used by JPL for planet orbit determination back in 2006?

Anyway, for a competition like GTOC the underlying data should be exactly defined and easy to compute, as it was the case for all GTOC competitions after GTOC1. We propose to use ESAs Astrotoolbox planet orbit approximation for further comparisons. Unfortunately the old results are not comparable to the new ones, but we can compare the planet sequences submitted by the best teams using the same new optimization method based on ESAs Astrotoolbox planet orbit approximation. 

Visualization of the planet sequences

The following pictures and videos visualize the best solutions we found for the given planet sequences.

JPL - EVEEEJSJA

Winning planet sequence submitted by JPL. See a video here https://youtu.be/JhJrPgbnrdE

DEIMOS - EVVEEVVEVEJSJA

Planet sequence which reached the 2nd place submitted by DEIMOS. See a video here https://youtu.be/9j8oLD1SSro 

JPL2 - EVVEEEEJSJA

Second JPL sequence from JPLs paper from 2007. See a video here https://youtu.be/JMPOHqOevWk 

GMV - EEVEEJSA

Sequence submitted by GMV which reached the 3rd place. See a video here https://youtu.be/9VokaMpBPRA 

JENA - EVVEVVEESJA

Sequence determined by Team Jena after the competition using backwards search. See a video here https://youtu.be/zk75TaJKG_8

Comparison of J-values using ESAs Astrotoolbox

Now we want to compare different planet sequences using ESAs Astrotoolbox planet orbit approximation. We want to show how J-values improve as we refine our optimization algorithm. No search will be used to determine the transfer timings, instead we use a global optimization algorithm followed by a continous thrust conversion of the resulting Lambert transfers. 

GTOP GTOC1

The simplest method is the one implemented as part of ESAs GTOP benchmark suite GTOP GTOC1 .

We modified the code to support also other planet sequences:

1. GMV - EEVEEJSA, J = 1708327
2. JPL - EVEEEJSJA, J = 1680247
3. JENA - EVVEVVEESJA, J = 1673184

These are the sequences submitted by JPL and GMV, together with the sequence we determined using backward search using impact velocity >= 53 km/s as constraint. If we finally find a timing with a very low overall dv budget, there are good chances this sequence can improve its current rank. GMV submitted a solution with only J= 1455000, it would be interesting to find out the reason for this. As we will see later it is not the missing low thrust conversion. GMV reported a low asteroid impact velocity of 48.5 km/s indicating some problem with their search algorithm GIOptImp. We find a first hint that JPL planet sequences don't like ESAs Astrotoolbox planet orbit approximation.

Using the new Lambert solver

We improved these result by replacing the used Lambert algorithm with PYKEP Lambert Solver thereby converting the problem into a mixed integer one:

1. GMV - EEVEEJSA, J = 1828322
2. JENA - EVVEVVEESJA, J = 1813613
3. JPL - EVEEEJSJA, J = 1810572

As optimization algorithm CMA-ES with retries was used. 

Using GCMA-ES 

Next we add the overall 30 years mission time limit and apply the new GCMA-ES optimization algorithm. We also introduce the DEIMOS planet sequence, where JPL already showed back in 2007 that it is superior to the sequence they submitted to the competition. 

1. JENA - EVVEVVEESJA, J = 1894013
2. Deimos - EVVEEVVEVEJSJA, J = 1887311
3. GMV - EEVEEJSA, J = 1846951
4. JPL - EVEEEJSJA, J = 1834679

This result is still based on a ballistic trajectory, no continuous thrust applied.
Specially the huge leap for the JENA sequence shows how much potential we missed before GCMA-ES was available. Improvements in optimization algorithms change the way how these kind of problems can be solved. Back in 2006 optimization was not powerful enough for GTOC1, so each team used search instead. 

Converting to a continuous thrust model

Continuous thrust conversion is difficult here because of the very low 0.04N max thrust limit. It succeeded only for the GMV and JPL impulse solutions. For Deimos and Jena an inferior impulse solution had to be chosen.  

1. Deimos - EVVEEVVEVEJSJA, J = 1861309
2. JENA - EVVEVVEESJA, J = 1847069
3. GMV - EEVEEJSA, J = 1795361
4. JPL - EVEEEJSJA, J = 1766122

We apply a simple conversion method: We keep the mission timing from the ballistic solution and adapt incoming and outgoing velocities from the Lambert arcs to fulfill the GA constraints (minimal rp / maximal angle). We divide each trajectory leg into 25 segments of equal length where we apply identical thrust vectors. Using CMA-ES these 25 thrust vectors were optimized to fulfill the max thrust constraint and to minimize propellant usage. ESAs Taylor propagator PYKEP Taylor is used for propagation. 

Using a Sims Flanagan based intermediary optimization phase

Now we improve the conversion method by introducing an intermediate optimization step. Very similar to the low thrust conversion optimization we use again 25 segments of equal length for each leg, but apply an impulse based approximation (Sims Flanagan) to speed up the computation. ESAs Lagrangian propagator PYKEP Lagrangian is used for orbit propagation, but we eliminated the repeated computation of sub-expressions like "cos(DE)" which leads to a 10% speedup using -O3 with the clang 6.0 compiler. In an outer CMA-ES optimization we vary the timings (+/- 10 days) and the incoming and outgoing velocities, where the inner Sims Flanagan optimization determines the propellant consumption. Advanced caching techniques are used to avoid recomputions of the inner optimization. Result of this nested optimization are optimal incoming/outgoing velocities fulfilling the GA constraints together with improved mission timings. As alternative to the expensive nested optimization we could combine both optimizations into one. CMA-ES could be replaced by GCMA-ES to solve the resulting hard optimization problem, something we  plan to try later. 

The optimized timings/velocities are then used as input to the (unchanged) low thrust conversion algorithm. We also introduce a 5th planet sequence called JPL2, the best sequence JPL found in (1), better than the sequence JPL submitted for the competition.

1. JENA - EVVEVVEESJA, J = 1898460
2. DEIMOS - EVVEEVVEVEJSJA, J = 1873294
3. JPL2 - EVVEEEEJSJA, J = 1847502
4. GMV - EEVEEJSA, J = 1819483
5. JPL - EVEEEJSJA, J = 1787113

Finally the JENA sequence was able to exploit its "impact = 53.18 km/s" potential. The JPL sequence suffers from problems with the starting EVE sequence which turned up to be "incompatible" with ESAs Astrotoolbox planet orbit approximation. 

Final results in detail:

JENA - EVVEVVEESJA, J = 1898460

v_inf = {-2,2097686852; 1,1226063577; 0,3179945906} = 2.499 km/s
mission time = 29.546 years
impact = 53.18 km/s
J = 1898460.6314493746
1: mass = 1499.999 fuel = 0.001 from = 2022-05-14T11:49:14.816 to = 2024-06-03T18:38:50.816 tf = 751.284 days
2: mass = 1499.998 fuel = 0.0 from = 2024-06-03T18:38:50.816 to = 2026-04-01T18:29:14.816 tf = 666.993 days
3: mass = 1499.998 fuel = 0.0 from = 2026-04-01T18:29:14.816 to = 2028-06-07T02:42:02.816 tf = 797.342 days
4: mass = 1499.996 fuel = 0.002 from = 2028-06-07T02:42:02.816 to = 2028-11-10T04:51:06.816 tf = 156.09 days
5: mass = 1497.061 fuel = 2.936 from = 2028-11-10T04:51:06.816 to = 2031-03-08T22:07:54.816 tf = 848.72 days
6: mass = 1492.971 fuel = 4.09 from = 2031-03-08T22:07:54.816 to = 2032-10-02T00:41:30.816 tf = 573.107 days
7: mass = 1477.806 fuel = 15.165 from = 2032-10-02T00:41:30.816 to = 2039-12-29T07:47:06.816 tf = 2644.296 days
8: mass = 1463.495 fuel = 14.311 from = 2039-12-29T07:47:06.816 to = 2042-04-14T02:56:58.816 tf = 836.799 days
9: mass = 1462.606 fuel = 0.888 from = 2042-04-14T02:56:58.816 to = 2050-06-23T06:15:22.816 tf = 2992.138 days
10: mass = 1462.045 fuel = 0.562 from = 2050-06-23T06:15:22.816 to = 2051-11-30T03:57:02.819 tf = 524.904 days

DEIMOS - EVVEEVVEVEJSJA, J = 1873294

v_inf = {-2,3484457158; -0,2427303926; -0,8221218154} = 2.5 km/s
mission time = 29.41 years
impact = 52.68 km/s
J = 1873294.0160389685
1: mass = 1466.274 fuel = 33.726 from = 2022-06-30T04:22:18.816 to = 2024-05-22T04:28:42.816 tf = 692.004 days
2: mass = 1462.898 fuel = 3.375 from = 2024-05-22T04:28:42.816 to = 2026-12-23T14:45:14.816 tf = 945.428 days
3: mass = 1460.858 fuel = 2.04 from = 2026-12-23T14:45:14.816 to = 2027-10-09T03:48:10.816 tf = 289.544 days
4: mass = 1460.857 fuel = 0.001 from = 2027-10-09T03:48:10.816 to = 2029-04-08T01:10:18.816 tf = 546.89 days
5: mass = 1460.856 fuel = 0.0 from = 2029-04-08T01:10:18.816 to = 2030-09-10T00:31:54.816 tf = 519.973 days
6: mass = 1460.856 fuel = 0.0 from = 2030-09-10T00:31:54.816 to = 2032-04-17T08:24:26.816 tf = 585.328 days
7: mass = 1460.855 fuel = 0.0 from = 2032-04-17T08:24:26.816 to = 2034-01-24T19:24:42.816 tf = 647.459 days
8: mass = 1460.855 fuel = 0.0 from = 2034-01-24T19:24:42.816 to = 2035-05-10T22:21:46.816 tf = 471.123 days
9: mass = 1460.855 fuel = 0.0 from = 2035-05-10T22:21:46.816 to = 2038-12-15T19:52:26.816 tf = 1314.896 days
10: mass = 1460.85 fuel = 0.005 from = 2038-12-15T19:52:26.816 to = 2040-02-29T04:52:10.816 tf = 440.375 days
11: mass = 1460.579 fuel = 0.271 from = 2040-02-29T04:52:10.816 to = 2041-06-23T19:05:30.816 tf = 480.593 days
12: mass = 1460.07 fuel = 0.509 from = 2041-06-23T19:05:30.816 to = 2050-05-31T03:41:46.816 tf = 3263.359 days
13: mass = 1459.096 fuel = 0.974 from = 2050-05-31T03:41:46.816 to = 2051-11-27T07:14:10.789 tf = 545.147 days

JPL2 - EVVEEEEJSJA, J = 1847502

v_inf = {-0,9957788351; 2,2829028434; -0,2162848027} = 2.5 km/s
mission time = 29.622 years
impact = 52.665 km/s
J = 1847502.2806930426
1: mass = 1447.263 fuel = 52.737 from = 2022-04-16T09:18:50.816 to = 2025-01-22T04:24:26.816 tf = 1011.796 days
2: mass = 1447.263 fuel = 0.0 from = 2025-01-22T04:24:26.816 to = 2027-05-06T23:54:34.816 tf = 834.813 days
3: mass = 1447.262 fuel = 0.0 from = 2027-05-06T23:54:34.816 to = 2028-06-08T23:17:14.816 tf = 398.974 days
4: mass = 1447.262 fuel = 0.0 from = 2028-06-08T23:17:14.816 to = 2032-01-09T21:18:50.816 tf = 1309.918 days
5: mass = 1445.575 fuel = 1.687 from = 2032-01-09T21:18:50.816 to = 2033-09-02T21:48:42.816 tf = 602.021 days
6: mass = 1441.931 fuel = 3.644 from = 2033-09-02T21:48:42.816 to = 2038-11-07T20:47:54.816 tf = 1891.958 days
7: mass = 1441.137 fuel = 0.793 from = 2038-11-07T20:47:54.816 to = 2040-02-22T23:40:42.816 tf = 472.12 days
8: mass = 1440.933 fuel = 0.204 from = 2040-02-22T23:40:42.816 to = 2041-06-19T05:56:10.816 tf = 482.261 days
9: mass = 1440.337 fuel = 0.596 from = 2041-06-19T05:56:10.816 to = 2050-06-03T18:42:02.816 tf = 3271.532 days
10: mass = 1439.432 fuel = 0.905 from = 2050-06-03T18:42:02.816 to = 2051-11-29T19:42:10.785 tf = 544.042 days

GMV - EEVEEJSA, J = 1819483

v_inf = {2,4996870315; -0,0372306106; -0,0125911107} = 2.5 km/s
mission time = 28.784 years
impact = 52.12 km/s
J = 1819483.3247717093
1: mass = 1459.431 fuel = 40.569 from = 2023-03-16T15:09:46.816 to = 2027-03-11T13:42:18.816 tf = 1455.939 days
2: mass = 1448.759 fuel = 10.672 from = 2027-03-11T13:42:18.816 to = 2030-08-02T16:47:54.816 tf = 1240.129 days
3: mass = 1447.57 fuel = 1.189 from = 2030-08-02T16:47:54.816 to = 2034-08-23T10:33:30.816 tf = 1481.74 days
4: mass = 1447.569 fuel = 0.0 from = 2034-08-23T10:33:30.816 to = 2038-11-18T22:51:38.816 tf = 1548.513 days
5: mass = 1447.108 fuel = 0.461 from = 2038-11-18T22:51:38.816 to = 2040-02-29T03:31:06.816 tf = 467.194 days
6: mass = 1446.681 fuel = 0.428 from = 2040-02-29T03:31:06.816 to = 2041-07-02T11:31:06.816 tf = 489.333 days
7: mass = 1446.33 fuel = 0.35 from = 2041-07-02T11:31:06.816 to = 2051-12-27T22:21:13.707 tf = 3830.451 days

JPL - EVEEEJSJA, J = 1787113

v_inf = {-0,3295540147; 2,4681388801; 0,2229002265} = 2.5 km/s
mission time = 27.675 years
impact = 52.282 km/s
J = 1787113.0262987711
1: mass = 1455.875 fuel = 44.125 from = 2024-04-04T10:11:06.816 to = 2026-04-09T06:26:02.816 tf = 734.844 days
2: mass = 1453.01 fuel = 2.865 from = 2026-04-09T06:26:02.816 to = 2029-08-26T18:17:30.816 tf = 1235.494 days
3: mass = 1452.718 fuel = 0.292 from = 2029-08-26T18:17:30.816 to = 2032-11-10T14:18:34.816 tf = 1171.834 days
4: mass = 1452.718 fuel = 0.0 from = 2032-11-10T14:18:34.816 to = 2037-09-21T14:40:58.816 tf = 1776.016 days
5: mass = 1449.678 fuel = 3.039 from = 2037-09-21T14:40:58.816 to = 2039-04-10T12:15:54.816 tf = 565.899 days
6: mass = 1431.437 fuel = 18.241 from = 2039-04-10T12:15:54.816 to = 2040-11-30T14:55:54.816 tf = 600.111 days
7: mass = 1405.893 fuel = 25.545 from = 2040-11-30T14:55:54.816 to = 2050-06-02T17:42:18.816 tf = 3471.116 days
8: mass = 1405.201 fuel = 0.692 from = 2050-06-02T17:42:18.816 to = 2051-12-07T13:20:18.513 tf = 552.818 days