Failed to get DSSTBatchLSModel/DSSTPropagator to produce correct solutions

I created a new configuration file and AbstractOrbitDetermination subclass, following the pattern for the Lageo2 test, but modified to use input ephemeris (PV) data.

I find that using the DSSTBatchLSModel against PV data does not seem to work correctly. Specifically,
It converges to a solution without difficulty, but a solution with large residuals. Even if I use only one data point a few seconds after the initial guess from the config file, the residual is unexpectedly large.

The DSSTPropagator is doing something I don’t understand, because with a PVT in the .in file used to initialize the estimated orbit, and a PVT data value a few seconds later, I would expect very small residuals for that data value (a few mm), but the residuals are surpringly large (100s of m and 10s of m/sec). This is true whether or not drag, radiation pression, sun/moon third body forces, etc., are turned on. Residuals are much larger (100s of km) when hours of data are used.

For example, I am providing all data in EME2000 coordinates and UTC (data comes from precision Grace ephemeris. As provided, the Grace ephemeris is in ITRF coordinates and modified GPS seconds - I am pretty sure that all PVT conversions to EME2000 and UTC are done correctly).

More notes:

mu is correct. Switching from central body degree/order 12 to 70 does not significantly change the results.

Using iers.conventions=2010. Switching from WGS-84 to CIO/2010 doesn’t change the results significantly.

Latest test using propagator.min.step = 0.1 propagator.max.step = 100.0 propagator.position.error = 0.0001 but using other step sizes does not siginicantly change the results.

Modifying DSSTBatchLSModel from

parallelizer.propagate(firstDate.shiftedBy(-1.0), lastDate.shiftedBy(+1.0));

to

parallelizer.propagate(firstDate.shiftedBy(-1.0e-3), lastDate.shiftedBy(+1.0e-3));

does not significantly change the results (but why was a shift of 1 sec chosen?? If the propagation region in this code line needs to strictly include the requested time range, shouldn’t the extensions be configurable? What if propagating backwards in time is impossible for some reason?).

Bottom line - something seems wrong with the physics model, but I am unable to determine what. I am aware that because I don’t have perfect knowledge of atmospheric drag, solar radiation forces, vehicle attitude, etc., the orbital predictions cannot be perfect. Nevertheless, I would expect accuracy in the mm range when I am propagating forward only 5 seconds. I am not achieving that. Suggestions please and thanks in advance!

Are you running DSST in mean or osculating mode?

Excellent question - thanks! I did not think of that. DId more tests. Here are results with 2 data points (for a total of three highly accurate ephemeris points at 5 sec intervals, including the initializer).

PropType=Mean StateType=Mean
Resids = 253,80,80,80 converges w/ large residual
w/either solver

PropType=Mean StateType=OSC
NPE immediately in DSSTTesseral$HansenObjects.computeHansenObjectsInitValues (see below)
w/either solver

PropType=Osc StateType=Mean
Resids = 28154,24149,16417,5916,3297,144,0.35,0.31,0.31 w/Levenberg-Marquart
Resids = 28154,329,0.46,0.31,0.31 w/Gauss-Newton

PropType=Osc StateType=Osc
Resids = 139,140,147,218,etc diverges slowly w/ Levenberg-Marquardt
Resids = 139,140,0.39,0.30,0.30 converges w/ Gauss-Newton

Emperically, it appears that I should use PropType=Osc, StateType=Osc. I presume that “osculating” in this context means “actual”. I am still surprised by

  1. Initial RMS error=139m. Where is that coming from? When my initial guess is highly accurate (2 cm, according to JPL), I wonder why propagating the correct initial data 5 seconds and 10 seconds results in errors of 139 m RMS. It makes me think that the initial PV is interpreted as mean, rather than actual. I don’t see a parameter in ParameterKey to specify that initial orbit PV parameters are actual rather than mean. I will dig into this further.

  2. After convergence, the RMS position error is 30cm, which is still much larger than JPL’s published error. Turning on friction, 3rd body gravity, solar pressure, and adjusting the corss-section terms made next to no difference, as I would expect for a scenario involving only 10 sec of motion.

  3. Why does the combination PropType=MEAN StateType=OSC blow up (NPE)?

Thanks in advance-

=====================================================================
java.lang.NullPointerException
at org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral$HansenObjects.computeHansenObjectsInitValues(DSSTTesseral.java:2798)
at org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral$UAnddU.(DSSTTesseral.java:2378)
at org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral.getMeanElementRate(DSSTTesseral.java:438)
at org.orekit.propagation.semianalytical.dsst.DSSTPropagator$Main.elementRates(DSSTPropagator.java:994)
at org.orekit.propagation.semianalytical.dsst.DSSTPropagator$Main.computeDerivatives(DSSTPropagator.java:973)
at org.orekit.propagation.integration.AbstractIntegratedPropagator$ConvertedMainStateEquations.computeDerivatives(AbstractIntegratedPropagator.java:677)
at org.hipparchus.ode.ExpandableODE.computeDerivatives(ExpandableODE.java:134)
at org.hipparchus.ode.AbstractIntegrator.computeDerivatives(AbstractIntegrator.java:265)
at org.hipparchus.ode.AbstractIntegrator.initIntegration(AbstractIntegrator.java:217)
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:196)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:468)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:410)
at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:141)
at org.orekit.estimation.leastsquares.DSSTBatchLSModel.value(DSSTBatchLSModel.java:264)
at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:616)
at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(GaussNewtonOptimizer.java:399)
at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:436)
at com.aac.common.AbstractOrbitDetermination.runBLS(AbstractOrbitDetermination.java:372)
at com.aac.common.AbstractOrbitDetermination.runBLS(AbstractOrbitDetermination.java:221)
at com.aac.Deter.testGrace(Deter.java:245)
at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
at java.base/jdk.internal.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
at java.base/java.lang.reflect.Method.invoke(Method.java:566)
at org.junit.runners.model.FrameworkMethod$1.runReflectiveCall(FrameworkMethod.java:50)
at org.junit.internal.runners.model.ReflectiveCallable.run(ReflectiveCallable.java:12)
at org.junit.runners.model.FrameworkMethod.invokeExplosively(FrameworkMethod.java:47)
at org.junit.internal.runners.statements.InvokeMethod.evaluate(InvokeMethod.java:17)
at org.junit.runners.ParentRunner.runLeaf(ParentRunner.java:325)
at org.junit.runners.BlockJUnit4ClassRunner.runChild(BlockJUnit4ClassRunner.java:78)
at org.junit.runners.BlockJUnit4ClassRunner.runChild(BlockJUnit4ClassRunner.java:57)
at org.junit.runners.ParentRunner$3.run(ParentRunner.java:290)
at org.junit.runners.ParentRunner$1.schedule(ParentRunner.java:71)
at org.junit.runners.ParentRunner.runChildren(ParentRunner.java:288)
at org.junit.runners.ParentRunner.access$000(ParentRunner.java:58)
at org.junit.runners.ParentRunner$2.evaluate(ParentRunner.java:268)
at org.junit.runners.ParentRunner.run(ParentRunner.java:363)
at org.junit.runner.JUnitCore.run(JUnitCore.java:137)
at com.intellij.junit4.JUnit4IdeaTestRunner.startRunnerWithArgs(JUnit4IdeaTestRunner.java:68)
at com.intellij.rt.junit.IdeaTestRunner$Repeater.startRunnerWithArgs(IdeaTestRunner.java:33)
at com.intellij.rt.junit.JUnitStarter.prepareStreamsAndStart(JUnitStarter.java:230)
at com.intellij.rt.junit.JUnitStarter.main(JUnitStarter.java:58)

Another data point - when I reduced the gravity model from 70/70 to 12/12, the initial RMS error for OSC/OSC was reduced to 1.05m. This is more like what I was expecting.

To get 70/70, I used the latest high order Grace-derived gravity model ITSG-Grace2018s.gfc. Why would using a high order gravity model worsen the DSST accuracy?

It may be related to the configuration of the short periodic terms. The relationshio between short period configuration and dynamic model is difficult (and I certainly don’t have this expertise).

Still have a few unresolved issues, but I am making progress. Sticking with the OSC/OSC case in DSST, i observed the following:

  1. When I fit posted accurate Grace ephemeris positions at 5 sec intervals to a spline (over periods of 10s of minutes), and then differentiate the splines, the resulting derivative matches the ephemeris velocity to better than 1E-5 m/sec.

  2. When I take the output of DSST and do the same thing, the numerical derivatives of DSST positions (computed at the same 5 sec intervals) differ from the DSST velocities by several cm RMS.

I conclude from this that the DSST outputs have positions and velocities that are inconsistent, well beyond what can be explained by the accelerations.

Is there a POC you can recommend who is an expert on DSST who might be able to help with this?

Thanks in advance-

Jim

That is an interesting trail. We already encountered a similar problem with another propagator using a mean model: Eckstein-Hechler. We had to implement an ugly workaround because the velocities were basically not the derivatives of the positions.

DSST works in equinoctial elements internally, so position/velocities are not its canonical output. I guess you extract them from the SpacecraftState. This object always hold the orbit just as the propagator provided it, i.e. as an EquinoctialOrbit in DSST case. So when one asks for position/velocity, as it is not what is already stored, a conversion is performed upon calling. Normally, it should take into account the first time derivatives of the elements, it does not only rely on Keplerian state, but when a model uses both mean elements and short periodics, perhaps something was forgotten there. Could you set up a minimalist junit test case that shows the discrepancy between positions derivatives and velocities and open a bug report with this test case?

Yes, @paulcefo who is the world best expert on DSST, which is is masterpiece. Paul did not do the Orekit port but helped us when we did it; he is also a member of the Orekit PMC (Project Management Committee).