Batch least square always exceeds the iterations

I am writing a code in Orekit to do orbit determination from the gcrf data of position and velocity data. I have some data that was recorded on orbit and I am also generating measurements in Orekit. however the measurements I generate when I use that data batchleastsquares is not able to converge. keeping everything the same and using the data from on orbit I am able to converge the solution. any help? I am attaching the code
orbit-det.py (10.0 KB)

I am beginner in Orekit. please let me know if I am doing something really wrong and there is better way to do it

Hi @flightdynamicsnoob , welcome

I just skimmed over your code, so my advice may be way off.
It seems to me you use very few and widely separated position-velocity points with huge noise. If I count well, you have less than 15 points, each one separated for its neighbor by 50 minutes which i more than one half orbit of your satellite, 100m of noise on position and 10m/s of noise on velocity (which is inconsistent with the noise on position, it should rather be 0.1m/s to be consistent, using the traditional 1/1000 factor between position and velocity accuracy data in SI units).

So I would suggest to reduce your T sampling rate from 3000s to something much lower (say 300s, or even 60s), to reduce position noise to 10m and to reduce velocity noise to 0.01m/s).

Also in order to generate your data, there are three other ways you could use. The first one would be to replace the loop with a single propagation, using a step handler to sample the dates. The second one would be to keep the ephemeris generator you are already using, but run it only once before the loop and use the loop only to sample points out of it. The third method is much more involved, it would require using the measurements generation feature, but it is also much more features rich as it would allow generating other types of measurements (range, range-rate, angular, turn-around, TDOA, FDOA…) taking into account all physical effects (signal tiume of flight,tropospheric effects…)

2 Likes

hello, thank you of your reply.

I tried with having the sampling rate of the 1 data point a minute and still the solution doesn’t converge.

is there an example of how to use measurements generator to get the PV data from an orbit

@flightdynamicsnoob Correct me if I’m mistaken, but if I’m reading your code correctly, the basic sequence is as follows:

  1. Generate a TLE orbit (in LEO, this is relevant)
  2. Use the TLE orbit to initialize a Keplerian orbit
  3. Create measurements using the Keplerian orbit
  4. Initialize your LS-estimator with a NumericalPropagator

My first guess as to why your estimator is not converging is that you’re feeding it bad data. LEO TLE orbits can deviate from the true orbit values by up to 40 km even when the orbit is stable (TLEs essentially trade away some short-term accuracy in exchange for long-term stability, especially in LEO where drag is such a strong factor). Additionally, drag can have a huge effect on LEO orbits even after only a few days that would cause it to deviate sharply from a pure keplerian orbit (this deviation can vary from 2-3 km for a 500-600 km altitude orbit to 20-30 km for a 300 km altitude orbit).

My recommendation would be to generate your measurements using the same type of propagator that you feed into the LS estimator. Since you are new to Orekit I would recommend you use a simple 2-body Keplerian propagator first, get that working, then work things out from there. For error add something very simple, like 3 m/s position in each direction and 0.05 m/s velocity in each direction. For frame, I would recommend you stick to EME2000 until and unless you are working with TLEs (TLEs are by default in the TEME frame). Also, bear in mind that the conversion between Keplerian elements and TLE values is not straightforward and you will definitely not want to generate a TLE orbit for propagation using the LS-estimator method in Orekit (which is not to say it can’t be done, just that you shouldn’t do it using the LS estimator).

P.S. If you’re not new to orbits I apologize if I over-explained the orbit details. I’d rather include them all and have you not need them than the converse.

1 Like

Also, I see where you set the solar radiation parameter to True in order to make sure the value is estimated by the solver, but did you also set the drag parameter to True? I don’t see that in there, and it would make a huge difference.

1 Like

Isn’t this the correct way to build propagator. I am adding Dragforce to the propagator

isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(CROSS_SECTION, CR)
solarRadiationPressure = SolarRadiationPressure(sun, wgs84ellipsoid, isotropicRadiationSingleCoeff)
propagatorBuilder.addForceModel(solarRadiationPressure)
for driver in solarRadiationPressure.getParametersDrivers():
    if driver.getName() == RadiationSensitive.REFLECTION_COEFFICIENT:
        driver.setSelected(True)
    propagatorBuilder.addForceModel(solarRadiationPressure)

# Atmospheric drag
msafe = MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
                                              MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
atmosphere = NRLMSISE00(msafe, sun, wgs84ellipsoid)
isotropicDrag = IsotropicDrag(CROSS_SECTION, CD)
dragForce = DragForce(atmosphere, isotropicDrag)
propagatorBuilder.addForceModel(dragForce)

yeah, I tried with actually data and works. it was most likely bad data generation. but I am wondering is there any way to get data fro convergence and additional data to get an idea if the estimate method will actually converge?
it tries to minimize the error any way to get data for each time step?

also your answer that I cannot build a TLE orbit from LS method of propagation. If I want to build a accurate orbit and then get TLE I guess this would be the best way to estimate the orbit with all the force models and then convert that orbit to TLE with something like this

estimated_propagator.propagate(FIRST_DATE, LAST_DATE)
tle_builder = TLEPropagatorBuilder(
    orekitTle,
    PositionAngleType.MEAN,
    100.0,
    LeastSquaresTleGenerationAlgorithm(),
)
fitter = FiniteDifferencePropagatorConverter(tle_builder, convergence_threshold, 1000)
fitted=fitter.convert(estimated_propagator, convergence_threshold, 60, "BSTAR")

#fitter.convert(tle_builder,50,26,'BSTAR')  # Setting BSTAR as free parameter
tle_propagator = TLEPropagator.cast_(fitted).getTLE()

It basically looks correct, but you need to do this bit

for driver in solarRadiationPressure.getParametersDrivers():
    if driver.getName() == RadiationSensitive.REFLECTION_COEFFICIENT:
        driver.setSelected(True)
    propagatorBuilder.addForceModel(solarRadiationPressure)

for drag as well as solar. Otherwise the LS estimator won’t realize it’s an estimated parameter and optimize for it, and that will make it harder for the estimator to fit the solution to the raw data.

I’m sorry, I don’t understand your question here. Do you mean get the difference between the LS solution and the measurement data for each time step?