Partial derivatives of a range measurement with respect to the initial PV coordinates

Orekit Python Wrapper usage
1

Hey there,

I am writing a Python code using your dedicated wrapper. I have created a KeplerianOrbit object in a standard way as

orbit = KeplerianOrbit(pvCoordinates, inertialFrame, initialDate, Constants.WGS84_EARTH_MU).

The aim of my work is to compute a range (in particular a one-way distance) measure between my satellite and a reference point, which is assumed to be fixed in time. This task is very easy and I already get my range measurements by simply considering the norm of r = r_sat - r_p, where r_sat is the 3D vector representing the satellite position in the inertialFrame and r_p is that of the reference point. So the code is able to collect a set of ranges r_1, r_2, ā€¦ r_n, where a single range is computed every s seconds, while the satellite orbit is propagated through the propagator defined as

propagator = KeplerianPropagator(orbit).

Now, I would like to estimate the derivatives of my measurements r_1, r_2, ā€¦ r_n with respect to the initial pvCoordinates I have chosen for the orbit. My problem here is that I donā€™t
know how to compute the derivatives of the satellite positions at the time instants t_1, t_2,ā€¦ t_n in which the original range measurements were calculated.

Can you please tell me which is the most convenient way to carry out this step? I must admit I am new to Orekit and I donā€™t know if there is already a method allowing it. This code is related to a Batch estimator and the partial derivatives of the ranges have to be taken with respect to an integration function G, which depends on the initial pvCoordinatesā€¦

Many thanks,
Samuele

2

Hi @samuele,

Welcome to the Orekit community!

First of all:

Are you aware that Orekit already provides classes and methods for batch least-squares orbit determination ?
If not, you should have a look the estimation package of the Orekit tutorials in Java. Starting for example with the NumericalOrbitDetermination tutorial.
For Python, here is a very good example written by @yzokras of Satellite Laser Ranging Orbit Determination.

That is indeed an effective way to do it but beware that you wonā€™t be able to have corrections to your measurements with this. Here you will miss the light time delay correction and latter on, if you want to add more complexity to your models, you wonā€™t be able to add tropospheric/ionospheric delays, station biases etc.
Note that Orekit also provides a measurement generation package. Youā€™ll find an example of how to use it in the MeasurementGenerator tutorial in the estimation/performance package of the Orekit tutorials, written by @bcazabonne.

For this, I advise you to have a look at the AbstractBatchLSModel.fetchEvaluatedMeasurements which performs the calculation of the Jacobian of the batch LS OD problem, i.e. the matrix of the derivatives of the measurements (range measurements in your case) with respect to the estimated parameters (the initial orbital parameters in your case).

Hope this will help you,
Maxime

3

To add to the measurement generation examples mentioned by Maxime, I wrote an example in Python here: orbit-determination-examples/00-generate-measurements-and-tle.ipynb at p/jonglez/measurement-generation Ā· GorgiAstro/orbit-determination-examples Ā· GitHub
(I plan to move this example to the official Python wrapper examples)

4

Thank you so much @MaximeJ and @yzokras for your helpful replies.

I am aware that Orekit provides many tools for batch least-squares orbit determination and the Laser Ranging Orbit Determination example written by @yzokras was my first reference in the development of my code. However, the Range (measurement) class requires a Ground Station and in my code I do not need one. In my code I create a set of measurements which represent the distance between the satellite and the Solar System Barycentre and it has no sense to generate a Ground Station at this point. Is there a way to use the Range class to deal with a fixed point such as the SSB instead of a Ground Station? I really just need to get the distance between the satellite and the SSB.

Starting from the observed range measurements r_1, ā€¦ r_n and the computed measurements ||r_sat - r_SSB|| + D_c, where D_c is a component due to the clock bias and the clock drift, I would like to determine the orbit (i.e. x, y, z, vx, vy and vz at time t), the clock bias and the clock drift of the satellite. What do you think could be the best way to proceed? I was thinking about re-writing the function estimate() of the class ā€œBatchLSEstimatorā€ including all of the computations done in AbstractBatchLSModel.fetchEvaluatedMeasurements but maybe thereā€™s a faster way to do thisā€¦

Many thanks,
Samuele

5

Hello,

at the end I decided to re-organize my code following the algorithm of the fetchEvaluatedMeasurements method. Now, the code is able to:

  • Compute the partial derivatives of the current Cartesian coordinates with respect to the current orbital state. I did it with:

aCY = JArray_double2D(6, 6)
reference_orbit.getJacobianWrtParameters(builder.getPositionAngle(), aCY)
dCdY = Array2DRowRealMatrix(aCY)

where reference_orbit is a KeplerianOrbit I defined as when I started this topic and the builder has been constructed with

integrator = DormandPrince853IntegratorBuilder(0.001, 1000.0, 1.0)
scale = 1.0
builder = NumericalPropagatorBuilder(reference_orbit, integrator, PositionAngle.MEAN, scale)

  • Compute the Jacobian of the measurement with respect to current orbital state dMdY. I did it by simply differentiating the range measurement formula (||R_sat - r_SSB||) with respect to the position (i.e. the P coordinates) of the satellite at the current state.

After these two steps, one has to compute the Jacobian of the measurement with respect to initial orbital state dMdY0. This can be done starting from the Jacobian of the current orbital state with respect to the initial orbital state dYdY0 i.e. leading to dMdY0 = dMdY*dYdY0.

To compute dYdY0, I followed what the method fetchEvaluatedMeasurements implements. If I understand this code correctly, the key method I have to use is getStateJacobian (SpacecraftState state, double[][] dYdY0) of the class JacobiansMapper. This method takes the current state as an argument and computes its partial derivatives with respect to the initial state. However, it is not too easy to me to understand whatā€™s going on.

I have created a JacobiansMapper object as

mapper = JacobiansMapper(ā€˜mapperā€™, None, OrbitType.KEPLERIAN, PositionAngle.MEAN)

where the None is because I do not need additional parameters (for the moment, clock error parameters will be added in a second stepā€¦) since I am interested only in the 6-dimensional state. However, this gives me a NotImplementedError and I am not able to continue. Maybe I am missing something obvious?

After this, I would like to get dYdY0 with

aYY0 = JArray_double2D(6, 6)
mapper.getStateJacobian(current_state, aYY0)
dYdY0 = Array2DRowRealMatrix(aYY0)

Thank you so much for your patience! I hope the first two points have been written correctly.

Have a nice evening,
Samuele

6

Hi @samuele,

To configure the JacobiansMapper you need to create it as done in BatchLSModel.configureDerivatives method.

Using the reference orbit here means you will always get the derivatives of the Cartesian parameters wrt to the Keplerian parameters for the initial reference orbit.
I donā€™t think this is what you want. You want the current orbit state derivatives (the orbit you are using when computing your range measurements)

7

Hi @MaximeJ,

thank you, your replies have helped me a lot. Now my code is almost done!

I am trying to configure the JacobiansMapper as done in BatchLSModel.configureDerivatives but I donā€™t understand how to create a PartialDerivativesEquations (equationName, (NumericalPropagator) propagator) object. So my questions are:

  1. In BatchLSModel.configureDerivatives the equationName is chosen to be BatchLSModel.class.getName() + ā€œ-derivativesā€. But what is the meaning of this? And how could I write it using your Python wrapper?

  2. In the second argument, there is a (NumericalPropagator) cast of the propagator. Do I have to create a new NumericalPropagator object or I can just get it from my NumericalPropagatorBuilder? I have seen that there is a method called buildPropagator(double[] normalizedParameters) in the NumericalPropagatorBuilder class, returning a NumericalPropagator. However, I donā€™t understand what I have to write as normalizedParameters. Normalized with respect to what?

Thank you again for your patience and sorry for all of these questions!

Best,
Samuele

8

Hi @samuele,

You are very welcome!

Method configureDerivatives adds what we call AdditionalEquations to the propagator, itā€™s a set of equations that will be integrated along with the orbit of the spacecraft during propagation.

In the case of PartialDerivativesEquations the equations provided are the so-called variational equations \frac{\partial \frac{\partial Y}{\partial Y_0}}{\partial t} where Y is the state vector (position/velocity) at t and Y_0 the state vector at t_0.
Integrating these gives you the matrix \frac{\partial Y}{\partial Y_0} that you retrieve in the code with JacobiansMapper.getStateJacobian

Regarding the name, I donā€™t think itā€™s important to you, you can just put it to ā€œderivativesā€.
Itā€™s useful in Orekit when you add several additional equations to a propagator. After propagation you can get each integrated additional state using their respective name.
Since you only need one additional equations in your case, you donā€™t really need this feature.
Plus youā€™re using the JacobiansMapper that already encapsulates the PartialDerivativesEquation so this name is really not important here.

Note that these interfaces are known to be cumbersome and @luc just changed them in the development branch of Orekit Java.
They will be deprecated in next Orekit version; and replaced by a more user-friendly API.

By the way, an example of an OD run using the propagator is given in AbstractBatchLSModel.value method.

Do:

propagator = propagator_builder.buildPropagator(propagator_builder..getSelectedNormalizedParameters())

Normalization was introduced to avoid numerical errors when manipulating matrix in the batch least-square (or Kalman) estimators.
If you donā€™t use Orekit OD features you will have to handle it by yourself.
Normalization for the batch LS is done in Orekit in AbstractBatchLSModel.fetchEvaluatedMeasurement.
Parameters are normalized using their scale attribute while measurements are normalized using their sigmas (.i.e. their theoretical standard deviations).
You can read more about normaliaztion (and about OD in Orekit) in the official documentation here.

Maxime

1 Like
9

Hi @MaximeJ,

all right, now I understand the meaning of this normalization. Will be surely helpful in the future since I will have to deal with some matrix manipulation in Batch LS estimators.

For what concerns the equationName, I am still having some problems. Running

    propagator = builder.buildPropagator(builder.getSelectedNormalizedParameters())
    partials = PartialDerivativesEquations('derivatives', propagator)
       
    rawState = propagator.getInitialState()       
    stateWithDerivatives = partials.setInitialJacobians(rawState)
    propagator.resetInitialState(stateWithDerivatives) 
    mapper = partials.getMapper()        
    
    aYY0 = JArray_double2D(6, 6)
    mapper.getStateJacobian(current_state, aYY0) 
    dYdY0 = Array2DRowRealMatrix(aYY0)

I get the error: org.orekit.errors.OrekitException: unknown additional state ā€œderivativesā€ at org.orekit.propagation.SpacecraftState.getAdditionalState(SpacecraftState.java:653). What could be the cause of it?

Thank you,
Samuele

10

Hi @samuele,

Are you calling the propagate method on the propagator on which you performed the call: propagator.resetInitialState(stateWithDerivatives) ?
It seems like your SpacecraftState object does not contain the additional state added by PartialDerivativesEquations.

I was thinking of something.
Instead of redoing the whole OD process you could probably use Orekitā€™s by writing a new measurement class for your range measurements.
For this you need to extend class AbstractMeasurement and provide an implementation for the method theoreticalEvaluation.
You could follow the example of Orekit Range class.
Replace the GroundStation attribute by a PVCoordinatesProvider which would be the Solar System Barycenter ephemeris that you will get with:

CelestialBody SSB = CelestialBodyFactory.getSolarSystemBarycenter()

Then in the theoreticalEvaluation method, provide the method to compute your range. Derivatives can be provided either analytically (as youā€™ve done it already) or by automatic differentiation using the MyClass<Gradient> pattern (a bit harder but usually better; this is what is done in the actual Orekit Range class).

Maxime

11

Hi @MaximeJ,

Here it is what I have done:

  1. I collect n measurements r_1, ā€¦ r_n at the instants t_1, ā€¦ t_n respectively.
  2. Once I have collected n measurements, I call my batch LS estimation function which takes as argument the reference_orbit (i.e. the predicted one).
    Starting from reference_orbit, I create a NumericalPropagatorBuilder as builder = NumericalPropagatorBuilder(reference_orbit, integrator, PositionAngle.MEAN, scale) and from the builder I get a NumericalPropagator with propagator = builder.buildPropagator(builder.getSelectedNormalizedParameters()).
  3. Then, I want to compute the derivatives (with respect to the initial state) of my estimated measurements e_1, ā€¦ e_n corresponding to the instants t_1, ā€¦ t_n. For this, I have created a for loop whose index i goes from 1 to n in which I call the propagate method on the propagator retrieved in the second step as current_state = propagator.propagate(t_i).
  4. I get the partial derivatives of the current Cartesian coordinates with respect to the current orbital state with
    aCY = JArray_double2D(6, 6)
    current_state.getOrbit().getJacobianWrtParameters(builder.getPositionAngle(), aCY)
    dCdY = Array2DRowRealMatrix(aCY)
  5. I get the Jacobian of the measurement with respect to the CURRENT orbital state by simply differentiating the distance between the predicted satellite position and the SSB. I carried it out through the sympy package tools.
  6. I (try toā€¦) compute the Jacobian of the measurement with respect to the INITIAL orbital state. For this, I have written:
    partials = PartialDerivativesEquations(ā€˜derivativesā€™, propagator)
    rawState = propagator.getInitialState()
    stateWithDerivatives = partials.setInitialJacobians(rawState)
    propagator.resetInitialState(stateWithDerivatives)
    mapper = partials.getMapper()
    aYY0 = JArray_double2D(6, 6)
    mapper.getStateJacobian(current_state, aYY0)
    dYdY0 = Array2DRowRealMatrix(aYY0)
    dMdY0 = dMdY.transpose().multiply(dYdY0)

At the next iteration, I re-compute current_state = propagator.propagate(t_i) and I repeat the total process.

In this way I would obtain, at the i-th iteration, the derivative of e_i with respect to the initial state. My idea is then to create a matrix whose i-th row corresponds to the derivative of e_i with respect to the initial state. Can you please tell me if I am doing something wrong?

For what concerns the alternative way you proposed me, it is surely a good way. But do you suggest me to change my approach at this point?

Have a nice weekend and thank you again @MaximeJ!
Samuele

12

Hi again @samuele,

Iā€™m not sure to understand. Do you do this:

in a loop (i.e. at each step of your propagation) or at the end of your process ?

You need to setup the propagator withā€¦

ā€¦before calling the propagate function.
Then use the mapper at each step of the propagation to get the derivatives; meaning, inside the for loop of step 3

Well, hard to say, youā€™ve been going very far in your design, far enough to say ā€œhey I donā€™t want to throw away all that workā€.
The risk in sticking to your initial way is that you may have a hard time setting it up or run into bugs that will be hard to correct.
The advantage of the way I propose is that you will be able to use Orekitā€™s already designed and validated OD.
I would probably give a try to the alternative method and see if it works. If it does, your work is finished. If it doesnā€™t, go back to your initial method.
I should have told you about this right at the beginning, Iā€™m sorryā€¦

Have a nice week-end too,
Maxime

14

Hi @MaximeJ ,

as you expected, I am running into several bugs following my original way. I should have followed the approach you suggested me since the beginningā€¦

My question is: how can I write a new measurement class if I am using the Python wrapper?

Best,
Samuele

15

Hi @samuele,

Unfortunately Iā€™m not very familiar with Pythonā€¦

I think you need to extend the PythonAbstractMeasurement class and write your own theoreticalEvaluation method in there.

See for example the AbstractDetectorTest on the Orekit Python Wrapper repo. There you will see how to extend PythonEventHandler and PythonEventDetector.

Best,
Maxime

16

Hello,

I have extended the PythonAbstractMeasurement class and I wrote my own theoreticalEvaluation method. I have plotted my theoretical evaluations (computed through the signalTimeOfFlight method) and they are correct but there is an orekit.InvalidArgsError in

estimated.setStateDerivatives(0, derivatives)

where estimated is an EstimatedMeasurement (org.orekit.estimation.measurements.EstimatedMeasurement@5a5f6a0c) and derivatives is

JArray <double> [0.3385106376097706, -0.8406770751239679, 0.4271676775157459, 0.0, 0.0, 0.0].

I guess the problem is related to the definition of derivatives as a JArray in some ways. I have solved several issues like this in the previous lines but this one is putting me in troubleā€¦ what is the problem here?

All the best,
Samuele

17

Hi @samuele,

This one might actually be a bit tricky.

The Java for the input derivatives of EstimatedMeasurement.setStateDerivatives is double[]... (notice the three dots).
So itā€™s actually expecting a double[][] (a matrix) or a list of double[].

In Python I think youā€™ll need a double-dimension JArray for this one.

18

Hi,

The EstimatedMeasurements.setStateDerivatives as I understand is this method:
public void setStateDerivatives(final int index, final double[]ā€¦ derivatives) {

In Python the variable number term is represented by a list like [stuff, stuff ]. The type of ech stuff should be lists of doubles.

so for your case it should likely be a python list of one or more JArray lists, I think :slight_smile:

Regards!

19

Hi @MaximeJ and @petrus.hyvonen,

thank you so much for your replies, now that line works! And I am finally able to run the entire code.

To test the code, I firstly set the estimated values in the theoreticalEvaluation method to be equal to the observed values. As we imagine, this returns the orbit we give in input. However, this happens only if the number of measurements is <= 4. In fact, if I consider > 4 measurements I get an org.orekit.errors.OrekitException: unable to solve: singular problem.

If instead I do not put the estimated values to be equal to the observed ones (i.e. considering the true estimated values) I get either the aforementioned exception or a non existing orbit (e.g. a non valid value for the true anomaly). What could be the problem here? I have followed the Satellite Laser Ranging Orbit Determination written by @yzokras so I guess the problem is in my extension of the of the PythonAbstractMeasurement classā€¦

Here it is the extended class: the difference between my class and Range is that the measurement consists in the distance between the satellite and the Solar System Barycentre (SSB) along a given unit vector n . Since I am considering an ICRF frame (centered on the SSB), I have defined my measurement as numpy.dot(sat_position, n) = n1*x + n2*y + n3*z i.e. the dot product between the satellite position and n. Notice that the vector of derivatives of the measurement with respect to x, y, z, vx, vy and vz is [n1, n2, n3, 0, 0, 0] i.e. the three velocity components do not impact on the ā€œistantaneousā€ range measurement.

    class NewRange(PythonAbstractMeasurement):
    
        def __init__(self, date, observed, sigma, base_weight, satellite):
    
            super().__init__(date, observed, sigma, base_weight, satellite)
    
    
        def theoreticalEvaluation(self, iteration, evaluation, states):
                 
            state = states[0]
            
            pvaDS = super().getCoordinates(state, 0, 6)
    
            n = [np.sin(90-22.0144)*np.cos(83.6292), np.sin(90-22.0144)*np.sin(83.6292), np.cos(90-22.0144)]
            
            sat_position = [state.getPVCoordinates().getPosition().getX(), state.getPVCoordinates().getPosition().getY(), state.getPVCoordinates().getPosition().getZ()]
            
            pSSB = FieldVector3D(Gradient.constant(6, 0.0), Gradient.constant(6, 0.0), Gradient.constant(6, 0.0))        
            vSSB = FieldVector3D(Gradient.constant(6, 0.0), Gradient.constant(6, 0.0), Gradient.constant(6, 0.0))     
            aSSB = FieldVector3D(Gradient.constant(6, 0.0), Gradient.constant(6, 0.0), Gradient.constant(6, 0.0))
            
            pvaSSB = TimeStampedFieldPVCoordinates(state.getDate(), pSSB, vSSB, aSSB)
    
    
            # CREATE THE ESTIMATEDMEASUREMENT OBJECT
            
            new_range = NewRange(self.getDate(), float(self.getObservedValue()[0]), float(self.getTheoreticalStandardDeviation()[0]), float(self.getBaseWeight()[0]), self.getSatellites())        
            estimated = EstimatedMeasurement(new_range, iteration, evaluation, [state], [pvaDS.toTimeStampedPVCoordinates(), pvaSSB.toTimeStampedPVCoordinates()])
    
    
            # SET THE ESTIMATED VALUE FOR THE MEASUREMENT  
            estimated.setEstimatedValue(float(np.dot(sat_position, n)))        
            
            # SET THE DERIVATIVE OF THE ESTIMATED MEASUREMENT WITH RESPECT TO THE CURRENT STATE (CARTESIAN)        
            derivative = [float(n[0]), float(n[1]), float(n[2]), 0.0, 0.0, 0.0]
            derivative_jarray = JArray('object')(1)
            derivative_jarray[0] = JArray('double')(derivative)
            
            estimated.setStateDerivatives(0, derivative_jarray)
            
    
            return estimated

Do you see something strange in my code? Another thing I was thinking about is that maybe I do not have enough measurements (the sigma here is very large, between 1 and 3 km) or maybe I have to divide the range_weight variable (used in the laser ranging example ) by the number of observations? I havenā€™t understood if it is done automatically or notā€¦

Thank you in advance.

Samuele

P.S. I am giving the true orbit as a guess orbit in these first experimentsā€¦

20

Hello,

I have run my code many times varying all of the parameters of interest: position scale, convergence threshold, number of iterations, number of evaluations and number of measurements. I am of the idea that my class NewRange works correctly and that the problem is that I do not have enough accurate measurements to determine the satellite orbit. My measurements have in fact an error of about 2 km and the semimajor axis here is very long (the period of the elliptical orbit is of 7 days). How do you suggest me to set the above parameters? I was thinking to put a quite high convergence threshold (about 1?). What do you think?

Many thanks as always,
Samuele

21

Hello @samuele

Have you tried debugging the code to see how the matrix and vector of the Gauss-Newton problem look ?

The actual errors of the measurements should not lead to a singular problem since they only impact the residuals computation.
The errors in input (the standard deviations you give for your measurements) do intervene in the computation of the Jacobian (derivatives matrix) but only as weights, so they should not lead to a singular problem neither.

The ā€œsingular problemā€ comes from the fact that your system cannot be decomposed with the QR decomposition in the Gauss-Newton algorithm.
It means that the pseudo-inverse of the Jacobian J of your problem cannot be computed.
Thereā€™s a threshold in the Hipparchus Gauss-Newton class to check if the problem is singular but it is fixed to 1e-11, so you cannot play around with that and see the effect.

The Jacobian J is such as: J_{i,j}= \frac{\partial m_i}{\partial y_{0j}}
Where m_i is the ith measurement and y_{0j} the jth component of the initial orbit y_0. For example if y_0 is in Cartesian coordinates: y_{01}=X_0, y_{02}=Y_0... and so on.
This matrix is decomposed as: J_{i,j}= \frac{\partial m_i}{\partial y_{0j}} = \sum_{k=1}^{k=6}{\frac{\partial m_i}{\partial y_{k}}.\frac{\partial y_k}{\partial y_{0j}}}
Where y_k is the state vector evaluated at measurementā€™s date.

In your case, y_k is states[0] in the theoreticalEvaluation method.
And the \frac{\partial m_i}{\partial y_{k}} are set in your derivative array.

The \frac{\partial y_k}{\partial y_{0j}} are the components of the so-called state transition matrix that is computed through integration along with the orbit.

In your case \frac{\partial m_i}{\partial y_{k}} seems to be a constant: \frac{\partial m_i}{\partial y_{k}} = [n_x, n_y, n_z, 0, 0, 0] = [\vec{n}, \vec{0}].

So Iā€™m afraid that if the measurements are too close in time and the orbit model vary slowly (which seems to be the case since the period is 7 days), then the components of the state transition matrix \frac{\partial y_k}{\partial y_{0j}} do not change much between two measurements, and, since the \frac{\partial m_i}{\partial y_{k}} are constant, the two lines of the Jacobian are almost identical, which is a classic cause of singular problem (non pseudo-inversible matrix) exception.

That being said, I may be wrong, itā€™s just a feeling since I donā€™t know the actual data.

Could you explain us a bit more what is behind your vector \vec{n} and why it can be considered constant in your case ?
Can you try using measurements that are further apart (say with 90Ā° anomaly separation from each other) to see what you get ?

I suggest you try to use the automatic differentiation feature of Orekit to compute the derivatives of your measurement. It may help having different derivatives for different measurements. And could also show that there may actually be a component on the velocity that wasnā€™t intuitively easy to apprehend.

Last question, I donā€™t see any light travel-time correction in your code, is that normal ? Iā€™m guessing it must already be corrected somewhere upstream otherwise you wouldnā€™t manage to get some orbit determinations to work :wink:

Hope this helps,
Maxime