Efficient single step propagation for UKF sim env

I am cross posting this from the python-wrapper issue tracker based on sage advice from Maxime…

I have been shopping around for an efficient numerical or SGP4 propagator for use with a python Gym environment I am working on. The environment will keep a true state for each space vehicle (SV) orbiting the earth as well as an Unscented Kalman Filter (UKF) propagation from the initial state. At each time-step, the an observation of one or more of the SVs will be made, updating their UKF state. I think Orekit might be the solution I am looking for, but I seem to be falling short in my implementation so far. What I am trying to accomplish is to have a method I can pass a current state (PV or Kep) to along with a time step and get a single step integration to the new step with the updated coordinates (PV or Kep - same as used for current state). My experience so far is that the propagation step is plenty fast (7e-5s per step), but the process of initializing each propagator slows the system down to more like 1.8e-1s per initialization and step. Code included below if you want to see how I am implementing Orekit currently to test the computation expense of each sigma point in the UKF.

I have started looking at Shiva-Iyer’s orbdetpy (/github.com/ut-astria/orbdetpy) based on feedback from Maxime, but at first pass their implementation is different than I would do myself. I had planned to keep the variables in memory vice writing them out to files (json), but I am still looking deeper into it as I am sure they had sound reasons for choosing that direction.

1.) Is Orekit even the right option for my propagator? I have run similar tests with PyEphem (SGP) (/rhodesmill.org/pyephem/) and the PyPi SGP4 (/pypi.org/project/sgp4/)
2.) Is there a faster / lighter wait way to initialize this?
3.) Does anyone have any good sample codes for working with Orekit and UKFs beyond orbdetpy?
4.) Any other advise?

//////////////////////// Begin Code Samples //////////////////////////
/////////////////////////////// Prep Orekit //////////////////////////////////

import orekit
orekit.initVM()

#%% Setup Orekit working directory

from orekit.pyhelpers import setup_orekit_curdir
setup_orekit_curdir()
orekit.pyhelpers.download_orekit_data_curdir()

#%% imports for orekit & datetime

from java.util import Arrays
from orekit import JArray_double
from org.orekit.bodies import  OneAxisEllipsoid, GeodeticPoint
from org.orekit.frames import  FramesFactory, TopocentricFrame
from org.orekit.data import DataProvidersManager, ZipJarCrawler
from org.orekit.time import TimeScalesFactory, AbsoluteDate, DateComponents, TimeComponents
from org.orekit.orbits import KeplerianOrbit, CartesianOrbit
from org.orekit.utils import Constants
from org.orekit.propagation.analytical import EcksteinHechlerPropagator, KeplerianPropagator
from org.orekit.propagation.analytical.tle import TLE, TLEPropagator
from org.orekit.propagation.conversion import FiniteDifferencePropagatorConverter
from org.orekit.propagation.conversion import TLEPropagatorBuilder
from datetime import datetime
from org.orekit.propagation import SpacecraftState
from org.orekit.orbits import OrbitType, PositionAngle
from org.orekit.propagation.numerical import NumericalPropagator
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.hipparchus.ode.nonstiff import ClassicalRungeKuttaIntegrator
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.utils import IERSConventions, PVCoordinates, Constants
from java.io import File

#%% imports for matplotlib, math, and numpy

import matplotlib.dates as mdates
from math import radians, pi
import numpy as np
import matplotlib.pyplot as plt
import time

#%% Some constants, ae, mu, utc
ae = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
mu = Constants.WGS84_EARTH_MU
utc = TimeScalesFactory.getUTC()

/////////////////////////////// Numerical Prop ///////////////////////////////

#%% Initial osculating orbits

numsat = 1000
states = []
# Inertial frame where the satellite is defined
inertialFrame = FramesFactory.getEME2000()

PV=True
if PV:
    for i in range(0,numsat):
        epochDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)
        
        # Orbit construction as Keplerian
        p=Vector3D(-1780159.5047552094 + np.random.normal(0,1000),
                   -1164835.539234647 + np.random.normal(0,1000), 
                   6842941.118987356 + np.random.normal(0,1000))
        v=Vector3D(-7205.19342727098 + np.random.normal(0,1),
                   -341.0428975779148 + np.random.normal(0,1),
                   -1924.0083602739878 + np.random.normal(0,1))
        initialOrbit = CartesianOrbit(PVCoordinates(p, v), 
                                      inertialFrame, 
                                      epochDate, mu)
        satellite_mass = 100.0 # kg
        initialState = SpacecraftState(initialOrbit, satellite_mass)
        states.append(initialState)
    print("states defined in cartesian coords")
else:
    for i in range(0,numsat):
        ra = 814 * 1000 + np.random.normal(0,1000)        # km Apogee!
        rp = 786 * 1000 + np.random.normal(0,1000)         # km Perigee!
        i = radians(98.55 + np.random.normal(0,0.5))      # inclination
        omega = radians(90.0 + np.random.normal(0,0.5))   # perigee argument
        raan = radians(5.1917 + np.random.normal(0,0.5))  # right ascension of ascending node
        lv = radians(0.0567634 + np.random.normal(0,0.5))    # True anomaly
        epochDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)
        
        a = (rp + ra + 2 * ae) / 2.0    # semi major axis in km
        e = 1.0 - (rp + ae) / a
        
        # Orbit construction as Keplerian
        initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,
                                      PositionAngle.TRUE,
                                      inertialFrame, epochDate, mu)
        satellite_mass = 100.0 # kg
        initialState = SpacecraftState(initialOrbit, satellite_mass)
        states.append(initialState)
    print("states defined in Keplerian elements")


#%% Set up propagator


integrator = ClassicalRungeKuttaIntegrator(30.0)
orbitType = OrbitType.CARTESIAN
propagator_num = NumericalPropagator(integrator)
propagator_num.setOrbitType(orbitType)

itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, True) # International Terrestrial Reference Frame, earth fixed
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                         Constants.WGS84_EARTH_FLATTENING,
                         itrf)
gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)
propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))


#%% Propagate one orbit numerically numsat steps forward 

## Create time vector
startDate = AbsoluteDate(2012, 1, 26, 11, 0, 00.000, utc)

step_time = 30

# Time array in orekit AbsoluteDate format
dt = startDate.shiftedBy(float(step_time))

states2 = []

start = time.time()
for initialState in states:
    propagator_num.setInitialState(initialState)
    state = propagator_num.propagate(dt)
    pv_org = state.getPVCoordinates()
    states2.append(state)
end = time.time()
print("total time: ", end-start)
print("per sat time: ", (end-start)/numsat)
# total time:  18.845316886901855
# per sat time:  0.018845316886901855

#%% Propagate numsat orbits numerically forward one step

TestState = states[0]
start = time.time()
propagator_num.setInitialState(TestState)
dt = startDate

propagator_num.setInitialState(TestState)
for i in range(1000):
    dt = dt.shiftedBy(float(step_time))
    state = propagator_num.propagate(dt)
    pv_org = state.getPVCoordinates()
end = time.time()
print("total time: ", end-start)
print("per sat time: ", (end-start)/numsat)

# total time:  0.07098388671875
# per sat time:  7.098388671875e-05

//////////////////////////// Keplarian Prop Only /////////////////////////////

#%%
inertialFrame = FramesFactory.getEME2000()

utc = TimeScalesFactory.getUTC()
initialDate = AbsoluteDate(2004, 1, 1, 23, 30, 00.000, utc) 

mu =  3.986004415e+14

numsat = 1000

keplers = []
for step in range(numsat):
    ra = 814 * 1000 + np.random.normal(0,1000)        # km Apogee!
    rp = 786 * 1000 + np.random.normal(0,1000)         # km Perigee!
    i = radians(98.55 + np.random.normal(0,0.5))      # inclination
    omega = radians(90.0 + np.random.normal(0,0.5))   # perigee argument
    raan = radians(5.1917 + np.random.normal(0,0.5))  # right ascension of ascending node
    lv = radians(0.0567634 + np.random.normal(0,0.5))    # True anomaly
    epochDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)
    
    a = (rp + ra + 2 * ae) / 2.0    # semi major axis in km
    e = 1.0 - (rp + ae) / a
    
    # Orbit construction as Keplerian
    initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,
                                  PositionAngle.TRUE,
                                  inertialFrame, epochDate, mu)

    kepler = KeplerianPropagator(initialOrbit)
    kepler.setSlaveMode()
    keplers.append(kepler)

stepT = 60.
start = time.time()
for step in range(numsat):
    #kepler = keplers[step]
    extrapDate = initialDate.shiftedBy(stepT*(1+step))
    currentState = kepler.propagate(extrapDate)
    #print("step :", step+1)
    #print(" time : ", currentState.getDate())
    #print(" ", currentState.getOrbit())
end = time.time()
print("total time: ", end-start)
print("per sat time: ", (end-start)/numsat)

# total time:  0.008999824523925781
# per sat time:  8.99982452392578e-06'''

Hi @Void

I am wondering if you could not just build as many propagators as you have states and loop on the propagators rather than on the states. This would factor the setInitialState() out of the loop.
I am not sure it will be sufficient if you tinker with the states in between loops as per the UKF as you would be forced to reset something anyway.

Hi @Void,

As we already discussed on the Python wrapper forge you may also look at FieldPropagation.
With a field propagator you can propagate the estimated state of the UKF once to get the predicted state and then use Taylor algebra to get the values of the sigma points around the predicted state.
I haven’t tried it so I don’t know if it’s faster than propagating each sigma point individually.

Maxine,

Thanks again for the initial point in the right direction! I don’t think the field method will be right for this situation. From the API documentation:

…depending on the number of derivatives, the payoff is still very important as soon as we evaluate a few hundreds of points. As Monte-Carlo analyses more often use several thousands of evaluations, the payoff is really interesting.

Since a UKF in a 6 element set will only be propagating 12+1 points for each object at each step, I don’t think I will see a performance boost, but I could be wrong. I still need to dive deeper on the JSON based UKF from UT that you recommended. I suspect their implementation will be quicker than FilterPy’s implementation that uses numpy for its binary back-end.

Luc,

Thanks for the feedback. I had tried that once before and found it slower; however I refactored my code a bit to remove anything from the loops that didn’t need to be there. I think what I did still results in valid propagations & is much closer to the runtime I was gunning for. As a side benefit of this implementation, it leaves open the possibility to multi-thread the process through Orekits own methods.

Do you see any issues using a common variable for the following:

inertialFrame = FramesFactory.getEME2000()
orbitType = OrbitType.CARTESIAN
itrf      = FramesFactory.getITRF(IERSConventions.IERS_2010, True) # International Terrestrial Reference Frame, earth fixed
earth     = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                             Constants.WGS84_EARTH_FLATTENING,
                             itrf)
gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)
satellite_mass = 100.0 # kg
startDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)
step_time = 30.0
dt = startDate.shiftedBy(float(step_time))

For those interested in the full refactored code loop:

//////////////////////// Numerical Multi-Prop Refactor ///////////////////////

#%% Initial osculating orbits

numsat = 240000
states = []

# Inertial frame where the satellite is defined

inertialFrame = FramesFactory.getEME2000()

# shared values for all propagators

orbitType = OrbitType.CARTESIAN
itrf      = FramesFactory.getITRF(IERSConventions.IERS_2010, True) # International Terrestrial Reference Frame, earth fixed
earth     = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                             Constants.WGS84_EARTH_FLATTENING,
                             itrf)
gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)
satellite_mass = 100.0 # kg
startDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)
step_time = 30.0
dt = startDate.shiftedBy(float(step_time))

# generate pool of states
states = []
for i in range(0,numsat):
     # Orbit construction as Keplerian
    p=Vector3D(-1780159.5047552094 + np.random.normal(0,1000),
               -1164835.539234647 + np.random.normal(0,1000), 
               6842941.118987356 + np.random.normal(0,1000))
    v=Vector3D(-7205.19342727098 + np.random.normal(0,1),
               -341.0428975779148 + np.random.normal(0,1),
               -1924.0083602739878 + np.random.normal(0,1))
    initialOrbit = CartesianOrbit(PVCoordinates(p, v), 
                                  inertialFrame, 
                                  epochDate, mu)
    initialState = SpacecraftState(initialOrbit, satellite_mass)
    states.append(initialState)


#%% Set up propagator

start = time.time()
for state in states:
    # generate Orekit Propagator
    propagator_num = NumericalPropagator(ClassicalRungeKuttaIntegrator(step_time))
    propagator_num.setOrbitType(orbitType)
    propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))
    propagator_num.setInitialState(state)
    state = propagator_num.propagate(dt)
end = time.time()
print("total time: ", end-start)
print("per sat time: ", (end-start)/numsat)
# total time:  15.953574895858765
# per sat time:  6.647322873274485e-05

I still welcome any further insights on how I can tighten up the speed here and will make sure to come back and update this as I progress.

Thanks!

Sharing the variables you selected seems good to me.
I even think you could share the HolmesFeatherstoneAttractionModel instance among all propagators, as long as you don’t change the internal parameter driver (your code above is safe in this respect). Sharing force models is in general not safe, but for this specific one, it should work: it does not preserve any state that depends on propagation.

Luc,

Thanks - pulling the force model out of the loop didn’t change the run time enough to notice.

In which situations would sharing the force models become unsafe?

Ash

Some models preserve internal state. The most prominent example is maneuver, with firing status up to 10.0, or detected on/off dates after 10.0 (to be released soon). Other less obvious examples are force model that have parameter drivers (all of them I think) if you happen to reset the value of the driver between propagators.

Maxine,

Thanks for advice to take a look at OrbDetPy. Shiva did some good work there and there is a lot I can learn from her code. That said, it looks like her run-time seems to be similar to what I started with on this thread. I am seeing about 0.012s per SV when propagating 1,000 SVs forward one 30 second time step. The UKF is a little more efficient with 0.042s per SV for one 30 second time step (as opposed to 0.012 x 12 points + matrix solution time), but still quite a bit more than the run times in the modified loops I worked out with Luc. I am going to take a crack at implementing a UKF using FilterPy (uses SciPy instead of Hipparchus as a back-end) and Orekit directly. I will update this thread when I knock that out so others can learn from lessons here. Thanks again for getting me pointed in some useful directions!

Ash