I’m trying to figure out the basic setup in OREKIT to do a maneuver calibration using ScaledConstantThrustPropulsionModel.

The scenario would be having sc state1 and state2 with the maneuver (or several) in between, with maneuver start/end times and durations known, and looking to determine the dv-rate (or total dv) of the maneuver.

" You could use orbit determination, by using two position-velocity measurements (one corresponding
to your state before and one corresponding to your state after), not estimating orbit itself (i.e. changing the orbital parameter drivers to not selected ), and estimate the maneuvers parameters. This corresponds in fact to maneuver calibration."

I’m having a hard time seeing how to implement this in Orekit without any examples, could someone expand upon this somewhat? I’ve seen this tutorial on orbit determination, but not clear at all how one would ‘invert’ this to do a maneuver calibration.

Also I’m working in python, and have implemented many examples of using pre-defined maneuvers, but not ones to be estimated or calibrated.

I am working right now (since yesterday), on something concerning partial derivatives of maneuvers.
This feature will probably be useful for you.
I expect some progress and primitive implementation (without any bells and whistles, and most probably with lots of rough edges) in the upcoming days.
Stay tuned, I will report about my progress on the forum.

Just a note, the second link you gave us points to an “Orekit Labs” page that is 2 years old and probably outdated.
You should better look into the official Orekit tutorials and more specifically the estimation package and the NumericalOrbitDetermination tutorial.

Now, to complement @luc’s old post that you linked, what you could try is:
(Preliminary note: the proposition below will work only if your maneuver date is set and you don’t want the OD to estimate it. If you want to estimate both date and thrust, then you will have to wait for Luc’s work to be finished (see his previous post))

Get started with the tutorial linked just above to see how to set up an OD with Orekit

Select the three ParameterDriver of your ScaledConstantThrustPropulsionModel (use ParameterDriver.setSelected(true) when looping on ScaledConstantThrustPropulsionModel.getParametersDrivers)

Add your Maneuver class with a ScaledConstantThrustPropulsionModel to the NumericalPropagatorBuilder that you set up for the OD

Deselect the orbital drivers in the propagator builder (use ParameterDriver.setSelected(false) when looping on NumericalPropagatorBuilder.getOrbitalParametersDrivers). This will prevent orbit estimation.

Use two PV measurements. The first one is your initial orbit, the second one your target orbit.

Run the OD

Tell us if it worked we are very interested in this too!

Thanks very much to you both! @MaximeJ this is exactly what I need to get the ball rolling - and to confirm, I’m trying to keep it simple initially and just do the ‘OD’ on the maneuver thrust, and will take the maneuver times as knowns in this equation.

Have been able to make some progress on a simple-as-possible example of this thrust estimation, but seem to have hit a wall where the estimator returns an ‘unable to solve: singular problem’ at the estimator.estimate() step.

Have attached the code sample in python, along with the output of a run, below.
The basic setup is as you’ve described above, and I’ve added a test to confirm which parameters in the estimator are left to be actually estimated, and they are just the maneuver thrust x, y, z scale factors.

To test a trivial case, here I’ve used actual pos/vel vectors for a satellite that has NOT undergone a maneuver at the stated time, to estimate a zero-thrust maneuver, but get the same ‘singular problem’ no matter what maneuver times I use.

It seems like I should be putting in an initial guess for the thrust, where the ScaledConstantThrustPropulsionModel scales this up or down to do the estimation - but could not find a place to fit that in.

I’ve also tried relaxing the convergence requirements, etc. to no avail. Any ideas?

The issue was including zeros in specifying the vector for the thrust. This also seems to be that ‘initial guess’ I was looking for, which I was specifying as a unit-vector in the +velocity direction at the time of firing, however it doesn’t like the zero values for the y and z coordinates. Setting these to a small, non-zero value provides estimated parameters!
eg. `thrust_vec = Vector3D(-1.0, 0.001, 0.001)

However this brought me to a new wall: This does not scale to multiple maneuvers between states. Is there a fundamental reason why the estimator could not handle this, ie estimating x, y, z scale factors for three maneuvers specified exactly as this single one is in the example?

The attempt brings me back to the singular problem at the estimation step. I have actual data for this so I don’t think that my inputs are too wild to try and fit within the estimator. But what I do not have is the measurements between each maneuver.

Nice that you solved it by yourself and shared the outcome!

I was about to tell you the same thing about the zeros. The problem with the ScaledConstantThrustPropulsionModel is that it cannot always explore new directions by itself as the scale factors are applied on the initial thrust vector. If you start with a thrust purely on X axis then it cannot try to add some thrust on the Y and Z axis.

The solution, as you found it, is either to estimate only the scale factor on X axis or use an initial thrust vector with components on the other axis.

Another thing, we may have mislead you on something with Luc… I think you should just add the target orbit as PV measurement, and not the initial orbit.
Indeed the initial orbit’s derivatives with respect to the different thrusts will always be zero.
So it’s basically useless and can cause a ‘Singular Matrix’ exception since the corresponding line of the Jacobian will always be a line of zeros.

That may come from adding the initial state but I’m not 100% sure.
This method is a bit of a hack, I’m not sure a maneuvers’ sequence should be optimized that way, i.e. using the orbit determination scheme.
The problem you may face is that you only have six “measurements” components (your target position/velocity vector) to optimize 3 x N unknowns, where N is the number of maneuvers you want to determine.
So as soon as you have more than 2 maneuvers your problem will be under-determined and you’ll step into a case where problems like ‘singular matrix’ can occur.
In some cases where the “minimum cost” of the problem is straightforward and the system well conditioned it will work; but it may fail in other cases.
You can add some conditioning to your system by setting min/max values of the thrust scales, and thus reducing the space of exploration of the algorithm.

Could you maybe share your new code ? So we can try it and see if we can make it work somehow.

Sure thing, attached is a working version where you can toggle estimation with a single burn, or three. The actual states and maneuver times I’ve put in correspond to an actual three-burn test, so if it’s going to work at all, it should work on this data. You can configure the 1 vs 3 burns in the scenario at line 170: test_single_burn = True/False

I’ve also removed the addition of the initial state as a ‘measurement’, as you suggest, and it gets rid of the singular error. Now the error we get with multiple maneuvers is:

org.orekit.errors.OrekitException:
unable to compute hyperbolic eccentric anomaly from the mean anomaly after 50 iterations

The number ‘50’ does not exist in the code so this smells like a deeper issue such as the under-determination. Also, I agree with the sentiment that this is a bit of a hack to estimate maneuver performance in this way. The issue is definitely lack of measurements between each burn, so looking to see how close we can come to something that works!

I’ve had another head-scratcher all along here to do with specifying the min/max deviation in the ScaledConstantThrustPropulsionModel construction:
How does one specify a min AND a max? The documentation just seems to list the input fields as a double for thrust, ISP, so it’s not clear what exactly this single (?) value is meant to represent (the max, min, absolute values, ratios, etc). Could you clarify?

Many thanks once again! Could not have made it this far without your help

Thanks for sharing your code. I won’t be able to look into it today but I’ll try to do that later.

This means that the propagated orbit is going hyperbolic… so there’s definitely an issue here.
It must try to “explore” the domain with huge maneuvers and pushes a “little bit too much”. So indeed it may come from the under-determination.
The number “50” is an internal maximum for the iterative loop that converts the mean anomaly to the hyperbolic eccentric anomaly. It’s hard-coded in Orekit…

I’ve had another head-scratcher all along here to do with specifying the min/max deviation in the ScaledConstantThrustPropulsionModel construction:
How does one specify a min AND a max? The documentation just seems to list the input fields as a double for thrust, ISP, so it’s not clear what exactly this single (?) value is meant to represent (the max, min, absolute values, ratios, etc). Could you clarify?

An ISP in seconds (note that this ISP will never change during your estimation process)

A thrust direction (unit vector) d = [u, v, w] along the maneuver’s attitude axis X, Y, Z

A name for the maneuver

It builds a PropulsionModel object that will contain an initial thrust vector: \overrightarrow{T_{0}} = T.\overrightarrow{d}
Or \overrightarrow{T_{0}} = [Tx_0, Ty_0, Tz_0] with Tx_0 = T.u, Ty_0 = T.v, Tz_0 = T.w

It also creates a set of three scale factors [sx, sy, sz] initialized to 1.
When calling the getThrustVector method it outputs: \overrightarrow{T} = [sx.Tx_0, sy.Ty_0, sz.Tz_0]

And then, during the estimation or optimization process, it will try to find the best [sx, sy,sz] that fit your problem.

So several things:

If you have Ty_0 = 0 or Tz_0 = 0 then sy and sz will never have an impact on the optimization process. This was the case when you set the direction of the thrust purely on the X axis

When you build the ScaledConstantThrustPropulsionModel you should use the initial guess of the thrust and the ISP of your thruster. I noted on your code that perhaps you think it is the thrust and ISP “deviations” instead. And note that with this model the ISP is always constant.

As said before, the thrust vector is \overrightarrow{T} = [sx.Tx_0, sy.Ty_0, sz.Tz_0].
The ParameterDriver's are the scale factors sx, sy and sz that will be optimized by the batch least-squares estimator.
They are the ones you can set a min/max value to to constrain your problem.
Example: for the 1st maneuver, if you want your thrust on X to be +/-50% of the initial thrust then do: man_param_lst[0].getParametersDrivers().get(0).setMinValue(0.5)
and man_param_lst[0].getParametersDrivers().get(0).setMaxValue(1.5)

One more thanks for the excellent explanation, @MaximeJ! It makes perfect sense.

It also shows my non-Java-native background, not having found those setMinValue() methods myself (which are members of the ParameterDriver class)!

In summary, after playing around with all of this for multiple maneuvers between measurements, and even adding more post-maneuver states to help convergence, it seems like it’s just an intractable problem - however does work fine for a single maneuver and not under-determining the problem

The behavior we get with multiple maneuvers in the most stable case is that the estimator will stop on the extremes of the thrust-range for the first burn, also the second, and not touch them for the third. As though it tries to force the solution via the first burn only, then ‘finds’ the later burns and doesn’t know what to do with them. Understandable.
I also tested just estimating the x-component only, similar results.

eg. for scale factors limits on x,y,z of thrust as follows, where x is the main thrust direction and the y,z are minor off pointings:

xs_min = 0.8
xs_max = 1.1
yzs_min = 0.0 # used for y and z
yzs_max = 2.0 # used for y and z

… for three burns we get estimated scale factors of: