Orbit Determination Sensitivity to Force Model Parametrization

[charlesc]

Hello all!

I was trying to get a rough sense of how sensitive covariances are to different force model parametrisations.

To do this I based my code on the tutorial that was kindly written by @yzokras (great tutorial btw thankyou!) and simply re-run the OD with more and more components in the force model.

With all the force models added (as per the tutorial) I get the same answer as in the tutorial for the position and velocity standard deviations:

Position std: cross-track 9.032e-02 m, along-track 3.718e-01 m, out-of-plane 2.267e-01 m
Velocity std: cross-track 3.581e-04 m/s, along-track 9.807e-05 m/s, out-of-plane 2.741e-04 m/s

Then I removed Drag and SRP and got something different.

Position std: cross-track 3.404e-01 m, along-track 1.784e-01 m, out-of-plane 2.239e-01 m
Velocity std: cross-track 1.483e-04 m/s, along-track 3.386e-04 m/s, out-of-plane 2.763e-04 m/s

So far, so good…

Except I cannot replicate this result. Now, whatever I do I always get the first result. Even if I use Earth Gravity only. This sounded to me like some variable is not being “reset” correctly. I tried restarting the notebook, restarting VSCode but I am completely unable to replicate this result.

UPDATE: As I am writing this I tried again (like 12 times…) with only Earth gravity and got a different answer once more, and then was unable to replicate it again…

Position std: cross-track 6.386e-02 m, along-track 2.629e-01 m, out-of-plane 1.603e-01 m
Velocity std: cross-track 2.532e-04 m/s, along-track 6.934e-05 m/s, out-of-plane 1.938e-04 m/s

I have two questions-

1. Could these “different” results have been some freak numerical error? In which case are OD fits really this insensitive to force model configurations?
2. Is there some kind of variable that is stored in orekit-data or some other semi-permanent variables that I should be manually emptying? I feel like this could be some kind of issue with the Java Virtual Machine.

Here is a link to a the notebook in question (note the first cell contains a lot of helper functions included for completeness): Orekit OD Sensitivity to Force Model · GitHub

Any help much appreciated!

[Serrof]

Hello there

I’m no Jupyter NB expert, but I’ve had experiences with them where some stuff get cached and the kernel needs restart. To remove this unknown completely, I’d advice you to create a .py file that you can run with an IDE or in a terminal. This way you should get reproductible results.

Best,
Romain.

[charlesc]

Hello Romain,

Thank you for your quick response!

I’ve implemented your suggestion of working in a .py file. Additionally, I switched to using Starlink ephemeris data, which helped in simplifying the process by eliminating laser ranging data for now.

To configure various propagatorBuilder options with different force model configurations, I created a function as shown below:

def configure_force_models(propagatorBuilder,cr, cd, cross_section, enable_gravity=True, enable_third_body=True,
                        enable_solar_radiation=True, enable_relativity=True, enable_atmospheric_drag=True):
    # Earth gravity field with degree 64 and order 64
    if enable_gravity:
        itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
        ecef = itrf
        gravityProvider = GravityFieldFactory.getNormalizedProvider(64, 64)
        gravityAttractionModel = HolmesFeatherstoneAttractionModel(ecef, gravityProvider)
        propagatorBuilder.addForceModel(gravityAttractionModel)

    # Moon and Sun perturbations
    if enable_third_body:
        moon = CelestialBodyFactory.getMoon()
        sun = CelestialBodyFactory.getSun()
        moon_3dbodyattraction = ThirdBodyAttraction(moon)
        sun_3dbodyattraction = ThirdBodyAttraction(sun)
        propagatorBuilder.addForceModel(moon_3dbodyattraction)
        propagatorBuilder.addForceModel(sun_3dbodyattraction)

    # Solar radiation pressure
    if enable_solar_radiation:
        itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
        ecef = itrf
        wgs84Ellipsoid = ReferenceEllipsoid.getWgs84(ecef)
        cross_section = float(cross_section)
        cr = float(cr)
        isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(cross_section, cr)
        solarRadiationPressure = SolarRadiationPressure(sun, wgs84Ellipsoid, isotropicRadiationSingleCoeff)
        propagatorBuilder.addForceModel(solarRadiationPressure)

    # Relativity
    if enable_relativity:
        relativity = Relativity(orekit_constants.EIGEN5C_EARTH_MU)
        propagatorBuilder.addForceModel(relativity)

    # Atmospheric drag
    if enable_atmospheric_drag:
        itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
        ecef = itrf
        wgs84Ellipsoid = ReferenceEllipsoid.getWgs84(ecef)
        msafe = MarshallSolarActivityFutureEstimation(
            MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
            MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
        atmosphere = DTM2000(msafe, sun, wgs84Ellipsoid)
        isotropicDrag = IsotropicDrag(float(cross_section), float(cd))
        dragForce = DragForce(atmosphere, isotropicDrag)
        propagatorBuilder.addForceModel(dragForce)

    return propagatorBuilder

Then, I followed more or less the same process as earlier (except using PV instead of Range and AngularAzEl):

def main():
    ephem_path = '/Users/charlesc/Documents/GitHub/ERP_tools/external/ephems/starlink/MEME_57632_STARLINK-30309_3530645_Operational_1387262760_UNCLASSIFIED.txt'
    spacex_ephem_dfwcov = spacex_ephem_to_df_w_cov(ephem_path)

    sat_list = {    
    'STARLINK-30309': {
        'norad_id': 57632,  # For Space-Track TLE queries
        'cospar_id': '2023-122A',  # For laser ranging data queries
        'sic_id': '000',  # For writing in CPF files
        'mass': 800.0, # kg; v2 mini
        'cross_section': 100.0, # m2; TODO: get proper value
        'cd': 1.5, # TODO: compute proper value
        'cr': 1.0  # TODO: compute proper value
    }
    }

    sc_name = 'STARLINK-30309'  # Change the name to select a different satellite in the dict

    odDate = datetime(2023, 12, 19, 6, 45, 42, 00000)
    collectionDuration = 1 * 1/24 * 1/60 * 5 # 5 minutes
    startCollectionDate = odDate + timedelta(days=-collectionDuration)

    #Get TLE for first guess
    # Space-Track
    identity_st = input('Enter SpaceTrack username')
    password_st = getpass.getpass(prompt='Enter SpaceTrack password for account {}'.format(identity_st))
    st = SpaceTrackClient(identity=identity_st, password=password_st)
    rawTle = st.tle(norad_cat_id=sat_list[sc_name]['norad_id'], epoch='<{}'.format(odDate), orderby='epoch desc', limit=1, format='tle')
    tleLine1 = rawTle.split('\n')[0]
    tleLine2 = rawTle.split('\n')[1]
    print(tleLine1)
    print(tleLine2)

    # Orbit propagator parameters
    prop_min_step = 0.001 # s
    prop_max_step = 300.0 # s
    prop_position_error = 5.0 # m

    # Estimator parameters
    estimator_position_scale = 1.0 # m
    estimator_convergence_thres = 1e-2
    estimator_max_iterations = 25
    estimator_max_evaluations = 35
    # tod = FramesFactory.getTOD(IERSConventions.IERS_2010, False) # Taking tidal effects into account when interpolating EOP parameters
    gcrf = FramesFactory.getGCRF()
    # Selecting frames to use for OD
    eci = gcrf
    # utc = TimeScalesFactory.getUTC()
    # mjd_utc_epoch = AbsoluteDate(1858, 11, 17, 0, 0, 0.0, utc)
    orekitTle = TLE(tleLine1, tleLine2)
    itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
    ecef = itrf
    wgs84Ellipsoid = ReferenceEllipsoid.getWgs84(ecef)
    nadirPointing = NadirPointing(eci, wgs84Ellipsoid)
    sgp4Propagator = SGP4(orekitTle, nadirPointing, sat_list[sc_name]['mass'])
    tleInitialState = sgp4Propagator.getInitialState()
    # tleEpoch = tleInitialState.getDate()
    tleOrbit_TEME = tleInitialState.getOrbit()
    tlePV_ECI = tleOrbit_TEME.getPVCoordinates(eci)
    tleOrbit_ECI = CartesianOrbit(tlePV_ECI, eci, wgs84Ellipsoid.getGM())
    integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)
    propagatorBuilder = NumericalPropagatorBuilder(tleOrbit_ECI,
                                                integratorBuilder, PositionAngleType.MEAN, estimator_position_scale)
    propagatorBuilder.setMass(sat_list[sc_name]['mass'])
    propagatorBuilder.setAttitudeProvider(nadirPointing)
    propagatorBuilders = []
    configurations = [
        {'enable_gravity': True, 'enable_third_body': False, 'enable_solar_radiation': False, 'enable_relativity': False, 'enable_atmospheric_drag': False},
        {'enable_gravity': True, 'enable_third_body': True, 'enable_solar_radiation': False, 'enable_relativity': False, 'enable_atmospheric_drag': False},
        {'enable_gravity': True, 'enable_third_body': True, 'enable_solar_radiation': True, 'enable_relativity': False, 'enable_atmospheric_drag': False},
        {'enable_gravity': True, 'enable_third_body': True, 'enable_solar_radiation': True, 'enable_relativity': True, 'enable_atmospheric_drag': False},
        {'enable_gravity': True, 'enable_third_body': True, 'enable_solar_radiation': True, 'enable_relativity': True, 'enable_atmospheric_drag': True},
    ]

    for config in configurations:
        propagatorBuilder = NumericalPropagatorBuilder(tleOrbit_ECI, integratorBuilder, PositionAngleType.MEAN, estimator_position_scale)
        propagatorBuilder.setMass(sat_list[sc_name]['mass'])
        propagatorBuilder.setAttitudeProvider(nadirPointing)
        configured_propagatorBuilder = configure_force_models(propagatorBuilder,sat_list[sc_name]['cr'], sat_list[sc_name]['cd'],
                                                              sat_list[sc_name]['cross_section'], **config)
        propagatorBuilders.append(configured_propagatorBuilder)

    results = {}
    # Outer loop over each propagatorBuilder (/force model configuration)
    for idx, propagatorBuilder in enumerate(propagatorBuilders):
        print(f"Running for configuration {idx}")
        print(f"propagatorBuilder: {propagatorBuilder}")
        print("allforcemodels: ", propagatorBuilder.getAllForceModels())
        matrixDecomposer = QRDecomposer(1e-7)
        optimizer = GaussNewtonOptimizer(matrixDecomposer, False)
        estimator = BatchLSEstimator(optimizer, propagatorBuilder)
        print(f"estimator: {estimator}")
        estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
        estimator.setMaxIterations(estimator_max_iterations)
        estimator.setMaxEvaluations(estimator_max_evaluations)

        # Try different number of points to use to use in OD process
        for points_to_use in range(1, 150, 45):
            print(f"Running for {points_to_use} points with configuration {idx}")

            observableSatellite = ObservableSatellite(0)

            for _, row in spacex_ephem_dfwcov.head(points_to_use).iterrows():
                date = datetime_to_absolutedate((row['UTC']).to_pydatetime())
                position = Vector3D(row['x']*1000, row['y']*1000, row['z']*1000)
                velocity = Vector3D(row['u']*1000, row['v']*1000, row['w']*1000)
                sigmaPosition = row['sigma_pos']
                sigmaVelocity = row['sigma_vel']
                baseWeight = 1.0
                orekitPV = PV(date, position, velocity, sigmaPosition, sigmaVelocity, baseWeight, observableSatellite)
                estimator.addMeasurement(orekitPV)

            estimatedPropagatorArray = estimator.estimate()

            date_start = datetime_to_absolutedate(startCollectionDate).shiftedBy(-86400.0)
            date_end = datetime_to_absolutedate(odDate).shiftedBy(86400.0)

            estimatedPropagator = estimatedPropagatorArray[0]
            estimatedInitialState = estimatedPropagator.getInitialState()
            actualOdDate = estimatedInitialState.getDate()
            estimatedPropagator.resetInitialState(estimatedInitialState)
            estimatedgenerator = estimatedPropagator.getEphemerisGenerator()
            estimatedPropagator.propagate(date_start, date_end)
            bounded_propagator = estimatedgenerator.getGeneratedEphemeris()

            lvlh = LocalOrbitalFrame(eci, LOFType.LVLH, bounded_propagator, 'LVLH')
            covMat_eci_java = estimator.getPhysicalCovariances(1.0e-10)
            eci2lvlh_frozen = eci.getTransformTo(lvlh, actualOdDate).freeze()
            jacobianDoubleArray = JArray_double2D(6, 6)
            eci2lvlh_frozen.getJacobian(CartesianDerivativesFilter.USE_PV, jacobianDoubleArray)
            jacobian = Array2DRowRealMatrix(jacobianDoubleArray)
            covMat_lvlh_java = jacobian.multiply(covMat_eci_java.multiply(jacobian.transpose()))

            # covarianceMat_eci = np.matrix([covMat_eci_java.getRow(iRow)
            #                             for iRow in range(0, covMat_eci_java.getRowDimension())])
            covarianceMat_lvlh = np.matrix([covMat_lvlh_java.getRow(iRow)
                                            for iRow in range(0, covMat_lvlh_java.getRowDimension())])

Here are the covariance matrices that I get when running the 5 different force models with 46 PV in the OD. As you can see there is some change- which is better than what I was getting earlier. I’m now just not sure whether this is a reasonable change to expect or if I may still have the same problem as earlier.

Thanks again for your help and expertise!

Best,
Charles

[charlesc]

Hello @Serrof,

Hope you are well and sorry to disturb you again (I completely understand if you are too busy to answer!).

I have made some progress… In each of my loops, I was not starting the OD process from a fresh state. I discovered this after creating an observer for my BatchLSEstimator and logging orbParams.getDrivers().

To tackle this, I modified my second for-loop to create fresh instances of the estimator each time:

for idx, configured_propagatorBuilder in enumerate(propagatorBuilders):
    for points_to_use in range(60, 130, 60):
        # Reset the initial state
        tleOrbit_ECI = CartesianOrbit(tlePV_ECI, eci, wgs84Ellipsoid.getGM())
        propagatorBuilder = NumericalPropagatorBuilder(
            tleOrbit_ECI,
            integratorBuilder, PositionAngleType.MEAN, estimator_position_scale
        )
        propagatorBuilder.setMass(sat_list[sc_name]['mass'])
        propagatorBuilder.setAttitudeProvider(nadirPointing)
        # Reapply force model configurations
        configured_propagatorBuilder = configure_force_models(
            propagatorBuilder, 
            sat_list[sc_name]['cr'], 
            sat_list[sc_name]['cd'],
            sat_list[sc_name]['cross_section'], 
            **configurations[idx]
        )

        # Reset and configure the estimator
        matrixDecomposer = QRDecomposer(1e-12)
        optimizer = GaussNewtonOptimizer(matrixDecomposer, False)
        estimator = BatchLSEstimator(optimizer, configured_propagatorBuilder)
        estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
        estimator.setMaxIterations(estimator_max_iterations)
        estimator.setMaxEvaluations(estimator_max_evaluations)
        estimator.setObserver(BLSObserver())

I can see that this has worked as the computation of the residuals gets better as I add force models (only adding first and last below for brevity).

These are the corresponding covariance matrices I get, and I feel like they are still too similar - although perhaps I just lack experience to see what is “acceptable”.

For context here is a gist of the code I am using: Investigating OD Sensitivity to Force Model Selection · GitHub

I am quite new to OD, so apologies if any of this is quite naive, and thanks again for any guidance!

Best,
Charles

[Serrof]

Hi there

What’s the difference between configuration 0 and 3 again?
Anyhow you have residuals of the order of 1m no so it’s not bad at all. Do you have info on the propagation model behind your reference ephemerides?

Just a remark, I don’t think your OD actually uses the attitude provider (you have a cannonball model right?) so you would gain in performance by defining an InertialProvider (or FrameAlignedProvider depending on the version of Orekit you use).

Cheers,
Romain.

[charlesc]

Hello

config/force model 0 = 64*64 gravity field
config/force model 3 = 64*64 gravity field + 3BP + SRP + Relativity + Drag

The reference ephemeris is the SpaceX Ephemeris (at https://www.space-track.org/#publicFiles)- I believe these coordinates are provided in MEME (a.k.a. MOD) reference frame as opposed to TEME which the TLEs are typically in. I just realized this detail today and converted the reference ephemeris from MOD to GCRF for consistency with the rest of the script.

            # Add measurements for the current number of points to use
            for _, row in spacex_ephem_dfwcov.head(points_to_use).iterrows():
                # existing code to add measurement
                date = datetime_to_absolutedate((row['UTC']).to_pydatetime())
                position = Vector3D(row['x']*1000, row['y']*1000, row['z']*1000)
                velocity = Vector3D(row['u']*1000, row['v']*1000, row['w']*1000)
                sigmaPosition = row['sigma_pos']
                sigmaVelocity = row['sigma_vel']
                baseWeight = 1.0
                observableSatellite = ObservableSatellite(0)
                orekitPV = PVCoordinates(position, velocity)
                #convert orekit PV from MOD to ECI
                MOD_frame = FramesFactory.getMOD(True)
                eci2mod = MOD_frame.getTransformTo(eci, date)
                eciPV = eci2mod.transformPVCoordinates(orekitPV)
                eci_position = eciPV.getPosition()
                eci_velocity = eciPV.getVelocity()
                PV_measurement = PV(date, eci_position, eci_velocity, sigmaPosition, sigmaVelocity, baseWeight, observableSatellite)
                estimator.addMeasurement(PV_measurement)

And this made my residuals about 10x worse… lol. It gets rid of the strange “jagged” edges I was getting before but degrades the overall fit dramatically (see below)

Thanks for the tip about the InertialProvider ! will give this a go once I manage to address the other issues. Do you have any intuition as to whether my covariance matrices should be evolving more strongly between force models?

Once again thank you very much for taking the time to reply! Truly appreciate it.

Best,
Charles

[charlesc]

Hello,

Just to let you know I ended up implementing my own Batch Least Squares method and the covariances I get from Orekit are of similar order to the ones that I get- so it seems this is the kind of change that one should expect! Still not sure what happened with the MEME to TEME conversion though- but I ended up just doing everything in MEME to get around this.

Anyways thanks for your help.

Best,
Charles

[Serrof]

Hi Charles,

It’s a shame you had to recode a least square process.
Covariance aside, I would expect the residuals to be orders of magnitude smaller with config 3 compared to 0. Is that what you obtained?

Cheers,
Romain.

[charlesc]

Hi Romain,

It was a good learning experience! It fortified my understanding of OD and least squares

The residuals ended up being very similar for all force models over a 20-min observation arc (1 observation per minute). However, config0 takes ~3x the number of iterations to converge than config3. Eventually the goodness of fit just plateaus for both at around the same values.

My impression now is that they will not get better than this because the SpaceX ephemerides seem to contain these discontinuities that I mentioned previously- I presume these are from their own OD fitting process.

However the covariance between config 0 and config 3 was ~1 order of magnitude better. Likewise the distribution of the residuals is much more gaussian in 3 than in 0.

Thanks again!
Charles

[charlesc]

Hello @Serrof !

Sorry to be coming back to this post… I just wanted to understand better what you meant about the InertialProvider. Do you mean that even for a Cannonball model there would be some benefit to including an AttitudeProvider? Or were you suggesting the use of something like the BoxAndSolarArraySpacecraft ?

Thanks again for your help.

Charles

[Serrof]

Hi Charles,

I saw you were using nadir pointing. If the attitude isn’t actually useful in your computation, you’re wasting time computing it. That’s all I meant.

Cheers,
Romain.

[charlesc]

Ah I see. got it. Thanks!

Charles