OutlierFilter with Batch LS Estimator

Orekit Python Wrapper usage
1

Hello Orekit Community,

I am trying to use the OutlierFilter with the batch LS estimator but not having much success. Does anyone have an example of how to implement the filter? Below is how I have already implemented it. However I reach an error once the code gets to the estimatedPropogatorArray. The main idea of the error seems to be “org.orekit.errors.OrekitException: unable to solve: singular problem”. I also checked the base weight of each data point before and after the filter had been implemented and the outliers still had a weight of 1.0.
###########################################
date = datetime_to_absolutedate(timestamp)
vector = Vector3D(pv_gcrf[[‘x’, ‘y’, ‘z’]].to_list())
orekit_position = Position(
date,
vector,
sigma_position,
1.0, # Base weight
observableSatellite
)
orekit_position.addModifier(rangeOutliersManager)
estimator.addMeasurement(orekit_position)
“”“Performing the orbit determination.”“”

estimatedPropagatorArray = estimator.estimate()
#########################################

Any advice or exmples would be welcome!
Thanks,
Olivia

2

Hi and welcome to this forum,

Just to be sure, your code only contains one measurement point for illustrative purposes? As in you removed the rest of the code that adds many points? Otherwise your system is underdetermined (more variables than equations) which could explain the singularity.

Anyhow, could you show us what’s inside rangeOutliersManager?

Best,
Romain.

1 Like
3

Hi Romain @Serrof ,

We are iterating through every data point individually, either adding both position and velocities measurements, or just position measurements if we don’t have velocities. The snippet I gave you was the part where we add each point that has only position measurements and not velocities. Here is rangeOutliersManager and also obsersvableSatellite.

rangeOutliersManager = OutlierFilter(2,5.0)
observableSatellite = ObservableSatellite(0)

Thanks,
Olivia

4

Hi @ojbBishop,

Welcome to the forum !

Maybe you’re editing too many measurements and the filter ends up with an underdetermined problem as Romain suggests.
I think of two potential reasons:

  • The theoretical standard deviations (sigma_position) you attribute to the measurements are too small
  • The measurements’ editing happens before the filter is somehow converged. In that case, try to increase the OutlierFilter.warmup iteration number

I think you got the idea of how it is used but here are some examples from the tutorials:

  • Creating a range outliers’ filter (values are parsed from an input file);
  • Adding it to Range measurements as an EstimationModifier

I strongly suggest you add a BatchLSObserver to your estimator. It will help you monitor what happens iteration after iteration and probably have a better understanding of the problem.
There’s an example of an implementation in the tutorials: OrbitDeterminationObserver

Hope this helps,
Maxime

5

Hi @MaximeJ,

Thanks for the feedback! I was able to get the filter to work, but now I’m wondering how to determine which points the filter rejects. I believe it has to do with the base weight of each point, but I’m not sure how to access the base weight of each point after the filter has finishing. Do you have any examples on how to tell which points the filter rejects?

Thanks!
Olivia

6

Hi @ojbBishop,

You’re welcome. I’m happy you managed to make it work!

It has to do with the value of the residual compared to the theoretical standard deviations (sigma_position) you attribute to the measurements.
if |residual| > Nsigma * sigma_position for each component of the measurement (in your case: pos_x, pos_y, pos_z), then the measurement is rejected.
Nsigma is the second value in the constructor of OutlierFilter.

Here’s the code inside OutlierFilter (the residual for each component is (observed[i] - theoretical[i]))

if (estimated.getIteration() > warmup) {

            // check if observed value is far to estimation
            final double[] observed    = estimated.getObservedMeasurement().getObservedValue();
            final double[] theoretical = estimated.getEstimatedValue();
            final double[] sigma       = estimated.getObservedMeasurement().getTheoreticalStandardDeviation();
            for (int i = 0; i < observed.length; ++i) {
                if (FastMath.abs(observed[i] - theoretical[i]) > maxSigma * sigma[i]) {
                    // observed value is too far, reject measurement
                    estimated.setStatus(EstimatedMeasurement.Status.REJECTED);
                }
            }
        }

From the previous piece of code, you’ll see that there is an enumerate EstimatedMeasurement.Status that can take the value PROCESSED or REJECTED.
So you can just check the status of your measurements at the end of the OD to know if they’ve been filtered or not.

Maxime

7

Hi @MaximeJ ,

Thank you again! I was stuck on trying to figure out which points were rejected so that’s very helpful.

I have one other question, when we add the modifier to our measurements we pass in a sigma_velocity as well as a sigma_position. Is the sigma_velocity used at all in rejecting the measurements? If not, what is sigma_velocitiy’s purpose?

Thanks,
Olivia

8

Hi Maxime @MaximeJ,

I am working with Olivia. We have a follow up question.

class BLSObserver(PythonBatchLSObserver):
    def evaluationPerformed(self, itCounts, 
           evCounts, orbits, orbParams, propParams, 
           measParams, provider, lspEval):
        print(f"iteration counts: {itCounts}")
        print(f"evaluation counts: {evCounts}")
        print(f"estimated orbital parameters: {orbParams.getDrivers()}")
        print(f"estimated propagator parameters: {propParams.getDrivers()}")
        print(f"estimated measurement parameters: {measParams.getDrivers()}")
        print(f"current evaluation of the underlying null: {lspEval.getCost()}")
        print(f"covariance matrix: {lspEval.getCovariances(.000001)}")
        # print(f"sigma: {lspEval.getSigma()}")
        for i in range(17):
            print(f"status of point: {provider.getEstimatedMeasurement(i).getStatus()}")
            print(f"time: {provider.getEstimatedMeasurement(i).getDate()}\n")
        listResiduals = lspEval.getResiduals()
        for i in range(0, listResiduals.getDimension(), 3):

# magnitudes of the residuals           
print(((listResiduals.getEntry(i)**2+listResiduals.getEntry(i+1)**2+listResiduals.getEntry(i+2)**2)**(1/2)))
        print(f"orbits: {orbits}")
        print()

We are using this code to determine the magnitudes of the residuals after each iteration of batch.

We are comparing these values from direct calculations of the residuals after the fit is completed:

bls_minus_measurements = np.linalg.norm(pos_bls - pos_measurements)

The values that we get using .getResiduals() don’t seem to agree with our direct calculations. We understand that they should disagree for the rejected measurements, but we don’t understand why the others don’t agree.

Are we misunderstanding what .getResiduals() is returning?

Thanks,
Ken Kiers, Olivia Bishop, and SeongHo Park

9

Hi @psh0078,

Sorry for the delay.

How much do they disagree?
Beware that getResiduals() will return the “weighted” residuals, i.e. the value vector from Batch LS model.
So there is a division by \sigma_{pos X,Y,Z} involved for each component. Is this what you observe?

Maxime