Skip to content
This repository has been archived by the owner on Oct 14, 2023. It is now read-only.

Lambert Izzo fails to converge under certain conditions, while Vallado works #1362

Closed
MLopez-Ibanez opened this issue Oct 23, 2021 · 15 comments

Comments

@MLopez-Ibanez
Copy link
Contributor

馃悶 Problem

A Lambert maneuver between two orbits fails to converge but it is unclear what is wrong.

The error seems very fickle and trying to manually recreate the orbits seems to avoid the issue, so I saved the orbits and the error to a pickle file. The error can be reproduced with this short code:

from poliastro.maneuver import Maneuver
from poliastro.bodies import Body
from poliastro.twobody import Orbit
from astropy import units as u
import pickle
import gzip
def from_pickle(filename):
    with gzip.open(filename, 'rb') as f:     
        return pickle.load(f)
        
db = from_pickle("two_shot_transfer_proposalj6y1by8g.pkl.gz")
from_orbit = db['from_orbit']
to_orbit = db['to_orbit']
e = db['e']

MU = 1.32712440018e11 # km^3/s^2

Sun = Body(
        parent=None,
        k = MU * u.km ** 3 / u.s ** 2,
        name="Sun")

man = Maneuver.lambert(from_orbit, to_orbit)

And this pickle file:

two_shot_transfer_proposalj6y1by8g.pkl.gz

Relevant packages:

# Paste your output here:
astropy==4.3.1
astroquery==0.4.3
numba==0.54.1
numpy==1.20.0
pandas==1.2.4
scipy==1.7.1
joblib==1.1.0
@jorgepiloto
Copy link
Member

jorgepiloto commented Oct 23, 2021

It looks like the maximum number of iterations has been reached. Try to modify this line as follows, @MLopez-Ibanez:

man = Maneuver.lambert(from_orbit, to_orbit, numiter=50)

which seems to converge 馃憤馃徑

@jorgepiloto
Copy link
Member

I am a bit surprised about Izzo's solver using more than 35 iterations, which is the default amount defined in poliastro. The transfer orbit is around ecc=0.30, not even a near-parabolic one... 馃

@MLopez-Ibanez
Copy link
Contributor Author

I am a bit surprised about Izzo's solver using more than 35 iterations, which is the default amount defined in poliastro. The transfer orbit is around ecc=0.30, not even a near-parabolic one... thinking

Well, this is what I want to understand. If for some reason it is some kind of unrealistic maneuver (e.g., huge velocity needed), then I'd like to bail out early rather than do more iterations.

@MLopez-Ibanez
Copy link
Contributor Author

I am a bit surprised about Izzo's solver using more than 35 iterations, which is the default amount defined in poliastro. The transfer orbit is around ecc=0.30, not even a near-parabolic one... thinking

I can easily generate dozens of examples of this error if you wish in case it helps debugging it.

@MLopez-Ibanez
Copy link
Contributor Author

This one fails even with numiter=1000:

two_shot_transfer_proposal1uklzh62.pkl.gz

I cannot see anything strange in the orbits:

In [2]: from_orbit
Out[2]: 2 x 4 AU x 13.1 deg orbit around Sun (None) at epoch 96884.93509314198 (UT1)
In [3]: to_orbit
Out[3]: 2 x 4 AU x 9.5 deg orbit around Sun (None) at epoch 97121.69041726712 (UT1)

@MLopez-Ibanez
Copy link
Contributor Author

I tried also without JIT to see if that could be the problem and no, it still fails.

@astrojuanlu
Copy link
Member

I confirm I can reproduce the error. The latest example doesn't converge even with 1_000_000 iterations. Working on it.

@astrojuanlu
Copy link
Member

Treating this as a bug, because if the transfer is impossible or some precondition is not met, it should be clearly communicated.

@astrojuanlu
Copy link
Member

tl;dr: @MLopez-Ibanez these are valid maneuvers with existing solutions, but our Izzo algorithm chokes on them. However, our Vallado algorithm works:

In [21]: from poliastro.iod.vallado import lambert as lambert_vallado

In [22]: Maneuver.lambert(from_orbit, to_orbit, numiter=1_000_000, method=lambert_vallado)
Out[22]: Number of impulses: 2, Total cost: 44.439604 km / s

More progress: the iteration starts with an initial guess of

-> x = _householder(x_0, T, ll, M, rtol, numiter)                                                                      
(Pdb) p x_0                                                                                                                                                                                                                                   
-0.48254756857234216 

and then it gets stuck in a loop like this

> /home/juanlu/Projects/LSF/poliastro/library/src/poliastro/core/iod.py(450)_householder()
-> if abs(p - p0) < tol:
(Pdb) p p, p0
(-3.6121719701854555, -0.5144239518586589)
(Pdb) c
> /home/juanlu/Projects/LSF/poliastro/library/src/poliastro/core/iod.py(450)_householder()
-> if abs(p - p0) < tol:
(Pdb) p p, p0
(-0.5144137836176683, -3.6121719701854555)
(Pdb) c
> /home/juanlu/Projects/LSF/poliastro/library/src/poliastro/core/iod.py(450)_householder()
-> if abs(p - p0) < tol:
(Pdb) p p, p0
(-3.6121652051680386, -0.5144137836176683)
(Pdb) c
> /home/juanlu/Projects/LSF/poliastro/library/src/poliastro/core/iod.py(450)_householder()
-> if abs(p - p0) < tol:
(Pdb) p p, p0
(-0.5144151669218648, -3.6121652051680386)

Jumping back and forth between the two same points.

Clearly there is a sign change in between:

(Pdb) !p0 = -0.51; _tof_equation_y(p0, _compute_y(p0, ll), T0, ll, M)
3.43690807694944
(Pdb) !p0 = -3.6; _tof_equation_y(p0, _compute_y(p0, ll), T0, ll, M)
-1.2575289342858709

But it turns out there's an asymptote at p0 = -1.0:

(Pdb) !p0 = -0.999; _tof_equation_y(p0, _compute_y(p0, ll), T0, ll, M)
35149.06051531792
(Pdb) !p0 = -0.9999; _tof_equation_y(p0, _compute_y(p0, ll), T0, ll, M)
1110802.6712689085
(Pdb) !p0 = -0.99999; _tof_equation_y(p0, _compute_y(p0, ll), T0, ll, M)
35124335.71514625

By doing some more magic and plotting:

tof

The solution for this case seems to be p = 0.941...:

In [40]: p0 = 0.942; _tof_equation_y(p0, _compute_y(p0, ll), T0, ll, M)
Out[40]: -0.0004182242581900475

In [41]: p0 = 0.941; _tof_equation_y(p0, _compute_y(p0, ll), T0, ll, M)
Out[41]: 0.0002703459920012641

which is outside of the interval being tested by the Householder iteration.

I have lots of questions because I don't remember the mathematical details. Hopefully @jorgepiloto can shed some light.

  • Why does the iteration jump to the other side of the asymptote? This requires more detailed debugging.
  • Why doesn't the initial guess lead the iteration to converge to p0 = 0.941...?
  • I don't remember the details of _tof_equation_y. Does p0 = 0.941... have physical meaning?
  • Is there anything special about this TOF curve because of the parameters of the problem? Or is it an unfortunate choice of initial guess?
  • In particular: is this asymptote always there?
  • Can the initial guess be improved so that this problem doesn't happen?
  • Is this a problem in our implementation? Or is it a limitation of the original algorithm?

To reproduce:

from poliastro.core.iod import _tof_equation_y, _compute_y
import numpy as np
ll = -0.9095948400461084
T0 = 1.2075613451614733
M = 0
p_array = np.linspace(-4, 0, num=10_000)
tof_array = np.array([_tof_equation_y(p0, _compute_y(p0, ll), T0, ll, M) for p0 in p_array])

One thing I know though: the Vallado method works on these ones!

In [22]: Maneuver.lambert(from_orbit, to_orbit, numiter=1_000_000, method=lambert_vallado)
Out[22]: Number of impulses: 2, Total cost: 44.439604 km / s

In [23]: from poliastro.iod.izzo import lambert as lambert_izzo

In [24]: Maneuver.lambert(from_orbit, to_orbit, numiter=1_000_000, method=lambert_izzo)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-24-2294576cfeb1> in <module>
----> 1 Maneuver.lambert(from_orbit, to_orbit, numiter=1_000_000, method=lambert_izzo)

~/Projects/LSF/poliastro/library/src/poliastro/maneuver.py in lambert(cls, orbit_i, orbit_f, method, short, **kwargs)
    188 
    189         # Compute all possible solutions to the Lambert transfer
--> 190         sols = list(method(k, r_i, r_f, tof, **kwargs))
    191 
    192         # Return short or long solution

~/Projects/LSF/poliastro/library/src/poliastro/iod/izzo.py in lambert(k, r0, r, tof, M, numiter, rtol)
     44     sols = izzo_fast(k_, r0_, r_, tof_, M, numiter, rtol)
     45 
---> 46     for v0, v in sols:
     47         yield v0 << kms, v << kms

~/Projects/LSF/poliastro/library/src/poliastro/core/iod.py in izzo(k, r1, r2, tof, M, numiter, rtol)
    229     sigma = np.sqrt(1 - rho ** 2)
    230 
--> 231     for x, y in xy:
    232         V_r1, V_r2, V_t1, V_t2 = _reconstruct(
    233             x, y, r1_norm, r2_norm, ll, gamma, rho, sigma

~/Projects/LSF/poliastro/library/src/poliastro/core/iod.py in _find_xy(ll, T, M, numiter, rtol)
    272     for x_0 in _initial_guess(T, ll, M):
    273         # Start Householder iterations from x_0 and find x, y
--> 274         x = _householder(x_0, T, ll, M, rtol, numiter)
    275         y = _compute_y(x, ll)
    276 

~/Projects/LSF/poliastro/library/src/poliastro/core/iod.py in _householder(p0, T0, ll, M, tol, maxiter)
    452         p0 = p
    453 
--> 454     raise RuntimeError("Failed to converge")

RuntimeError: Failed to converge
> /home/juanlu/Projects/LSF/poliastro/library/src/poliastro/core/iod.py(454)_householder()
    450         if abs(p - p0) < tol:
    451             return p
    452         p0 = p
    453 
--> 454     raise RuntimeError("Failed to converge")

Therefore, there is either

a) A problem with our implementation of the Izzo method, or
b) A limitation of the Izzo method for certain kinds of orbits, still to be determined.

@astrojuanlu astrojuanlu changed the title Lambert Izzo Failed to converge but unclear why Lambert Izzo fails to converge under certain conditions, while Vallado works Oct 29, 2021
@MLopez-Ibanez
Copy link
Contributor Author

@astrojuanlu Thanks for the detailed analysis. For the time being, do you think it would be correct to switch to Vallado within Izzo if Izzo fails to converge? I'm worried that we are missing true solutions to the problem due to some implementation detail.

@jorgepiloto
Copy link
Member

The bug is located in the initial guess! 馃悶

Equation which applies for when T1 < T < T0 uses a log2 base algorithm. I forced the natural logarithm log because of the transformation applied by Izzo in his paper and got the same results as all the solvers in 0.2dev0 version of lamberthub.

else:
# This is the real condition, which is not exactly equivalent
# elif T_1 < T < T_0
x_0 = (T_0 / T) ** (np.log2(T_1 / T_0)) - 1

Might this be an error in the article? I will review the maths involved tonight but the error is clearly located in previous lines 馃憖

@jorgepiloto
Copy link
Member

Solutions comparing the modified poliastro izzo2015 solver (using np.log instead of np.log2) against all lamberthub solvers:

poliastro modified solver izzo2015 computed
v1 = array([19.59117337, -7.22478863,  2.45951706])
v2 = array([19.92496512, -4.1345331 ,  2.49925548])

lamebrthub solver gauss1809 computed
v1 = array([19.59116491, -7.22478688,  2.459516  ])
v2 = array([19.93456908, -4.08699595,  2.50042674])

lamebrthub solver battin1984 computed
v1 = array([19.59019738, -7.22458745,  2.45939464])
v2 = array([19.92400524, -4.13418268,  2.49913497])

lamebrthub solver gooding1990 computed
v1 = array([19.59117337, -7.22478863,  2.45951706])
v2 = array([19.92496512, -4.1345331 ,  2.49925548])

lamebrthub solver avanzini2008 computed
v1 = array([19.59117337, -7.22478863,  2.45951706])
v2 = array([19.92496512, -4.1345331 ,  2.49925548])

lamebrthub solver arora2013 computed
v1 = array([19.59115074, -7.22478396,  2.45951422])
v2 = array([19.92494286, -4.13452498,  2.49925268])

lamebrthub solver vallado2013 computed
v1 = array([19.59117295, -7.22478854,  2.45951701])
v2 = array([19.9249647 , -4.13453295,  2.49925542])

lamberthub solver izzo2015 failed to converge! (Same implementation as poliastro 0.15 version)

@astrojuanlu
Copy link
Member

Amazing @jorgepiloto! Can you try the modified solver with the case from the first comment as well? #1362 (comment) If it passes the current tests plus these new weird cases, we have a winner.

@MLopez-Ibanez yes, for now it should be enough to switch to Vallado if Izzo doesn't converge.

@astrojuanlu
Copy link
Member

Closed in #1371!

@astrojuanlu
Copy link
Member

I confirm that the initial guess for the latest case is now x_0 = 0.9350..., which is much closer to the actual one! 馃挭馃徑

I'm sure this has improved the performance of the algorithm as well! Too bad we're not measuring it 馃槄

Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Projects
None yet
Development

No branches or pull requests

3 participants