Interplanetary transfer legs

class pykep.leg.sims_flanagan(rvs=[[1, 0, 0], [0, 1, 0]], ms=1., throttles=[0, 0, 0, 0, 0, 0], rvf=[[0, 1, 0], [-1, 0, 0]], mf=1., tof=pi / 2, max_thrust=1., veff=1., mu=1., cut=0.5)

This class represents an interplanetary low-thrust transfer between a starting and a final point in the augmented state-space \([\mathbf r, \mathbf v, m]\). The low-thrust transfer is described by a sequence of equally spaced impulses as described in:

Sims, J., Finlayson, P., Rinderle, E., Vavrina, M. and Kowalkowski, T., 2006, August. Implementation of a low-thrust trajectory optimization algorithm for preliminary design. In AIAA/AAS Astrodynamics specialist conference and exhibit (p. 6746).

The low-thrust transfer will be feasible is the state mismatch equality constraints and the throttle mismatch inequality constraints are satisfied.

Args:

rvs (2D array-like): Cartesian components of the initial position vector and velocity [[xs, ys, zs], [vxs, vys, vzs]]. Defaults to [[1,0,0], [0,1,0]].

ms (float): initial mass. Defaults to 1.

throttles (1D array-like): the Cartesan components of the throttle history [ux1, uy1, uz1, ux2, uy2, uz2, …..]. Defaults to a ballistic, two segments profile [0,0,0,0,0,0].

rvf (2D array-like): Cartesian components of the final position vector and velocity [[xf, yf, zf], [vxf, vyf, vzf]]. Defaults to [[0,1,0], [-1,0,0]].

mf (float): final mass. Defaults to 1.

tof (float): time of flight. Defaults to \(\frac{\pi}{2}\).

max_thrust (float): maximum level for the spacecraft thrust. Defaults to 1.

veff (float): effective velocity of the propulsion system. Defaults to 1.

mu (float): gravitational parameter. Defaults to 1.

cut (float): the leg cut, in [0,1]. It determines the number of forward and backward segments. Defaults to 0.5.

Note

Units need to be consistent.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> sf = pk.leg.sims_flanagan()

compute_mc_grad()

Computes the gradients of the mismatch constraints. Indicating the initial augmented state with \(\mathbf x_s = [\mathbf r_s, \mathbf v_s, m_s]\), the final augmented state with \(\mathbf x_f = [\mathbf r_f, \mathbf v_f, m_f]\), the total time of flight with \(T\) and introducing the throttle vector \(\mathbf u = [u_{x0}, u_{y0}, u_{z0}, u_{x1}, u_{y1}, u_{z1} ]\) and \(\mathbf {\tilde u} = [\mathbf u, T]\) (note the time of flight at the end), this method computes the following gradients:

\[\frac{\partial \mathbf {mc}}{\partial \mathbf x_s} \rightarrow (7\times7)\]

\[\frac{\partial \mathbf {mc}}{\partial \mathbf x_f} \rightarrow (7\times7)\]

\[\frac{\partial \mathbf {mc}}{\partial \mathbf {\tilde u}} \rightarrow (7\times(3\mathbf{nseg} + 1))\]

Returns:

tuple [numpy.ndarray, numpy.ndarray, numpy.ndarray]: The three gradients. sizes will be (7,7), (7,7) and (7, 3nseg + 1)

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> sf = pk.leg.sims_flanagan()
>>  sf.throttles = [0.8]*3
>>> sf.compute_mc_grad()

compute_mismatch_constraints()

In the Sims-Flanagan trajectory leg model, a forward propagation is performed from the starting state as well as a backward from the final state. The state values thus computed need to match in some middle control point. This is typically imposed as 7 independent constraints called mismatch-constraints computed by this method.

Returns:

list [float]: The seven mismatch constraints in the same units used to construct the leg.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> sf = pk.leg.sims_flanagan()
>>> sf.compute_mismatch_constraints()

compute_tc_grad()

Computes the gradients of the throttles constraints. Indicating the total time of flight with \(T\) and introducing the throttle vector \(\mathbf u = [u_{x0}, u_{y0}, u_{z0}, u_{x1}, u_{y1}, u_{z1} ]\), this method computes the following gradient:

\[\frac{\partial \mathbf {tc}}{\partial \mathbf u} \rightarrow (\mathbf{nseg} \times3\mathbf{nseg})\]

Returns:

tuple [numpy.ndarray]: The gradient. Size will be (nseg,nseg*3).

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> sf = pk.leg.sims_flanagan()
>>  sf.throttles = [0.8]*3
>>> sf.compute_tc_grad()

compute_throttle_constraints()

In the Sims-Flanagan trajectory leg model implemented in pykep, we introduce the concept of throttles. Each throttle is defined by three numbers \([u_x, u_y, u_z] \in [0,1]\) indicating that a certain component of the thrust vector has reached a fraction of its maximum allowed value. As a consequence, along the segment along which the throttle is applied, the constraint \(u_x ^2 + u_y ^2 + u_z^2 = 1\), called a throttle constraint, has to be met.

Returns:

list [float]: The throttle constraints.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> sf = pk.leg.sims_flanagan()
>>  sf.throttles = [0.8]*3
>>> sf.compute_throttle_constraints()

property cut

The leg cut: it determines the number of forward and backward segments.

property max_thrust

Maximum spacecraft thruet.

property mf

Final mass.

property ms

Initial mass.

property mu

Central body gravitational parameter.

property nseg

The total number of segments

property nseg_bck

The total number of backward segments

property nseg_fwd

The total number of forward segments

property rvf

The final position vector and velocity: [[xs, ys, zs], [vxs, vys, vzs]].

property rvs

The initial position vector and velocity: [[xs, ys, zs], [vxs, vys, vzs]].

property throttles

The Cartesan components of the throttle history [ux1, uy1, uz1, ux2, uy2, uz2, …..].

property tof

Time of flight.

property veff

Effective velocity of the propulsion system (Isp*G0 in the V units of the dynamics)

---

class pykep.leg.zoh(state0, controls, state1, tgrid, cut, tas)

This class implements an interplanetary low-thrust transfer between a starting and final state in the augmented state-space \([\mathbf{r}, \mathbf{v}, m]\). The transfer is modelled as a sequence of non-uniform segments along which a continuous and constant (zero-order hold) control acts. The time intervals defining these segments are also provided in tgrid.

The formulation generalises pykep.leg.sims_flanagan to arbitrary dynamics and non-uniform time grids. The dynamics are assumed to be zero-order hold and must be provided as compatible Taylor-adaptive integrators (tas).

Note

The requirements on the tas passed are: a) the first four heyoka parameters must be \(T, i_x, i_y, i_z\), b) the system dimension must be 7 c) for the variational integrator, variations on the state and the four parameters only are considered. These requirements are all fulfilled by pykep.ta.zoh_kep, pykep.ta.zoh_eq, pykep.ta.zoh_cr3bp and their variational versions.

A transfer is feasible when the state mismatch equality constraints are satisfied. In the intended usage, throttle equality constraints are also enforced to ensure a proper thrust representation as \(T \hat{\mathbf{i}}\) with \(|\hat{\mathbf{i}}| = 1\).

\[i_x^2 + i_y^2 + i_z^2 = 1, \quad \forall \text{segments}\]

Note

The variational integrator ta_var can be None or must have state dimension 84 (7 + 7x7 STM + 7x4 control sensitivity) and with the same dynamics as the nominal integrator ta. It is the user that must ensure the suitability of the integrators.

Args:

state0 (array-like): Initial state \([\mathbf{r}_0, \mathbf{v}_0, m_0]\) (length 7)

controls (array-like): Control parameters \([T, i_x, i_y, i_z] \times n_\text{seg}\)

state1 (array-like): Final state \([\mathbf{r}_1, \mathbf{v}_1, m_1]\) (length 7)

tgrid (array-like): Non-uniform time grid (length nseg+1)

cut (float): Forward/backward segment split ratio (0 ≤ cut ≤ 1)

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)

Raises:

ValueError: If state/parameter dimensions mismatch or input lengths are incompatible.

Examples:

import numpy as np
import heyoka as hy
# Define states, controls, time grid
state0 = np.array([1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0])
state1 = np.array([1.2, 0.1, 0.0, 0.0, 0.9, 0.1, 0.95])
controls = np.array([0.022, 0.7, 0.7, 0.1, 0.025, -0.3, 0.8, 0.4, 0.015, -0.2, 0.8, 0.4])
tgrid = np.array([0.0, 0.5, 1.0, 1.23])
# Get integrators
ta = pk.ta.get_zoh_eq(tol=1e-16)
ta_var = pk.ta.get_zoh_eq_var(tol=1e-16)
# Construct leg (50/50 split)
leg = zoh(state0, controls, state1, tgrid, cut=0.5, tas=(ta, ta_var))

compute_mc_grad()

Computes the gradients of the mismatch constraints. Indicating the initial augmented state with \(\mathbf x_s = [\mathbf r_s, \mathbf v_s, m_s]\), the final augmented state with \(\mathbf x_f = [\mathbf r_f, \mathbf v_f, m_f]\), the time grid as \(T_{grid}\) and the introducing the control vector \(\mathbf u = [T_0, i_{x0}, i_{y0}, i_{z0}, T_1, i_{x1}, i_{y1}, i_{z1}]\) (note the time of flight at the end), this method computes the following gradients:

\[\frac{\partial \mathbf {mc}}{\partial \mathbf x_s} \rightarrow (7\times7)\]

\[\frac{\partial \mathbf {mc}}{\partial \mathbf x_f} \rightarrow (7\times7)\]

\[\frac{\partial \mathbf {mc}}{\partial \mathbf u} \rightarrow (7\times(4\mathbf{nseg}))\]

\[\frac{\partial \mathbf {mc}}{\partial \mathbf T_{grid}} \rightarrow (7\times(\mathbf{nseg} + 1))\]

Returns:

tuple [numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]: The four gradients. sizes will be (7,7), (7,7) (7,4nseg) and (7,nseg+1)

compute_tc_grad()

Computes the gradients of the throttles constraints. Introducing the control vector as \(\mathbf u = [T_0, i_{x0}, i_{y0}, i_{z0}, T_1, i_{x1}, i_{y1}, i_{z1}, ...]\), this method computes the following gradient:

\[\frac{\partial \mathbf {tc}}{\partial \mathbf u} \rightarrow (\mathbf{nseg} \times 4\mathbf{nseg}) \]

Returns:

tuple [numpy.ndarray]: The gradient. Size will be (nseg,4nseg).

get_state_info(N=2)

This method returns state histories sampled along each ZOH segment, for both the forward and backward propagation parts of the leg. The sampling is performed by calling heyoka.taylor_adaptive.propagate_grid() on a uniformly-spaced grid of N points within each segment.

Args:

N (int): Number of sampling points per segment (including the segment endpoints). The default (N=2) returns only the segment endpoints.

Returns:

tuple: (state_fwd, state_bck) where:

• state_fwd (list): List of length nseg_fwd. Each entry contains the sampled 7D state history over the corresponding forward segment (from tgrid[i] to tgrid[i+1]).

• state_bck (list): List of length nseg_bck. Each entry contains the sampled 7D state history over the corresponding backward segment (from tgrid[-1-i] to tgrid[-2-i]).

Note

The backward propagation is carried out by integrating from the final time toward earlier times; depending on the backend conventions, the returned grids may thus be time-reversed with respect to a forward-time plot.

Note

This method uses the nominal integrator self.ta and overwrites its internal time, state and the first four parameters (T, i_x, i_y, i_z). If the integrator state must be preserved, call this method on a dedicated copy.

Examples:

ax = pk.plot.make_3Daxis()
fwd, bck = leg.get_state_info(N=100)
for segment in fwd:
    ax.scatter(segment[0,0], segment[0,1], segment[0,2], c='blue')
    ax.plot(segment[:,0], segment[:,1], segment[:,2], c='blue')
for segment in bck:
    ax.scatter(segment[0,0], segment[0,1], segment[0,2], c='darkorange')
    ax.plot(segment[:,0], segment[:,1], segment[:,2], c='darkorange')
ax.view_init(90, -90)

---

class pykep.leg.zoh_ss(state0, controls, state1, tgrid, cut, tas)

This class implements a simple solar sailing interplanetary leg between a starting and final state in the augmented state-space \([\mathbf{r}, \mathbf{v}]\). The transfer is modelled as a sequence of non-uniform segments along which the sail orientation is considered as fixed in the RTN frame (zero-order hold).

The class is thus very similar to pykep.leg.zoh but has different requirements on the dynamics it allows since it represents a solar sail trajectory (i.e. no mass equation and two controls only)

Note

The requirements on the tas passed are: a) the first two heyoka parameters must be \(\alpha, \beta\), b) the system dimension must be 6.

The time intervals defining these segments are also provided in tgrid.

Note

The tas requirements are fulfilled by pykep.ta.zoh_ss and pykep.ta.zoh_ss_var.

A transfer is feasible when the state mismatch equality constraints are satisfied.

Note

The variational integrator ta_var can be None or must have state dimension 54 (6 + 6x6 STM + 6x2 control sensitivity) and with the same dynamics as the nominal integrator ta. It is the user that must ensure the suitability of the integrators.

Args:

state0 (array-like): Initial state \([\mathbf{r}_0, \mathbf{v}_0]\) (length 6)

controls (array-like): Sail orientation controls \([\alpha, \beta] \times n_\text{seg}\)

state1 (array-like): Final state \([\mathbf{r}_1, \mathbf{v}_1]\) (length 6)

tgrid (array-like): Non-uniform time grid (length nseg+1)

cut (float): Forward/backward segment split ratio (0 ≤ cut ≤ 1)

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

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

• ta_var: Variational dynamics (state dim 54, same pars)

Raises:

ValueError: If state/parameter dimensions mismatch or input lengths are incompatible.

Examples:

import numpy as np
import heyoka as hy
# Define states, controls, time grid
state0 = np.array([1.0, 0.0, 0.0, 0.0, 1.0, 0.0])
state1 = np.array([1.2, 0.1, 0.0, 0.0, 0.9, 0.1])
controls = np.array([np.pi/2, 0.] * 3)
tgrid = np.array([0.0, 0.5, 1.0, 1.23])
# Get integrators
ta = pk.ta.get_zoh_ss(tol=1e-16)
ta_var = pk.ta.get_zoh_ss_var(tol=1e-16)
# Set their params
ta.pars[2] = 0.01
ta_var.pars[2] = 0.01
# Construct leg (50/50 split)
leg = zoh_ss(state0, controls, state1, tgrid, cut=0.5, tas=(ta, ta_var))

compute_mc_grad()

Computes the gradients of the mismatch constraints. Indicating the initial augmented state with \(\mathbf x_s = [\mathbf r_s, \mathbf v_s]\), the final augmented state with \(\mathbf x_f = [\mathbf r_f, \mathbf v_f]\), the time grid as \(T_{grid}\) and the introducing the control vector \(\mathbf u = [\alpha_1, \beta_1, \alpha_2, \beta_2, ...]]\) (note the time of flight at the end), this method computes the following gradients:

\[\frac{\partial \mathbf {mc}}{\partial \mathbf x_s} \rightarrow (6\times 6)\]

\[\frac{\partial \mathbf {mc}}{\partial \mathbf x_f} \rightarrow (6\times 6)\]

\[\frac{\partial \mathbf {mc}}{\partial \mathbf u} \rightarrow (6\times(2\mathbf{nseg}))\]

\[\frac{\partial \mathbf {mc}}{\partial \mathbf T_{grid}} \rightarrow (6\times(\mathbf{nseg} + 1))\]

Returns:

tuple [numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]: The four gradients. sizes will be (6,6), (6,6) (6,2nseg) and (6,nseg+1)

get_state_info(N=2)

This method returns state histories sampled along each ZOH segment, for both the forward and backward propagation parts of the leg. The sampling is performed by calling heyoka.taylor_adaptive.propagate_grid() on a uniformly-spaced grid of N points within each segment.

Args:

N (int): Number of sampling points per segment (including the segment endpoints). The default (N=2) returns only the segment endpoints.

Returns:

tuple: (state_fwd, state_bck) where:

• state_fwd (list): List of length nseg_fwd. Each entry contains the sampled 6D state history over the corresponding forward segment (from tgrid[i] to tgrid[i+1]).

• state_bck (list): List of length nseg_bck. Each entry contains the sampled 6D state history over the corresponding backward segment (from tgrid[-1-i] to tgrid[-2-i]).

Note

The backward propagation is carried out by integrating from the final time toward earlier times; depending on the backend conventions, the returned grids may thus be time-reversed with respect to a forward-time plot.

Note

This method uses the nominal integrator self.ta and overwrites its internal time, state and the first two parameters (\([\alpha, \beta]\)). If the integrator state must be preserved, call this method on a dedicated copy.

Examples:

ax = pk.plot.make_3Daxis()
fwd, bck = leg.get_state_info(N=100)
for segment in fwd:
    ax.scatter(segment[0,0], segment[0,1], segment[0,2], c='blue')
    ax.plot(segment[:,0], segment[:,1], segment[:,2], c='blue')
for segment in bck:
    ax.scatter(segment[0,0], segment[0,1], segment[0,2], c='darkorange')
    ax.plot(segment[:,0], segment[:,1], segment[:,2], c='darkorange')
ax.view_init(90, -90)