Halley convergence fixed for IodGooding

Orekit development
1

Hi everyone, I wanted to open this topic in order to get your opinion on the fixes I made to the IodGooding implementation (to fix the bug in the Halley convergence method). Basically, I followed the typical method and compared it with what was already implemented. I noticed several issues in the original code and made the following changes:

                final double Fc = P.dotProduct(C);
                final double Gc = EN.dotProduct(C);

I uncommented the G function and added it to the compute derivative function

    private void computeDerivatives(final Frame frame,
                                    final double x, final double y,
                                    final Vector3D lineOfSight1, final Vector3D lineOfSight3,
                                    final Vector3D Pin,
                                    final Vector3D Ein,
                                    final double F,
                                    final double G,                               
                                    final double T13, final double T12,
                                    final boolean withHalley,
                                    final int nRev,
                                    final boolean direction,
                                    final double[] FD, final double[] GD)

In the Halley method section, I found errors in the second-order derivative definitions.

            final double Fxx = (Fp1 + Fm1 - 2 * F) / hrho1Sq;
            final double Gxx = (Gp1 + Gm1 - 2 * F) / hrho1Sq;
            final double Fyy = (Fp3 + Fp3 - 2 * F) / hrho3Sq;
            final double Gyy = (Gm3 + Gm3 - 2 * F) / hrho3Sq;

After changes, it looks like this (replace F by G when computing Gxx and Gyy; in Fyy, change the second Fp3 to Fm3, and same for Gyy).

            final double Fxx = (Fp1 + Fm1 - 2 * F) / hrho1Sq;
            final double Gxx = (Gp1 + Gm1 - 2 * G) / hrho1Sq;                       
            final double Fyy = (Fp3 + Fm3 - 2 * F) / hrho3Sq;
            final double Gyy = (Gp3 + Gm3 - 2 * G) / hrho3Sq;

In the Fxy, and Gxy definitions we have :

final double Fxy = (Fp13 + Fm13) / (2 * dx * dy) - 0.5 * (Fxx * dx / dy + Fyy * dy / dx) - F / (dx * dy);
final double Gxy = (Gp13 + Gm13) / (2 * dx * dy) - 0.5 * (Gxx * dx / dy + Gyy * dy / dx) - F / (dx * dy);

I replaced the F function by G in Gxy :

final double Gxy = (Gp13 + Gm13) / (2 * dx * dy) - 0.5 * (Gxx * dx / dy + Gyy * dy / dx) - G / (dx * dy);

Finally in the terms of the jacobian we originally have :

            final double FxH = Fx + 0.5 * (Fxx * dx3NR + Fxy * dx1NR);
            final double FyH = Fy + 0.5 * (Fxy * dx3NR + Fxx * dx1NR);
            final double GxH = Gx + 0.5 * (Gxx * dx3NR + Gxy * dx1NR);
            final double GyH = Gy + 0.5 * (Gxy * dx3NR + Gyy * dx1NR);

But in the FyH definition we should have :
final double FyH = Fy + 0.5 * (Fxy * dx3NR + Fyy * dx1NR); (Change the Fxx for a Fyy)

After making those changes, I changed the boolean of the Halley method to be true at any iteration in order to only use Halley for convergence.

I also introduced an error in rho_init to have more than one or two iterations.

        // rho1 = rho1init;
        // rho3 = rho3init;
        rho1 = rho1init*6.0 ;
        rho3 = rho3init*6.0 ;

I gave three pre-computed RaDec (generated with an initial position of the target and the observer). The true position of the considered target is : x = -26649.115, y = 41027.775,z = -5294.460, x_dot = -2.993, y_dot = -1.943, z_dot = 0.267 (GCRF)

With only halley convergence and without any changes I obtained :

iter=0 rho1=4.6110494091158705 rho3=4.968694812279726 Fc=0.08130464752872525 stop=1.0E-13 detj=-1.8412824579164062E10 [H]
iter=1 rho1=4.609833976420186 rho3=4.970016431157314 Fc=0.08275745857588751 stop=0.016709747968573017 detj=4.5117283069085156E11 [H]

[...]

iter=99 rho1=4.638688189311865 rho3=4.938169143471605 Fc=0.08419265237286516 stop=0.017319291348148164 detj=9.270718313762109E11 [H]
=== 100 iter | crit=0.017312895545320228 | converge=false ===
SpacecraftState{orbit=Cartesian parameters: {P(-1.3157182332905114E8, 1.7937535613075933E8, -4.742736253553751E7), V(-108.16833248624134, 48.41570984076883, -188.23927712334591)}, attitude=org.orekit.attitudes.Attitude@78aab498, mass=1000.0, additional={}, additionalDot={}}

Basically the algorithm reaches 100 iterations (max iteration , so it does not converge) , after that i tested with only newton convergence (setting halley booleans to false):

iter=0 rho1=0.9199352246773724 rho3=1.0032745735207915 Fc=0.08130464752872525 stop=1.0E-13 detj=5.0865838066984505E-5 [N]
iter=1 rho1=0.9164790066979406 rho3=1.000035273457973 Fc=8.289453314661068E-5 stop=0.016709747968573017 detj=5.057591176035089E-5 [N]
iter=2 rho1=0.9164788018969453 rho3=1.0000350705179064 Fc=4.766519784284826E-9 stop=7.917761787949509E-5 detj=5.057217296105787E-5 [N]
iter=3 rho1=0.9164788018969446 rho3=1.0000350705179106 Fc=8.235447725389051E-17 stop=4.567270813927941E-9 detj=5.057217205322497E-5 [N]
=== 4 iter | crit=7.89119295009951E-17 | converge=true ===
SpacecraftState{orbit=Cartesian parameters: {P(-2.654762526549285E7, 4.088969679165287E7, -5266101.981770694), V(-2818.604309042942, -2176.784785455636, 361.48403445103736)}, attitude=org.orekit.attitudes.Attitude@78aab498, mass=1000.0, additional={}, additionalDot={}}

Four iteration and the algorithm converged to a solution which is not too far from the real position of the target, so the algorithm can converge.

And finally , I used only halley convergence with the fixes described above and obtained :

iter=0 rho1=0.9182210087747658 rho3=1.0806785155741272 Fc=0.08130464752872525 stop=1.0E-13 detj=5.19574749174509E-5 [H]
iter=1 rho1=0.9171477122553526 rho3=1.0018455500424541 Fc=0.0010129347671791087 stop=0.016709747968573017 detj=5.134398033016589E-5 [H]
iter=2 rho1=0.9164787920945492 rho3=1.0000350412456398 Fc=2.533557768076061E-6 stop=9.342238476293698E-4 detj=5.0571682332857805E-5 [H]
iter=3 rho1=0.9164788018969866 rho3=1.0000350705180212 Fc=5.1772430911083636E-11 stop=2.4248037492153602E-6 detj=5.0572175791752496E-5 [H]
iter=4 rho1=0.916478801896857 rho3=1.0000350705176695 Fc=4.906942494961125E-16 stop=4.960826211236196E-11 detj=5.057217540277937E-5 [H]
=== 5 iter | crit=4.701824516887309E-16 | converge=true ===
SpacecraftState{orbit=Cartesian parameters: {P(-2.654762526580604E7, 4.088969679140569E7, -5266101.981731054), V(-2818.6043085247147, -2176.784785039024, 361.4840343835512)}, attitude=org.orekit.attitudes.Attitude@78aab498, mass=1000.0, additional={}, additionalDot={}}

In this case, the algorithm converged around the same solution as Newton using only Halley. From what I saw, the code should be good, but I am not sure, and I would appreciate your opinions/tests on the question. I put the full modified Java code.

Hugo B

IodGooding.java (35.0 KB)

2

Hi Hugo,

Thanks you for your analysis.
It would make sense if Haley fails because of wrong second order derivatives.

@markrutten this will be of interest to you.

Cheers,
Romain.

3

Hi @HugoBocc2 thanks for the tag @Serrof, I’m definitely interested!

Those are really encouraging results @HugoBocc2, thank you very much for taking the time to correct the Halley equations.

The way that the Gooding algorithm was originally written it was intended to run for 50 iterations with standard Newton-Raphson, then transition to Halley after that. In my testing (you can see a bunch of new tests in my merge request), I couldn’t find any examples where the algorithm failed due to too many iterations - failure is always due to poor geometry and the Lambert solver fails (which is not in any way a criticism of the Lambert solver). I completely removed the Halley part of the algorithm because it was never being called.

Unless we can show its needed by one of the test examples I don’t think there’s a justification for a transition from Newton to Halley. Would you suggest just replacing Newton with your revised Halley algorithm?

4

Hi,

Thank you for your answer. I am currently investigating the convergence behaviour of Halley’s method versus Newton’s method across a wide range of scenarios, in order to gain a clearer understanding on the subject and also to verify a bit more my previous statement.

From my preliminary results, it appears that Newton’s method tends to be more resilient , whereas Halley’s method may offer faster convergence. I will keep you updated as soon as I finish this study.

Hugo B

5

For what it’s worth, Hipparchus provides a very robust and fast method with variable order to find roots without using derivatives. This is what we use for events detection for example, generally with maxOrder set to 5. It is implemented in the BracketingNthOrderBrentSolver class.
But of course, it has nothing to do with Halley or Newton-Raphson, so it is different from the original Gooding algorithm.

6

I made a small analysis over typical scenario for IoD ground based and I summarized it in the following document. I also checked the order of the convergence and it seems indeed to be cubic as expected. For the results I am a little bit concerned because Halley method converged everytime but with more wrong results than Newton’s. I don’t know if I’ll be able to look at what you proposed @luc but i’ll try if I can. If you have any questions or comments on the analysis feel free to ask (I used a keplerian propagator by the way). The code used is the following one :

def run_iod(orbit_true, dt_triplet, rho_factor):

    prop = KeplerianPropagator(orbit_true)
    date1 = date0.shiftedBy(dt_triplet[0])
    date2 = date0.shiftedBy(dt_triplet[1])
    date3 = date0.shiftedBy(dt_triplet[2])

    state1 = prop.propagate(date1)
    state2 = prop.propagate(date2)
    state3 = prop.propagate(date3)

    O1 = station_topo.getPVCoordinates(date1, gcrf).getPosition()
    O2 = station_topo.getPVCoordinates(date2, gcrf).getPosition()
    O3 = station_topo.getPVCoordinates(date3, gcrf).getPosition()

    los1 = state1.getPVCoordinates().getPosition().subtract(O1).normalize()
    los2 = state2.getPVCoordinates().getPosition().subtract(O2).normalize()
    los3 = state3.getPVCoordinates().getPosition().subtract(O3).normalize()

    rho1_true = state1.getPVCoordinates().getPosition().subtract(O1).getNorm()
    rho3_true = state3.getPVCoordinates().getPosition().subtract(O3).getNorm()

    iod = IodGooding(mu)

    try:

        orbit_est = iod.estimate(
            gcrf,
            O1, O2, O3,
            los1, date1,
            los2, date2,
            los3, date3,
            rho1_true * rho_factor,
            rho3_true * rho_factor

        )

small_analysis.pdf (261.0 KB)