Negative CD estimated in Drag force, during orbit determination

Orekit Python Wrapper usage
1

Hi all,
I have this problem with the CD estimate in the orbit determination (with GPS measurement).
In particular I got a negative value:

print(“Estimated Cd:”, driver.getValue())
Estimated Cd: -39.578737103873756

Can anyone help me solve this problem? I mean the CD should be in a range of 1.0 - 3.0

I set the CD estimate like this:

Atmospheric drag

msafe = MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES, MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)

atmosphere = NRLMSISE00(msafe, sun, wgs84Ellipsoid)
drag = DragForce(atmosphere, IsotropicDrag(cross_section, cd))

cdEstimated = True

for driver in drag.getParametersDrivers():
if driver.getName() == DragSensitive.DRAG_COEFFICIENT:
if cdEstimated:
driver.setSelected(True)

propagatorBuilder.addForceModel(drag)

Many thanks for your time and help.
Kind Regards.
Alessandro

2

Heya! You can modify your code to the following:

if cdEstimated:
    driver.setSelected(True)
    driver.setMinValue(1.0)
    driver.setMaxValue(3.0)

Which will constrain the Cd to 1 to 3. If you’re getting very large errors, it sounds like your atmospheric density isn’t matching - try the Strengthlevel.WEAK instead?

3

Some questions:

  • How many GPS measurements do you have? What length of time do they cover?
  • What mass and cross-sectional area (in m^2) have you put in for your spacecraft?
  • What sort of propagator are you using? If you’re using the NRLMSISE model, you should be using the numeric propagator.
  • What other forces have you added to your propagator?
  • What is the semi-major axis and eccentricity for your satellite (i.e. what orbital regime is it in)?

-39 is an absurd number - simply limiting the Cd range will stop it happening, but it won’t cure the real problem. Clearly something else is going on here.

4

Hi all,
thank you so much for your suggestions.
here are the answers:

  1. more or less 360 measurements in a 30 minute interval.
  2. 60-70 kg with a cross section close to about 0.8 m^2.
  3. the propagator is numerical.
  4. the forces considered are SRP, Sun/Moon body, relativity, NSG up to 64.64, tides, drag.
  5. the orbit is a SSO at about 600 km.

I noticed that I was missing propagatorBuilder.setMass(satellite_mass) (before it was the default value of 1000 kg).

Now I get values ​​that are about - 1.7 ; -2.2, so still negative, but more consistent (in absolute value).

I think it could be something about the STD of the measurement or the estimator setting: I am using those values:

STD

sigma_position = 5.0 # Noise (in terms of standard deviation of gaussian distribution) of input position data in [m]
sigma_velocity = 0.05 # Noise of input velocity data in [m/s]
weight_position = 1.0 #1/sigma_position**2 # Base weight

Estimator parameters

estimator_position_scale = 1.0 # m
estimator_convergence_thres = 1e-7
estimator_max_iterations = 200
estimator_max_evaluations = 200

Orbit estimator parameters

prop_min_step = 0.001 # s
prop_max_step = 60.0 # s
prop_position_error = 0.1 # m

tstep_OD = 1.0 # s
tprop_OD = 3600 * 1 # s

What do you think? Do you have any further suggestions?
Thank you for your availability and attention.
Regards

5

Most of this sounds good. The only things I can think of that might help are:

  1. initializing the NRLM_CSSI model with the SpaceWeather file (can be found in the orekit-data folder) and the proper switches:
CSSI_WEATHER_MODEL = CssiSpaceWeatherData('SpaceWeather-v1.2.txt')
NRLMSISE_CSSI_MODEL = NRLMSISE00(
            CSSI_WEATHER_MODEL,
            CelestialBodyFactory.getSun(),
            EARTH_BODY
        ).withSwitch(9, -1).withSwitch(1, 1)

You can find an explanation for what the switches mean here. Remember to run the update.sh executable first.

  1. Setting all the weights for the measurement points to 1. I’m not saying this is the reason for the odd Cd value, but generally I’ve always worked with data sets where all the points were equally weighted and it’s the only thing I can think of offhand that might make things a little wonky for you.

Try updating the atmosphere model first, then adjusting your weight values if that doesn’t help.

6

One other thing I can think of is really hard for me to know if it’s right or not because it’s dependent on how you’re flying your s/c. Essentially if you have a satellite that would have a significantly different cross-sectional area depending on its orientation, and you’re choosing an average value that is significantly different from what the answer would be for the particular data set you’re using, that might be the cause of the error.

7

Hi there

Your position error in your integrator is too big. You need 0.001 or even less.

Cheers,
Romain.

1 Like