Covariance propagation with modeled maneuvers

Hi all,

I’m trying to implement a covariance propagation model considering a modelled (known and error-free) maneuver.
The model is based on the following paper offering different methodologies to assert the covariance propagation with maneuvers.

On a first approach, as mentionned above, the maneuver would be considered as perfect, with no errors, and its effect can be reflected on the state transition matrix at the maneuver triggering time (see page 2 of above paper, method A).

Before proceeding with the implementation, I would like to make sure such applicative layer is not already available in Orekit.
Indeed, an unique impulsive maneuver, could be considered as an external force model influencing the covariance state transition matrix. Correct me if I’m wrong, but if so, then a propagator setup with a STM and maneuver could actually provide a Covariance matrix influenced by the maneuver, right ?

Similarly, could we consider the same approach for a finite maneuver scenario implementing a basic constant thrust propulsion model ?

In future development, with uncertainties added to the maneuver model, we might need to consider a custom implementation though…

Thanks in advance for your time and assistance.
Best regards,
Gueorguy Serafimov

Hi @Gueorguy,

Orekit already supports maneuvers’ estimation (thrust, orientation and dates), see for example the tutorial ManeuversEstimation by GC.
(The OD tutorials are a bit hard to understand because most of the work is done in AbstractOrbitDetermination with some customization for each type of propagator.)

To do so, there are some internal mechanisms computing the STM, so it’s actually already doable to propagate the covariance and have the STM or covariance of the maneuvers’ parameters.
Unfortunately (but that should change soon) you cannot directly use the StateCovarianceMatrixProvider since for now it only propagates the orbital covariance.
You’ll have to do it the “old way”, look for example how it’s done in method KalmanModel.getErrorStateTransitionMatrix.
You’ll need to:

  • First define the maneuver using for example a ConstantThrustManeuver.
  • “Select” the maneuver ParameterDrivers that you want to have in the STM:
    maneuver.getParameterDriver(driverName).setSelected(true)
  • Then setup your propagator with a matrices’ “harvester”
    MatricesHarvester harvester = propagator.setupMatricesComputation("stm", null, null); (assuming initial STM is null)
  • Then, during propagation, assuming you have a propagated state a time t, you can get:
    • the STM: RealMatrix dYdY0 = harvester.getStateTransitionMatrix(state);
    • the Jacobian of the propagation maneuvers: RealMatrix dYdPp = harvesters.getParametersJacobian(state);
      This last one contains the derivatives of the state a t with respect to the maneuver’s parameters selected in the second step. I’m assuming this is what you want ?

Yes, it’s done in the Kalman filter’s “prediction” step: P(t) = \Phi(t,t_0) * P(t_0) * \Phi^T(t,t_0) + Q(t)
With:

  • P(t_0), P(t) the covariances a t0 and t
  • \Phi(t, t_0) the STM + Jacobian of parameters between t0 and t
  • Q(t) the process noise matrix (can be hard to define)

That’s the one I actually provided. The ImpulseManeuver in Orekit is an EventDetector (not a ForceModel) so the mechanism above won’t work. There’s a work around though, which is to consider a very small constant thrust maneuver (1 second should work well).

For this I think you should try to define a proper process noise model (as mentioned above). There’s a basic class for this in Orekit called UnivariateProcessNoise. Maybe you will have to improve it.

Hope this helps.
Best regards,
Maxime

1 Like

Hi Maxime,

Thanks a lot for the detailed answer ! It helps a lot.

Just to make sure I fully understand the methodology on the Covariance computation:
The State Transition Matrix required to compute the Covariance matrix at the maneuver time tm is a sum of the propagation STM and the Jacobian matrix of the propagated matrix.
In this case scenario I still neglect the process noise matrix. We obtain the following equation:

And similarly, If I wish to compute the finale state Covariance, then I should apply the same formula at tf > tm, with the STM matrix and Jacobian paremeters extracted at tf.

A question came out from this:
Should we consider a loss in the Covariance computation model accuracy in case of long finite maneuvers ?
“Long” is a term yet to be defined. Thus, we might want to virtually discretize the maneuver into maximum duration thrust sections, perform intermediate Covariance computations, and reset the STM within the propagator at each step.

If we consider n sections, with 0 <i<n, we would obtain something like this at each computation steps:


And the propagation final state Covariance matrix would be obtained as followed:

If I’m not wring, this approach would also be necessary in case the process noise matrix is considered.

Thank you in advance for your time and support on this.
Best regards,
Gueorguy

Hi @Gueorguy,

You’re welcome !

Just to be clear it’s not a “sum” as in a matrix sum but more a matrix concatenation, see the comments at the beginning of KalmanModel.getErrorStateTransitionMatrix method for a proper equation.
Sorry if my first answer was misleading since I used the “+” sign.

Yes

I would tend to define this loss of accuracy with a process noise matrix instead and keep the linear model of the STM on the other hand.

I’m not sure to follow.
Know that the STM is a first order linearization and that you have:
\Phi(t_f,t_0) = \Phi(t_f, t_i) * \Phi(t_i, t_0) for any ti in [t0, tf]. So P_f = \Phi(t_f,t_0)*P_0*\Phi^T(t_f,t_0)=\Phi(t_f,t_i)*\Phi(t_i,t_0)*P_0*\Phi^T(t_i,t_0)*\Phi^T(t_f,t_i)=\Phi(t_f,t_i)*P_i*\Phi^T(t_f,t_i)
So stopping the propagation to get the covariance at each computation step won’t change the final value, unless you change the covariance at each ti to account for the loss of accuracy.
That being said, maybe I didn’t understand your point.

Best regards,
Maxime

Hi @MaximeJ

Thanks for the clarification on the STM computation.
In case several force models are implemented and assigned to the STM, I assume the harvester would provide the concatenated propagation jacobian parameters of each implemented force models.

Indeed you are right, this decomposition makes sense only when the process noise matrix is considered (covariance loss of accuracy).

Best regards,
Gueorguy

hi @MaximeJ,

I’m reopening this thread because of an issue encountered when reimplimentin the Error State Transition matrix computation model from the KalmanFilter class.

In line 547, the propagation parameters jacobian named “dYdPp” is retreived via the the harvester.getParametersJacobian(currentSpacecraftState) method. When implementing a simple numerical propagation scenario to test the feature, the retrieved matrix is null. The reason for it lies within the getJacobiansColumnsName() where an empty columnNames List object is created. If I understand it well, this object is supposed to list all the driver parameters names of the implemented force models.
This same object is present and populated at NumericalPropagator level but is unfortunatly not visible. The KalmanFilter class is retrieving the parameters from PropagatorBuilder instances used for orbit determination.

    // define initial state
    Frame inertialFrame = FramesFactory.getGCRF();
    AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
    double mu = Constants.IERS2010_EARTH_MU;
    Vector3D position = new Vector3D(7200000.0, 0.0, 7200000.0);
    Vector3D velocity = new Vector3D(0.0, -5763.39, 5763.39);
    PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
    CartesianOrbit initialOrbit = new CartesianOrbit(pvCoordinates, inertialFrame, date, mu);
    SpacecraftState state = new SpacecraftState(initialOrbit);
    
    // define propagator
    ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(1.0);
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setInitialState(state);

    // setup force models
    // Atmospheric Drag
    CelestialBody sun = CelestialBodyFactory.getSun();
    double ae = Constants.IERS2010_EARTH_EQUATORIAL_RADIUS;
    double f = Constants.IERS2010_EARTH_FLATTENING;
    Frame bodyFrame = FramesFactory.getGTOD(IERSConventions.IERS_2010, false);
    OneAxisEllipsoid earth = new OneAxisEllipsoid(ae, f, bodyFrame);
    Atmosphere atmosphere = new HarrisPriester(sun, earth);
    double crossSection = 2;
    double dragCoeff = 0.47; //for spherical shape. 
    IsotropicDrag scModel = new IsotropicDrag(crossSection, dragCoeff);
    ForceModel dragModel = new DragForce(atmosphere, scModel);
    propagator.addForceModel(dragModel);

    // SRP
    double dRef = Constants.JPL_SSD_ASTRONOMICAL_UNIT;
    double pRef = 4.56e-6;
    double cr = 0.5;
    RadiationSensitive spacecraft = new IsotropicRadiationSingleCoefficient(crossSection, cr);
    ForceModel srpModel = new SolarRadiationPressure(dRef, pRef, sun, earth, spacecraft);
    propagator.addForceModel(srpModel);

    // setup STM computation model
    MatricesHarvester harvester = (AbstractMatricesHarvester) propagator.setupMatricesComputation("stm", null, null);
    double[][] cartesianCovarianceMatrix = {{1.000000e+02, 1.000000e-02, 1.000000e-02, 1.000000e-04, 1.000000e-04, 1.000000e-04},
                                           {1.000000e-02, 1.000000e+02, 1.000000e-02, 1.000000e-04, 1.000000e-04, 1.000000e-04},
                                           {1.000000e-02, 1.000000e-02, 1.000000e+02, 1.000000e-04, 1.000000e-04, 1.000000e-04},
                                           {1.000000e-04, 1.000000e-04, 1.000000e-04, 1.000000e-04, 1.000000e-06, 1.000000e-06},
                                           {1.000000e-04, 1.000000e-04, 1.000000e-04, 1.000000e-06, 1.000000e-04, 1.000000e-06}, 
                                           {1.000000e-04, 1.000000e-04, 1.000000e-04, 1.000000e-06, 1.000000e-06, 1.000000e-042},};
    RealMatrix matrix = new Array2DRowRealMatrix(cartesianCovarianceMatrix);
    StateCovariance stateCovariance = new StateCovariance(matrix, date, inertialFrame, OrbitType.CARTESIAN, PositionAngleType.TRUE);
    StateCovarianceMatrixProvider stateCovarianceMatrixProvider = new StateCovarianceMatrixProvider("covariance", "stm", harvester, stateCovariance);
    propagator.addAdditionalStateProvider(stateCovarianceMatrixProvider);

    // propagate
    AbsoluteDate finalEpoch = date.shiftedBy(100);
    SpacecraftState finalState = propagator.propagate(finalEpoch);

    // retrieve STM
    RealMatrix dYdPp = harvester.getParametersJacobian(finalState);
    RealMatrix dYdY0 = harvester.getStateTransitionMatrix(finalState);
    System.out.println("Parameters Jacobian" + dYdPp);
    System.out.println("STM" + dYdY0);

Do you believe this same feature is available for orbit propagation ? If not available yet, I consider implementing a custom one relying on the getJacobiansColumnsName() script implemented within the NumericalPropagator class.

Thanks in advance for your time and assistance.

Best regards,
Gueorguy

Hi @Gueorguy,

Not sure I fully understand your question again, sorry…

Did you select the parameter drivers that you want to appear in the Jacobian? (ParameterDriver.setSelected(true)).
It is needed if you want a non-empty parameters’ Jacobian.
I see that you’re using the StateCovarianceMatrixProvider class. Beware that it only manages orbital parameters for now. There’s a feature issue open to extend it with propagation and measurement parameters (see here).
Contributions are welcome :wink: ! (unless @bcazabonne is already working on this).

Regards,
Maxime

Hi @MaximeJ,

Thanks for the quick feedback.

Indeed I’m trying to go around the current limitation with orbital covariance propagation. To do so, I need to implement a custom orbital state + propagation parameters covariance propagation.
The shared scenario was build to rapidly test the parameters jacobian generation feature.

Thanks for adding the missing brick in the code. It helps a lot.
Will also have a look at current opened issue on the topic

Best regards,
Gueorguy

1 Like