Orbit Determination with N-body attraction

1

Hello everyone!

I have embarked on a journey to perform Orbit Determination outside of the Earth’s Sphere of Influence.

In such a regime, we would need to consider the gravity of Earth, Moon, Sun, Jupiter, and an N number of bodies based on the use-case. An orbit in such a regime is generally non-Keplerian, non-planar, and not centered about a fixed barycenter (hinting to the inadecuacy of CR3BP limited realism for such a task in cislunar space, for example, to which @luc has eluded to here!.

The BLSEstimator I use requires a PropagatorBuilder. The underlying propagator that is built needs to have (inspired from the CR3BP tutorials):

  • numericalPropagator.setOrbitType(null), to clarify that the orbit is non-Keplerian
  • numericalPropagator.setIgnoreCentralAttraction(true), to clarify that there is no single central attraction point
  • numericalPropagator.addForceModel(...) with my choice of n-body accelerations from selected planetary bodies, similarly to @sebastien.herbiniere suggestion!

However, I do not have access to the first two methods from the NumericalPropagatorBuilder object. I noticed that the builder, inside the buildPropagator method, attempts to conditionally add a Newtonian attraction by itself from the input Orbit mu, if one such force is not present.

Adding the solar bodies as perturbation force models in an Earth-Centered initial orbit is probably not equivalent to what I would like to achieve, but yet again I am not sure now!

Question: What recommendations would you have to perform Orbit Determination within an n-body attraction environment, for a non-Keplerian orbit with known initial state guess in Cartesian coordinates (say in J2000?) ? I feel that if BatchLSEstimator had a NumericalPropagator argument instead of NumericalPropagatorBuilder, my problem would be solved, is that right or no? Alternatively, if the NumericalPropagatorBuilder could buildPropagator and keep a reference of the result in its instance as a member, I would also be able to reference it and use the setters on-the-fly (sounds too sketchy perhaps?).

Thank you!
Manny

P.S. Can I just say how happy I am for Orekit 12.0? :blush:

NewtonianAcceleration in BatchLSEtimator
2

Hi @Manny

That’s an interesting challenge :blush:

You can find a discussion on this topic here: Orbit determination for interplanetary missions

To summarize, I think that Orekit needs some additional improvements to be able to perform orbit determination for interplanetary missions.
For orbit propagation, everything is OK. But for orbit determination we need further developments that are summarized in the topic above.

In addition, if remember well, I think the NewtonianAttraction is added automatically to the NumericalPropagatorBuilder. So we need to be able to avoid it for orbit determination of interplanetary missions.
An easy hack is to remove the acceleration of the Newtonian attraction by adding a ForceModel performimg the inverse acceleration. That’s ugly, but needed with the current status of the orbit determination.

I’m very interested on that topic!

Best regards,
Bryan

1 Like
3

Bryan, Thank you so much for such a speedy reply! :pray: @bcazabonne

Your linked topic is excellent and will proceed to dig into it now!
The reverse force model is simple, genius and out-of-the-box thinking! I now have fresh weekend homework!

I’ll report any progress I make with this hack below. Thank you once again!

1 Like
4

Hi everyone,

I wanted to check if there have been any updates on the topic of interplanetary orbit determination. I couldn’t find any open issues related to it on the GitLab page, and I’m hoping it hasn’t gotten lost among other tasks.

Thanks in advance!

Best,
Davide

5

Hi @DDega!

There is no update on that topic.

6

Thanks for the update,
Davide

7

Hi everyone,

I am also stuck with an OD for interplanetary missions. As @Manny, I also had to fight a bit with the fact that the NumericalPropagatorBuilder requires an Orbit object.

It is not very elegant (sorry for this @bcazabonne…) I know, but I tried to put the central attraction coefficient equal to a very low number (0.00000001 m^3/s^2) and the propagated positions I get are very close to those I get with a NumericalPropagator with the same force models (cm-level differences, that for my case are ok). However, I am getting some strange results and so I started doubting about my “solution”.

First, I am considering three satellites around different Lagrangian points in the Solar System. The orbits (and times) of these three satellites are assumed to be known (I disabled the corresponding drivers). Then I have another satellite (the user) that needs to estimate its orbit from one-way signals sent by the three reference satellites. After I created the four NumericalPropagatorBuilder objects with my trick, I created the estimator as

builder_array = [propagator_builder_user, propagator_builder_se, propagator_builder_sm, propagator_builder_sv]
optimizer = LevenbergMarquardtOptimizer()
estimator = BatchLSEstimator(optimizer, builder_array)

and then I started creating the InterSatellitesRange measurements. For the first measurement within the first satellite, I wrote

meas_se = InterSatellitesRange(obs_user, obs_se, False, date, 0.0, float(sigma_range_se),      float(baseWeight)) 
theoretical_eval = meas_se.estimate(0, 0, [user_state, se_state])
range_orekit = theoretical_eval.getEstimatedValue()[0]
meas_se.setObservedValue(float(range_orekit  + np.random.normal(0, sigma_range_se))) 
estimator.addMeasurement(meas_se)

where user_state and se_state are the true states of user and reference satellite at the date of the reception of the one-way signal at the user.

If I do the same for many measurements with all the three satellites (I tested the code with a total of 20 measurements for a total time span of 10 hours) then if I use a BatchLSObserver I see that the residuals are already pretty high at the very beginning in the first iteration and evaluation (km-level residuals, while I assumed a perfect orbit knowledge for the user and I set its clock error to 0 s). As I was not sure using the estimate method in that way was reasonable, I implemented my own function to compute the value of the one-way range to be provided to the InterSatellitesRange object. My function does an iterative propagation to understand where the reference satellite was at the time of the signal transmission. I think it is well implemented and the residuals I get through the BatchLSObserver are much lower than those I get with the value I estimated with the estimate method of Orekit. In any case, from the second evaluation on, the residuals clearly diverge in both cases and there is no way I can have a good OD estimation.

What do you think about this? I have been struggling a lot and then I found this topic, so maybe we tackled similar problems..

All the best,
Samuele

P.S. I obtain the same even with a very low sigma_range_se value (like 1 mm).

1 Like
8

Hi,

There might be hidden side effects of using a small mu. You’d better add a force counteracting the Newtonian gravity like Bryan suggested.

Cheers,
Romain.

2 Likes
9

Hi!

Yes, I already experienced the anti-Newtoninan attraction for mission analyses and it works well. So, currently that’s the only solution without a deep refactoring of the OD…

Best regards,

Bryan

2 Likes
10

Just a side note here. This is going to change, see issue 1867 about this. It is however a huge modification (currently 491 files modified and 16000 lines changed in the library) and it is not ready for production, so the changes are only on a local branch on my computer. I hope to continue this work soon but am really busy these days.

3 Likes