EKF implementation Orekit_jpype

Orekit Python Wrapper usage
1

Hello all,

I’m trying to implement an EKF on the orekit_jpype version but I’m not sure if my poor results (position off by around 1000 km and velocity off by around 100 m/s) are due to bad covariance initialization or something else.

I have my observer set up as follows:

@JImplements(KalmanObserver)
class MyKalmanObserver:
    def __init__(self):
        self.step_count = 0

    @JOverride
    def evaluationPerformed(self, estimation):
        self.step_count += 1
        print(f"Estimation step {self.step_count} performed")
        
        # Get the corrected measurement+
        corrected_measurement = estimation.getCorrectedMeasurement()
        # print(f"Corrected Measurement value: {observed_value}")
        if corrected_measurement is not None:
            observed_value = corrected_measurement.getObservedValue()
            print(f"Measured (observed) value: {observed_value}")
        
        
        # Get the corrected spacecraft states
        corrected_states = estimation.getCorrectedSpacecraftStates()
        for state in corrected_states:
            estimated_position = state.getPVCoordinates().getPosition()
            estimated_velocity = state.getPVCoordinates().getVelocity()
            print(f"Estimated position: {estimated_position}")
            print(f"Estimated velocity: {estimated_velocity}")

        # Get the estimated state
        process_estimate = estimation.getEstimate()
        estimated_state = process_estimate.getState()
        
        print(f"Estimated state: {estimated_state}")

        # Covariance
        estimated_covariance = process_estimate.getCovariance()
        print(f"Estimated covariance: {estimated_covariance}")

And I have my EKF set up as follows:

initialOrbit = CartesianOrbit(pv_coordinates_eme, eme2000_frame, initialDate, mu)

min_step = 1.0e-6
max_step = 600.0
position_error = 1e4
velocity_error = 1e6

integrator_builder = DormandPrince853IntegratorBuilder(min_step, max_step, position_error)

propagator_builder = NumericalPropagatorBuilder(initialOrbit, integrator_builder, PositionAngleType.TRUE, 1.0)

coef_p = 1e3
coef_v = 1e6
covariance_matrix = BlockRealMatrix(6, 6)
covariance_matrix.setEntry(0, 0, coef_p)  
covariance_matrix.setEntry(1, 1, coef_p)  
covariance_matrix.setEntry(2, 2, coef_p) 
covariance_matrix.setEntry(3, 3, coef_v)  
covariance_matrix.setEntry(4, 4, coef_v)   
covariance_matrix.setEntry(5, 5, coef_v)   


constant_process_noise = ConstantProcessNoise(covariance_matrix, covariance_matrix)

kalman_builder = KalmanEstimatorBuilder()
kalman_builder.addPropagationConfiguration(propagator_builder, constant_process_noise)
kalman_builder.decomposer(QRDecomposer(1.0e-15))

kalman_estimator = kalman_builder.build()

my_observer = MyKalmanObserver()

kalman_estimator.setObserver(my_observer)

for measurement in measurements:
    kalman_estimator.estimationStep(measurement)
    print(kalman_estimator.getPhysicalEstimatedState())
    print(kalman_estimator.getPhysicalEstimatedCovarianceMatrix())

In my last estimation step I get the following output:

Estimation step 863 performed
Measured (observed) value: [-4971775.6363732, 4676129.498291034, -590578.2310049661, 5262.033526607344, -3849.086155134528, -75936.55217521796]
Estimated position: {-4,971,917.325422424; 4,675,930.244331793; -590,596.5369676582}
Estimated velocity: {258.547738749; -1,242.2032221552; -7,668.7775236856}
Estimated state: {924,926.1749803359; 1,210,874.9600476287; 366,566.2850783912; -3,312,354.312391164; -6,773,778.691221553; -457,916.6501600518}
Estimated covariance: BlockRealMatrix{{99992745.77583271,-5284.6998495584,-2887.2129397409,721878258.8615135,-724884839.2369523,-330410102.0415966},{-5284.6998495584,99987012.15671505,1889.7190272975,-724885369.0901498,-67109960.74374217,203119408.60165295},{-2887.2129397409,1889.7190272975,99996049.30283552,-330444847.30890274,203139907.68444553,37995838.938980006},{721878258.8615152,-724885369.0901481,-330444847.3089032,80423688589212.73,-1167534251.2945023,-1661587295.0364866},{-724884839.2369524,-67109960.74374244,203139907.68444544,-1167534251.3514972,80422454347701.64,816471892.3868656},{-330410102.04159635,203119408.60165343,37995838.93898024,-1661587295.1008139,816471892.4514334,80410992604994.03}}
Physical Estimated State {-4,971,917.325422424; 4,675,930.244331793; -590,596.5369676582; 258.547738749; -1,242.2032221552; -7,668.7775236856}
Physical Estimated Covariance Array2DRowRealMatrix{{99992745.77583271,-5284.6998495584,-2887.2129397409,80504.232526302,-80839.5279070426,-36847.5035191477},{-5284.6998495584,99987012.15671505,1889.7190272975,-80839.5869965421,-7484.1371356242,22651.9803026917},{-2887.2129397409,1889.7190272975,99996049.30283552,-36851.3783291231,22654.266370886,4237.3153858376},{80504.2325263022,-80839.5869965419,-36851.3783291231,1000215.2958931517,-14.5204185124,-20.664869482},{-80839.5279070426,-7484.1371356243,22654.266370886,-14.5204185131,1000199.9458481601,10.1543175867},{-36847.5035191477,22651.9803026918,4237.3153858376,-20.6648694828,10.1543175875,1000057.3981663161}}

The true final state is : (-5.852049e6,3.484293e6,-7.388627e5,515.842711,-721.091117,-7577.623566)

  • I can see that my estimated Pos/Vel match the Physical Estimated State but is the Estimated State way off?

  • Also, what’s the difference between the physical estimated covariance and the covariance matrix I get from (.getEstimate().getCovariance())?

I’m new to using Orekit and any help would be greatly appreciated.

Thanks in advance,

2

Hi @OrbitExplorer,

Welcome to the forum!

Physical covariance is the covariance in SI units whereas normalized covariance (the one you get with .getEstimate().getCovariance()) has been normalized using scales to avoid numerical errors during inversion.

I cannot tell you exactly why, it’s hard to debug a Kalman without being able to run the code.

  • What type of measurements are you using: Position/Velocity? If yes I think the velocity error is too large. In the example you gave the last component would be -75 936 m/s, which is a bit fast for a LEO :sweat_smile:
  • coef_p = 1e3
    coef_v = 1e6
    Same here, usually, it is advised to use a standard deviation \sigma for velocity that is 1000 times smaller than the position. Here \sigma_{vel}² = 1000 \sigma_{pos}², that may explain why your orbit is way off. 1 km/s of error is a very large error.

Hope this helps,
Maxime

1 Like
3

Hi there,

I believe the problem is the position error you pass to your integrator builder, it’s really big. So the propagation basically makes giant steps. Use 0.001 instead for example.

Cheers,
Romain.

1 Like
4

Thank you @MaximeJ and @Serrof I was able to make it work finally.

I just had to do some preprocessing on the measurement data before I fed it to the KMF.

Thanks again!

1 Like
5

Also how can I propagate the states in the future without measurement and without having to set it to zero? How can I access the estimated state at a random time (t) in between two measurements.

Is that even possible?

Thanks,

6

Hi @OrbitExplorer,

Method KalmanEstimator.estimationStep(observedMeasurement) returns an array of estimated propagators.
So you just have to propagate these propagators in the future to get what you want.

In Java, with one propagator, it looks like this:

Propagator[] estimatedPropagators = kalman_estimator.estimationStep(measurement);
SpacecraftState propagatedState = estimatedPropagators[0].propagate(myDate);

How can I access the estimated state at a random time (t) in between two measurements.

Same process overall but if you want to have access to the whole estimation interval at the end of your simulation then it’s a bit trickier.
For each measurement, you would have to:

  • Propagate the latest estimated propagator to the current measurement date and log it in an ephemeris that will fill an ephemeris list
  • Do the estimation and get the new estimated propagator

At the end of the loop, you can gather all the ephemeris in an AggregateBoundedPropagator. This last will then let you pick up a random state on the estimation interval.

It’s a bit overkill because you will be propagating twice, for the ephemeris and in the Kalman.
If you think that this feature is really important it can probably be added to the library (at least for the numerical EKF, I can’t tell for the others…)
If so, please open an issue on the forge. And, as usual, contributions are welcome!! :wink:

Cheers,
Maxime

1 Like
7

Hi Maxime,

Is there a similar method to propagate the covariance to a specified time? Or a way to propagate the two together?

Thank you,
CJM

8

Hi @cjm,

I don’t think the estimated propagator is set up for covariance computation, it’s a new object created on the fly using the propagator builder inside the Kalman.
So you will have to set it up by yourself using a StateCovarianceMatrixProvider (example linked) initialized with the latest estimated physical covariance of the Kalman.

Hope this helps,
Maxime

1 Like