I am a huge beginner in using Orekit and was looking for some help on understanding and getting orbit determination to work.
So far, I’ve used Orekit to simulate a satellite orbit and used an ElevationDetector to get azimuth and elevation (and range) of the satellite as it passes over a ground station. I am simulating FalconSAT3 using a TLE from N2YO.com. I plugged in the TLE into a “SGP4 Propagator” to generate the az/el points.
Now, I want to do a simple test to plug in those simulated azimuth and elevation points into an orbit determination estimator and see how close I can get to the simulated satellite parameters.
I’ve been following the very helpful guide of GorgiAstro on github(sorry I can’t link because I am new user)
and have made small modifications to make this work with AzEl instead of ranges, but when I run, I am getting an “unable to solve” run time error.
My “initial guess” is the same TLE that I used to simulate the orbit, so I expected the estimator to give me back the same/very close parameters.
I’ve attached a text file with the code snippet of my OD attempt. If it is easier, I can attach the beginning part with the code I used to get the simulated data too.
The variables: azArray, elArray, and timeStamps are arrays with the simulated data in
them that I used.
Some things I think might be wrong/I don’t understand
- To add the AngularAzEl measurement to the estimator, I needed to add a “ObservableSatellite” to the constructor. I gave this a propagator index of 0, but I don’t know what it is doing.
- Some objects needed a frame in the constructor. I have been using the variable “inertialFrame”, which is just FramesFactory.getEME2000(). Is it possible that I am using this improperly?
Let me know if you need more information!
odCodeSnippet.txt (3.7 KB)
Welcome to the Orekit community !
Observable satellite = 0 should be alright. It is an index used for multi-satellite orbit determination. If you only have one satellite, 0 is the value to use.
The way you’re defining the EME2000 frame seems good to me too.
I’ve quickly looked through your code and I can’t see any major error so far.
Just two remarks:
- You’re setting the standard deviation of the measurements to 1rad at the beginning of your code. This means you’re assuming a 1σ error of about 60deg on your az/el measurements. It’s a lot. For example in the OD tutorial the 1σ error for az/el measurements is set to 0.02rad.
Now if you’ve used 1rad for the 1σ error of your simulated measurements the error may come from there.
- You inverted sigma and weight in the signature of the AngularAzEl constructor (l. 75). But since since they’re both set to 1 the problem doesn’t come from here.
Without the whole code of your simulation it will be hard to find out where the error comes from. If you’re ok to share it we may be able to help.
Is your error “unable to solve: singular problem” ?
If yes it means the batch least square problem can’t be solved, the exception is most likely thrown at l.430 of the GaussNewtonOptimizer.
This indicates that the Jacobian matrix of the measurements with respect to the initial orbital state is probably ill conditioned.
The two most likely sources of this error are:
- An observability problem: You don’t have enough measurements to solve the problem. Or they’ve all been filtered out, but it doesn’t seem that you added any outlier filter. Check the number of measurements you have.
- A modelling difference between the measurements’ simulation and the orbit determination. It’s very common and can come from any configuration error.
To see if the problem may come from there you can use the debugger and put a break point in l.102 of the MeasurementHandler. From there check the difference between the Estimated measurements and the Observed measurements. If they are too large compared to the error you put in your measurements this indicates a mismatched configuration issue between simulation and OD.
I hope this will help you.
Happy late Thanksgiving! I’m sorry for the late response- I took some time off for the break. Thank you for the information you’ve provided- I appreciate your time. I need to check your second comment, about what is happening in the MeasurementHandler call.
Here is the full code I am using- I did not post it earlier because it was kind of messy/littered with comments.
The basic flow is that I:
- Define frames/earth model
- Create 3 ground stations
- Make ElevationDetector + BooleanDetector wrapper for all 3 stations
- create orbit from TLE (of falconsat-3)
- Propagate for 60,000 [seconds] (a bit arbitrary)
- Get the start and end times of a flyby that is visible by all 3 stations and then get the az and el in 30 second increments and put those in arrays.
Then, I try the OD using those measurements.
I think the simulated data collection be OK… unless I am misunderstanding something.
I changed the sigma to 0.02 rad and still got an “unable to compute hyperbolic eccentric anomaly from the mean anomaly after” error.
I’ll give your second comment a look now!
MatthewsimODfullCode.txt (16.0 KB)
Happy late Thanksgiving to you too!
I’ve just taken a quick look at your code and I think I found a mistake. It was actually in your previous code snippet and I missed it.
Line 285 of your last attached code file you convert the Az/El values in degrees before creating the “AngularAzEl” objects.
They should stay in radians. Orekit always uses SI units as input/output of its methods.
A note on measurements’ creation. You’re using TopocentricFrame.getAzimuth and getElevation methods to create your measurements.
It should be fine but you’ll miss the “light time”, the time the signal takes to go from the satellite to the ground station (and the other way around). For Az/El measurements on a LEO satellite it won’t change much your angles.
But it will introduce differences between simulation and estimation (OD takes light time into account).
If you want to improve that I suggest you take a look at the measurements generation package of Orekit. It takes care of all you need and is located in estimation.measurements.generation. You’'ll need:
You will find an example in the tests, in GeneratorTest.testIssue557.
Hope this will help you. Tell me if the degrees to radians conversion solves your initial issue.
Hmmm I fixed the degree/radian issue just now and took a look at the estimated vs observed measurement in debug mode (just looked at index 0).
For the first iteration, I get:
So these don’t match up at all, but I am a little confused… shouldn’t the observed value be equal to the 1st simulated azimuth value that I feed into the estimator, or at least close to it?
In this case, the first simulated az was: 4.302 rad
Yes it should, if your measurements are sorted in time. Have you checked the date of this measurement?
First measurement should be the first date (or last date) of your batch of measurements, depending on if the propagator is going forward or backward.
I don’t have much time these days so I don’t think I’ll be able to run your code and see what’s going on…
One thing I can tell you too is that you generate the measurements with a TLE propagation and you try to fit the orbit (through the OD) with a numerical propagator builder that only has Newtonian attraction as a force model.
Maybe you should add some force models, like spherical harmonics’ gravity and atmospheric drag to start with.
I’m not sure how good the OD can fit the orbit with AzEl measurements extracted from a TLE propagator…
Maybe someone else on the forum can tell us.
Hey @MaximeJ, no worries- thank you again for the advice you’ve been giving me now.
And yes, you are right- that measurement was the last date, it looks like the propagator is going backwards.
I see… I’ll look into making the model more realistic by adding these.
However, I’m trying to think if there’s a very simple model that’ll work just to verify that I’ve got the basic framework for OD down…
Do you know if there is a more compatible simulated orbit that I can give an “initial guess” to the LSEstimator?
Initially I tried to do this with simulated data from a KeplerianOrbit, but I could not figure out how to input an “initial guess” for this.
Propagator is going backwards because the direction of propagation is optimized in the BatchLSEstimator .
You’re first propagating up to final date to produce your measurements.
Then at line 239 of your code you’re extracting the current state of the propagator (ie at final date).
And this is the estimation date you give to the estimator. So the measurements are before the estimation date, that’s why it is going backward in time.
Produce your measurements from this KeplerianOrbit.
Then use this same KeplerianOrbit as input in the NumericalPropagatorBuilder. That way you should have a basic framework.
Afterwards you can displace the initial guess Keplerian orbit from its initial value to see how your OD manages to converge.
The test package can greatly help you. If you look at the BatchLSEstimatorTest class there is everything you need to set up a working context.
You can also have a look at the Numerical Orbit Determination Tutorial. It is more complete and complex, and uses real measurements data instead of simulated data.
Thank you! I will take a look at these. Hopefully will report back positively soon!
Happy New Years. Sorry for the really late response! I was able to get the basic framework by using the KeplerianOrbit as the input like you said. Thanks a lot, man!