Zero Order Hold: planet to planet low thrust transfer

Zero Order Hold: planet to planet low thrust transfer#

In this tutorial we demonstrate the use of the zoh_pl2pl class to find a low-thrust trajectory connecting two moving planets using the Zero-Order Hold (ZOH) direct method with free departure and arrival velocities.

Key Features#

The zoh_pl2pl class offers several features that make it a powerful tool for low-thrust trajectory optimization:

  1. Free departure and arrival velocities - The spacecraft can depart and arrive with excess velocity relative to the planets

  2. Non-dimensional formulation - Works internally with non-dimensional units for better numerical conditioning

  3. Taylor-adaptive integration - Uses heyoka as Taylor propagator for accurate trajectory propagation

  4. Analytical gradients - Optional analytical gradient computation for faster optimization

  5. Flexible time encoding - Supports both uniform and softmax time grid encodings

Decision Vector#

The decision vector for this class is:

\[ \mathbf{x} = [t_0, m_f, v^{\infty}_{\text{dep}}, \hat{i}_{\text{dep},x}, \hat{i}_{\text{dep},y}, \hat{i}_{\text{dep},z}, v^{\infty}_{\text{arr}}, \hat{i}_{\text{arr},x}, \hat{i}_{\text{arr},y}, \hat{i}_{\text{arr},z}, \text{controls}, T_{\text{tof}}] (+ [\mathbf{w}] \text{ for softmax}) \]

Where:

  • \(t_0\) is the departure epoch (MJD2000, in days)

  • \(m_f\) is the final spacecraft mass (non-dimensional)

  • \(v^{\infty}_{\text{dep/arr}}\) are the magnitudes of departure/arrival excess velocities (non-dimensional)

  • \(\hat{i}_{\text{dep/arr}}\) are unit direction vectors for the excess velocities

  • controls = \([T, \hat{i}_x, \hat{i}_y, \hat{i}_z] \times n_{\text{seg}}\) where \(T\) is thrust magnitude and \(\hat{i}\) is thrust direction

  • \(T_{\text{tof}}\) is the time of flight (non-dimensional)

  • \(\mathbf{w}\) are softmax weights (only for softmax encoding)

Note

This notebook can use the commercial solver SNOPT 7 or the open-source IPOPT solver. Both are demonstrated below.

Basic Imports#

import pykep as pk
import numpy as np
import time
import pygmo as pg
try:
    import pygmo_plugins_nonfree as ppnf
    has_snopt = True
except:
    has_snopt = False
    print("SNOPT not available, will use IPOPT instead")

from matplotlib import pyplot as plt

Problem Definition#

We define the problem parameters including planets, spacecraft characteristics, and non-dimensional scaling factors.

# Problem data
mu = pk.MU_SUN
max_thrust = 0.6
isp = 3000
veff = isp * pk.G0
# Initial mass
ms = 1500.0
# Number of segments
nseg = 30

# Non dimensional units
L = pk.AU
MU = mu  # (central body mu must be 1 in these units as a requirement of the ZOH integrator used)
TIME = np.sqrt(L**3 / MU)
V = L / TIME
ACC = V / TIME
MASS = ms
F = MASS * ACC
ms_nd = ms / MASS

# Source and destination planets
earth = pk.planet_to_keplerian(
    pk.planet(pk.udpla.jpl_lp(body="EARTH")), 
    when=pk.epoch(5000),
)
mars = pk.planet_to_keplerian(
    pk.planet(pk.udpla.jpl_lp(body="MARS")), 
    when=pk.epoch(5000),
)

# Low tolerances result in higher speed (the needed tolerance depends on the orbital regime)
tol = 1e-10
tol_var = 1e-6

# We instantiate ZOH Taylor integrators for Keplerian dynamics.
ta = pk.ta.get_zoh_kep(tol)
ta_var = pk.ta.get_zoh_kep_var(tol_var)

# We set the Taylor integrator parameters
veff_nd = veff / V

# Creating Taylor-Adaptive Integrators
# The ZOH method requires Taylor-adaptive integrators for high-precision trajectory propagation. We create two integrators:
# 1. Nominal dynamics integrator (for fitness evaluation)
# 2. Variational dynamics integrator (for gradient computation)
ta.pars[4] = 1.0 / veff_nd
ta_var.pars[4] = 1.0 / veff_nd
print(f"Nominal integrator state dimension: {len(ta.state)}")
print(f"Variational integrator state dimension: {len(ta_var.state)}")

Nominal integrator state dimension: 7
Variational integrator state dimension: 84

udp_uniform = pk.trajopt.zoh_pl2pl(
    pls=earth,
    plf=mars,
    ms=ms_nd,
    nseg=nseg,
    cut=0.6,
    t0_bounds=[7360, 8300.0],
    tof_bounds=[100.*pk.DAY2SEC/TIME, 350.*pk.DAY2SEC/TIME],
    mf_bounds=[2/3,1.],
    vinf_dep_bounds = [0., 0.],
    vinf_arr_bounds = [0., 0.],
    tas = (ta,None),
    max_thrust=max_thrust/F,
    # w_bounds_softmax=[-1., 1.],
    time_encoding='uniform',
)

udp_g_uniform = pk.trajopt.zoh_pl2pl(
    pls=earth,
    plf=mars,
    ms=ms_nd,
    nseg=nseg,
    cut=0.6,
    t0_bounds=[7360, 8300.0],
    tof_bounds=[100.*pk.DAY2SEC/TIME, 350.*pk.DAY2SEC/TIME],
    mf_bounds=[2/3,1.],
    vinf_dep_bounds = [0., 0.],
    vinf_arr_bounds = [0., 0.],
    tas = (ta,ta_var),
    max_thrust=max_thrust/F,
    # w_bounds_softmax=[-1., 1.],
    time_encoding='uniform',
    )

print("UDP instances created successfully")
print(f"With gradient: {udp_g_uniform.has_gradient()}")
print(f"Without gradient: {udp_uniform.has_gradient()}")

UDP instances created successfully
With gradient: True
Without gradient: False

Part 1: Gradient vs No-Gradient Comparison#

We first compare the performance of the UDP with and without analytical gradients. Both use uniform time encoding.

Gradient Performance Analysis#

Let’s compare the performance of analytical gradients versus numerical gradients. The analytical gradients use variational equations integrated alongside the dynamics, providing exact derivatives.

# Create a problem and population to test gradient
prob_g = pg.problem(udp_g_uniform)
pop_g = pg.population(prob_g, 1)
x_test = pop_g.champion_x

print("\n=== Gradient Performance Comparison ===")
print(f"\nDecision vector dimension: {len(x_test)}")
print(f"Fitness dimension: {len(prob_g.fitness(x_test))}")

print("\nAnalytical gradient (from variational equations + chain rule):")
%timeit udp_g_uniform.gradient(x_test)

print("\nNumerical gradient (simple finite difference):")
%timeit pg.estimate_gradient(udp_g_uniform.fitness, x_test)

print("\n✓ The analytical gradient is significantly faster!")

=== Gradient Performance Comparison ===

Decision vector dimension: 131
Fitness dimension: 40

Analytical gradient (from variational equations + chain rule):
310 μs ± 1.33 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Numerical gradient (simple finite difference):
12.2 ms ± 279 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

✓ The analytical gradient is significantly faster!

Solving with and without Gradients#

Now we solve the optimization problem using both approaches to compare convergence.

# Setup optimization algorithm
if has_snopt:
    # Use SNOPT if available
    library_snopt72 = "/Users/dario.izzo/opt/libsnopt7_c.dylib"
    uda = ppnf.snopt7(library=library_snopt72, minor_version=2, screen_output=False)
    uda.set_integer_option("Major iterations limit", 1000)
    uda.set_integer_option("Iterations limit", 10000)
    uda.set_numeric_option("Major optimality tolerance", 1e-3)
    uda.set_numeric_option("Major feasibility tolerance", 1e-9)
    algo_name = "SNOPT7"
else:
    # Use IPOPT as fallback
    uda = pg.ipopt()
    uda.set_numeric_option("tol", 1e-8)
    uda.set_numeric_option("constr_viol_tol", 1e-8)
    uda.set_integer_option("max_iter", 2000)
    algo_name = "IPOPT"

algo = pg.algorithm(uda)
print(f"Using {algo_name} optimizer")

Using SNOPT7 optimizer

Optimization with Gradients#

# Solve with gradients
prob_g_uniform = pg.problem(udp_g_uniform)
prob_g_uniform.c_tol = 1e-6

print("\n=== Optimization WITH Gradients (Uniform) ===")

masses = []
xs=[]
for i in range(10):
    pop_g = pg.population(prob_g_uniform, 1)
    pop_g = algo.evolve(pop_g)
    if(prob_g_uniform.feasibility_f(pop_g.champion_f)):
        print(".", end="")
        masses.append(pop_g.champion_x[1])
        xs.append(pop_g.champion_x)
    else:
        print("x", end ="")
print("\nBest mass is: ", np.max(masses)*MASS)
print("Worst mass is: ", np.min(masses)*MASS)
best_idx = np.argmax(masses)

=== Optimization WITH Gradients (Uniform) ===
..........
Best mass is:  1233.36286109
Worst mass is:  1207.65469155

Optimization without Gradients#

prob_uniform = pg.problem(udp_uniform)
prob_uniform.c_tol = 1e-6
print("\n=== Optimization WITHOUT Gradients (Uniform) ===")
masses_no_grad = []
for i in range(10):
    pop = pg.population(prob_uniform, 1)
    pop = algo.evolve(pop)
    if(prob_uniform.feasibility_f(pop.champion_f)):
        print(".", end="")
        masses_no_grad.append(pop.champion_x[1])
    else:
        print("x", end ="")
print("\nBest mass is: ", np.max(masses_no_grad)*MASS)
print("Worst mass is: ", np.min(masses_no_grad)*MASS)

=== Optimization WITHOUT Gradients (Uniform) ===
.x..xx.xx.
Best mass is:  1229.03971614
Worst mass is:  1178.63457089

if len(xs) > 0:
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')  
    t0 = udp_g_uniform.compute_t0(xs[best_idx])
    ax=pk.plot.add_planet_orbit(ax=ax, pla=earth, label='earth', units=1., color='green')
    ax=pk.plot.add_planet_orbit(ax=ax, pla=mars, label='mars', units=1., color='orange')
    ax = udp_g_uniform.plot(xs[best_idx], ax=ax, N=20, c='k', s=5)
    ax=pk.plot.add_planet(ax=ax, when=t0, pla=earth, units=1., color='darkgreen', s=20)
    ax=pk.plot.add_planet(ax=ax, when=t0+xs[best_idx][10+4*udp_g_uniform.nseg]*udp_g_uniform.TIME/pk.DAY2SEC, pla=mars, units=1., color='darkorange', s=20)
    ax.set_title("Best ZOH Transfer Trajectory (with gradients, uniform)")
    ax.legend()
../_images/18cd3081869308df9fd1a3ab6014df54f150aeb529b3c8d9d9e7c8490fc0e40a.png 
udp_g_uniform.plot_throttle(xs[best_idx])

<Axes: xlabel='time grid', ylabel='throttle value (nd)'>
../_images/1602d7c5cb249d8d4a9883b09918597f7e40fb9fb21134cbe89e813eb7ee3c3d.png

Part 2: Time Encoding Comparison (Uniform vs Softmax)#

Now we compare the two time encoding methods: uniform and softmax.

Time Encoding Explanation#

Uniform Encoding: Time grid nodes are evenly spaced from 0 to TOF: $\(t_i = \frac{i \cdot T_{\text{tof}}}{n_{\text{seg}}}, \quad i = 0, 1, ..., n_{\text{seg}}\)$

Softmax Encoding: Time grid nodes are determined by learnable weights \(\mathbf{w}\): $\(\Delta t_i = T_{\text{tof}} \cdot \text{softmax}(\mathbf{w})_i, \quad t_i = \sum_{j=0}^{i-1} \Delta t_j\)$

The softmax encoding adds \(n_{\text{seg}}\) additional decision variables but allows the optimizer to adaptively distribute time where it’s most needed. For a better discussion on the softmax encoding see also the tutorial: udp_zoh_point2point.ipynb.

Create UDP with Softmax Encoding#

udp_g_softmax = pk.trajopt.zoh_pl2pl(
    pls=earth,
    plf=mars,
    ms=ms_nd,
    nseg=nseg,
    cut=0.6,
    t0_bounds=[7360, 8300.0],
    tof_bounds=[100.*pk.DAY2SEC/TIME, 350.*pk.DAY2SEC/TIME],
    mf_bounds=[2/3,1.],
    vinf_dep_bounds = [0., 0.],
    vinf_arr_bounds = [0., 0.],
    tas = (ta,ta_var),
    max_thrust=max_thrust/F,
    w_bounds_softmax=[-1., 1.],
    time_encoding='softmax',
    )

print("UDP with softmax encoding created successfully")
print(f"Decision vector size:")
print(f"With softmax: {len(udp_g_softmax.get_bounds()[0])} (includes softmax weights)")
print(f"Without softmax: {len(udp_g_uniform.get_bounds()[0])}")

UDP with softmax encoding created successfully
Decision vector size:
With softmax: 161 (includes softmax weights)
Without softmax: 131

Solve with Softmax Encoding#

#solve with gradients:
prob_g_softmax = pg.problem(udp_g_softmax)
prob_g_softmax.c_tol = 1e-6
print("\n=== Optimization WITH Gradients (Softmax) ===")
masses_softmax = []
xs_softmax = []

for i in range(10):
    pop_g_sm = pg.population(prob_g_softmax, 1)
    pop_g_sm = algo.evolve(pop_g_sm)
    if(prob_g_softmax.feasibility_f(pop_g_sm.champion_f)):
        print(".", end="")
        masses_softmax.append(pop_g_sm.champion_x[1])
        xs_softmax.append(pop_g_sm.champion_x)
    else:
        print("x", end ="")
print("\nBest mass is: ", np.max(masses_softmax)*MASS)
print("Worst mass is: ", np.min(masses_softmax)*MASS)
best_idx_sm = np.argmax(masses_softmax)

=== Optimization WITH Gradients (Softmax) ===
..x.....x.
Best mass is:  1233.97957177
Worst mass is:  1215.06076073

Display Best Solution with Softmax#

if len(xs) > 0:
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')  
    t0 = udp_g_softmax.compute_t0(xs_softmax[best_idx_sm])
    ax=pk.plot.add_planet_orbit(ax=ax, pla=earth, label='earth', units=1., color='green')
    ax=pk.plot.add_planet_orbit(ax=ax, pla=mars, label='mars', units=1., color='orange')
    ax = udp_g_softmax.plot(xs_softmax[best_idx_sm], ax=ax, N=20, c='k', s=5)
    ax=pk.plot.add_planet(ax=ax, when=t0, pla=earth, units=1., color='darkgreen', s=20)
    ax=pk.plot.add_planet(ax=ax, when=t0+xs_softmax[best_idx_sm][10+4*udp_g_softmax.nseg]*udp_g_softmax.TIME/pk.DAY2SEC, pla=mars, units=1., color='darkorange', s=20)
    ax.set_title("Best ZOH Transfer Trajectory (with gradients, uniform)")
    ax.legend()
../_images/d5a6580d22513458df196ddacac98b12d6f4f085f0e09f6fa3b929dbd567bb90.png 
udp_g_softmax.pretty(xs_softmax[best_idx_sm])

Low-thrust ZOH transfer (free velocities, non-dimensional formulation)
Departure: earth(jpl_lp)(K)
Arrival: mars(jpl_lp)(K)

Launch epoch: 7452.11294 MJD2000, a.k.a. 2020-05-27T02:42:37.988699
Arrival epoch: 7802.11294 MJD2000, a.k.a. 2021-05-12T02:42:37.988699
Time of flight: 6.02073 (nd), 350.00000 days

Departure excess velocity:
  Magnitude: 0.000000 (nd), 0.000000 km/s
  Vector (nd): [-0.000000, -0.000000, 0.000000]
  Vector (km/s): [-0.000000, -0.000000, 0.000000]
  Direction: [-0.458249, -0.152246, 0.875688], norm: 1.000000

Arrival excess velocity:
  Magnitude: 0.000000 (nd), 0.000000 km/s
  Vector (nd): [-0.000000, 0.000000, 0.000000]
  Vector (km/s): [-0.000000, 0.000000, 0.000000]
  Direction: [-0.145760, 0.987435, 0.061043], norm: 1.000000

Final mass: 0.822653 (nd)

Scaling factors: L=149.598 Gm, V=29.785 km/s, TIME=58.132 days

Details on the ZOH leg:
Dynamics dimension: 7
Control dimension: 4
Number of segments: 30
Number of fwd segments: 18
Number of bck segments: 12
Cut parameter: 0.6
Variational integrator available: true
Maximum propagation steps: none

Initial state: [-0.41228838581227, -0.925570767758029, 2.8874712003856994e-05, 0.8973183813042093, -0.4107088381114756, 1.2812741965299295e-05, 1]
Final state: [-1.0344479955409878, 1.2809107326753661, 0.05223326215793006, -0.6021376946466708, -0.44184620554207943, 0.005524194056897131, 0.8226530478482634]
Time grid: [0, 0.1685490906988014, 0.33732802292945674, 0.503593188137969, 0.6671175546904856, 0.8302544649598823, 0.9933465232375999, 1.1569553904182452, 1.1988311512049292, 1.313949623571403, 1.4302554744611404, 1.6665154548116774, 1.8647987198274172, 2.1423042289218364, 2.447150077222255, 2.6559470918030947, 2.8757686231015436, 3.098569383413903, 3.3776407954400094, 3.563512676366831, 3.8607079525468526, 4.142042764683886, 4.345008366997576, 4.654430712635784, 4.8511154486701225, 5.08037077668507, 5.296174810752023, 5.463525174773438, 5.650075440410627, 5.835537225325473, 6.02073463248503]
Control values: [0.06745267560754033, 0.9108248640320349, -0.28513347431941516, -0.29849115720547403, 0.06745267560754033, 0.9548049454435172, -0.14201926610668633, -0.26110925770737275, 0.06745267560754033, 0.9760461288945445, 0.005177118466990535, -0.21750207353195697, 0.06745267560754033, 0.973740911554742, 0.15282661848975584, -0.16873844260517037, 0.06745267560754033, 0.9477550579707571, 0.2973243852611719, -0.11557925602030741, 0.06745267560754033, 0.8983118843217188, 0.4354129610071857, -0.05874788177724918, 0.06745267560754033, 0.8265801657065955, 0.5628187248628127, 0.0005600591667233601, 0, 0.5464471413223676, -0.7883826412699249, -0.28257447281180464, 0.06745267560754033, 0.7241542989164973, 0.6863533522427722, 0.06722820453371321, 0.06745267560754033, 0.6484537202463332, 0.7533739850938446, 0.10924933442412793, 0, 0.2562456545193487, 0.9033449912099412, 0.3439563800772957, 0, 0.07885949971211509, 0.8440556644412343, 0.5304255033742449, 0, 0.19203596493002237, 0.9386706458960573, 0.2863906540138387, 0, -0.24138615554125953, 0.4708729904059395, -0.8485348259318405, 0.06745267560754033, -0.22752952087966508, 0.8815545763154753, 0.4136325013739846, 0.06745267560754033, -0.36077440579069525, 0.8239262469475226, 0.4370210152423137, 0.06745267560754033, -0.4845542700275393, 0.7493280416183028, 0.45134759061880986, 0, -0.6131741386631627, 0.6789427432607051, 0.40379973631531735, 0, -0.7834205452466179, 0.4260920205107331, 0.4524354532334673, 0, -0.8532284854340219, 0.326182378872719, 0.40694742579008725, 0, -0.9185736936084746, 0.15421994713434195, 0.36392111413926753, 0, -0.9449751000935996, 0.021033457863548398, 0.3264653945725246, 0, -0.9608982013128101, -0.06713308268988197, 0.26864064458372855, 0, -0.9517989526589168, -0.1455184234778065, 0.270005818726047, 0, -0.9399170869818104, -0.2897698548830744, 0.18052506972902596, 0, -0.9072286308556379, -0.4006918908645311, 0.12799304649554832, 0, -0.8897223216075358, -0.44172327284520796, 0.11521606077379296, 0.06745267560754033, -0.8574394389752797, -0.5071784871419585, 0.0869919122561812, 0.06745267560754033, -0.8097016612572941, -0.5853825636332625, 0.04135789251476087, 0.06745267560754033, -0.7544142325999624, -0.6563904804467803, -0.0032719337987038937]

Mismatch constraints: [-3.4107927593396425e-10, -8.086553648922745e-11, 3.6761030851317766e-11, -2.948822297454967e-10, 1.4732881581380752e-11, 6.848102640510945e-11, -7.053346795515836e-11]

udp_g_softmax.plot_throttle(xs_softmax[best_idx_sm])

<Axes: xlabel='time grid', ylabel='throttle value (nd)'>
../_images/58121dec835036f50b39a2f03aab2d7e029fff0ccbade01ed1b1152010bbd91f.png