Numerical Propagation#

The backbone of numerical propagation in pykep is based on Lagrangian coefficients for Kepler’s dynamics and Taylor numerical integration, as implemented in the Heyoka [BI21] python package, for all other cases. The state transition matrix is also available and provided, in the case of numerical integration, seamlessly via variational equations.

The main routines are listed here:

Keplerian#

pykep.propagate_lagrangian(rv=[[1, 0, 0], [0, 1, 0]], tof=pi / 2, mu=1, stm=False)#

Propagates (Keplerian) the state for an assigned time and computes the State Transition Matrix (if requested) using the Lagrangian coefficients.

Args:

rv (2D array-like): Cartesian components of the initial position vector and velocity [[x0, y0, z0], [v0, vy0, vz0]]. Defaults to [[1,0,0], [0,1,0]].

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

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

stm (bool): requests the computations of the State Transition Matrix

Returns:

tuple (list, list): r and v, that is the final position and velocity after the propagation. (if stm is False) tuple (tuple (list, list), numpy.ndarray (6,6)): (r,v) and the STM. (if stm is True)

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> r0 = [1,0,0]
>>> v0 = [0,1,0]
>>> tof = pi/2
>>> mu = 1
>>> [r1,v1], stm = pk.propagate_lagrangian(rv=[r0,v0], tof = tof, mu = mu, stm = True)
>>> [r1,v1] = pk.propagate_lagrangian(rv=[r0,v0], tof = tof, mu = mu, stm = False)

pykep.propagate_lagrangian_grid(rv=[[1, 0, 0], [0, 1, 0]], tofs=[pi / 2], mu=1, stm=False)#

This function calls pykep.propagate_lagrangian() in a loop to provide the results on a time grid containing several epochs. The main advantage is in the API which is convenient and reminiscent of the heyoka.propagate_grid() interface.

Note that this function is not necessarily more efficient than calling pykep.propagate_lagrangian() in a loop, since there is no parallelization nor SIMD magic implemented atm. Nevertheless we offer this interface for cenvenience as it may allow more compact code.

Args:

rv (2D array-like): Cartesian components of the initial position vector and velocity [[x0, y0, z0], [v0, vy0, vz0]]. Defaults to [[1,0,0], [0,1,0]].

tof (1D array-like): time of flight. Defaults to [\(\frac{\pi}{2}\)].

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

stm (bool): requests the computations of the State Transition Matrix

Returns:

list [tuple ( list , list ) ]: For each time of flight: [r,v], that is the final position and velocity after the propagation and the flattened stm (if requested).

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> r0 = [1,0,0]
>>> v0 = [0,1,0]
>>> tofs = [pi/2, pi, 3/4pi]
>>> mu = 1
>>> res = pk.propagate_lagrangian_grid(rv = [r0, v0], tofs = tofs, mu = mu, stm = True)
>>> rs = [it[0][0] for it in res]
>>> vs = [it[0][1] for it in res]
>>> stms = [it[1] for it in res]

Non-Keplerian#

class pykep.stark_problem(mu=1., veff=1., tol=1e-16)#

Class representing the Stark problem. In pykep, abusing a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body and described by the equations:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot{\mathbf v} = -\frac{\mu}{r^3} \mathbf r + \frac{\mathbf T}{m} \\ \dot m = - \frac{|\mathbf T|}{I_{sp} g_0} \end{array}\right.\end{split}\]

Note

Similar and connected functionality is provided by the functions stark(), stark_var(): and stark_dyn():.

Args:

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

veff (float): propulsion system effective velocity (Isp g0). Defaults to 1.

tol (float): tolerance of the Taylor adaptive integration. Defaults to 1e-16.

Note

Units need to be consistent upon construction as well as later on when calling the propagate methods.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> mu = pk.MU_SUN
>>> veff = 3000. * pk.G0
>>> tol = 1e-14
>>> sp = pk.stark_problem(mu, veff, tol)
>>> sp.propagate(rvm_state = [1., 0., 0., 0., 1., 0., 1], thrust = [0., 0., 1e-8], tof = 7.32)
[0.5089647068650076, 0.8607873878989034, 0.0, -0.8607873878989032, 0.5089647068650074, 0.0, 1.0]

property mu#: The central body gravity parameter.

propagate(rvm_state, thrust, tof)#

Stark problem numerical propagation. In pykep, abusing a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.

The propagation will be singular for vanishing masses (infinite acceleration) and raise an exception.

Args:

rvm_state (list (7,)): position, velocity and mass flattened into a 7D list.

thrust (list (3,)): thrust flattened into a 3D list.

tof (float): time of flight.

Returns:

list (7,) : position, velocity and mass after the numerical propagation, flattened into a 7D list.

Note

Units need to be consistent with the ones used upon constructing the instance.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> mu = pk.MU_SUN
>>> veff = 3000. * pk.G0
>>> tol = 1e-14
>>> sp = pk.stark_problem(mu, veff, tol)
>>> sp.propagate(rvm_state = [1., 0., 0., 0., 1., 0., 1], thrust = [0., 0., 1e-8], 7.32)
[0.5089647068650076, 0.8607873878989034, 0.0, -0.8607873878989032, 0.5089647068650074, 0.0, 1.0]

propagate_var(rvm_state, thrust, tof)#

Stark problem numerical propagation via variational equations.

In pykep, abusing a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.

The propagation will be singular for vanishing masses (infinite acceleration) and raise an exception.

It also computes the system State Transition Matrix:

\[\mathbf M = \frac{d\mathbf x_f}{d\mathbf x_0}\]

as well as the gradients of the final states with respect to the thrust direction.

\[\mathbf U = \frac{d\mathbf x_f}{d\mathbf u}\]

Args:

rvm_state (list (7,)): position, velocity and mass flattened into a 7D list.

thrust (list (3,)): thrust flattened into a 3D list.

tof (float): time of flight.

Returns:

tuple (list (7,), numpy.ndarray (7,7), numpy.ndarray (7,3)): position, velocity and mass (after propagation) flattened into a 7D list, the state transition matrix and the gradient with respect to thrust.

Note

Units need to be consistent upon construction and later on when calling the propagate methods.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> mu = pk.MU_SUN
>>> veff = 3000. * pk.G0
>>> tol = 1e-14
>>> sp = pk.stark_problem(mu, veff, tol)
>>> sp.propagate(state = [1., 0., 0., 0., 1., 0., 1], thrust = [0., 0., 1e-8], 7.32)

property tol#: The Taylor integrator tolerance.

property veff#: The effective velocity (Isp g0)

Taylor adaptive propagators#

In addition to the above mentioned propagators, Taylor adaptive integrators are offered in pykep wrapping some of the functionalities of Heyoka [BI21] python package. Their variational version is also offered (at order one) as to be able to produce stms and, where needed, other useful quantities. Higher order variational equations can also be obtained directly using the available dynamics and using Heyoka [BI21] syntax.

Some of the Taylor adaptive integrators are associated to OCPs (Optimal Control Problems) of relevance to interplanetary flight. In particular they are born when applying Pontryagin principle to dynamics of interests and result in two point boundary value problems (TPBVP). In these cases, a number of auxiliary functiond describing the control, the Hamiltonian, the switching function, etc., are also provided.

Stark#

pykep.ta.get_stark(tol)#

Returns a Taylor adaptive propagator (Heyoka) for the Stark problem retreiving one from a global cache and making a copy.

In pykep, abusing a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.

If the requested propagator was never created this will create it, else it will return the one from the global cache, thus avoiding jitting.

The dynamics is that returned by stark_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_stark(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.,0.,0.,0.,1.,0.,1.]
>>> mu = 1.
>>> veff = 1.
>>> thrust = [0., 0., 0.]
>>> tof = 1.
>>> ta.pars[:] = [mu, veff] + thrust
>>> ta.propagate_until(tof)

pykep.ta.get_stark_var(tol)#

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the Stark problem retreiving one from a global cache and making a copy.

In pykep, abusing a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.

The dynamics is that returned by stark_dyn(): and also used in get_stark()

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_stark_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.,0.,0.,0.,1.,0.,1.]
>>> mu = 1.
>>> veff = 1.
>>> thrust = [0., 0., 0.]
>>> tof = 1.
>>> ta.pars[:5] = [mu, veff] + thrust
>>> ta.propagate_until(tof)

pykep.ta.stark_dyn()#

The dynamics of the Stark problem.

In pykep, abusing a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.

where \(\mu, v_{eff} = I_{sp}g_0\) and \(\mathbf T = [T_x, T_y, T_z]\) are parameters.

Returns:: list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

Circular Restricted Three Body Problem#

pykep.ta.get_cr3bp(tol)#

Returns a Taylor adaptive propagator (Heyoka) for the CR3BP problem retreiving one from a global cache and making a copy.

If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the CR3BP is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The specific dynamics used is that returned by cr3bp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_cr3bp(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = 0.01215058560962404
>>> tof = 5.7856656782589234
>>> ta.pars[0] = mu
>>> ta.propagate_until(tof)

pykep.ta.get_cr3bp_var(tol)#

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the CR3BP problem retreiving one from a global cache and making a copy.

If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the CR3BP is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The specific dynamics used is that returned by cr3bp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_cr3bp_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = 0.01215058560962404
>>> tof = 5.7856656782589234
>>> ta.pars[0] = mu
>>> ta.propagate_until(tof)

pykep.ta.cr3bp_dyn()#

The dynamics of the Circular Restricted Three Body Problem (CR3BP).

In pykep, the CR3BP is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The parameter \(\mu\) is defined as \(\frac{m_2}{m_1+m_2}\) where \(m_2\) is the mass of the secondary body (i.e. placed on the positive x axis).

The equations are non-dimensional with units \(L = r_{12}\) (distance between the primaries), \(M = m_1 + m_2\) (total system mass) and \(T = \sqrt{\frac{r_{12}^3}{m_1+m_2}}\) (period of rotation of the primaries).

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot v_x = 2v_y + x - (1 - \mu) \frac{x + \mu}{r_1^3} - \mu \frac{x + \mu - 1}{r_2^3} \\ \dot v_y = -2 v_x + y - (1 - \mu) \frac{y}{r_1^3} - \mu \frac{y}{r_2^3} \\ \dot v_z = -(1 - \mu) \frac{z}{r_1^3} - \mu \frac{z}{r_2^3} \end{array}\right.\end{split}\]

where \(\mu\) is the only parameter.

Returns:: list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

Low-thrust Pontryagin Cartesian TPBVP#

pykep.ta.get_pc(tol, optimality)#

Returns a Taylor adaptive propagator (Heyoka) for the TPBVP problem resulting from the application of

Pontryagin Maximum Principle (PMC) to the low-trhust problem (constant maximal thrust) in

Cartesian coordinates. If the requested propagator was never created, a call to this function will

trigger its compilation. Otherwise, it will return the one from a global cache, thus avoiding jitting.

The specific dynamics used is that returned by pc_dyn().

Both time optimal and mass optimal systems can be returned by setting the optimality parameter.

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

optimality (pykep.optimality_type): the optimality principle to be used.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_pc(tol = 1e-16, optimality = pk.optimality_type.TIME)
>>> ta.time = 0.
>>> # We set the initial conditions with some arbitrary values (all costates to 1.)
>>> ta.state[:14] = [1., 0., 0., 0., 1., 0., 10., 1., 1., 1., 1., 1., 1., 1.]
>>> ta.pars[:] = [1., 0.01, 1., 0.5, 1.]
>>> tof = 1.2345
>>> ta.propagate_until(tof)

pykep.ta.get_pc_var(tol, optimality)#

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the TPBVP problem resulting

from the application of Pontryagin Maximum Principle (PMP) to the low-thrust problem

(constant maximal thrust) in Cartesian coordinates. If the requested propagator was never created,

a call to this function will trigger its compilation. Otherwise, it will return the one from

a global cache, thus avoiding jitting.

The specific dynamics used is that returned by pc_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

optimality (pykep.optimality_type): the optimality principle to be used.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta_var = pk.ta.get_pc_var(tol = 1e-16)
>>> ta_var.time = 0.
>>> ta_var.state[:14] = [1., 0., 0., 0., 1., 0., 10., 1., 1., 1., 1., 1., 1., 1.]
>>> ta_var.pars[:] = [1., 0.01, 1., 0.5, 1.]
>>> tof = 1.2345
>>> ta_var.propagate_until(tof)

pykep.ta.pc_dyn()#

The augmented dynamics of the TPBVP originating when applying an indirect method to the low-thrust transfer problem in Cartesian coordinates.

The (non-augmented) dynamics is,

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf f_r = \mathbf v \\ \dot{\mathbf v} = \mathbf f_v = -\frac{mu}{r^3}\mathbf r + c_1 \frac um \hat{\mathbf i}\\ \dot{m} = f_m = - \frac{c_1}{c_2} u \end{array} \right.\end{split}\]

The state, containing the spacecraft state and the co-states, is:

\[\mathbf x = [\mathbf r, \mathbf v, m] = [x,y,z,v_x,v_y,v_z,l_x,l_y,l_z,l_{vx},l_{vy},l_{vz},l_m]\]

While the parameters:

\[\mathbf p = [\mu, c_1, c_2, \epsilon, \lambda_0]\]

describe, respectively, the gravitational parameter, the maximum thrust, the effective velocity (product of the specific impulse \(I_{sp}\) by \(g_0\)), an homotopy parameter and a factor multiplying the instantaneous cost (in theory this can be any positive number as the solution of the problem will not change)

The controls, representing magnitude and direction of the spacecraft thrust are:

\[\mathbf u = [u, \hat{\mathbf i}]\]

The final equations of motion to be possibly used in a shooting method, are derived from the Hamiltonian:

\[\mathcal H(\mathbf x, \mathbf \lambda, \mathbf u) = \mathbf \lambda_r \cdot \mathbf f_r + \mathbf \lambda_v \cdot \mathbf f_v + \lambda_m f_m + \lambda_0 \left(u + \epsilon\log(u(1-u))\right)\]

by taking the derivatives:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf x} = \frac{\partial \mathcal H}{\partial \mathbf \lambda} \\ \dot{\mathbf \lambda} = - \frac{\partial \mathcal H}{\partial \mathbf x} \\ \end{array}\right.\end{split}\]

and then made independent of the controls applying Pontryagin minimum principle:

\[\mathbf u^* = \argmin_{\mathbf u \in \mathcal U} \mathcal H(\mathbf x, \mathbf \lambda, \mathbf u)\]

and substituting in the equations the optimal controls found as function of the augmented system state.

Returns:: list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

Numerical Propagation

Contents

Numerical Propagation#

Keplerian#

Non-Keplerian#

Taylor adaptive propagators#

Stark#

Circular Restricted Three Body Problem#

Low-thrust Pontryagin Cartesian TPBVP#