Izzo's Algorithm and Multi-Revolution Transfers

L-similarity, the x-variable, Householder iterations, and the 2N+1 solutions of an N-revolution problem

Noordwijk, the Netherlands

Since 2005 the European Space Agency's Advanced Concepts Team has run the Global Trajectory Optimisation Competition, in which teams design tours that visit dozens of asteroids. A single search can require billions of Lambert solutions. At those numbers the cost of one solve stops being a detail. Dario Izzo, a member of that team, went back to a 1969 NASA report by Lancaster and Blanchard and asked whether there was a variable in which every Lambert problem's time-of-flight curve looks nearly the same — so that one good initial guess would work for all of them. There was. The resulting algorithm converges in about two iterations, handles multiple revolutions naturally, and now sits inside poliastro, pykep, Orekit and many flight-dynamics codes.

Watch first · presented by Shoaib Aziz

The numbers in this film are the labs' own: with λ = 0.5 the anchors are T(0) = 1.480210 and T(1) = 0.583333, and the multi-revolution minima are T_min = 4.4763 (N = 1) and 7.6777 (N = 2); Curtis Example 5.2 becomes λ = 0.354841, T = 0.830301 and converges from x₀ = 0.624048 to 0.6194516820 in three Householder steps; the printed middle case of the initial guess gives −0.71 at T = T₁ instead of 1; 5,000 random problems average 2.92 iterations (worst 5, one needing more than 4); and the 400 → 800 km rendezvous, 150° ahead in 5.5 h, has seven solutions costing 13.5, 10.9, 9.2, 7.8, 5.9 and 2.05 km/s and 257 m/s, against a Hohmann bound of 217 m/s.

Non-dimensionalisation (L-similarity). Lambert's theorem says the problem depends on

s

c

and

\Delta t

. Scale lengths by

s

and time by

\sqrt{s^3/2\mu}

and only two numbers are left:

\lambda^2 = 1 - \frac{c}{s}, \qquad T = \sqrt{\frac{2\mu}{s^3}}\,\Delta t

λ ∈ (−1, 1) carries the geometry; it is negative for long-way transfers. T is the non-dimensional flight time. Every Lambert problem with the same (λ, T) has the same solution up to scale.

The unknown is the Lancaster–Blanchard variable

x

, defined through

1 - x^2 = a_m/a

. It sweeps every conic on one axis:

-1 < x < 1

ellipses (

x = 0

is the minimum-energy ellipse),

x = 1

the parabola,

x > 1

hyperbolas. With

y = \sqrt{1 - \lambda^2(1 - x^2)}

the time-of-flight equation is

T(x) = \frac{1}{1 - x^2}\left[\frac{\psi + N\pi}{\sqrt{|1 - x^2|}} - x + \lambda y\right]

with cos ψ = xy + λ(1 − x²) for ellipses and cosh ψ = xy − λ(x² − 1) for hyperbolas, and N the number of complete revolutions. Near x = 1 the bracket is 0/0; there Battin's hypergeometric series 2T = η³·(4/3)·₂F₁(3, 1; 5/2; S₁) + 4λη is used instead.

💡Anchor values

Two points on the curve are known in closed form, for every

\lambda

T(0) = T_0 = \arccos\lambda + \lambda\sqrt{1-\lambda^2} + N\pi \qquad T(1) = T_1 = \tfrac{2}{3}(1 - \lambda^3)

T_0

is the minimum-energy time;

T_1

is the parabolic time (Euler's equation, non-dimensionalised). Between and beyond them the zero-revolution curve falls monotonically, with two oblique asymptotes on a log scale. Izzo inverts a piecewise approximation of those asymptotes to get the initial guess.

Loading interactive…

Initial guess (zero revolutions).

$T \ge T_0$ : $x_0 = (T_0/T)^{2/3} - 1$
$T < T_1$ : $x_0 = \tfrac{5}{2}\,T_1(T_1 - T)/[T(1 - \lambda^5)] + 1$
otherwise: $x_0 = 2^{\,\ln(T/T_0)/\ln(T_1/T_0)} - 1$

Iteration. With

f(x) = T(x) - T^*

and derivatives from Izzo's recurrences,

(1-x^2)T' = 3Tx - 2 + 2\lambda^3 x/y

and two more like it, a third-order Householder step

x_{n+1} = x_n - f\,\frac{f'^2 - f f''/2}{f'(f'^2 - f f'') + f''' f^2/6}

Cubic convergence from an initial guess that is already close: about three steps from guess to machine precision.

⚠A misprint in the paper

The middle case of the initial guess is printed in the paper as

(T_0/T)^{\log_2(T_1/T_0)} - 1

. Check the seams: at

T = T_0

it gives 0, as it should, but at

T = T_1

it does not give 1. Open-source implementations use the corrected form above, which interpolates

x

from 0 to 1 in log-time. In the lab, the printed version still converges on the textbook case, but in four iterations from the wrong side; the corrected one does it in three, and across 5,000 random problems never needs more than five. Reading a paper critically includes checking its formulas at the boundaries.

Velocities. With

\gamma = \sqrt{\mu s/2}

\rho = (r_1 - r_2)/c

and

\sigma = \sqrt{1 - \rho^2}

, the radial and transverse components come out directly:

V_{r1} = \gamma\frac{(\lambda y - x) - \rho(\lambda y + x)}{r_1},\quad V_{r2} = -\gamma\frac{(\lambda y - x) + \rho(\lambda y + x)}{r_2},\quad V_{t} = \gamma\sigma(y + \lambda x)

with

\vec v_1 = V_{r1}\hat r_1 + (V_t/r_1)\hat t_1

and

\vec v_2 = V_{r2}\hat r_2 + (V_t/r_2)\hat t_2

💡Multiple revolutions: 2N + 1 solutions

Allow the spacecraft to complete

N

full revolutions before arriving and the curve

T_N(x)

is no longer monotonic: it has a minimum

T_{\min}(N)

. Below it, no

N

-revolution transfer exists. Above it, the horizontal line

T = T^*

cuts the curve twice — a left branch (lower energy, longer period) and a right branch. So a problem with flight time long enough for

N

revolutions has up to

2N + 1

solutions: the zero-revolution arc plus a pair for each

N

T_{\min}(N)

is found by Halley's method on

T'(x) = 0

starting from

x = 0

; the two branch initial guesses are

x_{0,L} = \frac{q_L - 1}{q_L + 1},\ q_L = \Big(\frac{(N+1)\pi}{8T}\Big)^{2/3}; \qquad x_{0,R} = \frac{q_R - 1}{q_R + 1},\ q_R = \Big(\frac{8T}{N\pi}\Big)^{2/3}

📐Worked Example: seven ways to rendezvous

A chaser in a 400 km circular orbit must arrive at a point 800 km up and 150° ahead, in exactly 5.5 h (

\lambda = 0.132

). Minimum times: one revolution needs at least 2.28 h, two need 3.88 h, three need 5.47 h. At 5.5 h all three are feasible, so there are seven solutions. Total rendezvous

\Delta v

(departure plus matching burn): zero revolutions 13.5 km/s;

N=1

: 10.9 and 9.2 km/s;

N=2

: 7.8 and 5.9 km/s;

N=3

left 2.05 km/s,

N=3

right 257 m/s. The cheapest transfer is the one that circles the Earth three times on a nearly tangent orbit — close to the 217 m/s Hohmann bound. A zero-revolution solver would have reported 13.5 km/s and moved on.

💡Where multi-revolution Lambert matters

Low-Earth-orbit rendezvous (several orbits of phasing), Earth–Earth and Venus–Venus resonant legs in gravity-assist tours, asteroid tours, and low-thrust trajectory design (where multi-revolution Lambert arcs serve as initial guesses). Interplanetary transfers between different planets are mostly zero-revolution; Type III and IV transfers (one revolution plus) exist but are rarely competitive.

Method	Unknown	Strengths	Used in
Lagrange / bisection on a	a	Transparent; good teaching tool	Textbooks
Universal variables (BMW, Curtis)	z	One equation for all conics	Courses, much legacy software
Gooding (1990)	x (Lancaster–Blanchard)	Very robust, multi-rev, Halley	NASA/ESA tools, GMAT
Izzo (2015)	x	Fast (~2–3 iterations), multi-rev	poliastro, pykep, Orekit
Russell, ivLam (2019)	interpolated	Fastest; vectorised	Large-scale searches

↓But there is a problem

Three iterations, every conic, every revolution count. A Lambert call is now cheap enough to put inside a loop of a million. So put it inside a loop. Two spacecraft in Earth orbit: when should the chaser leave, and how long should it fly? Each choice is a Lambert problem, and each has a price. Mapping all of them turns targeting into reading a chart.