I was initially generating some test measurements using SPICE and trying to fit them using Orekit. I found that it wouldn’t converge, even when given the correct answer as a guess. As such, I went back and rewrote the propagation and measurement generation code to use Orekit also. Unfortunately, Orekit still fails to converge, even with the parameters given in most examples — 20 m and 0.001 m/s sigma on two-way range and range-rate, respectively.
The initial state is:
t0 = 2022-06-16T21:44:43.373
x0 = [-6.45306258e+06, -1.19390257e+06, -8.56858164e+04,
1.83609046e+03, -9.56878337e+03 -4.95077925e+03] in J2000
If I use that initial state, I get an exception: “org.orekit.errors.OrekitException: minimal step size (1.00E-03) reached, integration needs 7.43E-04”. Even lowering the min step, I get that exception. If I use the initial state for the W3B.aer example, it iterates forever without converging.
I’m only using a single newtonian attractor, with mu = 398600435436095.9.
Looking at your initial state, this behavior is expected. Indeed, the altitude of the perigee of your orbit is at 184 kilometers (with an eccentricity of 0.97). At this altitude, the dynamic is very strong and the numerical integrator needs a very low step to support this dynamic.
If I good understand, your are using the WB3 initial state (orbit and date) with your generated measurements ? If yes, that is normal if it iterates forever without converging because you are using measurements that does not correspond to the correct satellite. Moreover, you can encountered date issues because W3B data (measurements and initial state) are for the year 2010 and your data are given for the year 2022.
Looking at your initial orbit, I think you have to add some orbit perturbations. First, the drag because at the orbit’s perigee you are very low in the atmosphere. Moreover, at the orbit’s apogee you are at an altitude of 381052 kilometers (very very close to the Moon). Therefore, you have to also add the third body perturbation of the Moon. I don’t know if it will solve your problem because at the apogee, I think the influence of the Moon est higher than the influence of the Earth for your satellite. That’s a reason for what the convergence is difficult with your example because your it is difficult to model.
How did you determine the initial state? Maybe you can try to use the IODLaplace or IODGooding methods. These methods are built to give you an initial guess for the orbit determination by using only three measurements.
I generated the initial state using the patched conic strategy, and while it’s intended to eventually rendezvous with the moon, I’ve propagated/generated measurements with only a single Newtonian attractor — and thus, that’s the model I’m trying to fit. Yes, down the road it needs more complicated models, but for now I just want to get it working with a simple case.
Next I’ll add the third body attractor, the goal being to do a study of the effect of measurement error and measurement frequency on state knowledge when we arrive at the moon. I don’t expect that this study will be impacted much by leaving out drag.
Your note about the dynamics at low altitudes gave me the clue I needed to solve this. It looks like the easiest way to do it is to propagate through the low altitude portion of the trajectory with a large dP just once and start fitting measurements after that.
I was looking into your issue and I think I found the problem.
I used your initial state and got the same exception as you with a DP853 integrator.
Now try exactly the same thing but with a comma between the 2nd and 3rd component of the velocity
Replace: 1.83609046e+03, -9.56878337e+03 -4.95077925e+03
With: 1.83609046e+03, -9.56878337e+03, -4.95077925e+03
You’ll see it works much better, even with 1.00E-03 as minimum step size.
The problem I bumped in was that when building velocity without the comma like this: velocity = Vector3D(1.83609046e+03, -9.56878337e+03 -4.95077925e+03)
it gets into the constructor Vector3D(alpha, delta) instead of Vector3D(x, y, z); with a right ascension alpha = 1.83609046e+03 radians and delta = -9.56878337e+03 -4.95077925e+03 ~ -14519.56262 radians.
And this ultimately gives you a velocity vector of [0.1082544756, 0.6289622177, 0.7698619989]
In other words, you’re free-falling towards the Earth ground instead of going to the Moon
