Propagating kalman filter with drag/SRP

Note 1: I’m working in python, but I think this question applies to general Orekit useage
Note 2: Continuing this discussion from a prior thread because it’s a specific use case question.

@bcazabonne @Serrof

So yesterday I tried adding the drag and SRP paramter values back to my FKF (using the initial covariance from the LS optimization and dynamic process noise generation functions), and this time I’m getting what looks at first glance like good updates over time. However, I have one remaining issue.

I mentioned before that I am taking the FKF data and running it through a fixed-interval smoother that I wrote myself. One of the values I require for the FIS math is the STM. However, while the KalmanModel class provides the STM at each time step, I can’t always use it because it’s not accurate for timesteps that come after a rejected measurement.

Previously I was getting around this by obtaining this STM through an alternate method: I was taking the propagator from the getEstimatedPropagators() value from the KalmanModel parameter created at each time step and adding the corrected covariance at that time step, then backing out the correct (t2 - t1) STM from the propagator in the smoothing process. However, now that I have the drag and SRP parameters included in the FKF process I am having trouble adding the corrected covariance to the propagator since it’s now 8x8 instead of 6x6. I know the command for doing so is setupMatricesComputation() (here). My code currently reads:

p = model.getEstimatedPropagators()[0]
cov = self.correctedCovariance.getSubMatrix(0, 5, 0, 5) # trimming the drag/srp parameters
covInit = StateCovariance(cov, self.date, eme2000, OrbitType.CARTESIAN, PositionAngleType.TRUE)
stmInit = model.getPhysicalStateTransitionMatrix()

# Append covariance to propagator
harvester = p.setupMatricesComputation("stm", None, None)
                                               
provider = StateCovarianceMatrixProvider("covariance", "stm", harvester, covInit)
p.addAdditionalStateProvider(provider)

I have verified that this works for my FIS smoother when I just have the parameters of pos/vel. I know that in the line creating the harverster, the first None is to initialize the STM as an identity matrix and the second would be a DoubleArrayDictionary vaue that I can leave blank for the standard use case. So in summary, my questions are:

  1. How do I create a DoubleArrayDictionary value that would properly initialize the STM now that I have an 8x8 instead of a 6x6
  2. Is it possible to do so in such a way that it would correctly propagate the STM matrix given the additional drag and SRP parameters?

The STM from t1->t2->t3 where t2 is the time of the rejected measurement is just the product of the STM from t1->t2 (F2) with t2->t3 (F3) i.e. F3 * F2. That might make your code a bit simpler?

@markrutten I have been working in Python, so to double-check your assertion before I replied I moved over to the Java side to see what results I got there. On the Java side you seem to be correct; I wasn’t actually able to do the multiplication so I can’t be completely sure, but the results I got definitely support your statement. However, on the Python side this does not seem to be correct. Linking my code to prove it is unfortunately not possible because it’s all work-related, but here are some numeric results for the STM from one of my test cases.

First case: all data points accepted

Date: <AbsoluteDate: 2024-09-21T00:21:42.000Z>
STM: <RealMatrix: Array2DRowRealMatrix{
{1.0004559976,-0.0003226336,0.0007706352,30.0044704922,-0.0032102944,0.0077318896,0.0001148092,0.0},
{-0.0003226352,0.9995337127,-0.0002428801,-0.0032102504,29.9953357153,-0.0024459733,0.0000049138,0.0},
{0.0007706459,-0.0002428822,1.0000106129,0.0077319397,-0.0024459164,30.0001955994,-0.0001509391,0.0},
{0.0000298165,-0.0000214037,0.0000515481,1.000438194,-0.00031942,0.000775674,0.0000077098,0.0},
{-0.0000214038,-0.0000310893,-0.0000163079,-0.0003194178,0.9995334785,-0.0002463113,0.0000003077,0.0},
{0.0000515495,-0.0000163081,0.0000013159,0.000775683,-0.0002463087,1.0000286418,-0.0000100235,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0}}>

Date: <AbsoluteDate: 2024-09-21T00:22:12.000Z>
STM: <RealMatrix: Array2DRowRealMatrix{
{1.0004023049,-0.0003127656,0.0007845248,30.0039321985,-0.0031105494,0.007864653,0.0001188597,0.0},
{-0.0003127671,0.9995329456,-0.0002529564,-0.0031105028,29.9953276894,-0.0025456371,0.0000030995,0.0},
{0.0007845354,-0.0002529586,1.0000650725,0.0078647025,-0.0025455814,30.0007419193,-0.000146208,0.0},
{0.0000262277,-0.0000207386,0.0000524332,1.000384233,-0.0003093391,0.0007883343,0.0000079631,0.0},
{-0.0000207388,-0.0000311428,-0.0000169723,-0.0003093368,0.9995326405,-0.0002561667,0.0000001865,0.0},
{0.0000524346,-0.0000169726,0.0000049582,0.0007883432,-0.0002561642,1.0000834408,-0.0000096883,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0}}>

Second case: Point 21:42 rejected by filter, KF moves on to point 22:12

Date: <AbsoluteDate: 2024-09-21T00:22:12.000Z>
STM: Array2DRowRealMatrix{
{1.0004023047,-0.0003127657,0.0007845252,30.0039321967,-0.0031105506,0.0078646566,0.0001188612,0.0},
{-0.0003127673,0.9995329455,-0.0002529566,-0.0031105041,29.9953276878,-0.0025456392,0.0000030995,0.0},
{0.0007845358,-0.0002529588,1.0000650728,0.0078647062,-0.0025455835,30.0007419226,-0.0001462097,0.0},
{0.0000262276,-0.0000207386,0.0000524332,1.0003842328,-0.0003093393,0.0007883346,0.0000079632,0.0},
{-0.0000207388,-0.0000311428,-0.0000169723,-0.0003093369,0.9995326403,-0.000256167,0.0000001865,0.0},
{0.0000524346,-0.0000169726,0.0000049583,0.0007883436,-0.0002561644,1.0000834412,-0.0000096885,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0}}>

Note that the two matrices for 2024-09-21T00:22:12.000Z are identical to the point of numeric error.

Here is the product of the two STM matrices from the case where both data points are accepted:

self.prevStep.date
Out  [9]: <AbsoluteDate: 2024-09-21T00:21:42.000Z>
step.date
Out  [10]: <AbsoluteDate: 2024-09-21T00:22:12.000Z>

step.stm.multiply(self.prevStep.stm)
Out  [7]: <RealMatrix: Array2DRowRealMatrix{
{1.0017542745,-0.001277799,0.0031022626,60.0336353435,-0.0252899048,0.0624105765,0.0004648387,0.0},
{-0.0012778283,0.9981346121,-0.0009852899,-0.0252900019,59.9626637311,-0.0199723021,0.0000172457,0.0},
{0.0031024829,-0.0009853372,1.0001162781,0.0624128524,-0.0199725214,60.0037623929,-0.0005977217,0.0},
{0.0000561619,-0.0000421653,0.0001040331,1.0016107247,-0.0012512061,0.0031376656,0.0000156628,0.0},
{-0.0000421676,-0.0000621815,-0.0000332974,-0.0012512277,0.9981324673,-0.001011867,0.0000004943,0.0},
{0.0001040506,-0.0000333012,0.0000063637,0.0031378825,-0.0010119048,1.0002619545,-0.0000197015,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0}}>

Here is the STM from the case where I omit the 21:42 measurement from the data set altogether:


step.date
Out  [2]: <AbsoluteDate: 2024-09-21T00:22:12.000Z>
step.stm
Out  [1]: <RealMatrix: Array2DRowRealMatrix{
{1.0017542743,-0.0012777991,0.003102263,60.0336353362,-0.0252899098,0.0624105922,0.0004649186,0.0},
{-0.0012778285,0.9981346119,-0.0009852901,-0.0252900069,59.9626637242,-0.0199723106,0.0000172489,0.0},
{0.0031024834,-0.0009853374,1.0001162785,0.0624128682,-0.01997253,60.003762407,-0.0005978256,0.0},
{0.0000561619,-0.0000421653,0.0001040331,1.0016107241,-0.0012512064,0.0031376668,0.0000156648,0.0},
{-0.0000421676,-0.0000621815,-0.0000332974,-0.0012512281,0.9981324668,-0.0010118676,0.0000004944,0.0},
{0.0001040506,-0.0000333012,0.0000063637,0.0031378837,-0.0010119055,1.0002619556,-0.000019704,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0}}>

Again note that the two matrices are identical to the level of numeric error.

While this statement is not as good supporting evidence because you have nothing but my word, I will also say that I have been working on this code for months and that I encountered repeated issues with my backward kalman smoother covariance results (dependent on the STM) until I found an independent means to calculate the STM outside of using the FKF values directly.

As to why the STM calculations would work properly in Java but not in Python, I’m afraid I don’t have enough familiarity with the python binders to even venture a guess. Sorry I cannot be more help.

For completeness, here is the data I got from the test case I ran on the Java side.

Case with all steps accepted:

currentDate
AbsoluteDate@164 "2016-02-14T03:19:01.047Z"
model.getPhysicalStateTransitionMatrix()
Array2DRowRealMatrix@172 "Array2DRowRealMatrix{
{1.0004816265,-0.0011112766,-0.0001968388,84.0136878342,-0.0310528307,-0.0057665167},
{-0.0011112685,1.0002433815,0.0001767431,-0.0310527185,84.0065328201,0.0051372228},
{-0.0001968459,0.0001767508,0.9992755628,-0.0057666172,0.005137331,83.9797889378},
{0.0000116466,-0.0000264085,-0.0000049047,1.0004961823,-0.0011067497,-0.0002151153},
{-0.0000264081,0.0000055619,0.0000043677,-0.0011067417,1.0002233828,0.0001901085},
{-0.0000049051,0.0000043682,-0.0000171813,-0.0002151225,0.0001901162,0.9992810059}}"
currentDate
AbsoluteDate@175 "2016-02-14T03:21:17.446Z"
model.getPhysicalStateTransitionMatrix()
Array2DRowRealMatrix@183 "Array2DRowRealMatrix{
{1.0014058718,-0.002880407,-0.0006959119,136.4641914962,-0.1303372198,-0.0334570849},
{-0.0028803164,1.0004493379,0.0005890385,-0.1303351575,136.4174180407,0.0279535833},
{-0.0006959932,0.000589126,0.9981487636,-0.0334589323,0.0279555712,136.3154989191},
{0.0000210558,-0.0000420471,-0.0000107952,1.0014625127,-0.0028527412,-0.0007761737},
{-0.0000420438,0.0000059684,0.0000090119,-0.0028526505,1.0003618312,0.0006399018},
{-0.0000107981,0.0000090152,-0.0000269077,-0.000776255,0.0006399893,0.9981796317}}"
currentDate
AbsoluteDate@186 "2016-02-14T03:23:01.245Z"
model.getPhysicalStateTransitionMatrix()
Array2DRowRealMatrix@194 "Array2DRowRealMatrix{
{1.0003542161,-0.0013005187,-0.0005557631,103.8115277481,-0.0448060172,-0.0196774714},
{-0.0013005038,1.0003554652,0.0005557794,-0.0448057606,103.8106673841,0.019481725},
{-0.0005557741,0.0005557967,0.9992912133,-0.0196776617,0.0194820251,103.7748234353},
{0.0000069843,-0.0000249547,-0.0000109599,1.0003700357,-0.0012893238,-0.0005817048},
{-0.000024954,0.0000065066,0.0000108482,-0.0012893089,1.0003192093,0.0005701027},
{-0.0000109605,0.0000108491,-0.0000134564,-0.0005817158,0.0005701201,0.9993116488}}"
currentDate
AbsoluteDate@197 "2016-02-14T03:24:55.045Z"
model.getPhysicalStateTransitionMatrix()
Array2DRowRealMatrix@205 "Array2DRowRealMatrix{
{1.0008510988,-0.0017112484,-0.000811731,113.8326705851,-0.0644593434,-0.0315594316},
{-0.0017112181,1.000110771,0.0006545361,-0.064458769,113.8032437041,0.0251545388},
{-0.0008117561,0.0006545679,0.9990397437,-0.0315599066,0.0251551417,113.7641224261},
{0.0000151508,-0.00002987,-0.0000146248,1.0008716204,-0.0016871933,-0.0008522602},
{-0.0000298687,0.000001517,0.000011653,-0.001687163,1.00006075,0.0006713421},
{-0.0000146259,0.0000116544,-0.0000166111,-0.0008522852,0.0006713739,0.9990692417}}"
currentDate
AbsoluteDate@208 "2016-02-14T03:26:59.644Z"
model.getPhysicalStateTransitionMatrix()
Array2DRowRealMatrix@216 "Array2DRowRealMatrix{
{1.0010299225,-0.0020318995,-0.0010582046,124.6423068763,-0.0837513131,-0.0450643424},
{-0.0020318517,1.0000731722,0.0008374587,-0.0837503199,124.6006632554,0.0352124721},
{-0.0010582438,0.0008375095,0.9988992667,-0.045065157,0.0352135266,124.5540887478},
{0.0000167566,-0.0000323755,-0.0000174212,1.0010558302,-0.0020009511,-0.0011119695},
{-0.0000323736,0.0000006611,0.0000136076,-0.0020009033,1.0000076045,0.0008576805},
{-0.0000174228,0.0000136097,-0.0000173419,-0.0011120087,0.0008577313,0.9989389283}}"
currentDate
AbsoluteDate@219 "2016-02-14T03:28:41.244Z"
model.getPhysicalStateTransitionMatrix()
Array2DRowRealMatrix@227 "Array2DRowRealMatrix{
{1.0007523528,-0.0012695322,-0.0008569258,101.6256660591,-0.0426767385,-0.0295352951},
{-0.0012695147,0.9998705874,0.0006092941,-0.0426764414,101.5950132111,0.0207679168},
{-0.0008569407,0.0006093131,0.9993781275,-0.0295355465,0.0207682382,101.5793424563},
{0.0000149295,-0.0000248096,-0.0000171709,1.0007635097,-0.0012506632,-0.0008873564},
{-0.0000248087,-0.0000028893,0.000012071,-0.0012506456,0.9998352032,0.000616932},
{-0.0000171716,0.0000120719,-0.0000119982,-0.0008873713,0.000616951,0.9994023575}}"

Case with rejected data points:

currentDate
AbsoluteDate@171 "2016-02-14T03:19:01.047Z"
model.getPhysicalStateTransitionMatrix()
Array2DRowRealMatrix@173 "Array2DRowRealMatrix{
{1.0004816265,-0.0011112766,-0.0001968388,84.0136878342,-0.0310528307,-0.0057665167},
{-0.0011112685,1.0002433815,0.0001767431,-0.0310527185,84.0065328201,0.0051372228},
{-0.0001968459,0.0001767508,0.9992755628,-0.0057666172,0.005137331,83.9797889378},
{0.0000116466,-0.0000264085,-0.0000049047,1.0004961823,-0.0011067497,-0.0002151153},
{-0.0000264081,0.0000055619,0.0000043677,-0.0011067417,1.0002233828,0.0001901085},
{-0.0000049051,0.0000043682,-0.0000171813,-0.0002151225,0.0001901162,0.9992810059}}"
currentDate
AbsoluteDate@182 "2016-02-14T03:21:17.446Z"
model.getPhysicalStateTransitionMatrix()
null
currentDate
AbsoluteDate@191 "2016-02-14T03:23:01.245Z"
model.getPhysicalStateTransitionMatrix()
null
currentDate
AbsoluteDate@200 "2016-02-14T03:24:55.045Z"
model.getPhysicalStateTransitionMatrix()
null
currentDate
AbsoluteDate@209 "2016-02-14T03:26:59.644Z"
model.getPhysicalStateTransitionMatrix()
null
currentDate
AbsoluteDate@218 "2016-02-14T03:28:41.244Z"
model.getPhysicalStateTransitionMatrix()
Array2DRowRealMatrix@226 "Array2DRowRealMatrix{
{1.0010296923,-0.0013562177,-0.0008493595,101.6350595417,-0.0455201505,-0.0293237746},
{-0.0013561966,0.9997472828,0.0005376813,-0.0455197931,101.5908340474,0.0183480464},
{-0.0008493784,0.0005377017,0.9992242619,-0.0293240948,0.0183483909,101.5741315593},
{0.0000203921,-0.0000264635,-0.0000170479,1.0010409347,-0.0013319438,-0.0008824102},
{-0.0000264625,-0.0000053169,0.0000106641,-0.0013319227,0.9997117964,0.0005456185},
{-0.0000170488,0.0000106651,-0.0000150265,-0.0008824292,0.0005456388,0.9992485066}}"

Ah, so you’re missing the STMs when the measurements are rejected? Remember from our earlier discussion, the EKF processes a measurement at a time, regardless of whether the measurement is rejected … if the measurement is rejected then the updated state is just the predicted state. You need to extract the STM at all time-steps, but it looks like you’re getting nulls when the measurement was rejected? I’ll take a closer look at the orekit/hipparchus internals to see if that’s expected behaviour … we might need to make some changes.

1 Like

@markrutten In this post when you said

The STM from t1->t2->t3 where t2 is the time of the rejected measurement is just the product of the STM from t1->t2 (F2) with t2->t3 (F3) i.e. F3 * F2.

I interpreted that to mean that the output of the STM at t3 was supposed to equal F3*F2. In retrospect however, reading it again I can see how you were actually saying I should get the STM at t2 and t3 and use it to create F3*F2. Apologies for the confusion.

As for the actual STM value at the rejected measurement t2, no I cannot retrieve it. So perhaps ultimately the problem with the KF code is that the STM returns null for a rejected measurement.

Also, in retrospect the matrices on the Java side do match up… so long as you expect the numeric error to be quite bad. Or maybe it’s just that it’s dealing with a lot of large numbers. So I guess the python and java behaviors do match. So that’s good.

So it turns out that while the STM is not created directly in the KalmanModel for a rejected measurement, you can pull it indirectly by using the command:

model.getPredictedSpacecraftStates()[0]

The additional states for the core 6x6 dY/dY0 and additional columns dY/dPp for the added parameter values (drag, radiation, etc) are in there, and the STM can be put together from these values.

1 Like