Solving the PYKEP trajopt/gym Reference Problems
PYKEP provides a set of reference problems for testing optimization methods. Many of them replace the older GTOP problems. These newer implementations have several practical advantages:
-
They are written in Python, so they are easy to adapt.
-
They provide useful pretty-print and plotting functions for the computed trajectories.
-
They use a modern C++ backend. For example, the Lambert solver is more flexible.
The tradeoff is evaluation cost. These problems are about five times slower to evaluate than the old C++ versions. That makes it more important to reduce the number of objective calls than to reduce optimizer overhead.
Another limitation is the lack of reference solutions. The problems are implemented as PAGMO-UDPs, see Problem/UDP, so solving them with PYGMO/PAGMO is a natural starting point. But there is little guidance on how to configure a PAGMO-Archipelago well for these cases. Example 13 shows single-threaded runs, but not a parallel optimization workflow.
Most of these problems are single-objective and unconstrained, or can be made unconstrained with unconstrain. That makes them a good fit for fcmaes. They are also a useful way to show what coordinated retry can do on real trajectory optimization tasks.
Applying fcmaes coordinated retry to PYKEP trajopt/gym
See trajopt/gym for the full list of problems.
Tandem
Tandem models a rendezvous mission to Saturn as an MGA-1DSM transfer. The mission parameters are inspired by the TandEM mission. A launcher model, Atlas 501, is included, so the main objective is the final mass delivered at Saturn.
The problem includes several planet sequences. We start with sequence 6, [earth, venus, earth, earth, saturn], which later turns out to be the best one. Tandem has two variants. Here we use the time-constrained version, with total mission time limited to 10 years.
Single-objective PAGMO problems can be wrapped with fcmaes.optimizer.single_objective. We use pg.unconstrain to turn the time limit into a penalty added to the objective value.
import pygmo as pg
import pykep as pk
from fcmaes import advretry
from fcmaes.optimizer import single_objective, dtime, logger, de_cma
def test_problem(opt, problem, num_retries = 10000, num = 3, value_limit = math.inf, log = logger()):
log.info(problem.name + ' ' + opt.name)
for _ in range(num):
_ = advretry.minimize(problem.fun, problem.bounds, value_limit, num_retries, log, optimizer=opt)test_problem calls coordinated retry num times. Each call performs num_retries optimization runs.
We can now apply the fcmaes default optimizer de_cma, a sequence of differential evolution and CMA-ES, to Tandem:
tandem = pk.trajopt.gym.tandem(6)
uda = pg.unconstrain(tandem, method = "weighted", weights=[10]*tandem.get_nic())
opt = de_cma(de_max_evals = 750, cma_max_evals = 750)
pr = single_objective(pg.problem(uda))
test_problem(opt, pr, num_retries = 20000, value_limit = 0)We assign 750 function evaluations per run to both DE and CMA-ES. This budget grows every 100 runs, up to a factor of 50. Tandem is difficult, so we increase the number of retries to 20000. value_limit = 0 means poor runs are not stored and reused as seeds for new runs. For easier problems, the default setting is usually enough.
After about three hours on an AMD 3950X 16-core machine, we get:
solution = [7985.2005463509195, 0.9030011694679589, 0.3590118628421891, 3319.4873498469165, 0.013884406331823213, 169.39189888962179, -1.151661383078681, 1.050000000000008, 0.11740075310407928, 338.512464914102, -1.6796842963039957, 1.1003653467864567, 0.6019983639984013, 1178.9023407204934, -2.0265823052439456, 1.05000000000032, 0.05055146337411492, 1762.610851218632] fitness = -7.30677225 First Leg: earth to venus Departure: 2021-Nov-11 04:48:47.204719 (7985.2005463509195 mjd2000) Duration: 169.39189888962179days VINF: 3.3194873498469164 km/sec DSM after 2.351905953502622 days DSM magnitude: 8.139343635357746e-09m/s leg no. 2: venus to earth Duration: 338.512464914102days Fly-by epoch: 2022-Apr-29 14:13:07.268783 (8154.592445240542 mjd2000) Fly-by radius: 1.050000000000008 planetary radii DSM after 39.74161831603379 days DSM magnitude: 1.513243010013376e-09m/s leg no. 3: earth to earth Duration: 1178.9023407204934days Fly-by epoch: 2023-Apr-03 02:31:04.237361 (8493.104910154643 mjd2000) Fly-by radius: 1.1003653467864567 planetary radii DSM after 709.6972804276229 days DSM magnitude: 177.02562259814547m/s leg no. 4: earth to saturn Duration: 1762.610851218632days Fly-by epoch: 2026-Jun-25 00:10:26.475612 (9672.007250875136 mjd2000) Fly-by radius: 1.05000000000032 planetary radii DSM after 89.10255788819619 days DSM magnitude: 3.035843180184029e-09m/s Arrival at saturn Arrival epoch: 2031-Apr-22 14:50:04.020902 (11434.61810209377 mjd2000) Arrival Vinf: 6017.457155646168m/s Insertion DV: 696.6382951176456m/s Total mission time: 9.443990570137847 years (3449.417555742849 days) Initial mass: [2092.76404805] Final mass: [1490.35891505] Declination: 28.49999999999617
|
Note
|
The original tandem1.png result figure is not present in this repository.
|
In practice, a second retry run is often still needed.
Messenger
Messenger models a mission to Mercury as an MGA-1DSM transfer. The selected fly-by sequence, E-VVMeMeMe-Me, and the other parameters are inspired by the Messenger mission.
As already noted for the corresponding GTOP benchmark in Readme.adoc, this is a hard problem and needs many more optimization retries:
prob = single_objective(pg.problem(pk.trajopt.gym.messenger))
opt = de_cma(de_max_evals = 500, cma_max_evals = 500)
test_problem(opt, prob, num_retries = 100000, value_limit = 12000)Here we reduce the number of evaluations per retry and compensate with more retries. That makes the optimization slower, but it broadens the search and helps avoid local minima. Messenger is especially difficult because different Mercury resonances create poor local minima around fitness = 2400. The PAGMO team had already noted this in 2010 for the GTOP variant of the problem.
After at least five hours on an AMD 3950X 16-core machine, and often longer, we get:
solution = [2038.4736860390676, 0.6642050463374125, 0.634423807704011, 4049.999077261127, 0.4341666447516408, 451.63078577199883, -3.5111183225692955, 1.6980994436626364, 0.18066690678458555, 224.70119538674018, 1.5700355377783735, 1.1000011238527514, 0.695576590101669, 221.31897889904195, -3.6723464846600904, 1.0500107554228542, 0.7420710057520162, 265.99675406363514, -4.649955978159977, 1.0501275981581841, 0.8290240018499924, 357.96269678842157, 1.5807055451082577, 1.0503396683797719, 0.9029577050599082, 534.1293657177227] fitness = 1962.7194499524037 First Leg: earth to venus Departure: 2005-Jul-31 11:22:06.473775 (2038.4736860390676 mjd2000) Duration: 451.63078577199883days VINF: 4.049999077261127 km/sec DSM after 196.0830229251758 days DSM magnitude: 0.004660230945606301m/s leg no. 2: venus to venus Duration: 224.70119538674018days Fly-by epoch: 2006-Oct-26 02:30:26.364476 (2490.1044718110666 mjd2000) Fly-by radius: 1.6980994436626364 planetary radii DSM after 40.59606992132113 days DSM magnitude: 0.06026467291201602m/s leg no. 3: venus to mercury Duration: 221.31897889904195days Fly-by epoch: 2007-Jun-07 19:20:09.645890 (2714.8056671978065 mjd2000) Fly-by radius: 1.1000011238527514 planetary radii DSM after 153.94430066737883 days DSM magnitude: 602.3225699674765m/s leg no. 4: mercury to mercury Duration: 265.99675406363514days Fly-by epoch: 2008-Jan-15 02:59:29.422768 (2936.1246460968487 mjd2000) Fly-by radius: 1.0500107554228542 planetary radii DSM after 197.38847881477344 days DSM magnitude: 74.95713577573513m/s leg no. 5: mercury to mercury Duration: 357.96269678842157days Fly-by epoch: 2008-Oct-07 02:54:48.973866 (3202.121400160484 mjd2000) Fly-by radius: 1.0501275981581841 planetary radii DSM after 296.75966740455266 days DSM magnitude: 204.77715962193437m/s leg no. 6: mercury to mercury Duration: 534.1293657177227days Fly-by epoch: 2009-Sep-30 02:01:05.976385 (3560.0840969489054 mjd2000) Fly-by radius: 1.0503396683797719 planetary radii DSM after 482.29622627357935 days DSM magnitude: 172.40605634839977m/s Arrival at mercury Arrival epoch: 2011-Mar-18 05:07:23.174397 (4094.2134626666284 mjd2000) Arrival Vinf: 2281.269371304558m/s Insertion DV: 908.1916033350003m/s Total mission time: 5.628308765578537 years (2055.739776627561 days)
The remaining problems are much easier and can be solved in a single loop:
udas=[pk.trajopt.gym.cassini2,
pk.trajopt.gym.eve_mga1dsm,
pk.trajopt.gym.rosetta,
pk.trajopt.gym.juice,
pk.trajopt.gym.em5imp,
pk.trajopt.gym.em7imp]
opt = de_cma(de_max_evals = 750, cma_max_evals = 750)
for uda in udas:
test_problem(opt, single_objective(pg.problem(uda)), num_retries = 10000, value_limit = 100000)For these cases, the first try usually succeeds. Messenger often still needs a second or even a third try.
Juice
Juice models a rendezvous mission to Jupiter as an MGA-1DSM transfer. The selected fly-by sequence, E-EVEME-J, and the other parameters are inspired by the ESA Juice mission. A launcher model is included, namely an Ariane5 launch from Kourou. JUICE, JUpiter ICy moons Explorer, is the first large-class mission in ESA’s Cosmic Vision 2015-2025 program. Planned for launch in 2022 and arrival at Jupiter in 2029, it will spend at least three years making detailed observations of Jupiter and three of its largest moons: Ganymede, Callisto, and Europa.
Juice is a very easy optimization problem. On a 16-core machine it takes about three minutes to surpass -7.99, and about 12 minutes to reach the solution below. In fact, coordinated retry is only needed here if you want a nearly optimal solution. With a fixed 50000 evaluation budget for the DE → CMA-ES sequence, you will usually find a solution better than -7.99 in < 64 retries, which takes about 13 to 26 seconds. But that setup does not find a solution < -7.9994.
This matters when you are screening planet sequences before a final optimization pass. In that phase, regular retry (fcmaes.retry) can be faster. In one hour you can check about 150 sequences. Messenger is a reminder that this does not work equally well for every problem. If many resonances are present, more care is needed. You can also wrap other algorithms from PAGMO for use with fcmaes.retry, but some of them are hard to configure well.
solution = [8108.340864247543, 0.14265429780355834, 0.33964213005840077, 1000.0000080132397, 0.35092639600144443, 436.4274930596665, -4.052861631438546, 3.675356452528689, 0.10741066029466867, 149.0438085981111, -2.1631291025507076, 2.490048280168779, 0.010425812895849209, 315.36502811323714, 4.323135321315694, 1.3233114321446018, 0.03535422310799622, 164.26715796380074, 0.5490453386771064, 1.3808119312886746, 0.23609744554327003, 651.2673270168585, -1.346307881836456, 1.5235742458893275, 0.11717130589660309, 1107.410033856678] fitness = -7.999414363029065 First Leg: earth to earth Departure: 2022-Mar-14 08:10:50.670988 (8108.340864247543 mjd2000) Duration: 436.4274930596665days VINF: 1.0000000080132396 km/sec DSM after 153.15392725537419 days DSM magnitude: 858.1954903986241m/s leg no. 2: earth to venus Duration: 149.0438085981111days Fly-by epoch: 2023-May-24 18:26:26.071343 (8544.76835730721 mjd2000) Fly-by radius: 3.675356452528689 planetary radii DSM after 16.008893894355328 days DSM magnitude: 1.1778007011370638e-06m/s leg no. 3: venus to earth Duration: 315.36502811323714days Fly-by epoch: 2023-Oct-20 19:29:31.134220 (8693.812165905321 mjd2000) Fly-by radius: 2.490048280168779 planetary radii DSM after 3.287936777002836 days DSM magnitude: 3.1286927338502165e-06m/s leg no. 4: earth to mars Duration: 164.26715796380074days Fly-by epoch: 2024-Aug-31 04:15:09.563203 (9009.177194018557 mjd2000) Fly-by radius: 1.3233114321446018 planetary radii DSM after 5.8075377519686695 days DSM magnitude: 4.437463434131599e-06m/s leg no. 5: mars to earth Duration: 651.2673270168585days Fly-by epoch: 2025-Feb-11 10:39:52.011276 (9173.444351982358 mjd2000) Fly-by radius: 1.3808119312886746 planetary radii DSM after 153.7625522744738 days DSM magnitude: 2.233524200245386e-06m/s leg no. 6: earth to jupiter Duration: 1107.410033856678days Fly-by epoch: 2026-Nov-24 17:04:49.065532 (9824.711678999218 mjd2000) Fly-by radius: 1.5235742458893275 planetary radii DSM after 129.7566798299884 days DSM magnitude: 1.6416885164011999e-06m/s Arrival at jupiter Arrival epoch: 2029-Dec-06 02:55:15.990750 (10932.121712855896 mjd2000) Arrival Vinf: 5584.689535619848m/s Insertion DV: 1038.8238384424221m/s Total mission time: 7.7310906190509305 years (2823.7808486083522 days) Initial mass: [6058.99999482] Final mass: [2979.21273893] Declination: -5.8308136937633926e-08
Cassini2
Cassini2 is an MGA-1DSM mission to reach Saturn and be captured by its gravity. The fly-by sequence is Earth, Venus, Venus, Earth, Jupiter, Saturn, the same sequence used by the Cassini spacecraft. A score < 730 is usually reached after about three minutes. Unlike Juice, this problem does need coordinated retry.
solution = [-798.1255840834085, 0.9640536574861399, 0.3842218637483862, 4999.999999998122, 0.054074087669290435, 200.6367375026563, 4.499361877566994, 1.8874938091428237, 0.4971735916709705, 410.34433523595055, -1.9594016768524416, 1.0500000000006968, 0.1669972234371695, 53.27314941868789, 4.727848530339341, 1.3063321404420816, 0.8999999999402388, 588.1308355706944, 4.769745501055804, 70.0, 0.015517225867226283, 2199.999999999682] fitness = 729.0073774509763 First Leg: earth to venus Departure: 1997-Oct-24 20:59:09.535193 (-798.1255840834085 mjd2000) Duration: 200.6367375026563days VINF: 4.9999999999981215 km/sec DSM after 10.849248533399049 days DSM magnitude: 4.038300502236277e-10m/s leg no. 2: venus to venus Duration: 410.34433523595055days Fly-by epoch: 1998-May-13 12:16:03.655423 (-597.4888465807522 mjd2000) Fly-by radius: 1.8874938091428237 planetary radii DSM after 204.01236697109434 days DSM magnitude: 238.54041795196787m/s leg no. 3: venus to earth Duration: 53.27314941868789days Fly-by epoch: 1999-Jun-27 20:31:54.219809 (-187.14451134480169 mjd2000) Fly-by radius: 1.0500000000006968 planetary radii DSM after 8.896468036674339 days DSM magnitude: 8.213846331383754e-10m/s leg no. 4: earth to jupiter Duration: 588.1308355706944days Fly-by epoch: 1999-Aug-20 03:05:14.329584 (-133.87136192611376 mjd2000) Fly-by radius: 1.3063321404420816 planetary radii DSM after 529.3177519784776 days DSM magnitude: 16.41400147219909m/s leg no. 5: jupiter to saturn Duration: 2199.999999999682days Fly-by epoch: 2001-Mar-30 06:13:38.522892 (454.25947364458057 mjd2000) Fly-by radius: 70.0 planetary radii DSM after 34.13789690789289 days DSM magnitude: 5.749616809027544e-09m/s Arrival at saturn Arrival epoch: 2007-Apr-08 06:13:38.522864 (2654.2594736442625 mjd2000) Arrival Vinf: 4260.809617220865m/s Insertion DV: 474.0529580198345m/s Total mission time: 9.45211514778281 years (3452.385057727671 days)
Earth-Venus-Earth mga-1dsm
EVE is an MGA-1DSM mission to Venus and back to Earth. A good score is usually found in a few seconds.
solution = [946.5833262188481, 0.6183687175025218, 0.3663093221347987, 2499.999999999999, 0.14694954450113987, 123.97290946310139, -0.805792681834026, 1.1, 0.3272449823968612, 397.2667888747656] fitness = 4082.4897306048315 First Leg: earth to venus Departure: 2002-Aug-04 13:59:59.385308 (946.5833262188481 mjd2000) Duration: 123.97290946310139days VINF: 2.499999999999999 km/sec DSM after 18.2177625760838 days DSM magnitude: 302.2101312485047m/s leg no. 2: venus to earth Duration: 397.2667888747656days Fly-by epoch: 2002-Dec-06 13:20:58.762920 (1070.5562356819496 mjd2000) Fly-by radius: 1.1 planetary radii DSM after 130.00356333218025 days DSM magnitude: 2639.993724949317m/s Arrival at earth Arrival epoch: 2004-Jan-07 19:45:09.321700 (1467.8230245567152 mjd2000) Arrival Vinf: 1140.2858744070097m/s Total mission time: 1.4270765183788283 years (521.239698337867 days)
Rosetta
Rosetta is a rendezvous mission to comet 67P/Churyumov-Gerasimenko modeled as an MGA-1DSM transfer. The selected fly-by sequence, E-EMEE-C, is similar to the one planned for the Rosetta spacecraft. The objective is total mission delta V. No launcher model is used, and the final rendezvous with the comet is included in the delta V calculation.
A score < 1350 is usually reached after three to five minutes. It is possible to get such a score without coordinated retry, but it is not easy.
solution = [1543.2661408844722, 0.4934153600800337, 0.573351934777993, 4473.451666918325, 0.3729421882864208, 365.2561380887397, -1.6870315844374006, 3.0198930373784543, 0.8096906876626856, 707.9394006716375, 1.7884664717587317, 1.0499999999999998, 0.15802481556983658, 257.29067507087234, -1.7503278875941164, 3.1330604892167573, 0.643510171783615, 730.5078897991023, 4.332094684590686, 1.1005638208186208, 0.4358604536498662, 1849.999962373143] fitness = 1346.8640572197792 First Leg: earth to earth Departure: 2004-Mar-23 06:23:14.572418 (1543.2661408844722 mjd2000) Duration: 365.2561380887397days VINF: 4.473451666918325 km/sec DSM after 136.21942342386168 days DSM magnitude: 2.9724564566919556e-05m/s leg no. 2: earth to mars Duration: 707.9394006716375days Fly-by epoch: 2005-Mar-23 12:32:04.903286 (1908.522278973212 mjd2000) Fly-by radius: 3.0198930373784543 planetary radii DSM after 573.2119401533276 days DSM magnitude: 79.03578114684694m/s leg no. 3: mars to earth Duration: 257.29067507087234days Fly-by epoch: 2007-Mar-01 11:04:49.121315 (2616.461679644849 mjd2000) Fly-by radius: 1.0499999999999998 planetary radii DSM after 40.65831147591335 days DSM magnitude: 5.485266854903277e-05m/s leg no. 4: earth to earth Duration: 730.5078897991023days Fly-by epoch: 2007-Nov-13 18:03:23.447438 (2873.7523547157216 mjd2000) Fly-by radius: 3.1330604892167573 planetary radii DSM after 470.08925765390643 days DSM magnitude: 0.13130711687587276m/s leg no. 5: earth to Churyumov-Gerasimenko Duration: 1849.999962373143days Fly-by epoch: 2009-Nov-13 06:14:45.126081 (3604.260244514824 mjd2000) Fly-by radius: 1.1005638208186208 planetary radii DSM after 806.3418228521936 days DSM magnitude: 793.2628173165956m/s Arrival at Churyumov-Gerasimenko Arrival epoch: 2014-Dec-07 06:14:41.875120 (5454.260206887967 mjd2000) Arrival Vinf: 474.4340670622277m/s Total mission time: 10.707718182076645 years (3910.9940660034945 days)
|
Note
|
The original rosetta.png result figure is not present in this repository.
|
Earth-Mars 5 impulses
Earth-Mars 5 impulses is a 5-impulse deep-space maneuver mission to reach Mars. A good score is reached after a few seconds.
solution = [10359.962016407258, 700.0, 0.792397890631979, 1.0, 0.8009792273038916, 0.0, 0.386976567420959, 0.9334318605507639, 0.5661862201543786, 2433.127502980495, 0.0010315860640127835, 0.8853736291405103, 0.36888093872468103, 0.0, 0.0349218381975402] fitness = 5628.592098319032 Total DV (m/s): 5628.592098319032 Dvs (m/s): [0.0, 2433.127502980495, 0.0, 419.1083963569563, 2776.3561989815807] Tofs (days): [14.271327937417231, 58.22758851188319, 421.7557808109648, 205.74530273973488]
Earth-Mars 7 impulses
Earth-Mars 7 impulses is a 7-impulse deep-space maneuver mission to reach Mars. A good score is reached after a few seconds.
solution = [10359.96201621306, 699.99999999967, 0.8675858244378735, 0.980238773107442, 0.8654116976301882, 0.0, 0.5819397340477784, 0.9334318599048091, 0.5661862216940652, 2433.1275174375482, 0.969600805658548, 1.0, 0.9986677954486534, 0.0, 0.016922153243993637, 0.20270767737395895, 0.051295766859595894, 0.0, 0.8817674494223253, 0.9466891232243391, 0.490939062775394, 0.0, 0.12902376918710104] fitness = 5628.592098319043 Total DV (m/s): 5628.592098319043 Dvs (m/s): [0.0, 2433.1275174375482, 0.0, 0.0, 0.0, 419.1083950393422, 2776.3561858421526] Tofs (days): [14.271327774009604, 54.39513853881697, 3.1016979057927903, 409.8442782317333, 12.642260060314266, 205.74529748900304]
Multi Objective Optimization
We can also study a multi-objective variant of Juice. Our two objectives are:
-
Save propellant by increasing the final mass.
-
Reduce the total mission time.
Pygmo offers several algorithms for MO (multi-objective) optimization:
uda = pk.trajopt.gym.juice_mo
algo1 = pg.algorithm(pg.moead(generations))
algo2 = pg.algorithm(pg.nsga2(generations))
algo3 = pg.algorithm(pg.maco(generations))In practice, two issues appear immediately:
-
With the default configuration, all three algorithms have trouble finding good solutions.
-
It is not clear how to configure Pygmo’s parallel execution mechanism,
archipelago, well for this problem.
One possible setup is:
I = mp.cpu_count()
M = 1
P = 300
G = 100
algo = pg.algorithm(pg.nsga2(gen = G))
topo = pg.topology(pg.fully_connected(w=1.0))
rp = pg.r_policy(pg.fair_replace(M))
sp = pg.s_policy(pg.select_best(M))
archi = pg.archipelago(n=I, t=topo, algo = algo, prob = prob, pop_size=P, r_pol=rp, s_pol=sp)This fully uses 16 cores.
Another option is the weighted sum approach combined with fcmaes coordinated retry. Instead of using one fixed weight, we sweep over different weights and generate good solutions along the Pareto front. That gives us a practical benchmark for the dedicated MOO algorithms.
class two_objectives:
def __init__(self, pagmo_prob, weight):
self.pagmo_prob = pagmo_prob
self.name = pagmo_prob.get_name()
self.fun = self.fitness
lb, ub = pagmo_prob.get_bounds()
self.bounds = Bounds(lb, ub)
self.weight = weight
def fitness(self,X):
try:
val = self.pagmo_prob.fitness(X)
except Exception as ex:
return 1E12
return val[0] + self.weight*val[1]
def optim(weight):
uda = pk.trajopt.gym.juice_mo
pr = two_objectives(pg.problem(uda), weight)
logger().info(pr.name + ' weight = ' + str(weight))
advretry.minimize(pr.fun, pr.bounds, 10000, 10000, logger())
weight = 0.0001
while True:
optim(weight)
weight *= 2We weight the second objective relative to the first, starting with 0.0001, which strongly favors fuel or delta V. We then increase this weight exponentially in each iteration. The expected trend is simple: total mission time should go down while fuel or delta V gets worse. Execution time exceeds 10 hours on a 3950X 16-core CPU, but we do get reasonable results:
Juice (Trajectory Optimization Gym P13-14) weight = 0.0001 solution = [8110.116313009699, 0.15222656981841226, 0.33303598469012846, 1000.0000437885246, 0.3491143875020217, 0.4871833431602956, 2.229331840587739, 3.864027554505235, 0.383418611345456, 0.7835550811389311, -2.1690962755795145, 2.439639530064197, 0.074686310219425, 0.5940304368448553, 4.343738687598682, 1.305726623141431, 0.03188195412256487, 0.7650686059331681, 0.6255555639524654, 1.3305993533677483, 0.16379842866823457, 0.3407535577156361, -1.3951232209828814, 1.5679804383643228, 0.1672058375187888, 0.17031931310947931, 2789.5767102242344] fitness = (-7.997742270040018, 2789.576710224235) Total mission time: 7.637444791852799 years (2789.576710224235 days) Initial mass: [6058.99997739] Final mass: [2974.23538066] Juice (Trajectory Optimization Gym P13-14) weight = 0.0002 solution = [8178.893840594586, 0.6696744638172129, 0.6431354558268282, 2999.9991890745705, 0.24142402261301926, 0.48525432886155345, 1.0145512727192236, 1.2437380950571426, 0.2807995038706594, 0.7506038363731152, -2.157231592881495, 2.481065696638639, 0.0521871614241405, 0.5389611083471264, -1.9273581029317008, 1.3001669088367076, 0.06885787655660897, 0.7288991943110973, 0.6603320679485223, 1.3073296031259627, 0.14343312803516828, 0.27821737536586, -1.4117361081151207, 1.5799078891999832, 0.19077016051979012, 0.12523356744789776, 2709.2744355278537] fitness = (-7.987006391921133, 2709.2744355278537) Total mission time: 7.417589145866813 years (2709.2744355278537 days) Initial mass: [4731.79598124] Final mass: [2942.47514418] Juice (Trajectory Optimization Gym P13-14) weight = 0.0004 solution = [8179.443508482556, 0.6729509773841647, 0.6370156546367186, 3000.0130505441452, 0.24319652546458614, 0.4166617039422654, -5.280573254426214, 1.1324645750396043, 0.27462876304589007, 0.7092540727184542, -2.1950336354386852, 2.255189952426106, 0.03978925494280019, 0.4719587719252051, -1.8905822704179673, 1.2453108381383247, 0.12349063439035901, 0.6893539247353409, -5.469525789090203, 1.1543240184550114, 0.14142946154723224, 0.21140938606328202, -1.4812407219640134, 1.60655474357165, 0.04017554437749502, 0.09500591477807817, 2645.951340743001] fitness = (-7.968102442031603, 2645.9513407430013) Total mission time: 7.244219960966465 years (2645.9513407430013 days) Initial mass: [4728.98866135] Final mass: [2887.37320458] Juice (Trajectory Optimization Gym P13-14) weight = 0.0008 solution = [8186.179756871906, 0.20165084375795872, 0.26905080330409437, 3374.4146907104346, 0.23254985500140193, 0.47241555021382825, 2.0588168111700402, 6.739180591235458, 0.34794739269432007, 0.7432500909232481, -2.216423843954621, 2.1215147988023118, 0.05791330560377872, 0.5192854942075542, 4.414357837606136, 1.2063922900097728, 0.08125504813979208, 0.727712205494026, 0.8994676611088803, 1.050000040876132, 0.04285089314078716, 0.2577432506982371, -1.5142603957607905, 1.6008638622803881, 0.010523390650404502, 0.1400250711622686, 2601.7309912552946] fitness = (-7.946060835277903, 2601.730991255295) Total mission time: 7.123151242314292 years (2601.730991255295 days) Initial mass: [4397.26851686] Final mass: [2824.42712464] Juice (Trajectory Optimization Gym P13-14) weight = 0.0016 solution = [8187.696717353289, 0.20580116800841633, 0.26787059729909224, 3386.7416562244025, 0.22486520908302862, 0.39606425959350755, -4.253041811518965, 7.057814063074282, 0.4184571817165126, 0.6949302400904048, -2.207331400664618, 2.1533722260999126, 0.034130744023674366, 0.4452550938841818, -1.862656949932257, 1.2019648207938172, 0.06944571771256312, 0.6761285408299662, -5.3301227060661995, 1.0500000296263783, 0.8210782701673592, 0.18969985318261587, -1.5299955693533245, 1.6219442344127535, 0.023733219454938507, 0.13247032212681625, 2434.972454053947] fitness = (-7.751229257711435, 2434.972454053947) Total mission time: 6.666591249976583 years (2434.972454053947 days) Initial mass: [4386.34687932] Final mass: [2324.42798016] Juice (Trajectory Optimization Gym P13-14) weight = 0.0032 solution = [8378.693422586422, 0.9568081499212745, 0.0002686084118801893, 1999.9637038181368, 0.9379874262019485, 0.7041736922795353, 2.0292340441339287, 1.0500938586187636, 0.010048843483523123, 0.7511935654050826, 4.084548346708438, 2.1880766631468105, 0.04077348235789401, 0.5281892598024404, -1.8588571155044673, 1.200640632936553, 0.18289180035424735, 0.7345851701325535, 0.9846960017622043, 1.050230188317333, 0.8179188136864259, 0.2710583084636556, 4.741032108908588, 1.4637849008248247, 0.05503111350520763, 0.26346347776767015, 2109.962009208408] fitness = (-6.881259529430607, 2109.962009208408) Total mission time: 5.7767611477300695 years (2109.962009208408 days) Initial mass: [4044.04065172] Final mass: [973.85218314] Juice (Trajectory Optimization Gym P13-14) weight = 0.0064 solution = [8389.147244723335, 0.9201833254461056, 0.0002487612898780255, 1956.0353502490668, 0.8937743422244669, 0.6319118732854974, 2.114563855857773, 1.050030787279265, 0.010327707184346733, 0.6814842972922296, -2.1701389665942945, 2.3156941575389394, 0.012811607963323883, 0.41973744761815157, 4.432272455014974, 1.202716445919663, 0.016121147371226482, 0.6572308915905087, -5.275048116158234, 1.0500008005018495, 0.8227769638519822, 0.16979474039257636, -1.545538848627961, 1.203041913314907, 0.02316964539659546, 0.21177921119661475, 2000.0001268207916] fitness = (-6.412257537995305, 2000.0001268207918) Total mission time: 5.475701921480607 years (2000.0001268207918 days) Initial mass: [4093.24040772] Final mass: [609.26757437] Juice (Trajectory Optimization Gym P13-14) weight = 0.0128 solution = [8389.99246645576, 0.9310872905708278, 0.00023956109624324008, 1923.7145240133386, 0.8928089924310574, 0.6980785467270006, 2.1200442165161912, 1.050095711253613, 0.010022439877003894, 0.7402229683832497, -2.1675321412386, 2.328004999219424, 0.06651221131595277, 0.5056936685141469, -1.8502582115196264, 1.2030019613346181, 0.07384884550225786, 0.7191509814535274, 1.008687327515632, 1.0501070906961574, 0.8231781510283087, 0.24832655454479133, -1.54594248265311, 1.206481690018342, 0.05647684377465877, 0.2948982396357695, 2000.0310567382437] fitness = (-6.4123598926059335, 2000.0310567382435) Total mission time: 5.475786602979448 years (2000.0310567382435 days) Initial mass: [4129.43973311] Final mass: [609.32993891] Juice (Trajectory Optimization Gym P13-14) weight = 0.0256 solution = [8389.357022843442, 0.9226312928490936, 0.00024363029910228227, 1936.6928151536988, 0.8953886025771408, 0.5574895649768024, -4.165868167291542, 1.0509907137007959, 0.011190373445135158, 0.6137744810778374, 4.114054748168721, 2.3202702485380224, 0.012751492110167576, 0.33101487618584285, 4.432632210154407, 1.2026615518213288, 0.028664574239499935, 0.5859585191738412, 1.0094415826009, 1.0500174385146135, 0.8227607958750899, 0.10451845184815436, 4.737130508247658, 1.2036745840173237, 0.013232962157012154, 0.1383962744368863, 2000.003256297061] fitness = (-6.412206486795257, 2000.003256297061) Total mission time: 5.475710489519674 years (2000.003256297061 days) Initial mass: [4114.90404703] Final mass: [609.23647132]
The trend is clear. Below about 2400 days, the first objective degrades sharply.
Now assume management decides that 2400 days is the target mission time. We can adapt the weighted sum by turning the second objective into an equality constraint, (T == 2400), and use weight = 0.01 to prioritize that constraint:
...
def fitness(self,X):
try:
val = self.pagmo_prob.fitness(X)
except Exception as ex:
return 1E12
return val[0] + self.weight*abs(val[1] - 2400)
...
weight = 0.01
optim(weight)After about 80 minutes, we find:
solution = [8184.971833173281, 0.6785120362633554, 0.6399686066961614, 3384.5191657874484, 0.2177204427288376, 0.532791237252402, 0.9930523033160469, 1.0518066161610329, 0.5712277532411305, 0.7784361424802742, -2.2336337605711267, 2.043017304921004, 0.12767721565536183, 0.5762959058584155, 4.414415650075737, 1.1976599227921725, 0.18198049504078612, 0.7662669432203426, 0.9545084252936471, 1.0500426197011596, 0.8245192258310868, 0.3233439096671323, -1.534305327990903, 1.5981800484551762, 0.048586771966370075, 0.2696603273182286, 2400.0001374345097] fitness = (-7.6874721365395065, 2400.00013743451) Total mission time: 6.570842265392225 years (2400.00013743451 days) Initial mass: [4388.31145807] Final mass: [2180.85468555]
This weighted and penalized objective is much harder to optimize than the single-objective Juice problem. It would be interesting to see whether a dedicated MOO algorithm can reach a similar result.
Now add another engineering constraint. Suppose the team says: "We simply can not construct a space ship with such a low dry mass, 2180kg is not enough". Management agrees, and we are asked to target a first objective value of -7.9.
Our previous result with weight = 0.0008 was fitness = (-7.946060835277903, 2601.730991255295). Its dry mass was 2824 kg. Now let us see whether we can beat that dry mass while keeping the first objective value < -7.9.
We adapt fitness again:
...
def fitness(self,X):
try:
val = self.pagmo_prob.fitness(X)
except Exception as ex:
return 1E12
penalty = 10*max(0, val[0] - (-7.9))
return penalty + val[0] + self.weight*val[1]
...
weight = 0.0005
optim(weight)We lower the weight for the second objective and add a penalty whenever the first objective is worse than -7.9.
solution = [8179.473217982026, 0.6738225986673105, 0.6375260743595529, 2999.57061606059, 0.2453457164227346, 0.47258050731851226, 0.9976567917383342, 1.0892254502281873, 0.3784197076676981, 0.7460620377696593, -2.210741160866029, 2.162485225687267, 0.039709813601583197, 0.5253625125458127, 4.405422725768786, 1.221536455250864, 0.06893498354757457, 0.7299434421595803, 0.8667414770477118, 1.0907238132335606, 0.19067718225704694, 0.2640317495079226, 4.781028628270766, 1.6046636333696407, 0.19768474425666516, 0.14033398966827756, 2622.7076941048886] fitness = (-7.957616478970268, 2622.7076941048886) Total mission time: 7.180582324722487 years (2622.7076941048886 days) Initial mass: [4729.35404072] Final mass: [2857.25450358]
The mission is only 222 days longer, but dry mass rises to 2857kg.
First Leg: earth to earth Departure: 2022-May-24 11:21:26.033647 (8179.473217982026 mjd2000) Duration: 371.1683843156681days VINF: 2.99957061606059 km/sec DSM after 91.06457316339647 days DSM magnitude: 199.30118126994356m/s leg no. 2: earth to venus Duration: 145.06423724348917days Fly-by epoch: 2023-May-30 15:23:54.438521 (8550.641602297694 mjd2000) Fly-by radius: 1.0892254502281873 planetary radii DSM after 54.89516625071878 days DSM magnitude: 5.962198622316355e-08m/s leg no. 3: venus to earth Duration: 318.7374482405083days Fly-by epoch: 2023-Oct-22 16:56:24.536358 (8695.705839541182 mjd2000) Fly-by radius: 2.162485225687267 planetary radii DSM after 12.657004657474856 days DSM magnitude: 9.953831856214234e-09m/s leg no. 4: earth to mars Duration: 155.88003364375575days Fly-by epoch: 2024-Sep-05 10:38:20.064338 (9014.44328778169 mjd2000) Fly-by radius: 1.221536455250864 planetary radii DSM after 10.745587554627674 days DSM magnitude: 1.571277166941131e-08m/s leg no. 5: mars to earth Duration: 659.4377741312054days Fly-by epoch: 2025-Feb-08 07:45:34.971159 (9170.323321425447 mjd2000) Fly-by radius: 1.0907238132335606 planetary radii DSM after 125.73973664519721 days DSM magnitude: 5.455438517375732e-08m/s leg no. 6: earth to jupiter Duration: 972.4198165302618days Fly-by epoch: 2026-Nov-29 18:15:58.656095 (9829.761095556652 mjd2000) Fly-by radius: 1.6046636333696407 planetary radii DSM after 192.23256274089806 days DSM magnitude: 7.674593586671564e-08m/s Arrival at jupiter Arrival epoch: 2029-Jul-29 04:20:30.804309 (10802.180912086915 mjd2000) Arrival Vinf: 5668.335880969336m/s Insertion DV: 1067.5528671585835m/s Total mission time: 7.180582324722487 years (2622.7076941048886 days) Initial mass: [4729.35404072] Final mass: [2857.25450358] Declination: 4.9959251623060075
Tandem Planet Sequence
In the discussion of the Tandem problem above, we chose the 6th planet sequence, [earth, venus, earth, earth, saturn], because we suspected it was the best one. Now we can test that assumption. Instead of giving every sequence a full optimization budget up front, we use an iterative elimination approach. Each round applies parallel retry to all remaining sequences for a limited number of retries. We then remove the worst ones and continue with the rest. This concentrates more compute time on the better candidates.
Only the best remaining sequences get coordinated retry in the final stage. To make that work, we keep the evaluation store for each sequence across all iterations.
We define a class test_problem for one Tandem planet sequence together with its retry store. Here we use the Tandem problem without a time constraint, so mission time may exceed 10 years. This variant is slightly harder than the time-constrained one.
import pykep as pk
import numpy as np
from fcmaes import retry, advretry
from fcmaes.optimizer import single_objective, logger, de_cma
class test_problem:
def __init__(self, prob, num_evals, num_retries, id):
self.store = retry.Store(prob.bounds, logger = logger())
self.prob = prob
self.name = prob.name
self.fun = prob.fun
self.num_retries = num_retries
self.optimizer = de_cma(num_evals)
self.retries = 0
self.value = 0
self.id = id
self.ret = None
def retry(self):
retries = self.retries + self.num_retries
self.ret = retry.retry(self.fun, self.store, self.optimizer.minimize, retries)
self.retries = retries
self.value = self.store.get_y_best()
def minimize(self):
self.ret = advretry.minimize(self.fun, bounds=self.prob.bounds, num_retries = 20000,
value_limit = 0, logger = logger())
self.value = self.ret.funMethod retry performs the parallel retry for one iteration. Method minimize performs the coordinated retry at the end for the best sequences. Class test_problems represents the full set of sequences we want to compare.
class test_problems:
def __init__(self):
self.test_problems = []
def add(self, test_problem):
self.test_problems.append(test_problem)
def retry(self):
for tp in self.test_problems:
logger().info("retry for problem " + str(tp.id))
tp.retry()
def minimize(self):
for tp in self.test_problems:
logger().info("coordinated retry for problem " + str(tp.id))
tp.minimize()
def values(self):
return np.array([tp.value for tp in self.test_problems])
def remove_worst(self, n = 1):
for _ in range(n):
idx = self.values().argsort()
self.test_problems.pop(idx[-1])
def size(self):
return len(self.test_problems)
def dump(self):
idx = self.values().argsort()
for i in range(self.size()):
tp = self.test_problems[idx[i]]
logger().info(str(tp.id) + ' ' + str(tp.value))Method test_sequences initializes the test_problems instance with all Tandem variants, calls tps.retry() in a loop, and finally calls tps.minimize() to run coordinated retry for the best Tandem planet sequences.
def test_sequences():
tps = test_problems()
n = 24
delta = 4
for i in range(n):
tandem = pk.trajopt.gym.tandem(i+1)
prob = single_objective(pg.problem(tandem))
tps.add(test_problem(prob, 50000, 128, i+1))
while tps.size() > delta:
tps.retry()
tps.dump()
tps.remove_worst(delta)
tps.minimize()
tps.dump()We choose 128 as the parallel retry count for each increment. Make sure this is a multiple of your number of logical processors, mp.cpu_count(). As optimizer we use de_cma(50000), the default optimizer for parallel retry with an objective evaluation budget of 50000. The result is:
best planet sequences with their retry scores: 6 -7.234279427520673 20 -7.163017778254503 22 -7.0785510740228546 24 -7.042380192439136 2 -7.02035532005052 18 -6.985418830906635 8 -6.981895243731103 19 -6.9770645164161476
Sequence 6, [earth, venus, earth, earth, saturn], the one we used before, wins as expected. The GTOP result is worse, but the problems are not identical. At L116, the GTOP Tandem problem uses the "wrong direction" for Vinfz, and this is corrected in pykep. A different Lambert solver is also used.
It is still an open question how good the final result below really is. But fitness = -7.4338 is at least a strong target for other algorithms on the unconstrained Tandem problem with sequence 6.
solution = [8985.7529255691, 0.6288683258483504, 0.5114559852888231, 3069.6887063429026, 0.9257586246693946, 1508.947519175164, -1.1855127646915973, 1.0557970110155512, 0.8773779112867278, 2354.6044622468585, 4.703963836863605, 1.2164967307344021, 0.9043728224239871, 2481.2229452637835, -1.3794795580387544, 1.3379481557564792, 0.08143963789806034, 2459.0975767600025] fitness = (-7.4338140059429945, 5151.372503445809) First Leg: earth to venus Departure: 2024-Aug-07 18:04:12.769170 (8985.7529255691 mjd2000) Duration: 1508.947519175164days VINF: 3.0696887063429026 km/sec DSM after 1396.9211800498947 days DSM magnitude: 45.55784305233069m/s leg no. 2: venus to earth Duration: 2354.6044622468585days Fly-by epoch: 2028-Sep-24 16:48:38.425904 (10494.700444744265 mjd2000) Fly-by radius: 1.0557970110155512 planetary radii DSM after 2065.877944992558 days DSM magnitude: 0.0010049610627437088m/s leg no. 3: earth to earth Duration: 2481.2229452637835days Fly-by epoch: 2035-Mar-07 07:19:03.964033 (12849.304906991123 mjd2000) Fly-by radius: 1.2164967307344021 planetary radii DSM after 2243.9505980713657 days DSM magnitude: 0.007327435540476821m/s leg no. 4: earth to saturn Duration: 2459.0975767600025days Fly-by epoch: 2041-Dec-21 12:40:06.434824 (15330.527852254907 mjd2000) Fly-by radius: 1.3379481557564792 planetary radii DSM after 200.26801620733227 days DSM magnitude: 0.010057156024015437m/s Arrival at saturn Arrival epoch: 2048-Sep-14 15:00:37.066888 (17789.625429014908 mjd2000) Arrival Vinf: 5225.071650761917m/s Insertion DV: 554.072387224358m/s Total mission time: 24.103689263369773 years (8803.87250344581 days) Initial mass: [2172.69961397] Final mass: [1692.24953049] Declination: 15.480672242471464
|
Note
|
The original tandem2.png result figure is not present in this repository.
|