Regarding Kalman Estimator Usage

Orekit Python Wrapper usage
1

Greetings,

From the repository of Clément Jonglez (a.k.a GorgiAstro), I have managed to determine orbits by using BLS estimator with ILRS data.

As a next step towards understanding orbit determination and propagation, I’ve dived into OD with Kalman Estimator and retraced the following Java codes and Orekit forum posts:

So basically, I have tried to implement KalmanEstimator instead of BLS estimator in the code of GorgiAstro as following:

from org.orekit.estimation.sequential import PythonKalmanObserver
class MyObserver(PythonKalmanObserver):
    def __init__(self):
        super(MyObserver, self).__init__()
    def evaluationPerformed(self, estimation):
        print(estimation.getCorrectedMeasurement().getObservedValue()) 
        
# The propagator builder to use in the Kalman filter.  initial covariance matrix
initP = MatrixUtils.createRealDiagonalMatrix([.1e4, 4.e3, 1., 5.e-3, 6.e-5, 1.e-4]) 
# The process noise matrices provider to use, consistent with the builder. process noise matrix
initQ = MatrixUtils.createRealDiagonalMatrix([1.e-4, 1.e-4, 1.e-4, 1.e-10, 1.e-10, 1.e-10])

# Setting up the estimator    
kalman_estimator_builder = KalmanEstimatorBuilder()
kalman_estimator_builder.addPropagationConfiguration(propagatorBuilder,ConstantProcessNoise(initP, initQ)) 
kalman_estimator = kalman_estimator_builder.build()

#Adding the observer
kalman_observer = MyObserver()
kalman_estimator.setObserver(kalman_observer)

In the above piece of code I am confused with the output of the line:

print(estimation.getCorrectedMeasurement().getObservedValue())

It gives me bunch of range values but I am somewhat lost here. If I understand the implementation correctly, estimationStep method updates the observation data at each iteration but where are these printed range values coming from and how should associate them?

Another issue I am facing is in the following code. Aforementioned estimationStep method is utilized here to feed the observation range data to the KalmanEstimator for each observation. However, I am not sure if I should do it this way.

observableSatellite = ObservableSatellite(0) # Propagator index = 0

for receiveTime, slrData in slrDataFrame.iterrows():
    if slrData['station-id'] in stationData.index: # Checking if station exists in the SLRF2014_POS+VEL, because it might not be up-to-date
        if not np.isnan(slrData['range']):  # If this data point contains a valid range measurement
            orekitRange = Range(stationData.loc[slrData['station-id'], 'OrekitGroundStation'], 
                                True, # Two-way measurement
                                datetime_to_absolutedate(receiveTime),
                                slrData['range'],
                                range_sigma,
                                range_weight,
                                observableSatellite)

            newPropagator = kalman_estimator.estimationStep(orekitRange) # Range data is now processed

In the above Java codes and posts, processMeasurements method is often used, so I thought maybe I should also utilize this method instead of estimationStep. Can processMeasurements method be used in the loop instead of estimationStep method? I also tried it, however, it does not accept orekitRange as input argument, in fact it does not even accept any kind of list argument for my case.

I am somewhat new to Python and I have no knowledge in Java so these might be a bit messy, sorry for that in advance. I can further try to clarify any vague part, if there is.

Thank you for the discussions!
Baris

2

Hi @Echulion

Basically, you can consider they are three different measurements:

  1. The observed measurement (i.e., the one printed in your example). The observed measurement is the raw measurement from your sensor. Please note that there is no different by calling estimation.getCorrectedMeasurement().getObservedValue() or estimation.getPredictedMeasurement().getObservedValue()

  2. The predicted measurement. It is the estimated measurement computed using the predicted state of the Extended Kalman Filter.

  3. The corrected measurement. It is the estimated measurement computed using the corrected state of the Extended Kalman Filter.

The estimationStep() method performs all the steps of the estimation (i.e., prediction and correction).

There is only one difference between estimationStep() and processMeasurements(). processMeasurements() process all the measurements in one call whereas estimationStep() process only one measurement per call. Therefore, there is no difference in the results. Calling estimationStep() in a loop will work perfectly.
For creating a List, did you tried List.cast_(...)?

It looks good. Based on my last comment, maybe you can also try:

  1. Read all the Range measurements
  2. Store them in a list
  3. Sort the measurements in chronological order You could try:
sorter = ChronologicalComparator()
measurements.sort(sorter)
  1. Call the processMeasurements() method.

Best regards,
Bryan

3

Hello @bcazabonne

Thank you so much for the detailed explanation! Since what I’ve been doing was not wrong, I decided to continue with estimationStep method.

Now, I have managed to save all the range measurements coming from the sensors. However, I was unable to save the corrected range values from the last iteration that is obtained from the KalmanEstimator. I tried the below code but it doesn’t really save the range values. In BLS there was getLastEstimations() method to obtain the timestamped range values, but apparently it is not the case with KalmanEstimator.

Is there a practical way to obtain the corrected range values of the last iteration of KalmanEstimator (compatible with the timestamps of observation data, so that I can calculate residuals in the end).

estimation_list = []

#My Kalman Observer class
from org.orekit.estimation.sequential import PythonKalmanObserver
class MyObserver(PythonKalmanObserver):
    def __init__(self):
        super(MyObserver, self).__init__()
    def evaluationPerformed(self, estimation):
        estimation_list = np.append(estimation_list, estimation.getCorrectedMeasurement().getEstimatedValue())
        print(estimation_list)
        print(estimation.getPredictedMeasurement().getEstimatedValue()) 
        print(estimation.getCorrectedMeasurement().getDate())

Another question came up as I am writing this reply, is this KalmanEstimator automatically update the covariance matrix at each iteration? Or should I update it manually?

Thank you for your time and reply,
Baris

4

Hi,

Please note that the Kalman FIlter is a sequential algorithm. Therefore, you can’t save directly all the corrected range in one method like the getLastEstimations() method. For a batch least squares, one iteration is a processing of all the measurements of your sample, whereas for the Kalman Filter one iteration is the processing of the measurements of your sample one by one (because it is sequential). So, you will have to save them one by one for each processed measurement. This is done in the observer, and what you do looks good. We do a close thing in the Observer of the ExtendedSemianalyticalKalmanFilter.java tutorial.

Measurements can have a dimension different to 1. Indeed, the dimension of an angular measurement, like azimuth-elevation, is 2. The dimension of a position measurement is 3 (i.e., X, Y, Z). Therefore, the return parameter of .getEstimatedValue() (or .getObservedValue()) is an array.
Because you works with range measurement, for which the dimension is 1, you could use instead estimation.getCorrectedMeasurement().getEstimatedValue()[0]

Also, to save the observed measurement and compute the residuals maybe you could add a observed_list = [] and fill it using estimation.getCorrectedMeasurement().getObservedValue()[0]

Yes , it does. The prediction step at iteration k+1 uses the corrected values at iteration k.
You can access the covariance matrix for each processed measurement in the Observer using estimation.getPhysicalEstimatedCovarianceMatrix(). Because the observer is called after the correction step, the returned matrix is the corrected covariance matrix.

Best regards,
Bryan