Partial derivatives w.r.t Maneuver objects

Hi Orekit community,

I am searching for the best way to compute the differential effect of the thrust parameter on a given propagation.
When it comes to compute the partial derivatives of a fieldSpacecraftState propagated through a fieldPropagator starting from an initial fieldKeplerianOrbit, the DerivativesStructure methodology is just fine.
However, the thrust maneuver case seems more complicated.
First, I thought that using parameterDriver object can be the relevant tool, but I am not sure to understand how such tool permits to compute the partial derivatives of propagated fieldSpacecraftState w.r.t. the parameters of the thrust.
So, my question is how can be computed the partial derivatives of a propagated fieldSpacecraftState w.r.t a Maneuver-like object (for instance a ContantThrustManeuver object).

Christophe

Hi @ChristopheL,

The way the maneuver models are designed right now, it does not seem possible to do what you want easily.

That is if you are really bound to use a FieldPropagator.
In that case I think you will have to design your own “Fielded” maneuver model.
I’m talking about this a bit more below.

If it is ok for you to use a classical numerical propagator, there are two options:

  1. Propagate the state transition matrix and extract the Jacobian of the maneuver’s parameter drivers. For this you can use a method close to the one used in the batch least-square process, see for example the method AbstractBatchLSModel.fetchEvaluatedMeasurement#409 (note that you don’t need the measurement derivatives part)
  2. Use a SmallManeuverAnalyticalModel (if your maneuver is small) and extract the jacobians directly from the model.

Creating a fielded version of a Maneuver object:

This is a bit tricky since a Maneuver object is a force model and that a ForceModel already uses some fielded methods (addContribution and acceleration).
You’ll have to be very careful on how you handle the fields.
And the down side of the new (since 10.2) maneuvers package is that you will need to field several classes (maneuver, propulsion model and maneuver triggers) instead of just the maneuver class.

But maybe in your specific case you could just create a fielded version of a PropulsionModel, for example FieldConstantThrustPropulsionModel.

This model will implements the PropulsionModel and should have as attributes a RealFieldElement thrust, flowRate and isp and a FieldVector3D direction.

The problem with this is that, when you implement for example the fielded method getAcceleration, there is a conflict between the fielded parameters from the method signature and your fielded attributes.
To avoid this problem each of your attributes should be a Map<Field, RealFieldElement> and when you want to get your thrust attribute in the getAcceleration method you will have to do something like:

// Field from the methode signature (s = FieldSpacecraftState in input of the method)
final Field<T> field = s.getDate().getField();

// Thrust attribute, note that the field used in the propagator and in the model must be the same (which is usually the case)
final T fieldThrust = (T) thrust.get(field);

That way you don’t need to define two different type of fields in your class and you just let the propagator handle the different fields you want to propagate.

Example:

  • Define class, thrust attribute and field used
public class FieldConstantThrustPropulsionModel implements PropulsionModel {

    /** Thrust value (N). */
    private final Map<Field<?>, RealFieldElement<?>> thrust;

[...]

  /** Field used to define all the previous attributes.
     *  By construction it is unique.
     */
    private final Field<?> fieldUsed;
  • Constructor
public <T extends RealFieldElement<T>> FieldConstantThrustPropulsionModel(final T thrust,
                                                                              final T isp,
                                                                              final FieldVector3D<T> direction) {

        // Initialize the fields' maps
        this.thrust    = new HashMap<>();
        this.flowRate  = new HashMap<>();
        this.isp       = new HashMap<>();
        this.direction = new HashMap<>();

        // Extract the field
        this.fieldUsed = thrust.getField();

        this.thrust.put(fieldUsed, thrust);
        this.isp.put(fieldUsed, isp);
        this.flowRate.put(fieldUsed, thrust.negate().divide(isp.multiply(Constants.G0_STANDARD_GRAVITY)));
        this.direction.put(fieldUsed, direction.normalize());
    }
  • Fielded method getAcceleration
    @Override
    public <T extends RealFieldElement<T>> FieldVector3D<T>
        getAcceleration(final FieldSpacecraftState<T> s, final FieldAttitude<T> maneuverAttitude, final T[] parameters) {
        // Field from the method, i.e. from the propagator
        final Field<T> field = s.getDate().getField();

        // compute thrust acceleration in inertial frame
        final T fieldThrust = (T) thrust.get(field);
        final FieldVector3D<T> fieldDirection = (FieldVector3D<T>) direction.get(field);
        return new FieldVector3D<>(s.getMass().reciprocal().multiply(fieldThrust),
                                   maneuverAttitude.getRotation().applyInverseTo(fieldDirection));
    }
  • Then, build a Maneuver force model, assuming you have previously defined a DSFactory (see the FieldPropagation tutorial for more):
            // Define thrust derivative structure
            final double thrust = 10.;
            // Index 0 here means the thrust will be your first derivative parameter
            // factory is a DSFactory you would have built earlier
            final DerivativeStructure thrustDs = factory.variable(0, thrust);
           
            [...]

            // Get the field of the factory
            final Field<DerivativeStructure> field = thrustDs.getField();

            // Fielded propulsion model, only the thrust is derivative parameter here
            final double isp = 300.;
            final Vector3D direction = new Vector3D(1., 0., 0.);
            final double duration = 10.;
            final FieldConstantThrustPropulsionModel propModel =
                            new FieldConstantThrustPropulsionModel(thrustDs, field.getZero().add(isp), new FieldVector3D<>(field, direction));

           // Fielded maneuver, manDate is the maneuver date, frame is the propagation frame
           final AbsoluteDate manDate = ...
           final Frame frame = ...
           final Maneuver maneuver = new Maneuver(new LofOffset(frame, LOFType.TNW), new DateBasedManeuverTriggers(manDate, duration), propModel);
            [...]

            // Then, use this maneuver in a FieldNumericalPropagator

I hope this will work and will avoid you “fielding” the whole maneuver package.
Keep me posted on the results, I’m really interested!

Regards,
Maxime

3 Likes

Hi Maxime,
Your help has been very precious! I managed to compute a propagated fieldOrbit that depend on the free variables so I can compute partial derivatives.
Thanks a lot!
Christophe

You’re very welcome Christophe !!
I’m glad it worked.
Do you think the FieldConstantThrustPropulsionModel would be worth contributing to Orekit?
If yes, would you like to do it?

It will be a very nice feature indeed. Your proposition is a elegant method to insert fielded options in maneuver objects. For instance, I adapted and used it to “field” AttitudeProvider object.
Anyway, when the time is right, it can be cool to contribute.

@MaximeJ and @bcazabonne Thank you for the example you have provided. I am following the steps that you have posted above to implement a FieldConstantThrustPropulsionModel.

When I code the PropulsionModel interface I am required to implementt three more methods, in addition to the fielded version of getAcceleration. I have used the quick fix facility (in Visual Studio Code) to generate generic method stubs. For example,

@Override
public Vector3D getAcceleration(SpacecraftState s, Attitude maneuverAttitude, 
                                double[] parameters) {
    // TODO Auto-generated method stub - do I leave it generic???
    return null;
}

@Override
public double getMassDerivatives(SpacecraftState s, double[] parameters) {
    // TODO Auto-generated method stub - do I leave it generic???
    return 0;
}

@Override
public <T extends CalculusFieldElement<T>> T 
    getMassDerivatives(FieldSpacecraftState<T> s, T[] parameters) {
    // TODO Auto-generated method stub - do I leave it generic???
    return null;
}

So I have a question related to the implementation. Should I leave these stubs as they are or should they implement the proper methods? I haven’t been able to figure it out by myself - I’m still trying to set up the field propagation.

Thank you.

Bogdan

Hi @bogdan.udrea

Tell me if I’m right, but your question was initially for the drag models, not the maneuvers?
In this case, you have to implement fielded versions of DragForce and IsotropicDrag classes.

Bryan

Hello @bcazabonne. My original plan was to start with drag model uncertainty propagation. However, I decided to work through @MaximeJ 's maneuver uncertainty example and I unashamedly copied and pasted the code snippets in this thread. That was the reason why I asked about completing the implementation of the FieldConstantThrustPropulsionModel. I will proceed with the drag model after I am done with the maneuver.

I have answered my own question. Of course, the fielded version of getMassDerivatives is needed. This is my implementation:

    // "fielded" version of getMassDerivatives
    @Override
    public <T extends CalculusFieldElement<T>> T 
        getMassDerivatives(final FieldSpacecraftState<T> s, 
                           final T[] parameters) 
    {
        
        // extracts the field from the propagator
        final Field<T> field = s.getDate().getField();

        // and the flow rate
        final T fieldFlowRate = (T) flowRate.get(field);

        // returns the flow rate that is equal to 

        return fieldFlowRate;

    }

Seems to work good. I’ll update as I progress through the other uncertainties of interest.

Great! Do not hesitate if you have new questions.

Bryan

1 Like

@bcazabonne I have spent some quality time with the IsotropicDrag and DragForce methods and I have concluded that the drag coefficient(dragCoeff) is the only(?) parameter I can consider uncertain and use in field propagation. However, the dragParametersDrivers confuses me a bit. At the first look it seems that I do not need it in a fielded version of IsotropicDrag. However, if I implement my FieldIsotropicDrag as

public class FieldIsotropicDrag implements DragSensitive { ...

I am required to provide a method that returns it, such as

@Override
public List<ParameterDriver> getDragParametersDrivers() { ... 

So I wonder if dragParametersDrivers is needed and I need to field the getDragParameterDrivers method.

Thank you.

Bogdan

Hi,

In the case of IsotropicDrag, the drag coefficient is the only coefficient that you can consider. It exists another parameter (not applicable in you case) related to the drag effect that can be estimated. It is the lift ratio. The lift ratio is the proportion of atmosphere modecules that will experience specular reflection when hitting spacecraft instead of experiencing diffuse reflection. The lift ratio is only used for complex DragSensitive model like BoxAndSolarArraySpacecraft. However, in your case you only have to consider the classical drag coefficient.

Yes it is needed by the DragForce to know the dynamic parameters (i.e., the drag coefficient in your example) of the model. I think you can initialize dragParameterDrivers in your class constructor like that:

this.dragParametersDrivers = new ParameterDriver(DragSensitive.DRAG_COEFFICIENT,
                                                                                   dragCoeff.getValue(), SCALE,
                                                                                  dragCoeffMin, dragCoeffMax);

Thanks you that, the implementation of the getDragParameterDrivers method would be:

/** {@inheritDoc} */
@Override
public List<ParameterDriver> getDragParametersDrivers() {
    return Collections.singletonList(dragParametersDrivers);
}

Could you try and tell me if it works?

Best regards,
Bryan

Thank you @bcazabonne . I am getting back to this during the long weekend and I’ll post updates and for sure additional questions.
Bogdan

Hello @bcazabonne. Unfortunately I haven’t been able to get my own implementation of FieldIsotropicDrag to work.

I have a couple of observations:

  1. You suggested using dragCoeff.getValue() in the call to ParameterDriver. However, getValue() is not available so I have changed to getReal()
  2. If I use dragCoeff = parameters[0], then element 0 of the derivative structure stores the nominal value and all the others are zero. If I use final T fieldDragCoeff = (T) dragCoeff.get(field), then element 0 stores the nominal value and element 36 stores a value of 1.

The simulation with FieldIsotropicDrag gives the same results as the simulation with IsotropicDrag and that makes me believe that the uncertainty in the dragCoeff does not contribute to the propagation.

I’m pasting below most of the code. I would appreciate some advice - I think I’m missing something fundamental in organizing fielded versions of force models.

Thank you and best regards.

Bogdan

public class FieldIsotropicDrag implements DragSensitive
{

    /** Parameters scaling factor.
     * <p>
     * We use a power of 2 to avoid numeric noise introduction
     * in the multiplications/divisions sequences.
     * </p>
     */
    private final double SCALE = FastMath.scalb(1.0, -3);

    /** Drivers for drag coefficient parameter */
    private final ParameterDriver dragParametersDrivers;

    /** NB fielded parameters */
    // private final Map<Field<?>, CalculusFieldElement<?>> fieldedHalf;
    private final double crossSection;   // m2
    private final Map<Field<?>, CalculusFieldElement<?>> dragCoeff;

    // this field is used to define all the fielded attributes
    private final Field<?> fieldUsed;

    /** Constructor with drag coefficient min/max set to ±∞.
     * @param crossSection Surface (m²)
     * @param dragCoeff drag coefficient
     */
    public <T extends CalculusFieldElement<T>>
        FieldIsotropicDrag(final double crossSection, final T dragCoeff) 
    {

        // initializes the maps of each field
        this.dragCoeff = new HashMap<>();

        // extracts the field in current use
        this.fieldUsed = dragCoeff.getField();
            
        // assigns values to the fielded parameters
        this.dragCoeff.put(fieldUsed, dragCoeff);
        this.crossSection = crossSection;

        // initializes the drag parameters driver
        // bcazabonne suggested using dragCoeff.getValue() but there exists no
        // such method for (generic) type T
        this.dragParametersDrivers = new
            ParameterDriver(DragSensitive.DRAG_COEFFICIENT, dragCoeff.getReal(), 
            SCALE, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);

    }

    
    /**
     * @param date
     * @param frame
     * @param position
     * @param rotation
     * @param mass
     * @param density
     * @param relativeVelocity
     * @param parameters
     * @return FieldVector3D<T>
     */
    @Override
    public <T extends CalculusFieldElement<T>> FieldVector3D<T> 
        dragAcceleration(final FieldAbsoluteDate<T> date, final Frame frame,
        final FieldVector3D<T> position, final FieldRotation<T> rotation,
        final T mass, final T density, final FieldVector3D<T> relativeVelocity, 
        T[] parameters) 
    {

        // extracts the field from the date
        final Field<T> field = date.getField();

        /** debug debug */
        // drag coefficient manipulation
        // final T fieldDragCoeff = parameters[0];
        final T fieldDragCoeff = (T) dragCoeff.get(field);

        FieldVector3D<T> dragAcc = new 
            FieldVector3D<>(relativeVelocity.getNorm().
                                multiply(density.
                                multiply(fieldDragCoeff).
                                multiply(crossSection / 2)).
                                divide(mass),
                            relativeVelocity);

        return dragAcc;

    }

    /** 
     * @param date
     * @param frame
     * @param position
     * @param rotation
     * @param mass
     * @param density
     * @param relativeVelocity
     * @param parameters
     * @return Vector3D dragAcceleration
    */
    @Override
    public Vector3D dragAcceleration(final AbsoluteDate date, 
        final Frame frame, final Vector3D position, 
        final Rotation rotation, final double mass,
        final double density, final Vector3D relativeVelocity,
        double[] parameters) 
    {

        // TODO Auto-generated method stub
        return null;

    }


    
    /** 
     * @return List<ParameterDriver>
     */
    @Override
    public List<ParameterDriver> getDragParametersDrivers() {
        return Collections.singletonList(dragParametersDrivers);
    }
    
}

Hi @bogdan.udrea

My apologies

final T fieldDragCoeff = (T) dragCoeff.get(field) is the good solution.

I tired your class in a gradient propagation context. It gives interesting results. Indeed, the partial derivatives of the orbital parameters with respect to the grad coefficient are not equals to zero. I used a very low orbit (i.e., about 450.0 km of altitude). So the drag effect has a significant impact. Please find attached the example.


public class GradientPropagation {

    /** Private constructor for utility class. */
    private GradientPropagation() {
        // empty
    }

    /** Program entry point.
     * @param args program arguments (unused here)
     */
    public static void main(final String[] args) {

        // configure Orekit
        final File home       = new File(System.getProperty("user.home"));
        final File orekitData = new File(home, "orekit-data");
        if (!orekitData.exists()) {
            System.err.format(Locale.US, "Failed to find %s folder%n",
                              orekitData.getAbsolutePath());
            System.err.format(Locale.US, "You need to download %s from %s, unzip it in %s and rename it 'orekit-data' for this tutorial to work%n",
                              "orekit-data-master.zip", "https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip",
                              home.getAbsolutePath());
            System.exit(1);
        }
        final DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
        manager.addProvider(new DirectoryCrawler(orekitData));

        //nominal values of the Orbital parameters
        final double a    = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 450.0e3;
        final double e    = 1e-3;
        final double i    = FastMath.toRadians(98.3);
        final double pa   = 0.5 * FastMath.PI;
        final double raan = 0.0;
        final double ni   = 0.0;

        //gradient
        final int freeParameters = 7;
        final Gradient aG    = Gradient.variable(freeParameters, 0, a);
        final Gradient eG    = Gradient.variable(freeParameters, 1, e);
        final Gradient iG    = Gradient.variable(freeParameters, 2, i);
        final Gradient paG   = Gradient.variable(freeParameters, 3, pa);
        final Gradient raanG = Gradient.variable(freeParameters, 4, raan);
        final Gradient niG   = Gradient.variable(freeParameters, 5, ni);
        final Gradient cdG   = Gradient.variable(freeParameters, 6, 1.0);

        //sometimes we will need the field of the Gradient to build new instances
        final Field<Gradient> field = aG.getField();

        //initializing the FieldAbsoluteDate with only the field it will generate the day J2000
        final FieldAbsoluteDate<Gradient> date = new FieldAbsoluteDate<>(field);

        //initialize a basic frame
        final Frame frame = FramesFactory.getEME2000();

        //initialize the orbit
        final Gradient mu = Gradient.constant(freeParameters, 3.9860047e14);

        final FieldKeplerianOrbit<Gradient> KO = new FieldKeplerianOrbit<>(aG, eG, iG, paG, raanG, niG, PositionAngle.ECCENTRIC, frame, date, mu);

        //force models
        final ForceModel fModel_Sun  = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
        final ForceModel fModel_Moon = new ThirdBodyAttraction(CelestialBodyFactory.getMoon());
        final ForceModel fModel_HFAM =
                        new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
                                                              GravityFieldFactory.getNormalizedProvider(18, 18));
        final ForceModel fModel_Drag = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
                                                                        new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                                                                                             Constants.WGS84_EARTH_FLATTENING,
                                                                                             FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
                                                     new FieldIsotropicDrag(1.0, cdG));

        //setting an hipparchus field integrator
        final OrbitType type = OrbitType.CARTESIAN;
        final double[][] tolerance = NumericalPropagator.tolerances(0.001, KO.toOrbit(), type);
        final AdaptiveStepsizeFieldIntegrator<Gradient> integrator =
                        new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);

        integrator.setInitialStepSize(60);

        final FieldSpacecraftState<Gradient> state0 = new FieldSpacecraftState<>(KO);
        final FieldNumericalPropagator<Gradient> numProp = new FieldNumericalPropagator<>(field, integrator);

        numProp.setOrbitType(type);
        numProp.setInitialState(state0);
        numProp.addForceModel(fModel_Sun);
        numProp.addForceModel(fModel_Moon);
        numProp.addForceModel(fModel_HFAM);
        numProp.addForceModel(fModel_Drag);

        final FieldSpacecraftState<Gradient> finalState = numProp.propagate(date.shiftedBy(3600.0));

        System.out.println(finalState.getDate());
    }

}

Thank you again @bcazabonne .I like the example you provided. My orbit analysis is for a lower orbit than in your example so I definitely see the influence of the uncertainty in the drag coefficient.

After a few hours of debugging =D I reached the conclusion that
final T fieldDragCoeff = (T) dragCoeff.get(field)
is the way to go. So you have confirmed my guess. I believe passing the dragCoeff.getReal() through parameters strips away the information in the DerivativeStructure.

Regards,

Bogdan

Because the .getReal() method is only used in the class constructor to initialize the reference value of the ParameterDriver, I think there is no problem. The important point is to never use the .getReal() method with the dragCoeff parameter inside the dragAcceleration() method because you will loose the derivatives information.

Bryan

1 Like