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

MaximeJ · December 2, 2021, 4:32pm

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

samuele · December 3, 2021, 11:54am

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

MaximeJ · December 3, 2021, 1:08pm

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

samuele · December 3, 2021, 6:01pm

Hi @MaximeJ,

Here it is what I have done:

I collect n measurements r_1, … r_n at the instants t_1, … t_n respectively.
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()).
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).
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)
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.
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

MaximeJ · December 3, 2021, 7:05pm

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

samuele · January 31, 2022, 12:02pm

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

MaximeJ · January 31, 2022, 1:07pm

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

samuele · February 4, 2022, 4:41pm

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

MaximeJ · February 4, 2022, 5:03pm

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.

petrus.hyvonen · February 5, 2022, 12:35pm

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

Regards!

samuele · February 7, 2022, 5:08pm

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…

samuele · February 9, 2022, 8:51am

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

MaximeJ · February 9, 2022, 2:52pm

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

Hope this helps,
Maxime

samuele · February 14, 2022, 2:33pm

Hi @MaximeJ,

thank you a lot for your helpful reply. Now I will try to explain you a bit more what my code is about.

It is a code about pulsar navigation, a field which is quite complex to describe in a few lines: our aim is to exploit the periodic signals coming from pulsar stars to perform an orbit determination. Pulsars are very far away and their line-of-sight directions can be easily considered to be constant in the Solar System.
The phase (e.g. that of the signal maximum) of the received signal depends on the position of the satellite along the line-of-sight of the source and refers to a specific point on the satellite orbit. All of the successive received periods (each about 10 ms let’s say) must be adjusted during the observation shifting the phase of the successive received periods to match the phase of the signal received at the specific point, in which you estimate the distance between this specific point and the SSB. WIthout such “folding” (or sum) procedure, the SNR of the received signal is not sufficient. Then, you make your SSB-SAT distance estimation at the end of your signal integration, referring to the beginning of this orbit segment (i.e. at the specific point). Thus, the theoretical value of this distance is always a function of the sole specific point and does not depend on the satellite velocity. Of course, the derivative of the current state y with respect to the initial state y0 will depend on the initial velocity! But as far as I understand this is already considered in AbstractBatchLSModel.fetchEvaluatedMeasurements that I had the opportunity to study in deep while re-writing the whole OD process.

Indeed, a pulsar-based navigation system does not pretend to achieve a precision of 10 m. The best result achieved so far is in fact a positioning accuracy of about 5 km. But someone like me is looking to boost such performance . Such low-accuracy measurements make my parameters setting a bit harder. I have put the estimation_position_scale = 1000 m, estimator_convergence_threshold = 1.0, and in the creation of the measurement I give a sigma = 2000.0 and a weight which is equal to 1/m, where m is the number of measurements. Do you think these values are reasonable?

For what concerns the Gauss-Newton problem, I do not find anymore a singular problem if I choose the QRDecomposer singularity threshold to be very small (about 1e-20) and the formNormalEquations parameter of the GaussNewtonOptimizer to be True (leading thus to the use of the pseudoinverse of the Jacobian). Are those too strong assumption in your opinion?
Setting those values, I performed an OD covering half of the orbital period (180° anomaly separation) in which I collected 50 measurements, with a sigma of about 2000 m each. The estimated orbit is very similar to the true one and if I make a plot of the two orbits together they are almost aligned. This “almost”, however, means a large number of meters. These are the distances (3D norm) between the real initial position and the initial position of the estimated orbit in the successive least squares iterations:

iteration #6 = 84 km
iteration #7 = 46 km
iteration #8 = 106768 km (the inclination drastically changed here…)
iteration #11 = 12879 km
iteration #12 = 12979 km
iteration #20 = 13244 km
iteration #26 = 12794 km

It seems in fact that the estimator is not able to converge… my class NewRange seems to be correct to me, it is very simple. Maybe that’s something wrong with the setting of the estimator parameters? Some hints would be very appreciated at this point

Let me know if there is something else I can report you that would help!

Best,
Samuele

MaximeJ · February 14, 2022, 4:33pm

Hi @samuele,

Yes you are right.

Maybe you could increase the estimation_position_scale to 5 to 10km since it’s probably the accuracy you will achieve at the end (according to " a positioning accuracy of about [5 km]").
I would set the weights to 1. This parameter is a bit misleading (and not very useful) in Orekit and we actually always set it to 1. The problem is already weighted by the sigmas of the measurements.
That being said, since your measurements all have the same weights it probably won’t change your results.

I think it is ok, normal equations should be equivalent to straight Jacobian pseudo-inverse and a too-high singularity threshold was probably the reason why you were locked with a “singular problem”.

Indeed it does not converge… Could you share your code ? It’s hard to find out the issue without being able to actually debug the code.

Best,
Maxime

samuele · February 16, 2022, 8:49am

Hi @MaximeJ,

I have tried to use the LevenbergMarquardtOptimizer() instead of the GaussNewton one and the results are better (about a order of magnitude). But the positioning accuracy still does not converge and the positioning error increases quite fast… this is the code I call every l measurements along the orbit. kepler is the propagator associated to my guess orbit, Oj is the array containing the observed measurements and dates_array contains the AbsoluteDates associated to the measurements. NewRange is the class I reported in this topic some days ago i.e. it is my extension of the class PythonAbstractMeasurement.

prop_min_step = 0.001 # s
prop_max_step = 0.03 # s
prop_position_error = 1.0 # m

# Orbit determinator parameters
estimator_position_scale = 10000.0 
range_weight = 1.0  

initialDate = dates_array[0]

# get the satellite state vector at the initial point 
x_guess = kepler.propagate(initialDate).getPVCoordinates().getPosition().getX()
y_guess = kepler.propagate(initialDate).getPVCoordinates().getPosition().getY()
z_guess = kepler.propagate(initialDate).getPVCoordinates().getPosition().getZ()
vx_guess = kepler.propagate(initialDate).getPVCoordinates().getVelocity().getX()
vy_guess = kepler.propagate(initialDate).getPVCoordinates().getVelocity().getY()
vz_guess = kepler.propagate(initialDate).getPVCoordinates().getVelocity().getZ()    

position = Vector3D(x_guess, y_guess, z_guess) # first entry of the state vector for the initial position (at the initialDate)
velocity = Vector3D(vx_guess, vy_guess, vz_guess) # first entry of the state vector for the initial velocity (at the initialDate)
pvCoordinates = PVCoordinates(position, velocity)  

reference_orbit = KeplerianOrbit(pvCoordinates, inertialFrame, initialDate, Constants.WGS84_EARTH_MU) # orbit definition line

integrator = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error) # args: min_step (s), max_step (s), position_error (m)

propagatorBuilder = NumericalPropagatorBuilder(reference_orbit, integrator, PositionAngle.MEAN, estimator_position_scale)

# Construct the estimator
optimizer = LevenbergMarquardtOptimizer()
estimator = BatchLSEstimator(optimizer, propagatorBuilder)

estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)

observableSatellite = ObservableSatellite(0)

for i in range(0, l):

    sat = ArrayList()
    sat.add(observableSatellite)
    
    psr_range = NewRange(dates_array[i], float(Oj[i]), float(sigma), float(range_weight), sat)
    estimator.addMeasurement(psr_range)


# Estimate the orbit (could take some time depending on l and on the orbit determinator parameters)
print('\n\n-----------------ESTIMATING -------------------')
estimatedPropagatorArray = estimator.estimate()
print('\n\n---------------- FINISHED --------------------\n\n')

estimatedPropagator = estimatedPropagatorArray[0]
estimatedInitialState = estimatedPropagator.getInitialState()
actualOdDate = estimatedInitialState.getDate()

estimatedOrbit_init = estimatedInitialState.getOrbit()

I can see that the propagator returned by the estimate() method is always associated with a minimization of the computed-measurement coefficients and maybe it is the accuracy I can achieve with such noisy measurements… I will try to play a bit with the number of the measurements etc to see what happens. However, let me know if you notice some errors in my code, it would help me a lot.

Thanks as always!
Samuele

MaximeJ · February 21, 2022, 10:44am

Hi @samuele,

Glad to see that you had more luck with Levenberg-Marquardt instead of Gauss-Newton.
This is bit weird, though, since LM mainly brings a bigger convergence radius with respect to GN.
Since you start of from the desired solution you shouldn’t need the bigger convergence radius that LM provides.

I don’t see any obvious error in the code you provided.
Still, if you start from the solution your OD should converge nicely.
You are using simulated measurements right? Have you tried simulating with smaller errors?

Maxime

samuele · February 28, 2022, 12:32pm

Hi @MaximeJ ,

I tried using simulated measurements with smaller errors and it was definitely the cause of this unwanted behaviour… so, thank you so much! This would have been very difficult without your help. I consider the topic’s question to be answered.

I have a last question for you but I don’t know if it’s better to ask it here or to create another topic. I put it here before selecting the “solution” in your replies because I don’t know if it’s still possible to reply on a closed topic. If you think it’s better to open another topic, I will create a new one

My question is: I would like to add the satellite clock’s bias and drift to my problem. Hence I would like to estimate these two parameters along with the 6 of the state. Mathematically speaking it is quite easy since the effects of these parameters lead to a redefinition of my range measurement in which I have to add a (bias + drift*t)*c component to the range measurement. But what is the easiest way to add them to my NewRange class I reported in this topic (21 days ago)? I imagine I will have to manage those parameters through the ParameterDriver class but I haven’t studied too much this class yet so maybe there is a more convenient way to add the clock parameters to my class…

All the best,
Samuele

MaximeJ · February 28, 2022, 2:58pm

Hi @samuele,

You’re very welcome. I should have started with this simple solution!

I think it is possible (never actually tried it though).

If you look into the Orekit Range class constructor you will see that the clock bias (called clock offset in Orekit) is added from the ObservableSatellite.
The clock drift works the same way.
So you will have to do the same in your NewRange class.
Then in the theoreticalEvaluation method you will have to set the derivatives of your clock bias and drift manually.
In Orekit Range.theoreticalEvaluation it is done within this loop.
There is no explicit computation of the derivatives since they are computed by automatic differentiation (here’s an article about automatic differentiation in Orekit if you’re interested).

All the best,
Maxime

samuele · February 28, 2022, 3:05pm

Hi @MaximeJ ,

OK, got it. I will follow what it is done for the clock offset in the Range class constructor!

Best,
Samuele