Multiple impulses: planet to planet transfer

In this tutorial, we use the pykep.trajopt.pl2pl_N_impulses to study a single direct transfer between planets.

The encoding is original with pykep and makes use of a decision vector (chromosome) in the form:

\[ \mathbf x = [t_0,T] + \underbrace{[\alpha_0,u_0,v_0,\Delta V_0] + [\alpha_1,u_1,v_1,\Delta V_1] + ...}_{\text{propagations}} + \underbrace{[\alpha_{N-1}]}_{\text{Lambert}} + \underbrace{([t_f])}_{\text{Optional phase term}} \]

to define an \(N\) impulse trajectory starting at \(t_0\) at the departure planet with a \(\Delta V_0\) and arriving at the target planet with a \(\Delta V_{N-1}\) and applying the remaining \(\Delta V_i\) at epochs encoded using the \(\alpha\)-encoding for times (see Multiple Gravity Assist (MGA)).

This optimization problem is a nice tool to visualize and study the famous “How many impulses?” question which has a very complicated answer in general.

We start, as often, with some fundamental imports:

import pykep as pk
import numpy as np
import pygmo as pg

We define the problem as an Earth-Venus transfer. We set the transfer time as fixed as well as the departure epoch.

start = pk.planet(pk.udpla.jpl_lp("earth"))
target = pk.planet(pk.udpla.jpl_lp("venus"))

tof_bounds = [360.0, 360.0]
DV_max_bounds = [0.0, 4]
phase_free = False
multi_objective = False
t0_bounds = [140, 140.]

udp3 = pk.trajopt.pl2pl_N_impulses(
    start=start,
    target=target,
    N_max=3,
    tof_bounds=tof_bounds,
    phase_free=phase_free,
    multi_objective=multi_objective,
    t0_bounds=t0_bounds,
    DV_max_bounds=DV_max_bounds
)

def solve(udp, N=20):
    # We use CMA-ES 
    uda = pg.cmaes(4500, force_bounds=True, sigma0=0.5, ftol=1e-4)
    # But if you prefer a self adaptive version of differential evolution thats also OK.
    #uda = pg.sade(2500, ftol=1e-4, xtol=1e-4)
    algo = pg.algorithm(uda)

    print("Multi-start:")

    res = list()
    for i in range(N):
        pop = pg.population(udp, 20)
        pop = algo.evolve(pop)
        res.append([pop.champion_f, pop.champion_x])
        print(i, pop.champion_f[0], end= '\r')
        

    best_x = sorted(res, key =  lambda x: x[0][0])[0][1]

    print(f"\nThe best solution found has a DV of {udp.fitness(best_x)[0]/1000:.5e} km/s")
    return best_x

best_x_3 = solve(udp3)

Multi-start:
0 6015.33207581
1 5986.38946014
2 5986.3894582
3 5986.38946074
4 6015.3320758
5 6015.33207581
6 6291.37907204

7 5986.38946093
8 5986.38945962
9 5986.3894605
10 5986.38945796
11 5986.3894637
12 5986.38945986
13 6015.3320758
14 5986.38945803

15 5986.38945895
16 6015.33207584
17 5986.38945803
18 6015.33207582
19 6015.33207581
The best solution found has a DV of 5.98639e+00 km/s

We can easily visualize the solution:

def visualise(udp, best_x):
    import matplotlib.pyplot as plt
    import numpy as np

    # Create the figure and define the grid
    fig = plt.figure(figsize=(10, 10))
    gs = fig.add_gridspec(
        2, 2, width_ratios=[2, 1], height_ratios=[1, 1], wspace=0.01, hspace=0.01
    )

    # Top view (spans rows and columns)
    ax1 = fig.add_subplot(gs[:, 0], projection="3d")
    udp.plot(best_x, ax=ax1)
    ax1.view_init(90, 0)
    ax1.axis("off")

    # Ecliptic view 1
    ax1 = fig.add_subplot(gs[0, 1], projection="3d")
    udp.plot(best_x, ax=ax1)
    ax1.view_init(0, 0)
    ax1.axis("off")

    # Ecliptic view 2
    ax1 = fig.add_subplot(gs[1, 1], projection="3d")
    udp.plot(best_x, ax=ax1)
    ax1.view_init(0, 90)
    ax1.axis("off")
    return fig


visualise(udp3, best_x_3);

… and quickly investigate a human readable format for the solution found.

udp3.pretty(best_x_3)

Total DV (m/s):  5986.389457960815
Dvs (m/s):  [219.5038165281643, 2899.3602556143114, 2867.525385818339]
Total DT (days):  360.0
Tofs (days):  [205.94310275292256, 154.05689724707744]

Are more impulses needed?

In the previous design we allowed maximum 3 impulses, but would more impulse be needed? We can quickly answer this question for this specific case by looking at the primer vector.

ax, (grid, p) = udp3.plot_primer_vector(best_x_3)

… the primer vector goes above one this trajevtory IS NOT OPTIMAL. We need more impulses. Lets add one.

udp4 = pk.trajopt.pl2pl_N_impulses(
    start=start,
    target=target,
    N_max=4,
    tof_bounds=tof_bounds,
    phase_free=phase_free,
    multi_objective=multi_objective,
    t0_bounds=t0_bounds,
    DV_max_bounds=DV_max_bounds
)

best_x_4 = solve(udp4, N=20)

Multi-start:

0 5986.38945626

1 5986.44812541

2 5986.44295199
3 5986.38945625

4 6015.3865159
5 6015.332076

6 5986.44299703
7 5986.38945626

8 5937.92739105
9 6015.33207595
10 5986.38945626

11 5986.38945625

12 5986.44298194
13 5937.92738758

14 6312.58212984
15 6015.3320759
16 6015.33207589

17 5980.12401477
18 6015.33207592

19 5986.38945626
The best solution found has a DV of 5.93793e+00 km/s

visualise(udp4, best_x_4);

udp4.pretty(best_x_4)

Total DV (m/s):  5937.927387579902
Dvs (m/s):  [197.6814973083379, 2732.318921327866, 231.0582773627758, 2776.8686915809217]
Total DT (days):  360.0
Tofs (days):  [209.6654151022237, 102.53307439075782, 47.80151050701844]

We found a better solution! Four impulses are better than 4 in this case! This is quite rare and interesting as a case and its here found as the optimization problem was defined to have a fixed duration and starting date. Typically, when freeing these degrees of freedom, three impulses are enough.

Try to find other interesting cases where the number of impulses matter, and do not be fooled by artefacts. Some times impulses with zero \(\Delta V\) are added to move around the starting date outside of the bounds. Whats the maximum number of impulses you can find to be optimal?

We conclude showing the primer vector for this case:

ax, (grid, p) = udp4.plot_primer_vector(best_x_4)

The very same trajectory is also studied in the notebook dedicated to the primer vector where the whole theory is also developed and coded as to show its inner workings.