Jacobian computation and ephemeris generator

Hello,

I try to combine the use of STM and jacobian computation with the EphemerisGenerator during a numerical propagation.
With only a STM, it seems to work well but if I add parameters and the computation of the jacobian of the final orbital parameters w.r.t. these parameters, an exception is raised when I use the BoundedPropagator built during the propagation.

The code is called from Matlab but I don’t think it’s the origin of the problem.

Java exception occurred:
java.lang.NullPointerException

at org.orekit.time.AbsoluteDate.isEqualTo(AbsoluteDate.java:1237)

at org.orekit.time.AbsoluteDate.isAfterOrEqualTo(AbsoluteDate.java:1288)

at org.orekit.forces.maneuvers.jacobians.TriggerDate.init(TriggerDate.java:204)

at org.orekit.propagation.AbstractPropagator.initializeAdditionalStates(AbstractPropagator.java:248)

at
org.orekit.propagation.analytical.AbstractAnalyticalPropagator.propagate(AbstractAnalyticalPropagator.java:126)

at org.orekit.propagation.AbstractPropagator.propagate(AbstractPropagator.java:276)

In my propagation, I add constant thrust maneuvers and I want the jacobian w.r.t. the median date and the duration of the maneuvers. This jacobian matrix is used for an optimization but I also want to use this function to compute the covariance at any time of the propagation (by using the STM and the jacobian) and that’s why I use the EphemerisGenerator.
Here is an extract of the code where I tested it (the exception is raised on the last line):

for idx = 1:numel(manMedianDate)
name = “man” + idx + “_”;
fireDate = epoch.shiftedBy((manMedianDate(idx)-0.5*manDuration(idx)) * 86400);
duration = manDuration(idx) * 86400;
constantThrustPropulsionModel = BasicConstantThrustPropulsionModel(simData.thrust, simData.isp, manDirection(idx), name);
maneuver = ConstantThrustManeuver( …
propagator.getAttitudeProvider(), …
DateBasedManeuverTriggers(name, fireDate, duration), …
constantThrustPropulsionModel);
maneuver.getParameterDriver(name+“_MEDIAN”).setSelected(true);
maneuver.getParameterDriver(name+“_DURATION”).setSelected(true);
propagator.addForceModel(maneuver);
end

% add initial state
state = SpacecraftState(EquinoctialOrbit(initialOrbit), simData.satelliteMass);
propagator.setInitialState(state);

% propagate
harvester = propagator.setupMatricesComputation(“STM”, , );
generator = propagator.getEphemerisGenerator();

finalState = propagator.propagate(startDate, endDate);

ephemeris = generator.getGeneratedEphemeris();
intermediateState = ephemeris.propagate(midDate);

Did anyone already try to combine the EphemerisGenerator with the matrices computation and particularly the jacobian ?

Thank you!

Chris

Hi Chris,

I believe your problem is similar to what is discussed here.
In your case, a workaround is to use the Field version of propagation, more precisely with the Gradient class. You’ll be able to retrieve in your ephemerides the first order partial derivatives w.r.t. the parameters you defined.

Romain.

1 Like

Hi @Christophe,

I have good and bad news :sweat_smile: .

The good news is that we are aware of this issue (issue 949) and the bad news is that it still need to be adressed…

Cheers,
Vincent

You beat me to it @Serrof !

Hello Romain,

Thanks for the advice.
But in fact, I don’t see how to declare a parameter of a force model (for example the thrust of a ConstantThrustManeuver) as part of the gradient I want to compute.
For the initial orbital parameters, I know how to buiId my FieldOrbit but I don’t see how to declare the thrust for example as the 7th free parameter of my gradient, in order to be able to build the 6x7 jacobian.
Is there a way to declare a parameter of a force model as a Gradient object or do I misunderstand something in the usage of the Gradient object ?

Thanks!

Christophe

Hi @Christophe,

Force models parameters aren’t “fielded”, unfortunately.
You will have the same issue with the drag or SRP coefficient.

This can be resolved in two different manners:

  1. The hard way: we manage to find an architecture allowing ParameterDriver to be fielded
  2. The simpler but verbose way: we duplicate some models (implementations of PropulsionModel, DragSensitive, RadiationSensitive…) by adding a “fielded” version of each of them.

Given that the ParameterDriver handling architecture is already quite difficult (all the more since the TimeSpanMap logic was added for v12), I would tend to prefer solution (2).
This solution also has the advantage that we can add the models one by one, on users’ demand.

In your case, you could open an issue asking for a fielded version of the PropulsionModel?
If we opt to implement solution (1), then your issue will be solved in the process.

Maxime

There should be an ugly workaround to do locally, by overwritting the Field version of the acceleration method in a new force class extending ConstantThrustManeuver. In order to redefine there the independent variables for the automatic differentiation. This would imply an unchecked cast to Gradient, but would work for that particular case.
Obviously not something we want for an official release of Orekit, but that would do the trick for Christophe I think.

Edit: since you’re using Python, the type is handled a bit under the hood and I’m not even sure an IDE would complain about safety.

Cheers,
Romain.