Trajectory Optimization

Trajectory Optimization#

In pykep both direct and indirect optimization transcriptions as well as evolutionary encodings are provided to perform spacecraft trajectory optimization. Most direct techniques provided are variations to the “Sims” transcription [SF97], while indirect techniques are based on some version of Pontryagin Maximum Principle. The evolutionary encoding are mostly based on the work performed at ESA’ Advanced Concepts Team ([Izz10], [ISM+13]).

Most of the classes in this Trajectory Optimization module are provided as optimization problems (NLPs) in the for of User Defined Problems (UDP) compatible with the pygmo [BI20] python package.

Direct#

When solving optimal control problems (OCPs), the unknown is a function (i.e. the control), hence one is dealing with an infinite dimensional search space. A common approach to overcome this difficulty is to use direct methods. These are based on the idea of discretizing the problem introducing some time grid and thus transforming the OCP into a NLP (Non Linear Programming) problem to then solve it using numerical solvers available for the task.

Sims-Flanagan#

One of pykep’s historical favourite transcription methods is the Sims [SF97] transcription which is called Sims-Flanagan as per the authorship of the first brief conference paper introducing it (anecdotically the work was, though, entirely done by Sims).

class pykep.trajopt.sf_point2point(rvs=[_np.array([1, 0.1, -0.1]) * _pk.AU, _np.array([0.2, 1, -0.2]) * _pk.EARTH_VELOCITY], rvf=[_np.array([-1.2, -0.1, 0.1]) * _pk.AU, _np.array([0.2, -1.023, 0.44]) * _pk.EARTH_VELOCITY], ms=1000., mu=_pk.MU_SUN, max_thrust=0.12, veff=3000 * _pk.G0, tof_bounds=[80., 400.], mf_bounds=[200.0, 1000.0], nseg=10, cut=0.6, mass_scaling=1000, r_scaling=_pk.AU, v_scaling=_pk.EARTH_VELOCITY, with_gradient=True)#

Represents the optimal low-thrust transfer between two fixed points using the Sims-Flanagan (direct) method.

This problem works internally using the sims_flanagan and manipulates its transfer time T, final mass mf and the controls as to link two fixed points in space with a low-thrust trajectory.

It can be used to better profile and understand performances of optimizers on this type of direct approach, but has a limited use in the design of interplanetary trajectories as per the fixed point limitation.

The decision vector is:

z = [mf, throttles, tof]

where throttles is a vector of throttles structures as [u0x, u0y,u0z, …]. By throttles we intend non dimensiona thrust levels in [0,1].

Initializes the sf_point2point instance with given parameters.

Args:

rvs (list): Initial position and velocity vectors. Defaults to two vectors scaled by AU and Earth’s velocity.

rvf (list): Final position and velocity vectors. Defaults to two vectors scaled by AU and Earth’s velocity.

ms (float): Initial spacecraft mass in kg. Defaults to 1000 kg.

mu (float): Gravitational parameter, default is for the Sun (MU_SUN).

max_thrust (float): Maximum thrust in Newtons. Defaults to 0.12 N.

isp (float): Specific impulse in seconds. Defaults to 3000 s.

tof_bounds (list): Bounds for time of flight in days. Defaults to [0, 400] days.

mf_bounds (list): Bounds for final mass in kg. Defaults to [200.0, 1000.0] kg.

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the sims_flanagan. Defaults to 0.6.

mass_scaling (float): Scaling factor for mass (used to scale constraints). Defaults to 1000.

r_scaling (float): Scaling factor for distance, (used to scale constraints). Defaults AU (AU).

v_scaling (float): Scaling factor for velocity (used to scale constraints). Defaults the Earth’s velocity (EARTH_VELOCITY).

with_gradient (bool): Indicates if gradient information should be used. Defaults True.

plot(x, ax=None, units=_pk.AU, show_midpoints=False, show_gridpoints=False, show_throttles=False, length=0.1, arrow_length_ratio=0.05, **kwargs)#

Plots the trajectory leg 3D axes.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to _pk.AU.

show_midpoints (bool, optional): Whether to show midpoints on the trajectory. Defaults to False.

show_gridpoints (bool, optional): Whether to show grid points on the trajectory. Defaults to False.

show_throttles (bool, optional): Whether to show throttle vectors. Defaults to False.

length (float, optional): Length of the throttle vectors. Defaults to 0.1.

arrow_length_ratio (float, optional): Arrow length ratio for the throttle vectors. Defaults to 0.05.

**kwargs: Additional keyword arguments for the plot.

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis with the plotted trajectory.

pretty(x)#

Prints a detailed representation of the Point to point problem.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

class pykep.trajopt.sf_pl2pl(pls=_pk.planet(_pk.udpla.jpl_lp(body='EARTH')), plf=_pk.planet(_pk.udpla.jpl_lp(body='MARS')), ms=1500, mu=_pk.MU_SUN, max_thrust=0.12, veff=3000 * _pk.G0, t0_bounds=[6700.0, 6800.0], tof_bounds=[200.0, 300.0], mf_bounds=[1300.0, 1500.0], vinfs=3.0, vinff=0.0, nseg=10, cut=0.6, mass_scaling=1500, r_scaling=_pk.AU, v_scaling=_pk.EARTH_VELOCITY, with_gradient=True)#

Represents the optimal low-thrust transfer between two planet using the Sims-Flanagan (direct) method.

This problem works internally using the sims_flanagan and manipulates its initial and final states, as well as its transfer time T, final mass mf and the controls as to link the two planets with a low-thrust trajectory.

The particular transcription used is suitable only for few revolutions, after which convergence will start to be problematic.

The decision vector is:

z = [t0, mf, Vsx, Vsy, Vsz, Vfx, Vfy, Vfz, throttles, tof] - all in S.I. units except t0 and tof in days

where throttles is a vector of throttles structured as [u0x, u0y,u0z, …]. By throttles we intend non dimensiona thrust levels in [0,1].

Args:

pls (planet): Initial planet. Defaults to jpl_lp Earth.

plf (planet): Final planet. Defaults to jpl_lp Mars.

ms (float): Initial spacecraft mass in kg. Defaults to 1000 kg.

mu (float): Gravitational parameter, default is for the Sun (MU_SUN).

max_thrust (float): Maximum thrust in Newtons. Defaults to 0.12 N.

veff (float): Exhaust velocity (Specific impulse * _pk.G0) in m/s. Defaults to 3000 s * _pk.G0.

t0_bounds (list): Bounds for departure epoch in MJD2000. Defaults to [6700.0, 6800.0].

tof_bounds (list): Bounds for time of flight in days. Defaults to [200, 300] days.

mf_bounds (list): Bounds for final mass in kg. Defaults to [1300.0, 1500.0] kg.

vinfs (float): Allowed magnitude for the departure’s relative velocity in km/s. Defaults to 3.

vinff (float): Allowed magnitude for the arrival’s relative velocity in km/s. Defaults to 0.

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the sims_flanagan. Defaults to 0.6.

mass_scaling (float): Scaling factor for mass (used to scale constraints). Defaults to 1500.

r_scaling (float): Scaling factor for distance, (used to scale constraints). Defaults AU (AU).

v_scaling (float): Scaling factor for velocity (used to scale constraints). Defaults the Earth’s velocity (EARTH_VELOCITY).

with_gradient (bool): Indicates if gradient information should be used. Defaults True.

plot(x, ax=None, units=_pk.AU, show_midpoints=False, show_gridpoints=False, show_throttles=False, length=0.1, arrow_length_ratio=0.05, **kwargs)#

Plots the trajectory leg 3D axes.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to _pk.AU.

show_midpoints (bool, optional): Whether to show midpoints on the trajectory. Defaults to False.

show_gridpoints (bool, optional): Whether to show grid points on the trajectory. Defaults to False.

show_throttles (bool, optional): Whether to show throttle vectors. Defaults to False.

length (float, optional): Length of the throttle vectors. Defaults to 0.1.

arrow_length_ratio (float, optional): Arrow length ratio for the throttle vectors. Defaults to 0.05.

**kwargs: Additional keyword arguments for the plot.

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis with the plotted trajectory.

pretty(x)#

Prints a human readable representation of the transfer.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

class pykep.trajopt.sf_pl2pl_alpha(pls=_pk.planet(_pk.udpla.jpl_lp(body='EARTH')), plf=_pk.planet(_pk.udpla.jpl_lp(body='MARS')), ms=1500, mu=_pk.MU_SUN, max_thrust=0.12, veff=3000 * _pk.G0, t0_bounds=[6700.0, 6800.0], tof_bounds=[200.0, 300.0], mf_bounds=[1300.0, 1500.0], vinfs=3.0, vinff=0.0, nseg=10, cut=0.6, mass_scaling=1500, r_scaling=_pk.AU, v_scaling=_pk.EARTH_VELOCITY)#

Represents the optimal low-thrust transfer between two planet using the Sims-Flanagan (direct) method.

This problem works internally using the sims_flanagan and manipulates its initial and final states, as well as its transfer time T, final mass mf and the controls as to link the two planets with a low-thrust trajectory.

The particular transcription used is suitable only for few revolutions, after which convergence will start to be problematic.

The decision vector is:

z = [t0, mf, Vsx, Vsy, Vsz, Vfx, Vfy, Vfz, alphas, throttles, tof] - all in S.I. units except t0 and tof in days

where throttles is a vector of throttles structured as [u0x, u0y,u0z, …]. By throttles we intend non dimensiona thrust levels in [0,1].

Args:

pls (planet): Initial planet. Defaults to jpl_lp Earth.

plf (planet): Final planet. Defaults to jpl_lp Mars.

ms (float): Initial spacecraft mass in kg. Defaults to 1000 kg.

mu (float): Gravitational parameter, default is for the Sun (MU_SUN).

max_thrust (float): Maximum thrust in Newtons. Defaults to 0.12 N.

veff (float): Exhaust velocity (Specific impulse * _pk.G0) in m/s. Defaults to 3000 s * _pk.G0.

t0_bounds (list): Bounds for departure epoch in MJD2000. Defaults to [6700.0, 6800.0].

tof_bounds (list): Bounds for time of flight in days. Defaults to [200, 300] days.

mf_bounds (list): Bounds for final mass in kg. Defaults to [1300.0, 1500.0] kg.

vinfs (float): Allowed magnitude for the departure’s relative velocity in km/s. Defaults to 3.

vinff (float): Allowed magnitude for the arrival’s relative velocity in km/s. Defaults to 0.

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the sims_flanagan. Defaults to 0.6.

mass_scaling (float): Scaling factor for mass (used to scale constraints). Defaults to 1500.

r_scaling (float): Scaling factor for distance, (used to scale constraints). Defaults AU (AU).

v_scaling (float): Scaling factor for velocity (used to scale constraints). Defaults the Earth’s velocity (EARTH_VELOCITY).

plot(x, ax=None, units=_pk.AU, show_midpoints=False, show_gridpoints=False, show_throttles=False, length=0.1, arrow_length_ratio=0.05, **kwargs)#

Plots the trajectory leg 3D axes.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to _pk.AU.

show_midpoints (bool, optional): Whether to show midpoints on the trajectory. Defaults to False.

show_gridpoints (bool, optional): Whether to show grid points on the trajectory. Defaults to False.

show_throttles (bool, optional): Whether to show throttle vectors. Defaults to False.

length (float, optional): Length of the throttle vectors. Defaults to 0.1.

arrow_length_ratio (float, optional): Arrow length ratio for the throttle vectors. Defaults to 0.05.

**kwargs: Additional keyword arguments for the plot.

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis with the plotted trajectory.

pretty(x)#

Prints a human readable representation of the transfer.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

Zero-Order-Hold#

The pykep original extension to the original Sims-Flanagan model is called the Zero-Order-Hold transcription method as it substitutes impulses with constant control segments. It has been developed and constantly refined throughout years starting from the early appearence in 2010 in the work from Yam et al. [YIB10] who reported it as a high-fidelity version of the original transcription.

class pykep.trajopt.zoh_point2point(states=None, statef=None, ms=1.0, max_thrust=0.22, tof_bounds=None, mf_bounds=None, nseg=10, cut=0.6, tas=(_pk.ta.get_zoh_kep(1e-10), None), time_encoding='uniform', w_bounds_softmax=None, inequalities_for_tc=False, max_steps=None, state2cart=None)#

Represents the optimal low-thrust transfer between two fixed points using pykep’s Zero Order Hold (direct) trajectory legs.

This problem works internally using the zoh and manipulates its transfer time T, final mass mf and the controls as to link two fixed points in space with a low-thrust trajectory.

It can be used to better profile and understand performances of optimizers on this type of direct approach, but has a limited use in the design of interplanetary trajectories as per the fixed point limitation.

The decision vector is:

x = [mf] + controls + tof (+ [weights for softmax])

where controls is a vector of control parameters \([T, i_x, i_y, i_z] \times n_\text{seg}\) representing magnitude and direction of the thrust applied in each segment.

Initializes the zoh_point2point instance with given parameters.

Args:

states (list): Initial state (only the first 6 states, i.e. no mass). Units as expected by the numerical integrator, defaults to [1.2, 0.0, -0.01, 0.01, 1.0, -0.01].

statef (list): Final state (only the first 6 states, i.e. no mass). Units as expected by the numerical integrator, defaults to [1.0, 0.0, -0.0, 0.01, 1.1, -0.0].

ms (float): Initial mass. Units as expected by the numerical integrator, defaults to 1.

max_thrust (float): Maximum thrust. Units as expected by the numerical integrator, defaults to 0.22.

tof_bounds (list): Bounds for time of flight. Units as expected by the numerical integrator, defaults to [3.4, 8.6].

mf_bounds (list): Bounds for final mass. Units as expected by the numerical integrator, defaults to [0.2, 1].

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the zoh. Defaults to 0.6.

w_bounds_softmax (list): Bounds for the softmax weights (only used if time_encoding is ‘softmax’). Defaults to [-1.0, 1.0].

state2cart (callable): Optional callable for converting dynamics states back to Cartesian.

It maps a 6D state (in the dynamics’ representation) to 6D Cartesian positions and velocities (in integrator units). Used by plot() and pretty() methods for visualization.

inequalities_for_tc (bool): If True, the throttle constraints are formulated as inequalities (|i_u|**2 - 1 <= 0), otherwise they are formulated as equalities (|i_u| = 1). Defaults to False (equality formulation).

tas (tuple): (ta, ta_var) Taylor-adaptive integrators

  • ta: Nominal dynamics (state dim 7, pars ≥ 4)

  • ta_var: Variational dynamics (state dim 84, same pars). When None, no gradients will be used.

max_steps (int or None): Maximum number of Taylor integrator steps per propagation call. When None, uses the default integrator behavior.

plot(x, ax=None, N=30, mark_segments=True, mark_mismatch=True, orbit_color=None, tof=None, **kwargs)#

Plots the trajectory of the zero order hold point to point problem.

Args:

x (list): The decision vector containing final mass, thrust direction, time of flight, and (if time encoding is softmax) the softmax weights.

ax (matplotlib.axes.Axes): The matplotlib axes to plot on. If None, a new figure and axes will be created.

N (int): The number of points to plot along the trajectory.

mark_segments (bool): adds markers at each segment edge

mark_mismatch (bool): highlights mismatch constraints

orbit_color (str): when not None is the color of the border trajectories propagated for the same length as the trajectory encoded by x.

tof (float): time of flight for the border trajectories (only used if orbit_color is not None). If None, uses the tof encoded in x.

**kwargs (dict): Additional keyword arguments to pass to the plotting functions (e.g., marker size, color, etc.)

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The modified Axes object with the Lambert’s problem trajectory added.

pretty(x)#

Prints a detailed representation of the zero order hold point to point problem.

Args:

x (list): The decision vector containing final mass, thrust direction, time of flight, and (if time encoding is softmax) the softmax weights.

class pykep.trajopt.zoh_pl2pl(pls=None, plf=None, ms=1.0, max_thrust=0.22, t0_bounds=None, tof_bounds=None, mf_bounds=None, vinf_dep_bounds=None, vinf_arr_bounds=None, nseg=10, cut=0.6, tas=None, cart2state=None, state2cart=None, inequalities_for_tc=False, time_encoding='uniform', w_bounds_softmax=None, max_steps=None, nrevs=None, L=_pk.AU, V=_pk.EARTH_VELOCITY)#

Represents the optimal low-thrust transfer between two planet using the Zero Order Hold (direct) method with free departure and arrival velocities and under a generic (low-thrust) dynamics.

This problem works internally using the zoh and manipulates its initial and final states, as well as its transfer time T, final mass mf and the controls as to link two orbits described via planets to a low-thrust trajectory. The spacecraft can also be allowed to depart and arrive with residual relative velocities.

This class manages two distinct and independent unit systems:

  1. Ephemeris Units (from planet)

    • Time input: MJD2000 days

    • Position/velocity output: user-defined (e.g., SI: meters, meters/second)

    • Critical contract: velocities and accelerations must be per second

  2. Integrator Units (from the Taylor adaptive integrator tas)

    • These are free, typically assuming μ = 1 as well as a length scale

The transition between systems is controlled by three user-supplied parameters:

  • L (default: AU): length scale mapping ephemeris position to Integrator Units

  • V (default: EARTH_VELOCITY): velocity scale mapping ephemeris velocity to Integrator Units

  • TIME = L / V: derived time scale (in SI seconds) converting seconds to Integrator Units

  • ACC = V² / L: derived acceleration scale converting ephemeris acceleration to Integrator Units

As a consequence, the class applies these standardized conversions to the units used in the backbone numerical integration:

  • Position: $r_{tas} = r_{pla} / L$

  • Velocity: $v_{tas} = v_{pla} / V$

  • Acceleration: $a_{tas} = a_{pla} / text{ACC}$

  • Time: $t_{tas} = t_{pla} cdot text{DAY2SEC} / text{TIME}$

where $t_{pla}$ is in days.

In a nuthsell: choose L and V to match the units returned by pls/plf, but those ephemeris rates must be per second (not per day or something else).

The decision vector is:

x = [t0_fractional, mf, vinf_dep, idep_x, idep_y, idep_z, vinf_arr, iarr_x, iarr_y, iarr_z] + controls + [tof] (+ [weights] for softmax)

where:

  • t0_fractional: non-dimensional departure time offset in integrator time units (TIME); maps to MJD2000 via t0 = t0_bounds[0] + x[0] * TIME * SEC2DAY

  • mf: non-dimensional final mass

  • vinf_dep, vinf_arr: non-dimensional excess velocity magnitudes (scaled by V)

  • idep, iarr: unit Cartesian direction vectors of the relative velocities.

  • controls: [T, i_x, i_y, i_z] * nseg where T is non-dimensional throttle.

    The interpretation of i_x, i_y, i_z is dynamics-dependent: any frame is admissible, as long as they are constrained compatibly with throttle constraints.

  • tof: non-dimensional time of flight (scaled by TIME)

  • weights: softmax weights (only if time_encoding='softmax')

Note

The API differs from the point2point problem: unit scales (L, V) must be explicitly provided to correctly bridge ephemeris and integrator unit systems. Mismatched scales will silently produce incorrect trajectories and gradients.

Initializes a planet-to-planet ZOH low-thrust transfer problem with free departure and arrival excess velocities.

Args:
pls (planet): Departure planet or ephemeris provider.

It must implement eph(t_mjd2000) and, if gradients are requested, also acc(t_mjd2000). The returned position, velocity and acceleration units must be consistent with L and V. Defaults to JPL low-precision Earth.

plf (planet): Arrival planet or ephemeris provider.

It follows the same contract as pls. Defaults to JPL low-precision Mars.

ms (float): Initial spacecraft mass in integrator units.

This is the mass used in self.leg.state0 and therefore must match the mass normalization assumed by the chosen dynamics in tas. Defaults to 1.0.

max_thrust (float): Maximum thrust magnitude in integrator units.

The throttle component stored in each segment of the decision vector is scaled by this value before being assigned to the ZOH leg. Defaults to 0.22.

t0_bounds (list): Lower and upper bounds for the departure epoch,

expressed in MJD2000 days. The first decision variable is a non-dimensional time in integrator units (TIME), bounded by [0, (t0_bounds[1]-t0_bounds[0]) * DAY2SEC / TIME] and mapped to MJD2000 via t0 = t0_bounds[0] + x[0] * TIME * SEC2DAY. Defaults to [6700.0, 6800.0].

tof_bounds (list): Bounds for the time of flight in integrator time units.

Conversion to days is performed internally using TIME / DAY2SEC when querying the arrival ephemeris. Defaults to [3.4, 8.6].

mf_bounds (list): Bounds for the final spacecraft mass in integrator units.

The optimization objective is -mf. Defaults to [0.2, 1.0].

vinf_dep_bounds (list): Bounds for the magnitude of the departure excess

velocity in integrator velocity units, that is, after division by V. Defaults to [0.0, 0.2].

vinf_arr_bounds (list): Bounds for the magnitude of the arrival excess

velocity in integrator velocity units, that is, after division by V. Defaults to [0.0, 0.2].

nseg (int): Number of constant-control segments used by the ZOH transcription.

Each segment contributes four decision variables: throttle magnitude and a direction vector. Defaults to 10.

cut (float): Cut parameter passed to zoh.

It defines the internal forward/backward split used by the transcription. Defaults to 0.6.

inequalities_for_tc (bool): If True, the throttle constraints are

formulated as inequalities (|i_u|**2 - 1 <= 0), otherwise they are formulated as equalities (|i_u| = 1). Defaults to False (equality formulation).

tas (tuple): (ta, ta_var) Taylor-adaptive integrators

  • ta: Nominal dynamics (state dim 7, pars ≥ 4)

  • ta_var: Variational dynamics (state dim 84, same pars). When None, no gradients will be used.

cart2state (list): Optional list [transform, jacobian] for non-Cartesian dynamics.

transform is a callable that maps a 6D Cartesian state (in integrator units) to the 6D state expected by the dynamics. jacobian is a callable that returns the 6×6 derivative of transform w.r.t. the Cartesian state, as a flattened 36-element array (row-major order). When provided, the class applies this transformation to boundary states before propagation, and uses the Jacobian for gradient computation. Both functions receive the state as a 1D array. Defaults to None (Cartesian dynamics).

state2cart (callable): Optional callable for converting dynamics states back to Cartesian.

It maps a 6D state (in the dynamics’ representation) to 6D Cartesian positions and velocities (in integrator units). Used by plot() and pretty() methods for visualization. Ignored if cart2state is not provided. Defaults to None.

time_encoding (str): Time-grid encoding scheme.

'uniform' uses equally spaced segment boundaries over the total time of flight. 'softmax' adds nseg weights to the decision vector and maps them to positive segment durations summing to the total time of flight. Defaults to 'uniform'.

nrevs (int): Use only in case the dynamics is in mean orbital elements (MEE) form, this
parameter specifies the number of revolutions to be completed. It is used to add a mean longitude

offset to the final state to enforce the desired number of revolutions.

w_bounds_softmax (list): Lower and upper bounds for the softmax weights.

These bounds are used only when time_encoding='softmax'. Defaults to [-1.0, 1.0].

max_steps (int or None): Maximum number of internal integrator steps

allowed by zoh. If None, the integrator default is used.

L (float): Length scale used to map ephemeris positions to integrator units.

If pls and plf return positions in meters, a typical choice is pykep.AU or another length unit consistent with the integrator dynamics.

V (float): Velocity scale used to map ephemeris velocities to integrator units.

If ephemeris velocities are returned in meters per second, V must also be expressed in meters per second.

plot(x, ax=None, N=30, mark_segments=True, mark_mismatch=True, orbit_color=None, tof=None, **kwargs)#

Plots the trajectory of the zero order hold point to point problem.

Args:

x (list): The decision vector. Layout: [t0_frac, mf, vinf_dep, idep_x, idep_y, idep_z, vinf_arr, iarr_x, iarr_y, iarr_z] + controls + [tof] (+ [weights]) in non-dimensional units.

ax (matplotlib.axes.Axes): The matplotlib axes to plot on. If None, a new figure and axes will be created.

N (int): Number of points per segment used when sampling the propagated trajectory.

mark_segments (bool): Adds markers at each segment boundary.

mark_mismatch (bool): If True, highlights the mismatch point between forward and backward integrations.

orbit_color (str or None): If provided, plots the departure and arrival reference orbits (propagated without thrust) using this color.

tof (float or None): If orbit_color is provided, this overrides the tof extracted from x for plotting the border orbits.

**kwargs (dict): Additional keyword arguments passed to the underlying matplotlib plotting calls (e.g., marker size, color for markers).

Returns:

The matplotlib Axes object with the trajectory added.

pretty(x)#

Prints a detailed representation of the zero order hold planet to planet problem with free velocities.

class pykep.trajopt.zoh_ss_point2point(states=None, statef=None, tof_bounds=None, nseg=10, cut=0.6, tas=(_pk.ta.get_zoh_ss(1e-10), None), time_encoding='uniform', w_bounds_softmax=None, max_steps=None)#

Represents the time optimal Solar Sail transfer between two fixed points using Zero Order Hold (direct) trajectory legs.

This problem works internally using zoh with dim_dynamics=6 and dim_controls=2 and manipulates its transfer time tof and the sail clock and cone angles (controls) as to link two fixed points in space with a Solar Sail trajectory.

It can be used to better profile and understand performances of optimizers on this type of direct approach, but has a limited use in the design of interplanetary trajectories as per the fixed point limitation.

The decision vector is:

x = controls + [tof] (+ [weights for softmax])

where controls is a vector of control parameters \([\alpha, \beta] \times n_\text{seg}\) representing the sail cone and clock angles, assumed fixed in the RTN frame along each segment.

Initializes the zoh_ss_point2point instance with given parameters.

Args:

states (list): Initial state. Units as expected by the numerical integrator, defaults to [1.2, 0.0, -0.01, 0.01, 1.0, -0.01].

statef (list): Final state. Units as expected by the numerical integrator, defaults to [1.0, 0.0, -0.0, 0.01, 1.1, -0.0].

tof_bounds (list): Bounds for time of flight. Units as expected by the numerical integrator, defaults to [3.4, 8.6].

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for zoh. Defaults to 0.6.

tas (tuple): (ta, ta_var) Taylor-adaptive integrators

  • ta: Nominal dynamics (state dim 6, pars ≥ 2)

  • ta_var: Variational dynamics (54, same pars). When None, no gradients will be used.

time_encoding (str): Time-grid encoding scheme.

'uniform' uses equally spaced segment boundaries over the total time of flight. 'softmax' adds nseg weights to the decision vector and maps them to positive segment durations summing to the total time of flight. Defaults to 'uniform'.

w_bounds_softmax (list): Lower and upper bounds for the softmax weights.

These bounds are used only when time_encoding='softmax'. Defaults to [-1.0, 1.0].

max_steps (int or None): Maximum number of Taylor integrator steps per propagation call. When None, uses the default integrator behavior.

plot(x, ax=None, N=30, to_cartesian=lambda state: ..., mark_segments=True, mark_mismatch=True, plot_sail=True, sail_size=0.05, **kwargs)#

Plots the trajectory of the zero order hold point to point problem.

Args:
x (list): The decision vector containing controls,

time of flight, and (if time encoding is softmax) the softmax weights.

ax (matplotlib.axes.Axes): The matplotlib axes to plot on.

If None, a new figure and axes will be created.

N (int): The number of points to plot along the trajectory.

to_cartesian (Callable): A function that converts whatever

state is used in the internal Taylor integrator to Cartesian (r,v).

mark_segments (bool): Adds markers at each segment edge.

mark_mismatch (bool): Marks the terminal mismatch point.

plot_sail (bool): Adds a visualization of a rectangular sail.

sail_size (float): Half-side length of the rendered sail square,

in the same units as the trajectory positions.

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The modified Axes object with the trajectory added.

Notes:

If propagation fails on one or more segments, this method still plots the successfully propagated forward/backward segments. This is intentional so users can visually inspect where and how the trajectory integration breaks.

pretty(x)#

Prints a detailed representation of the zero order hold point to point problem.

Args:

x (list): The decision vector containing controls, time of flight, and (if time encoding is softmax) the softmax weights.

Indirect#

Indirect methods are based on the Pontryagin Maximum Principle (PMP), which provides necessary conditions for optimality. These methods involve deriving the optimal control laws and the corresponding state trajectories by solving a two point boundary value problem (TPBVP) derived from the PMP.

class pykep.trajopt.pontryagin_cartesian_mass(posvel0=_posvel0, posvelf=_posvelf, tof=250, mu=_pk.MU_SUN, lambda0=1.0, eps=1e-4, T_max=0.6, Isp=3000, m0=1500, L=_pk.AU, MU=_pk.MU_SUN, MASS=1500, with_gradient=False, taylor_tolerance=1e-16, taylor_tolerance_var=1e-8)#

This class is a pygmo (http://esa.github.io/pygmo/) UDP representing a single-shooting indirect method for the minimum mass, fixed time, optimization of a low-thrust trajectory.

The decision vector is:

x = [lx, ly, lz, lvx, lvy, lvz, lm, l0]

Note

Gradients are also provided via the variational version of the augmented dynamics.

pykep.trajopt.pontryagin_cartesian(start=default, target=default, tof=250, mu=1.32712440018e+20, eps=1e-4, T_max=0.6, Isp=3000, m0=1500, L=1.495978707e+11, TIME=31557600.0, MASS=1500, with_gradient=False)

Args:

posvel0 (list): the initial position and velocity of the spacecraft.

posvelf (list): the final position and velocity of the spacecraft.

tof (float): the time of flight of the trajectory (in days).

mu (float): the gravitational parameter of the central body.

lambda0 (float or None): multiplicative factor for the objective. If None, lambda0 is added to the decision vector and the constraint ||lambda||=1 is added so that all costates must be in [-1,1].

eps (float): the accuracy of the numerical integrator.

T_max (float): the maximum thrust of the spacecraft.

Isp (float): the specific impulse of the spacecraft.

m0 (float): the initial mass of the spacecraft.

L (float): units for length. Default is the astronomical unit (AU).

MU (float): units for the gravitational parameter of the central body. Default is the gravitational parameter of the Sun.

MASS (float): units for mass. Default is 1500 kg.

with_gradient (bool): whether to use the gradient of the constraints or not.

taylor_tolerance (float): the tolerance for the Taylor integrator.

taylor_tolerance_var (float): the tolerance for the variational Taylor integrator.

plot(x, N=100, ax3D=None)#

This function plots the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot. ax3D (matplotlib.axes._axes.Axes): the axis to use for the plot. If None, a new axis is created.

Returns:

ax3D (matplotlib.axes._axes.Axes): the axis of the plot.

plot_misc(x, N=100)#

This function plots the throttle, thrust direction, switching function, mass costate and Hamiltonian of the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot.

Returns:

axs (list): the list of axis of the plots.

class pykep.trajopt.pontryagin_cartesian_time(start=default, target=default, tof=250, mu=1.32712440018e+20, eps=1e-4, T_max=0.6, Isp=3000, m0=1500, L=1.495978707e+11, TIME=31557600.0, MASS=1500, with_gradient=False)#

This class is a pygmo (http://esa.github.io/pygmo/) UDP representing a single-shooting indirect method for the minimum time, optimization of a low-thrust trajectory.

The decision vector is:

x = [lx, ly, lz, lvx, lvy, lvz, lm, lJ, tof]

Note

Gradients are also provided via the variational version of the augmented dynamics.

Args:

source (list): the initial planet.

target (list): the final planet.

t0 (epoch): the initial epoch.

tof_guess (float): a guess for the time of flight. Bound will be defined as tof_guess/2 and 2*tof_guess.

lambda0 (float or None): multiplicative factor for the objective. If None, lambda0 is added to the decision vector and the constraint ||lambda||=1 is added so that all costates must be in [-1,1].

T_max (float): the maximum thrust of the spacecraft.

Isp (float): the specific impulse of the spacecraft.

m0 (float): the initial mass of the spacecraft.

L (float): units for length. Default is the astronomical unit (AU).

MU (float): units for the gravitational parameter of the central body. Default is the gravitational parameter of the Sun.

MASS (float): units for mass. Default is 1500 kg.

with_gradient (bool): whether to use the gradient of the constraints or not.

taylor_tolerance (float): the tolerance for the Taylor integrator.

taylor_tolerance_var (float): the tolerance for the variational Taylor integrator.

plot(x, ax3D=None, N=100, **kwargs)#

This function plots the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. ax3D (matplotlib.axes._axes.Axes): the axis to use for the plot. If None, a new axis is created. N (int): the number of points to use in the plot. **kwargs: Additional keyword arguments to pass to the trajectory plot.

Returns:

ax3D (matplotlib.axes._axes.Axes): the axis of the plot.

plot_misc(x, N=100, **kwargs)#

This function plots the throttle, switching function and Hamiltonian of the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot. **kwargs: Additional keyword arguments to pass to the plt.subplots function.

Returns:

axs (list): the list of axis of the plots.

class pykep.trajopt.pontryagin_equinoctial_mass(start=default, target=default, tof=250, mu=1.32712440018e+20, eps=1e-4, T_max=0.6, Isp=3000, m0=1500, L=1.495978707e+11, TIME=31557600.0, MASS=1500, with_gradient=False)#

This class is a pygmo (http://esa.github.io/pygmo/) UDP representing a single-shooting indirect method for the minimum mass, fixed time, optimization of a low-thrust trajectory. The coordinates used for the dynamics and hence the TPBVP are the equinoctial elements.

The decision vector is:

x = [lp, lf, lg, lh, lk, lL, lm, l0]

Note

Gradients are also provided via the variational version of the augmented dynamics.

Args:

posvel0 (list): the initial position and velocity of the spacecraft.

posvelf (list): the final position and velocity of the spacecraft.

tof (float): the time of flight of the trajectory (in days).

mu (float): the gravitational parameter of the central body.

eps (float): the accuracy of the numerical integrator.

lambda0 (float or None): multiplicative factor for the objective. If None, lambda0 is added to the decision vector and the constraint ||lambda||=1 is added so that all costates must be in [-1,1].

T_max (float): the maximum thrust of the spacecraft.

Isp (float): the specific impulse of the spacecraft.

m0 (float): the initial mass of the spacecraft.

L (float): units for length. Default is the astronomical unit (AU).

n_rev (int): number of revolutions. If negative any number will be allowed.

MU (float): units for the gravitational parameter of the central body. Default is the gravitational parameter of the Sun.

MASS (float): units for mass. Default is 1500 kg.

with_gradient (bool): whether to use the gradient of the constraints or not.

taylor_tolerance (float): the tolerance for the Taylor integrator.

taylor_tolerance_var (float): the tolerance for the variational Taylor integrator.

plot(x, N=100, ax3D=None)#

This function plots the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot. ax3D (matplotlib.axes._axes.Axes): the axis to use for the plot. If None, a new axis is created.

Returns:

ax3D (matplotlib.axes._axes.Axes): the axis of the plot.

plot_misc(x, N=100)#

This function plots the throttle, thrust direction, switching function, mass costate and Hamiltonian of the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot.

Returns:

axs (list): the list of axis of the plots.

class pykep.trajopt.pontryagin_equinoctial_time(source=_pl0, target=_plf, t0=_pk.epoch(0), tof_guess=250, lambda0=1.0, T_max=0.6, veff=3000 * _pk.G0, m0=1500, n_rev=-1, L=_pk.AU, MU=_pk.MU_SUN, MASS=1500, with_gradient=False, taylor_tolerance=1e-16, taylor_tolerance_var=1e-8)#

This class is a pygmo (http://esa.github.io/pygmo/) UDP representing a single-shooting indirect method for the minimum time optimization of a low-thrust trajectory. The coordinates used for the dynamics and hence the TPBVP are the modified equinoctial elements.

The decision vector is:

x = [lp, lf, lg, lh, lk, lL, lm, (l0), tf]

Note

Gradients are also provided via the variational version of the augmented dynamics.

pykep.trajopt.pontryagin_cartesian_time(start=default, target=default, tof=250, mu=1.32712440018e+20, eps=1e-4, T_max=0.6, Isp=3000, m0=1500, L=1.495978707e+11, TIME=31557600.0, MASS=1500, with_gradient=False)

Args:

source (list): the initial planet.

target (list): the final planet.

t0 (epoch): the initial epoch.

tof_guess (float): a guess for the time of flight. Bound will be defined as tof_guess/2 and 2*tof_guess.

lambda0 (float or None): multiplicative factor for the objective. If None, lambda0 is added to the decision vector and the constraint ||lambda||=1 is added so that all costates must be in [-1,1].

T_max (float): the maximum thrust of the spacecraft.

veff (float): the effective velocity of the spacecraft propulsion system.

m0 (float): the initial mass of the spacecraft.

n_rev (int): number of revolutions. If negative any number will be allowed.

L (float): units for length. Default is the astronomical unit (AU).

MU (float): units for the gravitational parameter of the central body. Default is the gravitational parameter of the Sun.

MASS (float): units for mass. Default is 1500 kg.

with_gradient (bool): whether to use the gradient of the constraints or not.

taylor_tolerance (float): the tolerance for the Taylor integrator.

taylor_tolerance_var (float): the tolerance for the variational Taylor integrator.

plot(x, N=100, ax3D=None)#

This function plots the trajectory encoded in the decision vector x. x = [lx, ly, lz, lvx, lvy, lvz, lm, (l0), tof]

Args:

x (list): the decision vector. N (int): the number of points to use in the plot. ax3D (matplotlib.axes._axes.Axes): the axis to use for the plot. If None, a new axis is created.

Returns:

ax3D (matplotlib.axes._axes.Axes): the axis of the plot.

plot_misc(x, N=100)#

This function plots the throttle, thrust direction, switching function, mass costate and Hamiltonian of the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot.

Returns:

axs (list): the list of axis of the plots.

Evolutionary Encodings#

Some type of interplanetary trajectories can be evolved: shocking?. This AI approach to trajectory design resulted in two silver ([TAF+05], [PGJ+18]) and one gold ([ISM+13]) Humies award for human-competitive results that were produced by any form of genetic and evolutionary computation. In pykep we offer some classes that help instantiating these type of transfers amenable to evolutionary techniques.

class pykep.trajopt.mga(seq, t0, tof, vinf, multi_objective=False, tof_encoding="direct"", orbit_insertion=False, e_target=None, rp_target=None)#

The Multiple Gravity Assist (MGA) encoding of an interplanetary trajectory.

This class may be used as a User Defined Problem (UDP) for the pygmo (http://esa.github.io/pygmo/) optimisation suite.

The decision vector (chromosome) is:

direct encoding: z = [t0, T1, T2 ... ] in [mjd2000, days, days ... ]
alpha encoding: z = [t0, T, a1, a2 ...] in [mjd2000, days, nd, nd ... ]
eta encoding: z = [t0, n1, n2, n3 ...] in [mjd2000, nd, nd ...]

Note

The time of flights of a MGA trajectory (and in general) can be encoded in different ways. When they are directly present in the decision vector, we talk about a direct encoding. This is the most ‘evolvable’ encoding but also the one that requires the most problem knowledge (e.g. to define the bounds on each leg) and is not very flexible in dealing with constraints on the total time of flight. The alpha and eta encodings, instead, allow to only specify bounds on the time of flight of the entire trajectory, and not on the single legs: a property that is attractive for multi-objective optimization, for example.

In the alpha encoding each leg time-of-flight is decoded as follows,

\(T_i = T \log\alpha_i / \sum_n\log\alpha_n\).

In the eta encoding each leg time-of-flight is decoded as follows,

\(T_i = (T_{max} - \sum_{j=0}^{(i-1)}T_j) \eta_i\)

The chromosome dimension for the direct and eta encoding is the same, while the alpha encoding requires one more gene.

Note

The resulting problem is box-bounded (unconstrained).

Args:

seq (list [planet]): sequence of planetary encounters including the departure body.

t0 (list [float or epoch]): lower and upper bounds for the launch epoch. When floats are used MJD2000 is assumed.

tof (list or float): defines the bounds on the time of flight. If tof_encoding is ‘direct’, this contains a list of 2D lists defining the upper and lower bounds on each leg. If tof_encoding is ‘alpha’, this contains a 2D list with the lower and upper bounds on the total time-of-flight. If tof_encoding is ‘eta’ tof is a float defining an upper bound for the time-of-flight.

vinf (float): the vinf provided at launch for free

multi_objective (bool): when True constructs a multiobjective problem (dv, T). In this case, ‘alpha’ or eta encodings are recommended

tof_encoding (str): one of ‘direct’, ‘alpha’ or ‘eta’. Selects the encoding for the time of flights

orbit_insertion (bool): when True the arrival dv is computed as that required to acquire a target orbit defined by e_target and rp_target

e_target (float): if orbit_insertion is True this defines the target orbit eccentricity around the final planet

rp_target (float): if orbit_insertion is True this defines the target orbit pericenter around the final planet (in m)

max_revs (int): maximal number of revolutions for lambert transfer

plot(x, ax=None, units=_pk.AU, N=60, c_orbit='dimgray', c_lambert='indianred', leg_ids=[], figsize=(5, 5), **kwargs)#

Plots the trajectory encoded into x in 3D axes.

Args:

x (list): The decision vector in the correct tof encoding.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to pk.AU.

N (int, optional): The number of points to use when plotting the trajectory. Defaults to 60.

c_orbit (str, optional): The color of the planet orbits. Defaults to ‘dimgray’.

c (str, optional): The color of the trajectory. Defaults to ‘indianred’.

figsize (tuple): The figure size (only used if a*ax* is None and axis have to be created.), Defaults to (5, 5).

leg_ids (list): selects the legs to plot. Optional, defaults to all legs.

**kwargs: Additional keyword arguments to pass to the trajectory plot (all Lambert arcs)

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis where the trajectory was plotted.

pretty(x)#

Prints a human readable representation of the transfer.

Args:

x (list): The decision vector in the correct tof encoding.

to_planet(x: List[float])#

Returns a planet representing the spacecraft trajectory encoded into x.

Args:

x (list): The decision vector in the correct tof encoding.

Example:

sc = mga_udp.to_planet(population.champion_x) r, v = sc.eph(7000.)

class pykep.trajopt.mga_1dsm(seq, t0, tof, vinf_bounds=[0.5, 2.5], add_vinf_dep=False, add_vinf_arr=True, tof_encoding='direct', multi_objective=False, orbit_insertion=False, e_target=None, rp_target=None, eta_bounds=[0.01, 0.99], rp_ub=30)#

This class transcribes a Multiple Gravity Assist trajectory with one deep space maneuver per leg (MGA-1DSM) into an optimisation problem. It may be used as a User Defined Problem (UDP) for the pygmo (http://esa.github.io/pygmo/) optimisation suite.

  • Izzo, Dario. “Global optimization and space pruning for spacecraft trajectory design.” Spacecraft Trajectory Optimization 1 (2010): 178-200.

The decision vector (chromosome) is:

direct encoding: [t0] + [u, v, Vinf, eta1, T1] + [beta, rp/rV, eta2, T2] + ...
alpha encoding:  [t0] + [u, v, Vinf, eta1, a1] + [beta, rp/rV, eta2, a2] + ... + [T]
eta encoding:    [t0] + [u, v, Vinf, eta1, n1] + [beta, rp/rV, eta2, n2] + ...

where t0 is a mjd2000, Vinf is in km/s, T in days, beta in radians and the rest non dimensional.

Note

The time of flights of a MGA-1DSM trajectory (and in general) can be encoded in different ways. When they are directly present in the decision vector, we have a direct encoding. This is the most ‘evolvable’ encoding but also the one that requires the most problem knowledge (e.g. to define the bounds on each leg) and is not very flexible in dealing with constraints on the total time of flight. The alpha and eta encodings, instead, allow to only specify bounds on the time of flight of the entire trajectory, and not on the single legs: a property that is attractive for multi-objective optimization, for example.

In the alpha encoding each leg time-of-flight is decoded as follows, T_i = T log(alpha_i) / sum_n(log(alpha_n)). In the eta encoding each leg time-of-flight is decoded as follows, T_i = (tof_max - sum_0^(i-1)(T_j)) * eta_i

The chromosome dimension for the direct and eta encoding is the same, while the alpha encoding requires one more gene.

Note

The resulting problem is box-bounded (unconstrained).

seq (list [planet]): sequence of planetary encounters including the departure body.

t0 (list [float or epoch]): lower and upper bounds for the launch epoch. When floats are used MJD2000 is assumed.

tof (list or float): defines the bounds on the time of flight. If tof_encoding is ‘direct’, this contains a list of 2D lists defining the upper and lower bounds on each leg. If tof_encoding is ‘alpha’, this contains a 2D list with the lower and upper bounds on the total time-of-flight. If tof_encoding is ‘eta’, tof is a float defining an upper bound for the time-of-flight.

vinf_bounds (list [float]): the minimum and maximum allowed initial hyperbolic velocity (at launch), in km/sec.

add_vinf_dep (bool): when True the computed Dv includes the initial hyperbolic velocity (at launch).

add_vinf_arr (bool): when True the computed Dv includes the final hyperbolic velocity (at the last planet).

tof_encoding (str): one of ‘direct’, ‘alpha’ or ‘eta’. Selects the encoding for the time of flights.

multi_objective (bool): when True constructs a multiobjective problem (dv, T).

orbit_insertion (bool): when True the arrival dv is computed as that required to acquire a target orbit defined by e_target and rp_target.

e_target (float): if orbit_insertion is True this defines the target orbit eccentricity around the final planet.

rp_target (float): if orbit_insertion is True this defines the target orbit pericenter around the final planet (in m).

eta_bounds (list [float]): lower and upper bounds for the eta variable (defining the position of the DSM).

rp_ub (float): upper bound for the pericenter radius of the flybys (in planetary radii).

plot(x, ax=None, units=_pk.AU, N=60, c_orbit='dimgray', c_lambert='indianred', c_ballistic='royalblue', leg_ids=[], figsize=(5, 5), **kwargs)#

Plots the trajectory encoded into x in 3D axes.

Args:

x (list): The decision vector in the correct tof encoding.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): the 3D axis. If None a new one is created.

units (float, optional): length units to be used. Defaults to pk.AU.

N (int, optional): number of points to be used in the plots. Defaults to 60.

c_orbit (str, optional): color of the planetary orbits. Defaults to ‘dimgray’.

c (str, optional): color of the trajectory. Defaults to ‘indianred’.

figsize (tuple, optional): size of the figure. Defaults to (5, 5).

**kwargs: Additional keyword arguments to pass to the trajectory plotting functions.

pretty(x)#

  • x: encoded trajectory

Prints human readable information on the trajectory represented by the decision vector x

Example:

print(prob.pretty(x))
class pykep.trajopt.pl2pl_N_impulses(start='earth', target='venus', N_max=3, tof=[20., 400.], vinf=[0., 4.], phase_free=True, multi_objective=False, t0=None)#

This class is a pygmo (http://esa.github.io/pygmo/) UDP representing a single leg transfer between two planets allowing up to a maximum number of impulsive Deep Space Maneuvers.

The decision vector is:

[t0,T] + [alpha,u,v,V_inf] * (N-2) + [alpha] + ([tf])

… in the units: [mjd2000, days] + [nd, nd, m/sec, nd] + [nd] + [mjd2000]

Each time-of-flight can be decoded as follows, T_n = T log(alpha_n) / sum_i(log(alpha_i))

Note

The resulting problem is box-bounded (unconstrained). The resulting trajectory is time-bounded.

Args:

start (planet): initial body.

target (planet): target body.

N_max (int): maximum number of impulses.

tof_bounds (list): the box bounds [lower,upper] for the total time of flight (days).

DV_max_bounds (list): the box bounds [lower,upper] for each DV magnitude (km/sec).

phase_free (bool): when True, no rendezvous condition are enforced and start and arrival anomalies will be free.

multi_objective (bool): when True, a multi-objective problem is constructed with DV and time of flight as objectives.

t0_bounds (list): the box bounds on the launch window containing two pykep.epoch. This is not needed if phase_free is True.

plot(x, ax=None, units=_pk.AU, N=60, c_orbit='dimgray', c_segments=['royalblue', 'indianred'], figsize=(5, 5), **kwargs)#

Plots the trajectory encoded into x in 3D axes.

Args:

x (list): The decision vector in the correct tof encoding.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to pk.AU.

N (int, optional): The number of points to use when plotting the trajectory. Defaults to 60.

c_orbit (str, optional): The color of the planet orbits. Defaults to ‘dimgray’.

c_segments (list, optional): The colors to alternate the various trajectory segments (inbetween DSMs). Defaults to [“royalblue”, “indianred”].

figsize (tuple): The figure size (only used if a*ax* is None and axis have to be created.), Defaults to (5, 5).

**kwargs: Additional keyword arguments to pass to the trajectory plot (common to Lambert arcs and ballistic arcs)

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis where the trajectory was plotted.

plot_primer_vector(x, N=200, ax=None)#

Plots the primer vector magnitude along the trajectory encoded in x.

Args:

x (list): The decision vector in the correct tof encoding.

N (int, optional): The number of points to use when plotting the primer vector. Defaults to 200.

ax (matplotlib.axes.Axes, optional): The axis to plot on. Defaults to None.

Returns:

matplotlib.axes.Axes: The axis where the primer vector was plotted. tuple: A tuple containing the grid and the primer vector magnitude.

Utilities#

In order to facilitate the use of the classes in this module, some utilities are provided.

pykep.trajopt.primer_vector(DVi, DVj, Mji, Mjk)#

This function computes the primer vector in a point k, relative to finite impulses in i and j.

Args:

DVi (ndarray - (3,)): the impulse at point i.

DVj (ndarray - (3,)): the impulse at point j.

Mji (ndarray - (6,6)): the state transition matrix from i to j (dxj = Mji dxi).

Mjk (ndarray - (6,6)): the state transition matrix from k to j (dxj = Mji dxi).

Returns:

tuple: The primer vector, the Aik matrix, the Ajk matrix.

Note:

The impulse transfer matrix Anm is defined as those matrices that allow to compute the variation of the impulse at point n given the variation of the impulse at point m. In formal terms, dDVn = Anm dDVm. All variations are such that the terminal state is kept fixed.

pykep.trajopt.primer_vector_surrogate(DVk, Mki, Mkj)#

This function computes the surrogate primer vector at two points i and j, corresponding to a single finite impulse at point k.

Args:

DVk (ndarray - (3,)): the impulse at point k.

Mki (ndarray - (6,6)): the state transition matrix from i to k (dxk = Mki dxi).

Mkj (ndarray - (6,6)): the state transition matrix from j to k (dxk = Mkj dxj).

Returns:

tuple: The surrogate primer vector, the Aij matrix, the Akj matrix.

Note:

The impulse transfer matrix Anm is defined as that matrix that allows to compute (in this surrogate case) the variation of the impulse at point n given the variation of the impulse at point m. In formal terms, dDVn = Anm dDVm. All variations are such that the terminal state is kept fixed.

class pykep.trajopt._launchers#

This class contains an API to access the tabulated data on the mass launchers can deliver to a certain declination / vinf.

Note

In pykep the object pykep.trajopt.launchers is an instance of this class created upon import and is used to access launchers data.

Examples:
>>> import pykep as pk
>>> mass = pk.trajopt.launchers.soyuzf(4.5, 33.21)
>>> mass2 = pk.trajopt.launchers.atlas501(4.5, 33.21)
ariane5(vinfs, decls)#

Estimates the mass that the Ariane5 launcher can deliver to a certain vinf and declination, assuming a launch from Kourou. Data provided to ESOC by Arianespace when Ariane launch for ExoMars was an option. If the inputs are arrays, then a mesh is constructed and the mass is returned on all points of the mesh.

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

decls (float or numpy.ndarray (N)):: the declinations in degrees (in the ECF frame).

Returns:

numpy.ndarray (N, N): masses delivered to escape with said declinations and magnitudes.

ariane64(vinfs, decls)#

Estimates the mass that the Ariane 64 launcher can deliver to a certain vinf and declination, assuming a launch from Kourou. Data based on the Ariane 6 User’s Manual (Issue 2 Rev. 0, Feb. 2021) and ESA mission call technical annexes. If the inputs are arrays, a mesh is constructed and the mass is returned on all points of the mesh.

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

decls (float or numpy.ndarray (N)): the declinations in degrees (in the ECF frame).

Returns:

numpy.ndarray (N, N): masses delivered to escape in kg.

atlas501(vinfs, decls)#

Estimates the mass that the Atlas 501 launcher can deliver to a certain vinf and declination. If the inputs are arrays, then a mesh is constructed and the mass is returned on all points of the mesh.

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

decls (float or numpy.ndarray (N)):: the declinations in degrees (in the ECF frame).

Returns:

numpy.ndarray (N, N): masses delivered to escape with said declinations and magnitudes.

atlas551(vinfs)#

Estimates the mass that the Atlas 551 launcher can deliver to a certain vinf

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

Returns:

numpy.ndarray (N): masses delivered to escape with said vinf magnitudes.

falcon9(vinfs)#

Estimates the mass that a Falcon 9 (Block 5, expendable) can deliver to a certain vinf for interplanetary escape. Launch from Cape Canaveral (28.5 deg N). The declination dependence is not officially published and is not modeled here. Data from SpaceX Falcon 9 Payload User’s Guide (2025) and NASA LSP.

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

Returns:

numpy.ndarray (N): masses delivered to escape in kg.

soyuzf(vinfs, decls)#

Estimates the mass that the Soyutz-Fregat launcher can deliver to a certain vinf and declination. If the inputs are arrays, then a mesh is constructed and the mass is returned on all points of the mesh.

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

decls (float or numpy.ndarray (N)):: the declinations in degrees (in the ECF frame).

Returns:

numpy.ndarray (N, N): masses delivered to escape with said declinations and magnitudes.