Kalman model/estimator usage suggestions

muld00n · October 29, 2020, 8:27pm

Hi Orekit developers,

I am currently trying to implement a Kalman estimator for orbit determination. I tried to follow the tutorials in Gitlab, but I’m still having some trouble getting the parameter drivers to work. When I run the estimator using .processMeasurements() I end up with this error:

   at org.orekit.estimation.sequential.KalmanModel.checkDimension(Int32 , ParameterDriversList , ParameterDriversList , ParameterDriversList )
   at org.orekit.estimation.sequential.KalmanModel.getEvolution(Double previousTime, RealVector previousState, MeasurementDecorator measurement)
   at org.orekit.estimation.sequential.KalmanModel.getEvolution(Double d, RealVector rv, Measurement m)
   at org.hipparchus.filtering.kalman.extended.ExtendedKalmanFilter.estimationStep(Measurement measurement)
   at org.orekit.estimation.sequential.KalmanEstimator.estimationStep(ObservedMeasurement observedMeasurement)
   at org.orekit.estimation.sequential.KalmanEstimator.processMeasurements(Iterable observedMeasurements)

My code is as follows:

            ODEIntegratorBuilder integrator = new DormandPrince853IntegratorBuilder(0.001, 1000, 1);
        PositionAngle angle = PositionAngle.MEAN;
        double scale = 1; // scaling factor used for orbital parameters normalization (typically set to the expected standard deviation of the position)
        NumericalPropagatorBuilder propagatorBuilder = new NumericalPropagatorBuilder(initialGuess, integrator, angle, scale);

        // Define the univariate functions for the standard deviations      
        UnivariateFunction[] lofCartesianOrbitalParametersEvolution = new UnivariateFunction[6];

        // Evolution for position error
        lofCartesianOrbitalParametersEvolution[0] = new PolynomialFunction(new double[] { 100, 0, 1e-4 });
        lofCartesianOrbitalParametersEvolution[1] = new PolynomialFunction(new double[] { 100, 1e-1, 0 });
        lofCartesianOrbitalParametersEvolution[2] = new PolynomialFunction(new double[] { 100, 0, 0 });

        // Evolution for velocity error
        lofCartesianOrbitalParametersEvolution[3] = new PolynomialFunction(new double[] { 1, 0, 1e-6 });
        lofCartesianOrbitalParametersEvolution[4] = new PolynomialFunction(new double[] { 1, 1e-3, 0 });
        lofCartesianOrbitalParametersEvolution[5] = new PolynomialFunction(new double[] { 1, 0, 0 });

        UnivariateFunction[] propagationParametersEvolution = new UnivariateFunction[] 
        { 
            new PolynomialFunction(new double[] { 10, 1, 1e-4 }), 
            new PolynomialFunction(new double[] { 1000, 0, 0 }),
        };

        double threshold = 1e-10; // threshold to identify matrix singularity
        RealMatrix cov = blse.getPhysicalCovariances(threshold);
        UnivariateProcessNoise covarianceMatricesProvider = new UnivariateProcessNoise(cov, LOFType.TNW, angle, lofCartesianOrbitalParametersEvolution, propagationParametersEvolution);
        
        // Get estimatedMeasurementParameters
        ParameterDriversList estimatedMeasurementParameters = new ParameterDriversList();

        foreach (ObservedMeasurement ob in obs)
            foreach (ParameterDriver driver in ob.getParametersDrivers().toArray())
                estimatedMeasurementParameters.add(driver);

        //foreach (DelegatingDriver delegating in propagatorBuilder.getOrbitalParametersDrivers().getDrivers().toArray())
        //    foreach (ParameterDriver driver in delegating.getRawDrivers().toArray())
        //        estimatedMeasurementParameters.add(driver);
        //foreach (DelegatingDriver delegating in propagatorBuilder.getPropagationParametersDrivers().getDrivers().toArray())
        //    foreach (ParameterDriver driver in delegating.getRawDrivers().toArray())
        //        estimatedMeasurementParameters.add(driver);

        KalmanEstimator kalman = new KalmanEstimatorBuilder()
            .addPropagationConfiguration(propagatorBuilder, covarianceMatricesProvider)
            .estimatedMeasurementsParameters(estimatedMeasurementParameters)
            .build();

I’ve tried adding the parameter drivers from the measurements themselves and the numerical propagator (propagation parameters and/or orbital parameters). It seems that there may be some issue either with the way I am inputting the measurements (I’m adding the measurements using the Position class with a sigma = 1) or with the drivers themselves. It may also be possible that I am using the Kalman estimator entirely wrong.

Any help is appreciated, as this is my first time building a Kalman estimator for orbit determination.

Thanks,

bcazabonne · October 30, 2020, 7:38am

Hi @muld00n

Please note that there is an issue opened with UnivariateProcessNoise class since it does not consider measurement parameters in the covariance and process noise matrices. An ugly workaround is to set the measurement parameter covariance and noise in the propagationParametersEvolution object.
That 's probably the reason why you have an exception when the dimension of the process noise matrix is verified.

Furthermore, I have some questions:

Do you add force models to your numerical orbit propagator ?
Are you estimating other parameters than the orbit ?

I ask the second question because I see you are adding estimated measurement parameters in the following lines

    foreach (ObservedMeasurement ob in obs)
        foreach (ParameterDriver driver in ob.getParametersDrivers().toArray())
            estimatedMeasurementParameters.add(driver);

However, any parameters is set to estimated. In Orekit, when you want to estimate a parameter which is optionally estimable (all propagation and measurement parameters are optionally estimable), you have to specify it by calling the method .setSelected(true). Please note that the orbital parameters are by default estimated.

Best regards,
Bryan

muld00n · October 30, 2020, 6:04pm

I added some force models and selected them, but I think this may be increasing the requiredDimension too much in checkDimension(). This leads me to another question: If for example I have 6 parameter drivers from the orbital parameters only, would I need to adjust the sizes of lofCartesianOrbitalParametersEvolution and propagationParametersEvolution accordingly by adding 6 polynomial functions to each? And if I added more parameter drivers from propagation and measurements, would I need to add more dimensions to lofCartesianOrbitalParametersEvolution and propagationParametersEvolution ?

Due to the bug with UnivariateProcessNoise I tried to use ConstantProcessNoise instead. I made the initial covariance matrix equal to the covariance from a previous batch least squares estimation, but I am still getting the dimension error from before, even with a 6x6 covariance matrix per 6 delegated parameter drivers. What sort of initial covariance matrix or set of parameter drivers could I use to avoid this error?

bcazabonne · November 2, 2020, 10:23am

You have to define one UnivariateFunction per orbital parameter. Therefore, the size of lofCartesianOrbitalParametersEvolution is 6. That’s what you’ve done in your example.

After that, you have to define one UnivariateFunction per estimated propagation parameter. Therefore, the size of propagationParametersEvolution is the same as (propagatorBuilder.getNbSelected() - 6).

Since issue #720, you can’t define measurement parameters covariance and noise matrices using UnivariateProcessNoise class. A workaround is to define the estimated measurement parameter covariance inside the propagationParametersEvolution object. However, it is just a workaround. The number of estimated measurement parameters is define by estimatedMeasurementParameters.size();

At the end, you must verify that:

lofCartesianOrbitalParametersEvolution.length = 6
propagationParametersEvolution.length = propagatorBuilder.getNbSelected() - 6 (optionally + estimatedMeasurementParameters.size() for the workaround of estimated measurement parameters)
lofCartesianOrbitalParametersEvolution.length + propagationParametersEvolution.length = propagatorBuilder.getNbSelected() (optionally + estimatedMeasurementParameters.size() for the workaround of estimated measurement parameters)

In your example, I see two possible errors:

propagationParametersEvolution size is equal to 2. Therefore, it implies that you have two estimated propagation parameters. However, I don’t see any estimated propagation parameters. Make sure that propagatorBuilder.getNbSelected() - 6 is equal to 2.
You add drivers in the estimatedMeasurementParameters list. This means that you want to estimate that parameters but they are not set as selected. If you add drivers in estimatedMeasurementParameters list, you have to define UnivariateFunction for that parameters too. Since issue #720, that UnivariateFunction must be in the propagationParametersEvolution array.

muld00n · November 2, 2020, 9:03pm

I wasn’t able to get the estimator to build using UnivariateProcessNoise even when I check the sizes and dimensions. Do you mean to add the estimatedMeasurementParameters to the propagationParametersEvolution such that propagationParametersEvolution.length = propagatorBuilder.getNbSelected() - 6 + estimatedMeasurementParameters.size()?

i.e. if:
lofCartesianOrbitalParametersEvolution.length = 6
propagationParametersEvolution.length = 5
estimatedMeasurementParameters = 6 + 5 = 11 (and active = true)
Then for the workaround would I need to add the estimatedMeasurementParameters covariances to propagationParametersEvolution such that the size becomes 16 instead of 5? Or would I just replace the propagationParametersEvolution functions with estimated measurement parameter covariances so that the size is 11?

I tried building with ConstantProcessNoise instead which allowed the estimator to build, but still with an error when using processMeasurements():

at org.orekit.orbits.KeplerianOrbit.checkParameterRangeInclusive(String , Double , Double , Double )
at org.orekit.orbits.KeplerianOrbit..ctor(Double a, Double e, Double i, Double pa, Double raan, Double anomaly, Double aDot, Double eDot, Double iDot, Double paDot, Double raanDot, Double anomalyDot, PositionAngle type, Frame frame, AbsoluteDate date, Double mu)
at org.orekit.orbits.KeplerianOrbit..ctor(Double a, Double e, Double i, Double pa, Double raan, Double anomaly, PositionAngle type, Frame frame, AbsoluteDate date, Double mu)
at org.orekit.orbits.OrbitType.4.mapArrayToOrbit(Double[] , Double[] , PositionAngle , AbsoluteDate , Double , Frame )
at org.orekit.orbits.OrbitType.4.mapArrayToOrbit(Double[] , Double[] , PositionAngle , AbsoluteDate , Double , Frame )
at org.orekit.propagation.numerical.NumericalPropagator.OsculatingMapper.mapArrayToState(AbsoluteDate , Double[] , Double[] , PropagationType )
at org.orekit.propagation.integration.StateMapper.mapArrayToState(Double t, Double[] y, Double[] yDot, PropagationType type)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.ConvertedMainStateEquations.computeDerivatives(Double , Double[] )
at org.hipparchus.ode.ExpandableODE.computeDerivatives(Double t, Double[] y)
at org.hipparchus.ode.AbstractIntegrator.computeDerivatives(Double t, Double[] y)
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(ExpandableODE equations, ODEState initialState, Double finalTime)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbsoluteDate tEnd, Boolean activateHandlers)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbsoluteDate tStart, AbsoluteDate tEnd)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbsoluteDate target)
at org.orekit.estimation.sequential.KalmanModel.predictState(AbsoluteDate )
at org.orekit.estimation.sequential.KalmanModel.getEvolution(Double previousTime, RealVector previousState, MeasurementDecorator measurement)
at org.orekit.estimation.sequential.KalmanModel.getEvolution(Double d, RealVector rv, Measurement m)
at org.hipparchus.filtering.kalman.extended.ExtendedKalmanFilter.estimationStep(Measurement measurement)
at org.orekit.estimation.sequential.KalmanEstimator.estimationStep(ObservedMeasurement observedMeasurement)
at org.orekit.estimation.sequential.KalmanEstimator.processMeasurements(Iterable observedMeasurements)

The error seems to stem from the eccentricity of the processed orbit being negative, but not really sure why that would be the case. The measurements I used are approximately 2 hours apart in LEO, which may heavily affect the result of each given measurement process. If I were to use more closely spaced measurements, would that prevent this error from occurring?

bcazabonne · November 3, 2020, 7:30am

I meant to add the covariance evolution of the estimatedMeasurementParameters to the propagationParametersEvolution.

To consider the workaround yes.

I think I explained it badly. The size of your covariance and process noise matricies must correspond to the number of estimated parameters. For instance, if you estimate 6 orbital parameters, 2 propagation parameters and 3 measurement parameters, the size of your matrices must by 11*11. In order to know how many parameters do you estimate, could you tell me how many propagation parameters do you estimate (i.e. propagatorBuilder.getNbSelected() - 6) and how many measurement parameters do you estimate (i.e. estimatedMeasurementParameters.size())?

Does this error occur at the first measurement or randomly?

Could you tel me the epoch of the first measurement and the epoch of the last measurement? Furthermore, could you show me the initialization of your first point (i.e. epoch, orbital parameters and frame)?

muld00n · November 3, 2020, 4:49pm

I’m using 6 orbital, 5 propagation, and 0 measurement parameters in my current code. If measurement parameters are created from the drivers in getOrbitalParametersDrivers() and getPropagationParametersDrivers() do they also need to be added manually from the measurements themselves?

My measurements are here:

     Position ob1 = new Position(Resources.StringToAbsoluteDate("2020-10-27T08:18:12.184000Z"), new Vector3D(-3890.837837541537, -1657.937908623541, 5457.215954646945).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob2 = new Position(Resources.StringToAbsoluteDate("2020-10-27T10:22:30.184000Z"), new Vector3D(-3908.473081715485, -147.4737420180856, -5710.632147870466).scalarMultiply(1000), 1, 1, new ObservableSatellite(0)); 
     Position ob3 = new Position(Resources.StringToAbsoluteDate("2020-10-27T12:18:07.184000Z"), new Vector3D(4394.221158196921, 1735.496425283955, -5057.250383657628).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob4 = new Position(Resources.StringToAbsoluteDate("2020-10-27T14:18:05.184000Z"), new Vector3D(4868.153971365797, 502.8529674994862, 4872.016815083606).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob5 = new Position(Resources.StringToAbsoluteDate("2020-10-27T16:18:23.184000Z"), new Vector3D(-4984.26223209633, -1805.7834897521839, 4427.000752317712).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob6 = new Position(Resources.StringToAbsoluteDate("2020-10-27T18:17:58.184000Z"), new Vector3D(-4390.673710586525, -337.6068293981518, -5339.332525781674).scalarMultiply(1000), 1, 1, new ObservableSatellite(0)); 
     Position ob7 = new Position(Resources.StringToAbsoluteDate("2020-10-27T20:17:57.184000Z"), new Vector3D(5293.174196465872, 1846.096430649118, -4057.081811101456).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob8 = new Position(Resources.StringToAbsoluteDate("2020-10-27T22:18:02.184000Z"), new Vector3D(3823.924116933546, 155.3430260544334, 5745.599263588878).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob9 = new Position(Resources.StringToAbsoluteDate("2020-10-28T01:27:17.184000Z"), new Vector3D(4345.79947191626, 347.08517819695584, 5353.933545615597).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob10 = new Position(Resources.StringToAbsoluteDate("2020-10-28T02:18:28.184000Z"), new Vector3D(-3052.854266306893, 88.91045573243765, -6211.629490449016).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob11 = new Position(Resources.StringToAbsoluteDate("2020-10-28T03:26:30.184000Z"), new Vector3D(-5124.192342687291, -1858.568843612926, 4242.178718336686).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob12 = new Position(Resources.StringToAbsoluteDate("2020-10-28T04:18:14.184000Z"), new Vector3D(6076.519695985133, 1891.556386303349, -2713.227648271093).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob13 = new Position(Resources.StringToAbsoluteDate("2020-10-28T06:18:27.184000Z"), new Vector3D(2431.9295418249253, -269.35408227960994, 6453.496045991255).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob14 = new Position(Resources.StringToAbsoluteDate("2020-10-28T08:18:13.184000Z"), new Vector3D(-6338.606456362155, -1881.2255902869852, 2013.973342019202).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob15 = new Position(Resources.StringToAbsoluteDate("2020-10-28T10:24:43.184000Z"), new Vector3D(967.9950921507909, 1170.331701032592, -6753.5918691485795).scalarMultiply(1000), 1, 1, new ObservableSatellite(0)); 
     Position ob16 = new Position(Resources.StringToAbsoluteDate("2020-10-28T12:18:27.184000Z"), new Vector3D(6543.713545110223, 1850.9063470958629, -1261.89292247667).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));
     Position ob17 = new Position(Resources.StringToAbsoluteDate("2020-10-28T14:18:45.184000Z"), new Vector3D(962.7672502344093, -681.5692854313795, 6799.435255443729).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));

I tried using both the first 3 and last 3 to do an initial OD, and both trigger the error upon the Kalman/BLSE estimation. When I add propagation parameters to the least squares estimator the same error occurs, which tells me that the error must be with my measurements.

bcazabonne · November 4, 2020, 8:10am

Measurement parameters are not created from orbital and propagation parameter drivers. They are independent parameters. Orbital parameters are the 6 orbital parameters corresponding to the orbit of your satellite, propagation parameters are force model parameters (i.e. attraction coefficient, drag coefficient, etc.) and measurement parameters are parameters link directly to the measurements (i.e. clock offsets, station bias, etc.). As you already done, I think you should use the ConstantProcessNoise class.

Are your position values in the same frame as the initial spacecraft state?
Sometimes, position measurements are given in terrestrial frame but initial spacecraft state is initialized in inertial frame (e.g. EME2000). Due to the frame difference, problems occur during the orbit determination.
If the frame is different, please transform the position vectors defining your measurements in position vectors defined in the same frame as the initial spacecraft state.

If it doesn’t work, could you provide an executable code to debug your application?

Bryan

muld00n · November 17, 2020, 8:50pm

Hi Bryan, sorry for the delay.

I’ve created a crude console app in C# .NET framework that shows the gist of what I’m trying to do. You will need to create an orekit 10.2 DLL using IKVM and reference those DLLs when building. Let me know if you have any issues running it.

OrbitDeterminationScripts.zip (9.9 KB)

Thanks,
Connor

bcazabonne · November 18, 2020, 3:07pm

Hi,

I never used C#.NET. Therefore I am not going to try to run the code as it would take too much time to install. However, I can try to provide some comments just by looking at your code.

I saw that you initialize the measurement epochs using a StringToAbsoluteDate(String epoch). Please note that Orekit is able to initialize AbsoluteDate directly with the string representation of the epoch. Below an example using your first observation;

TimeScale utc = TimeScalesFactory.getUTC();
Position ob1 = new Position(new AbsoluteDate("2020-10-27T08:18:12.184", utc), new Vector3D(-3890.837837541537, -1657.937908623541, 5457.215954646945).scalarMultiply(1000), 1, 1, new ObservableSatellite(0));

Furthermore, could you tell me in which frame is expressed the position value when you initialize the measurement ? It is very important to initialize the position value in the same frame as the spacecraft state. Using different frames is prone to error.

Bryan

muld00n · November 18, 2020, 5:25pm

Thanks for the suggestion. The frame I’m using is J2000 throughout.