I am trying to implement a smoother to be used after orbit determination with a Kalman filter and have a question on how to extract the relevant information from Orekit and Hipparchus.
For the smoother, it is necessary to access the normalised corrected state estimates, normalised transition matrices and normalised predicted covariances at every step.
My current approach was to store all the necessary information using a KalmanObserver, and then apply the smoother afterwards, but some of the necessary components seem hard to retrieve.
The algorithm I use is derived Section 4.15 from Tapley, Schutz & Born (2004) and corresponds to the smoother proposed by Rauch, Tung & Striebel (1965). The main equations are
where k = \ell - 1 and all of the elements should be normalised based on the scale of the parameters.
\bm P_k^k is the covariance matrix at t_k after the correction step of the Kalman filter. This matrix is available to the KalmanObserver, but only in its physical (not normalised) form.
\bm \Phi_{l, k} is the normalised state transition matrix from t_k to t_\ell and it is also only available in unnormalised form
\bm P_\ell^k = \bm \Phi_{l,k}\bm P_k^k \bm \Phi_{l,k}^T + \bm\Gamma_{l, k}\bm Q_{k} \bm\Gamma_{l,k} is the predicted covariance at time t_\ell, including state noise compensation. This matrix is computed inside the AbstractKalmanFilter::predict method in Hipparchus, and added to the AbstractKalmanFilter::predicted field. Inside Orekit, KalmanEstimator::estimationStep passes only the KalmanModel to the KalmanObserver, which has no way to access the prediction step of the filter.
\xhat_k^k is the normalised corrected state vector at time t_k. This vector is available, but only in its unnormalised form.
It is very likely that I’m missing something or that there is something wrong with the approach in general. I am wondering if there have been previous efforts to implement a smoother or if you have any suggestions on how to obtain the desired, normalised vectors and matrices without recomputing and normalising them myself.
Thank you for considering this question.
Best,
Simon
I think you are on the right track, using an observer seems like the way to go.
If you have everything as unnormalized, then you coudl just do your smoothing directly with these quantities no? (if you’re afraid of numerical noise, you can do your own scaling).
Anyway it does look like the predicted covariance is not available as is. However you could recompute it, given that you still have access to the noise provider no?
A smoother for the EKF would be a great addition to Orekit, would you consider contributing it?
I tag @markrutten as he is quite knowledgeable on this.
Hi @svanhulle. You’re on the right track and you’re right … some of the quantities you need can’t be extracted out of the current orekit/hipparchus classes, because (for good reason) you don’t have access to the underlying Kalman filter. You could try and reconstruct the prediction as @Serrof suggested; the predicted covariance should be OK, but the predicted mean will be more challenging.
Unfortunately (because it would be a fair effort and require understanding the internals of the Kalman model implementation) my recommendation would be to extend the current classes and add the features you need to KalmanModel (e.g. methods to extract predicted mean and covariance). Unless your scenario is numerically demanding, maybe start with the physical quantities rather than the normalised ones? The normalised mean needs to be treated very carefully.
I agree with @Serrof that this would be a good orekit feature for the Kalman and unscented filters.
Thank you for your answers. I will try recreating the predicted covariance myself, and otherwise attempt to extend the KalmanModel.
Just out of curiosity, I was wondering if you (@markrutten) could elaborate on the good reason that we cannot access the underlying Kalman filter. Furthermore, is there a specific reason that KalmanEstimationCommon has a field correctedEstimate, but not predictedEstimate? It seems like such a field (along with the necessary getters) might already solve most problems, but maybe there is a specific reason why this is not desirable?
Thanks for your help; I will let you know if I make any progress.
Best,
Simon
I don’t think there’s a specific reason that the predicted quantities aren’t available … just no one has needed them until now! I agree that adding them should satisfy your requirements for a smoother. If you don’t mind, please add an issue to the orekit gitlab to record your idea. My answer above was assuming that you needed the functionality immediately, but waiting for an official release with that functionality is also an option! You can contribute your changes if you’re comfortable doing that.
I suspect that the underlying Kalman filter is hidden from a user because the output doesn’t follow a standard Kalman recursion. The predicted mean isn’t related to the previous posterior estimate in the usual way. Giving access would require carefully documenting the non-standard behaviour, which is complex and apparently hasn’t been needed before.
Thank you for the very clear explanation. I will explore the code a bit further and see if I can contributed the changes. Otherwise, I’ll create an issue.
Long-term design question. There’s more than 1 type of backward smoother. Would it be possible to write into the code the ability to choose between different types? I have some experience using the fixed-interval smoother; I would be willing to contribute this code.
I’m sorry, I come a bit late on the topic. Do you plan to implement a new recursive filter? Like fo example the BSEKF (Backward Smoothing Extended Kalman Filter)?
In that case, my opinion is to first contribute the mathematical part in Hipparchus first (in hipparchus-filtering) and after that to implement the “flight dynamics” part in Orekit.
Usually you’d run a filter forward, collect what you need on the forward pass, then smooth backward using the data you’ve collected. We’d need to think carefully about how we structure that in code.
It makes sense to have an implementation of the backward smoothing pass in Hipparchus. I suspect we’d have to give an option to collect the smoothing data in the filter. And we’d have a parallel astrodynamics version of the smoother in orekit (in a similar way to the parallel with the EKF).
As a bonus, if we’re careful about what we collect in the forward pass, we should be able to use exactly the same backwards smoothing recursion with output from the unscented filter.