Hohmann Transfer via Constrained Optimization
Li Xie∗
State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources
North China Electric Power University, Beijing 102206, P.R. China

arXiv:1712.01512v1 [] 5 Dec 2017

Yiqun Zhang† and Junyan Xu‡
Beijing Institute of Electronic Systems Engineering, Beijing 100854, P.R. China

In the first part of this paper, inspired by the geometric method of Jean-Pierre Marec, we
consider the two-impulse Hohmann transfer problem between two coplanar circular orbits
as a constrained nonlinear programming problem. By using the Kuhn-Tucker theorem, we
analytically prove the global optimality of the Hohmann transfer. Two sets of feasible solutions
are found, one of which corresponding to the Hohmann transfer is the global minimum, and the
other is a local minimum. In the second part, we formulate the Hohmann transfer problem as
two-point and multi-point boundary-value problems by using the calculus of variations. With
the help of the Matlab solver bvp4c, two numerical examples are solved successfully, which
verifies that the Hohmann transfer is indeed the solution of these boundary-value problems.
Via static and dynamic constrained optimization, the solution to the orbit transfer problem
proposed by W. Hohmann ninety-two years ago and its global optimality are re-discovered.

I. Introduction

I

n 1925, Dr. Hohmann, a civil engineer published his seminal book [1] in which he first described the well-known
optimal orbit transfer between two circular coplanar space orbits by numerical examples, and this transfer is now

generally called the Hohmann transfer. Hohmann claimed that the minimum-fuel impulsive transfer orbit is an elliptic
orbit tangent to both the initial and final circular orbits. However a mathematical proof to its optimality was not
addressed until 1960s. The first proof targeted to the global optimality was presented by Barrar [2] in 1963, where the
Whittaker theorem in classical analytical dynamics was introduced and the components of the velocity at any point on
an elliptic orbit, perpendicular to its radius vector and the axis of the conic respectively, were used to be coordinates
in a plane; see also Prussing and Conway’s book [3, eq. (3.13)] in detail. Before that, Ting in [4] obtained the local
optimality of the Hohmann transfer. By using a variational method, Lawden investigated the optimal control problem
of a spacecraft in an inverse square law field, and invented the prime vector methodology [5]. According to different
thrusts, the trajectory of a spacecraft is divided into the arcs of three types: (1) null thrust arc; (2) maximum thrust arc;
∗ Professor,

School of Control and Computer Engineering; the correspondence author, lixie@ncepu.edu.cn
Research Scientist, yiqunzhang@hotmail.com
‡ Associate Research Scientist, junyan_Xu@sina.cn
† Senior

(3) intermediate arc. If we approximate the maximum thrust by an impulse thrust, then the orbit transfer problem can
be studied by the the prime vector theory and the Hohmann transfer is a special case, which can also be seen in the
second part of this paper. Lawden derived all necessary conditions that the primer vector must satisfy. A systemic
design instruction can be found in [6].
In 2009, Pontani in [7] reviewed the literature concerning the optimality of the Hohmann transfer during 1960s
and 1990s, for example [8–14]. Among of them, Moyer verified the techniques devised respectively by Breakwell
and Constensou via variational methods for general space orbit transfers. The dynamical equations involved were
established for orbital elements. Battin and Marec used a Lagrange multiplier method and a geometric method in light of
a hodograph plane, respectively. It is Marec’s method that enlightens us to consider the Hohmann transfer in a different
way. Based on Green’s theorem, Hazelrigg established the global optimal impulsive transfers. Palmore provided an
elemental proof with the help of the gradient of the characteristic velocity, and Prussing simplified Palmore’s method by
utilizing the partial derivatives of the characteristic velocity, and a similar argument also appeared in [15] which was
summarized in the book [16]; see more work of Prussing in [3, Chapter 6]. Yuan and Matsushima carefully made use of
two lower bounds of velocity changes for orbit transfers and showed that the Hohmann transfer is optimal both in total
velocity change and each velocity change.
Recently, Gurfil and Seidelmann in their new book [17, Chapter 15] consider the effects of the Earth’s oblateness J2
on the Hohmann transfer. Then the velocity of a circular orbit is calculated by
v
u
t
v=

2
3J2 req
µ
1+
r
2r 2

!

By using Lagrange multipliers, the optimal impulsive transfer is derived. This transfer is referred as to the extended
Hohmann transfer which degenerates to the standard Hohmann transfer as J2 = 0. Avendaño et al. in [18] present a pure
algebraic approach to the minimum-cost multi-impulse orbit-transfer problem. By using Lagrange multipliers, as a
particular example, the optimality of the Hohmann transfer is also provided by this algebraic approach. These authors
are all devoted to proving the global optimality of the Hohmann transfer.
In the first part of this paper, we present a different method to prove the global optimality of the Hohmann transfer.
Inspired by the geometric method of Marec in [10, pp. 21-32], we transform the Hohmann transfer problem into a
constrained nonlinear programming problem, and then by using the results in nonlinear programming such as the
well-known Kuhn-Tucker theorem, we analytically prove the global optimality of the Hohmann transfer. Here by the
global optimality, we mean that the Hohmann transfer is optimum among all possible two-impulse coplanar transfers.
Two sets of feasible solutions are found, one of which corresponding to the Hohmann transfer is the global minimum,
and the other is a local minimum. In the second part of the paper, we consider the Hohmann transfer problem as
dynamic optimization problems. By using variational method, we first present all necessary conditions for two-point and
2

multi-point boundary-value problems related to such dynamic optimization problems, and we then solve two numerical
examples with the help of Matlab solver bvp4c, which verifies that the Hohmann transfer is indeed the solution of these
constrained dynamic optimization problems. By formulating the Hohmann transfer problem as constrained optimization
problems, the solution to the orbit transfer problem proposed by W. Hohmann 92 years ago and its global optimality are
re-discovered.

II. The optimality of the Hohmann Transfer
The Hohmann transfer is a two-impulse orbital transfer from one circular orbit to another; for the background see,
e.g., [3, 10, 19, 20]. We use the same notation as Marec in [10, pp. 21-32]. Let O0 denote the initial circular orbit
of a spacecraft, and the radius and the velocity of the initial circular orbit are r0 and V0 respectively. Let O f denote
the final circular orbit, and its radius and velocity are r f and Vf respectively. For simplicity, let r f be greater than r0 .
We use boldface to denote vectors. The center of these two coplanar circular orbits is at the origin of the Cartesian
inertia reference coordinate system. Suppose all transfer orbits are coplanar to O0 . Hence we need only to consider the
x − y plane of the reference coordinate system. Assume that the spacecraft is initially located at the point (r0, 0) of the
initial circular orbit. At the initial time t0 , the first velocity impulse vector ∆V0 is applied, and its components in the
direction of the radius and the direction perpendicular to the radius are ∆X0, ∆Y0 respectively. Similarly, at the final time
t f , the second velocity impulse ∆V f occurs, and its components are ∆X f , ∆Yf . Hence at the time t0+ just after t0 , the
components of the velocity of the spacecraft are (∆X0, V0 + ∆Y0 ). Then in order to enter the final circular orbit, the
components of the velocity of the spacecraft must be (−∆X f , Vf − ∆Yf ) at the time t −f just before t f .
During the period from t0+ to t −f , the spacecraft is on a transfer orbit (conic) and the angular momentum and the
energy (per unit mass) are conserved, hence we have
h = r0 (V0 + ∆Y0 ) = r f (Vf − ∆Yf )
i µ
 µ
1
1h
E=
(V0 + ∆Y0 )2 + ∆X02 −
=
(Vf − ∆Yf )2 + ∆X f2 −
2
r0 2
rf

(1)

µ
, µ is gravitational constant. It follow from (1)
r
that we can express ∆X0, ∆Y0 in terms of ∆X f , ∆Yf and vice versa. In order to simplify equations, at first, using the
Notice that for a circular orbit, its radius and velocity satisfy V =

r

initial orbit radius and velocity as reference values, we define new variables:
r̄ f =

rf
,
r0

vf =

Vf
,
V0

y0 =

∆Y0
,
V0

x0 =

3

∆X0
,
V0

yf =

∆Yf
,
V0

xf =

∆X f
V0

(2)

Substituting (2) into (1) yields the following non-dimensional angular momentum and energy equalities
h = (1 + y0 ) = r̄ f (v f − y f )
i 1

1h
1
(1 + y0 )2 + x02 − 1 =
(v f − y f )2 + x 2f −
E=
2
2
r̄ f

(3)

see [10, p.22]. Then we express x f , y f in terms of x0, y0
y f = v f − (1 + y0 )r̄ f−1 = r̄ f−1/2 − (1 + y0 )r̄ f−1

(4)

x 2f = x02 + (1 + y0 )2 − 2(1 − r̄ f−1 ) − (v f − y f )2
= x02 + (1 + y0 )2 (1 − r̄ f−2 ) − 2(1 − r̄ f−1 )

(5)

from which we have

∆v 2f



Vf
=
v0

2

= x 2f + y 2f


2
= x02 + (1 + y0 )2 (1 − r̄ f−2 ) − 2(1 − r̄ f−1 ) + r̄ f−1/2 − (1 + y0 )r̄ f−1

2
= x02 + y0 + 1 − r̄ f−3/2 − (r̄ f − 1)(2r̄ f2 − r̄ f − 1)r̄ f−3

2
= x02 + y0 + 1 − r̄ f−3/2 − (r̄ f − 1)2 (2r̄ f + 1)r̄ f−3
∆v02 = x02 + y02
where we define ∆v0 =

q

x02 + y02 > 0 and ∆v f =

q

x 2f + y 2f > 0. Thus the cost functional (i.e., the characteristic

velocity) can be written as
∆v(x0, y0 ) = ∆v0 + ∆v f
=

q

r
x02

+

y02

+


2
x02 + y0 + 1 − r̄ f−3/2 − (r̄ f − 1)2 (2r̄ f + 1)r̄ f−3

Observing (5), it is noted that we must make the following constraint on the first impulse
x02 + (1 + y0 )2 (1 − r̄ f−2 ) − 2(1 − r̄ f−1 ) ≥ 0

(6)

That is, when we use the energy conservation to calculate the non-dimensional component of the second impulse x f , the
inequality (6) must hold such that a non-negative number is assigned to x 2f , which actually requires that the transfer
orbit intersects the inner circle and the outer circle.

4

The above background can be found in [10, Section 2.2] which is specified here. Marec used the independent
variables x0, y0 as coordinates in a hodograph plan and in geometric language, an elegant and simple proof was given
to show the optimality of the Hohmann transfer. Marec’s geometric method enlightens us to consider the Hohmann
transfer in a different way. Then the global optimality of the Hohmann transfer analytically appears.
We now formulate the Hohmann transfer problem as a nonlinear programming problem subject to the inequality
constraint (6) in which x0 and y0 are independent variables.
Theorem II.1 The classical Hohmann transfer is the solution of the following constrained optimization problem

min ∆v(x0, y0 )

(7)

x0,y0

s.t.

x02 + (1 + y0 )2 (1 − r̄ f−2 ) − 2(1 − r̄ f−1 ) ≥ 0

Also this solution is the global minimum.
Proof. We say that the Hohmann transfer is the solution to the above optimization problem if the pair (x0, y0 )
corresponding to the Hohmann transfer is feasible to (7). As usual, we use Lagrange multiplier method and define the
Lagrangian function as follows


F(x0, y0 ) = ∆v(x0, y0 ) + λ −x02 − (1 + y0 )2 (1 − r̄ f−2 ) + 2(1 − r̄ f−1 )

(8)

In view of the Kuhn-Tucker theorem, a local optimum of ∆v(x0, y0 ) must satisfy the following necessary conditions
∂F
x0
x0
=
+
− 2λx0 = 0
∂ x0 ∆v0 ∆v f

(9)

y0 + 1 − r̄ f−3/2
∂F
y0
=
+
− 2λ(1 + y0 )(1 − r̄ f−2 ) = 0
∂ y0 ∆v0
∆v f


λ −x02 − (1 + y0 )2 (1 − r̄ f−2 ) + 2(1 − r̄ f−1 ) = 0

(11)

− x02 − (1 + y0 )2 (1 − r̄ f−2 ) + 2(1 − r̄ f−1 ) ≤ 0

(12)

λ≥0

(13)

(10)

Notice that we have assumed that ∆v0, ∆v f > 0. In order to find a feasible solution, we divide the proof of the first part
of this theorem into two cases.
(i) If λ = 0, then due to the equations (9) and (10), we have x0 = 0 and
y0 + 1 − r̄ f−3/2
y0
+
=0
∆v0
∆v f
5

(14)

respectively. We claim that the equality (14) does not hold because the equality x0 = 0 implies that
y0
y0
= ±1
=
∆v0
|y0 |
which contradicts

1<

y0 + 1 − r̄ f−3/2
∆v f

y0 + 1 − r̄ f−3/2
=r
,

2
−3/2
−3
2
y0 + 1 − r̄ f
− (r̄ f − 1) (2r̄ f + 1)r̄ f

or

y0 + 1 − r̄ f−3/2
∆v f

< −1

Therefore the equation (14) has no solution under the assumption λ = 0.
(ii) We now consider the case λ > 0. Then the equation (11) yields
−x02 − (1 + y0 )2 (1 − r̄ f−2 ) + 2(1 − r̄ f−1 ) = 0

(15)

which together with (5) leads to x 2f = 0. Thus
∆v f =

q

x 2f + y 2f =

q

y 2f = y f

(16)

and also
x02 = −(1 + y0 )2 (1 − r̄ f−2 ) + 2(1 − r̄ f−1 )

(17)

∆v02 = x02 + y02 = r̄ f−2 (1 + y0 )2 − (1 + 2y0 ) + 2(1 − r̄ f−1 )

2
= r̄ f−1 (1 + y0 ) − r̄ f + 2(1 + y0 ) − r̄ f2 − (1 + 2y0 ) + 2(1 − r̄ f−1 )

2
= r̄ f−1 (1 + y0 ) − r̄ f + 3 − r̄ f2 − 2r̄ f−1

(18)

Then we obtain ∆v02 as follows

Thanks to (9)
x0
x0
+
− 2λx0 = 0
∆v0 ∆v f

(19)

We now show that x0 = 0 by contradiction. If x0 , 0, then the equality (19) implies that
2λ =

1
1
+
∆v0 ∆v f
6

(20)

Substituting (20) into (10) gives


y0 + 1 − r̄ f−3/2
y0
1
1
−2
+
= (1 + y0 )(1 − r̄ f )
+
∆v0
∆v f
∆v0 ∆v f
Arranging the above equation and using (4) and (16), we get
y0 − (1 + y0 )(1 − r̄ f−2 )
∆v0

r̄ f−3/2 − (y0 + 1)r̄ f−2

=

∆v f

=

r̄ f−1/2 − (y0 + 1)r̄ f−1
∆v f

r̄ f−1 =

yf
yf

r̄ f−1

where y f , 0 since we have assumed v f > 0 and just concluded x f = 0. By further rearranging the above equation, we
have
r̄ f−1 (1 + y0 ) − r̄ f
∆v0

=

yf
yf

=





 1,


if y f > 0




 −1,


if y f < 0

(21)

It is noted that (18) gives a strict inequality

2

2
∆v02 = r̄ f−1 (1 + y0 ) − r̄ f + 3 − r̄ f2 − 2r̄ f−1 < r̄ f−1 (1 + y0 ) − r̄ f

(22)

where the function 3 − r̄ f2 − 2r̄ f−1 < 0 since it is decreasing with respect to r̄ f and r̄ f > 1. Hence (21) contradicts to
(22) and does not hold, which implies that the assumption λ > 0 and the equality x0 , 0 do not hold simultaneously.
Therefore λ > 0 leads to x f = 0 and x0 = 0.
Summarizing the two cases above, we now conclude that a feasible solution to (7) must have λ∗ > 0, x ∗f = 0 and
x0∗ = 0. Then the corresponding normal component y0 can be solved from (15)
s
y0∗ =

s

2r̄ f
− 1,
1 + r̄ f

ŷ0∗ = −

2r̄ f
−1
1 + r̄ f

(23)

Substituting them into (4) and (10) yields
s

y ∗f

=

r̄ f−1/2

ŷ ∗f = r̄ f−1/2

!
2
1−
,
1 + r̄ f
s
!
2
1+
,
1 + r̄ f

y0∗ + 1 − r̄ f−3/2 ª
1
©
λ =
­1 +
®>0
y ∗f
2(1 + y0∗ )(1 − r̄ f−2 )
«
¬
∗ + 1 − r̄ −3/2
ŷ
1
0
f
©
ª
λ̂∗ =
­1 +
®>0
∗
∗
−2
ŷ
2(1 + ŷ0 )(1 − r̄ f )
f
«
¬
∗

(24)

One can see that there exist two sets of feasible solutions, one of which (x0∗, y0∗ ) is corresponding to the Hohmann
transfer, and its the Lagrange multiplier is λ∗ .

7

We are now in a position to show that the Hohmann transfer is the global minimum. We first show that (x0∗, y0∗ ) is a
strict local minimum. Theorem 3.11 in [21] gives a second order sufficient condition for a strict local minimum. To
apply it, we need to calculate the Hessian matrix of the Lagrangian function F defined by (8) at (x0∗, y0∗ )
∂2 F
©
­ ∂ x2
­
0
∇2 F(x0∗, y0∗, λ∗ ) = ­
­ ∂2 F
­
« ∂ y0 ∂ x0

∂2 F
ª
∂ x0 ∂ y0 ®®
®
∂ 2 F ®®
∂ y02 ¬(x ∗,y ∗ )
0 0

(25)

After a straightforward calculation, by the equations (9) and (10) and using the fact x0∗ = 0, we have
∂2 F
∂2 F
(x0∗, y0∗ ) =
(x ∗, y ∗ ) = 0
∂ x0 ∂ y0
∂ y0 ∂ x0 0 0
and also
1 + r̄ f−1 (1 + y0∗ )
1
1
∂2 F ∗ ∗
∗
(x
,
y
)
=
>0
+
−
2λ
=
y0∗ y ∗f
∂ x02 0 0
y0∗ (1 + y0∗ )(1 + r̄ f−1 )

(26)

∂2 F ∗ ∗
(x , y ) = 
∂ y02 0 0

(27)

−b
(y0∗

+

a)2

−b

∗
−2
 3/2 − 2λ (1 − r̄ f ) < 0

where a = 1 − r̄ f−3/2, b = (r̄ f − 1)2 (2r̄ f + 1)r̄ f−3 > 0. Then one can see that the Hessian matrix is an indefinite matrix.
Thus we cannot use the positive definiteness of the Hessian matrix, that is,
zT ∇2 F(x0∗, y0∗, λ∗ )z > 0,

∀z , 0

as a sufficient condition to justify the local minimum of (x0∗, y0∗ ) as usual. Fortunately, Theorem 3.11 in [21] tells us that
for the case λ∗ > 0, that is, the inequality constraint is active, when we use the positive definiteness of the Hessian
matrix to justify the local minimum, we need only to consider the positive definiteness of the first block of the Hessian
matrix with the non-zero vector z , 0 defined by
z ∈ Z(x0∗, y0∗ ) = {z : z T ∇g(x0∗, y0∗ ) = 0}
where g(x0, y0 ) is the constraint function
g(x0, y0 ) = −x02 − (1 + y0 )2 (1 − r̄ f−2 ) + 2(1 − r̄ f−1 )

8

(28)

With this kind of z, if zT ∇2 F(x0∗, y0∗, λ∗ )z > 0, then (x0∗, y0∗ ) is a strict local minimum. Specifically, the vector z defined
by the set (28) satisfies



z T ∇g(x0∗, y0∗ ) = z1





−2x0




=0
z2 

−2(1 + y )(1 − r̄ −2 )


0
f

 (x0∗,y0∗ )

(29)

Notice that x0∗ = 0, −2(1 + y0∗ )(1 − r̄ f−2 ) , 0. Hence (29) implies that the components z1 , 0 and z2 = 0. Obviously it
follows from (26) that



 
z 
 1
∂2 F ∗ ∗
2
∗ ∗ ∗  
0 ∇ F(x0, y0, λ )   = z1 2 (x0, y0 )z1 > 0
∂ x0
0
 
 


z1

(30)

Therefore in light of Theorem 3.11 in [21, p.48], we can conclude that (x0∗, y0∗ ) is a strict local minimum. The same
argument can be used to show that ( x̂0∗, ŷ0∗ ) is also a strict local minimum in view of (30). Meanwhile a straightforward
calculation gives
∆v(x0∗, y0∗ ) < ∆v( x̂0∗, ŷ0∗ )
Hence the pair (x0∗, y0∗ ) is a global minimum, and further it is the global minimum due to the uniqueness; see [22, p.194]
for the definition of a global minimum. This completes the proof of the theorem.

III. The Hohmann transfer as dynamic optimization problems
In this section, the orbit transfer problem is formulated as two optimal control problems of a spacecraft in an inverse
square law field, driven by velocity impulses, with boundary and interior point constraints. The calculus of variations is
used to solve the resulting two-point and multi-point boundary-value problems (BCs).

A. Problem formulation
Consider the motion of a spacecraft in the inverse square gravitational field, and the state equation is
rÛ = v
vÛ = −

µ
r
r3

(31)

where r(t) is the spacecraft position vector and v(t) is its velocity vector. The state vector consists of r(t) and v(t). We
use (31) to describe the state of the transfer orbit, which defines a conic under consideration; see, e.g., [19, Chapter 2].
Problem III.1 Given the initial position and velocity vectors of a spacecraft on the initial circular orbit, r(t0 ), v(t0 ).
9

The terminal time t1 is not specified. Let t0+ signify just after t0 and t1− signify just before t1 .∗ During the period from t0+
to t1− , the state evolves over time according to the equation (31). To guarantee at the time t1+ , the spacecraft enters the
finite circular orbit, we impose the following equality constraints on the finial state
gr1 (r(t1+ )) = r(t1+ ) − r f = 0
gv1 (v(t1+ )) = v(t1+ ) − v f = 0
g2 (r(t1+ ), v(t1+ )) = r(t1+ ) · v(t1+ ) = 0
where v f is the orbit velocity of the final circular orbit v f =

p

(32)

µ/r f . Suppose that there are velocity impulses at time

instants t0 and t1
v(ti+ ) = v(ti− ) + ∆vi,

i = 0, 1

(33)

where v(t0− ) = v(t0 ). The position vector r(t) is continuous at these instants. The optimal control problem is to design
∆vi that minimize the cost functional
J = |∆v0 | + |∆v1 |

(34)

subject to the constraints (32).
In control theory, such an optimal control problem is called impulse control problems in which there are state or
control jumps. Historically an optimal problem with the cost functional (34) is also referred to as the minimum-fuel
problem. By Problem III.1,† the orbit transfer problem has been formulated as a dynamic optimization problem instead
of a static one considered in Section II. The advantage of this formulation is that the coplanar assumption of the transfer
orbit to the initial orbit is removed, but the computation complexity follows.
In the next problem we consider the orbit transfer problem as a dynamic optimization problem with interior point
constraints.
Problem III.2 Consider a similar situation as in Problem III.1. Let t HT be the Hohmann transfer time. Instead of the
unspecified terminal time, here the terminal time instant t f > t HT is given and the time instant t1 now is an unspecified
interior time instant. The conditions (32) become a set of interior boundary conditions. The optimal control problem is
to design ∆vi that minimize the cost functional (34) subject to the interior boundary conditions (32).
∗ In

side.

mathematical language, t0+ represents the limit of t0 approached from the right side and t1− represents the limit of t1 approached from the left

† A Matlab script for the Hohmann transfer was given as a numerical example in the free version of a Matlab-based software GPOPS by an direct
method. The constrained (32) was also used to describe the terminal conditions. We here use the variational method (i.e., indirect method) to optimal
control problems

10

B. Two-point boundary conditions for Problem III.1
Problem III.1 is a constrained optimization problem subject to static and dynamic constraints. We use Lagrange
multipliers to convert it into an unconstrained one. Define the augmented cost functional
J˜ : = |∆v0 | + |∆v1 |
 +

 +

T
T
+ qr1
r(t0 ) − r(t0 ) + qr2
r(t1 ) − r(t1− )




+ qTv1 v(t0+ ) − v(t0 ) − ∆v0 + qTv2 v(t1+ ) − v(t1− ) − ∆v1
+ γr1 gr1 (r(t1+ )) + γv1 gv1 (v(t1+ )) + γ2 gr2 (r(t1+ ), v(t1+ ))
∫ t−
 µ

1
+
prT (v − rÛ ) + pTv − 3 r − vÛ dt
r
t0+
where Lagrange multipliers pr , pv are also called costate vectors; in particular, Lawden [5] termed −pv the primer
vector. Introducing Hamiltonian function

H(r, v, p) :=

prT v

−

µ
pTv 3 r,
r

T



p = pr

pv

(35)

By taking into account of all perturbations, the first variation of the augmented cost functional is
δ J˜ =

∆vT0
|∆v0 |

δv0 +

∆vT1
δv1
|∆v1 |





T
T
+ qr1
dr(t0+ ) − dr(t0− ) + qr2
dr(t1+ ) − dr(t1− )




+ qTv1 dv(t0+ ) − dv(t0− ) − δv0 + qTv2 dv(t1+ ) − dv(t1− ) − δv1
+ drT (t1+ )

∂gv1 (v(t1+ ))
∂gr1 (r(t1+ ))
T +
γ
+
dv
(t
)
γv1
r1
1
∂r(t1+ ) ∗
∂v(t1+ ) ∗

∂g2 (r(t1+ ), v(t1+ ))
∂g2 (r(t1+ ), v(t1+ ))
T +
γ
+
dv
(t
)
γ2
2
1
∗
∗
∂r(t1+ )
∂v(t1+ )


(∗) + H∗ − prT rÛ − pTv vÛ −∗ δt1
+ drT (t1+ )

t1

(∗∗) + prT (t0+ )δr(t0+ ) − prT (t1−∗ )δr(t1−∗ ) + pTv (t0+ )δv(t0+ ) − pTv (t1−∗ )δv(t1−∗ )
T

T #
∫ t −∗ " 
1
∂H(r, v, p)
∂H(r, v, p)
+
+ pÛ r (t) δr +
+ pÛ v (t) δv dt
∂r
∂v
t0+
∗
∗

(36)

where we use d(·) to denote the difference between the varied path and the optimal path taking into account the
differential change in a time instant (i.e., differential in x), for example,
dv(t1+ ) = v(t1+ ) − v∗ (t1+∗ )

11

and δ(·) is the variation, for example, δv(t1∗ ) is the variation of v as an independent variable at t1∗ . Notice that
dr(t0− ) = 0, dv(t0− ) = 0 since t0 is fixed. The parts of δ J˜ in (36) marked with asterisks are respectively due to the linear
term of
∫

t1−∗ +δt1

t1−∗




H(r, v, p) − prT rÛ − pTv vÛ dt

and the first term in the right-hand side of the following equation
∫

t1−∗

t0+

t −∗

prT δÛrdt = prT δr t1+ −
0

∫

t1−∗

t0+

pÛ rT δrdt

obtained by integrating by parts. In order to derive boundary conditions, we next use the following relation
dv(t1− ) = δv(t1−∗ ) + vÛ (t1−∗ )δt1

(37)

see [23, Section 3.5] and [20, Section 3.3].
In view of the necessary condition δ J˜ = 0 and the fundamental lemma, we have the costate equations
pÛ r (t) = −

∂H(r, v, p)
,
∂r

pÛ v (t) = −

∂H(r, v, p)
∂v

(38)

Based on the definition of the Hamiltonian function in (35), the costate equation (38) can be rewritten as
pÛ r =

∂ µ 
µ 3
r pv = − 3 ( 2 rrT − I3×3 )pv,
3
∂r r
r r

pÛ v = −pr

(39)

where I3×3 is the 3 × 3 identity matrix.
Using (37) and regrouping terms in (36) yields the following terms or equalities for ∆v1, ∆v2, dv


∆v0
− qv1
|∆v0 |

T
∆v0

T
∆v1
− qv2 ∆v1
|∆v1 |


qTv1 dv(t0+ ) + pTv (t0+ )δv(t0+ ) = qTv1 + pTv (t0+ ) dv(t0+ )



−qTv2 dv(t1− ) − pTv (t1− ) δv(t1−∗ ) + vÛ (t1−∗ )δt1 = − qTv2 + pTv (t1− ) dv(t1− )


∂gv1 (v(t1+ ))
∂g2 (r(t1+ ), v(t1+ ))
qTv2 + γv1
+
γ
dv(t1+ )
2
∗
∂vT (t1+ ) ∗
∂vT (t1+ )


(40)

As usual, in order to assure δ J˜ = 0, we choose Lagrange multipliers to make the coefficients of ∆v0 , ∆v1 , dv(t0+ ), dv(t1− ),

12

dv(t1+ ) in (40) vanish respectively
qv1 −

∆v0
= 0,
|∆v0 |

pv (t0+ ) + qv1 = 0,
qv2 + γv1

qv2 −

∆v1
=0
|∆v1 |

pv (t1− ) + qv2 = 0

∂gv1 (v(t1+ ))
∂g2 (r(t1+ ), v(t1+ ))
+
γ
2
∂v(t1+ ) ∗
∂v(t1+ )

∗

=0

(41)

Applying the similar argument as above for pr gives
pr (t0+ ) + qr1 = 0,
qr2 + γr1

pr (t1− ) + qr2 = 0

∂g2 (r(t1+ ), v(t1+ ))
∂gr1 (r(t1+ ))
+
γ
2
∂r(t1+ ) ∗
∂r(t1+ )

∗

=0

(42)

We now choose H(t1− ) to cause the coefficient of δt1 in (36) to vanish
H(t1− ) = prT (t1− )v(t1− ) − pTv (t1− )

µ
r 3 (t1− )

r(t1− ) = 0

Finally, by solving the first and second component equations of the last vector equation in (41), we obtain the
Lagrangian multipliers γv1, γ2 , and then substituting them into the first component equations of the last vector equation
in (42) yields the Lagrangian multiplier γr1 . With these Lagrangian multipliers and rearranging the second and third
component equations of the last vector equation in (42), we obtain a boundary value equation denoted by

g r(t1+ ), v(t1+ ); pr (t1− ), pv (t1− ) = 0

III.1 Two-point boundary conditions for Problem III.1

(1) r(t0 ) − r0 = 0,
(2) p M v (t0+ ) +

v(t0+ ) − v(t0− ) − ∆v0 = 0

∆v0
= 0,
|∆v0 |

p Mv (t1− ) +

∆v1
=0
|∆v1 |

(3) H(t1− ) = 0




r(t1+ ) − r f = 0






(4) v(t + ) − v f = 0

1






 r(t1+ ) · v(t1+ ) = 0



(5) g r(t1+ ), v(t1+ ); pr (t1− ), pv (t1− ) = 0

In summary, a complete list including 19 BCs is shown in List III.1 where r(t1+ ) = r(t1− ), v(t1+ ) = v(t1− ) + ∆v1 .
13

C. Multi-point boundary conditions for Problem III.2
Define the augmented cost functional
J˜ : = |∆v0 | + |∆v1 |
 +

 +

T
T
+ qr1
r(t0 ) − r(t0 ) + qr2
r(t1 ) − r(t1− )




+ qTv1 v(t0+ ) − v(t0 ) − ∆v0 + qTv2 v(t1+ ) − v(t1− ) − ∆v1
+ γr1 gr1 (r(t1+ )) + γv1 gv1 (v(t1+ )) + γ2 gr2 (r(t1+ ), v(t1+ ))
∫ tf
∫ t−

 µ

 µ
1
T
T
Û
Û
prT (v − rÛ ) + pTv − 3 r − vÛ dt
+
pr (v − r) + pv − 3 r − v dt +
r
r
t1+
t0+
We introduce the Hamiltonian functions for two time sub-intervals [t0+, t1− ] ∪ [t1+, t f ]
T
Hi (r, v, pi ) := pri
v − pTvi

µ
r,
r3


pi = pri

T
pvi

,

i = 1, 2

(43)

A tedious and similar argument as used in Problem III.1 can be applied to Problem III.2 to derive the costate
equations and boundary conditions. The two-phase costate equations are given by
pÛ ri (t) = −

∂Hi (r, v, pi )
,
∂r

pÛ vi (t) = −

∂Hi (r, v, pi )
,
∂v

i = 1, 2

A complete list for boundary conditions is given in List III.2.
III.2 31 boundary conditions for Problem III.2

(1)





 r(t0 ) − r0 = 0,


v(t0− ) − v(t0+ ) − ∆v0 = 0




 r(t1− ) − r(t1+ ) = 0, v(t1− ) − v(t1+ ) − ∆v1 = 0



∆v1
∆v0

+

= 0, pv1 (t1− ) +
=0
 pv1 (t0 ) +

|∆v
|
|∆v
0
1|
(2)



 pv2 (t f ) = 0, pr2 (t f ) = 0

T (t − )∆v = 0
(3) H1 (t1− ) − H2 (t1+ ) = 0 or − pr1
1
1



+


 r(t1 ) − r f = 0




(4) v(t + ) − v f = 0

1






 r(t1+ ) · v(t1+ ) = 0



(5) g r(t1+ ), v(t1+ ); pr1 (t1− ), pr1 (t1+ ), pv1 (t1− ), pv1 (t1+ ) = 0

14

IV. Numerical Examples
7

The magnitude of the primer vector

x 10

1
4
3

0.95

2

ry m

1

0.9

0
0.85

−1
−2

0.8

−3
−4
−6

−4

−2

0

2

rx m

Fig. 1

4

0.75

6

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

t secon d

7

x 10

1.8

2
4

x 10

The Hohmann transfer and the magnitude of the primer vector

In order to use Matlab solver bvp4c or bvp5c to solve two-point and multi-point boundary value problems with
unspecified switching time instants in Section III, two time changes must be introduced. For the time change, we refer
to, e.g., [20, Appendix A] and [24, Section 3] for details. The solver bvp4c or bvp5c accepts boundary value problems
with unknown parameters; see [25] and [26].
Example IV.1 This example in [10, p.25] is used to illustrate that the Hohmann transfer is the solution to Problem
III.1. The altitude of the initial circular orbit is 300km, and the desired final orbit is geostationary, that is, its radius is
42164km. By the formulas of the Hohmann transfer, the magnitudes of two velocity impulses are given by
∆v1 = 2.425726280326563e + 03,

∆v2 = 1.466822833675619e + 03

m/s

(44)

Setting the tolerance of bvp4c equal to 1e-6, we have the solution to the minimum-fuel problem III.1 with an unspecified
terminal time given by bvp4c
dv1 = 1.0e+03*[0.000000000000804;

2.425726280326426;

dv2 = 1.0e+03*[0.000000000000209;

-1.466822833675464;

0]
0]

Compared with (44), the accuracy of the numerical solution is found to be satisfactory since only the last three digits of
fifteen digits after decimal place are different. It is desirable to speed up the computational convergence by scaling
the state variables though bvp4c already has a scale procedure. The computation time is about 180 seconds on Intel
Core i5 (2.4 GHz, 2 GB). Figure 1 shows the Hohmann transfer and the magnitude of the primer vector. The solution

15

provided by bvp4c depends upon initial values. Figure 2 shows the transfer orbit and the magnitude of the primer vector
corresponding to the local minimum given by (24) in the proof of Theorem II.1.
7

The magnitude of the primer vector

x 10

1
4
0.9

3
2

0.8

ry m

1
0.7

0
−1

0.6
−2
−3

0.5

−4
−6

−4

−2

0

2

4

Fig. 2

0.4

6

rx m

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

t secon d

7

x 10

1.8

2
4

x 10

The orbit transfer and the magnitude of the primer vector for the local minimum

6

The magnitude of the primer vector

x 10
8

1
6

0.8

4

2

ry m

0.6
0
0.4
−2

0.2

−4

−6
0
−8
−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

rx m

Fig. 3

0

1

500

1000

1500

2000

2500

3000

t secon d

7

x 10

The Hohmann transfer and the magnitude of the primer vector

Example IV.2 Consider the Hohmann transfer as the multi-point boundary value problem (Problem III.2). The altitudes
of the initial and final circular orbits are 200km and 400km respectively. By the formulas of the Hohmann transfer, the
magnitudes of two velocity impulses are given by
∆v1 = 58.064987253967857,

16

∆v2 = 57.631827424189602

(45)

The Hohmann transfer time t HT = 2.715594949192177e + 03 seconds, hence we choose the terminal time t f = 2800
seconds. The constants, the initial values of the unknown parameters, and the solver and its tolerance are specified in
Table 1.

Table 1
Constants

Constants and initial values
Gravitational constant µ = 3.986e + 14

Earth’s radius Re = 6378145

r0 m

v0 m/s

The initial values of the state
6578145

0

0

0

7.7843e+03

pv

0

pr

The initial values of the costate
-0.0012

0

0

0

-0.9

∆v1

0

∆v2

The initial values of the velocity impulses
0

102

0

The initial value of the scaled time

0

102

0

0.9

Solver bvp4c

Tolerance=1e-7

By trial and error, this example is successfully solved by using bvp4c. We obtain the following solution message
The solution was obtained on a mesh of 11505 points.
The maximum residual is

4.887e-08.

There were 2.2369e+06 calls to the ODE function.
There were 1442 calls to the BC function.
Elapsed time is 406.320776 seconds.
The first velocity impulse vector dv1 = [0.000000000000001,

58.064987253970472, 0]

The second velocity impulse vector dv2 = [0.000000000000924, -57.631827424187151, 0]
The scaled instant of the second velocity impulse 0.969855338997211
The time instant of the second velocity impulse 2.715594949192190e+03
The maximal error of boundary conditions 4.263256e-13
Figure 3 shows the Hohmann transfer and the magnitude of the primer vector. The magnitude of the velocity and the
velocity vector on the vx − vy plane are shown in Figure 4, where the symbol + corresponds to the time instant at which
the velocity impulse occurs.
We now assume the initial values of the velocity impulses
dv10=[0, 45, 0],

dv20=[0, -60, 0]
17

The magnitude of the velocity vector

The velocity vector
8000

7850

6000
7800

4000
2000
vy m/s

7750

7700

0
−2000
−4000

7650
−6000

0

500

1000

Fig. 4

1500
time second

2000

−8000
−8000

2500

−7000

−6000

−5000

−4000 −3000
vx m/s

−2000

−1000

0

1000

The magnitude of the velocity and the velocity vector on the vx − vy plane

In stead of bvp4c, we use the other solver bvp5c. Setting Tolerance=1e-6, the solution corresponding to the local
minimum is found.
The solution was obtained on a mesh of 2859 points.
The maximum error is

5.291e-14.

There were 809723 calls to the ODE function.
There were 1382 calls to the BC function.
Elapsed time is 190.920474 seconds.
The first velocity impulse vector
The second velocity impulse vector

dv1 = 1.0e+4*[0.000000000000001, -1.562657038966310, 0]
dv2 = 1.0e+4*[0.000000000000011, -1.527946697280544, 0]

The scaled instant of the second velocity impulse 0.969855338997205
The time instant of the second velocity impulse 2.715594949192175e+03
The maximal error of boundary conditions 4.403455e-10
Figure 5 shows the orbit transfer and the magnitude of the primer vector corresponding to the local minimum, where the
spacecraft starts from the initial point on the initial orbit, moves clockwise along the transfer orbit until the second pulse
point, and then travels counterclockwise along the final circle orbit until the terminal point marked with the symbol ×.

V. Conclusion
In this paper, by a static constrained optimization, we study the global optimality of the Hohmann transfer using a
nonlinear programming method. Specifically, an inequality presented by Marec in [10, pp. 21-32] is used to define
an inequality constraint, then we formulate the Hohmann transfer problem as a constrained nonlinear programming
18

6

The magnitude of the primer vector

x 10
8

1
6

0.8

4

2

ry m

0.6
0
0.4
−2

0.2

−4

−6
0
−8
−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

rx m

Fig. 5

0.6

0.8

0

1

500

1000

1500

2000

2500

3000

t secon d

7

x 10

The orbit transfer and the magnitude of the primer vector corresponding to the local minimum

problem. A natural application of the well-known results in nonlinear programming such as the Kuhn-Tucker theorem
clearly shows the the global optimality of the Hohmann transfer. In the second part of the paper, we introduce two
optimal control problems with two-point and multi-point boundary value constraints respectively. With the help of
Matlab solver bvp4c, the Hohmann transfer is solved successfully.

Acknowledgments
The first author is supported by the National Natural Science Foundation of China (no. 61374084).

References
[1] Hohmann, W., The Attainability of Heavenly Bodies (1925), NASA Technical Translation F-44, 1960.
[2] Barrar, R. B., “An Analytic Proof that the Hohmann-Type Transfer is the True Minimum Two-Impulse Transfer,” Astronautica
Acta, Vol. 9, No. 1, 1963, pp. 1–11.
[3] Prussing, J. E., and Conway, B. A., Orbital Mechanics, 1st ed., Oxford University Press, 1993.
[4] Ting, L., “Optimum Orbital Transfer by Impulses,” ARS Journal, Vol. 30, 1960, pp. 1013–1018.
[5] Lawden, D. F., Optimal Trajectories for Space Navigation, Butter Worths, London, 1963.
[6] Prussing, J. E., “Primer Vector Theory and Applications,” Spacecraft Trajectory Optimization, edited by B. A. Conway,
Cambridge, 2010, Chap. 2, pp. 16–36.
[7] Pontani, M., “Simple Method to Determine Globally Optimal Orbital Transfers,” Journal of Guidance, Control, and Dynamics,
Vol. 32, No. 3, 2009, pp. 899–914.

19

[8] Moyer, H. G., “Minimum Impulse Coplanar Circle-Ellipse Transfer,” AIAA Journal, Vol. 3, No. 4, 1965, pp. 723–726.
[9] Battin, R. H., An Introduction to the Mathematics and Methods of Astrodynamics, AIAA Education Series, AIAA, New York,
1987.
[10] Marec, J.-P., Optimal Space Trajectories, Elsevier, New York, 1979.
[11] Hazelrigg, G. A., “Globally Optimal Impulsive Transfers via Green’s Theorem,” Journal of Guidance, Control, and Dynamics,
Vol. 7, No. 4, 1984, pp. 462–470.
[12] Palmore, J. I., “An Elementary Proof the Optimality of Hohmann Transfers,” Journal of Guidance, Control, and Dynamics,
Vol. 7, No. 5, 1984, pp. 629–630.
[13] Prussing, J. E., “Simple Proof the Global Optimality of the Hohmann Transfer,” Journal of Guidance, Control, and Dynamics,
Vol. 15, No. 4, 1992, pp. 1037–1038.
[14] Yuan, F., and Matsushima, K., “Strong Hohmann Transfer Theorem,” Journal of Guidance, Control, and Dynamics, Vol. 18,
No. 2, 1995, pp. 371–373.
[15] Vertregt, M., “Interplanetary orbits,” Journal of the British Interplanetary Society, Vol. 16, 1958, pp. 326–354.
[16] Cornelisse, J. W., Schöyer, H. F. R., and Walker, K. F., Rocket Propulsion and Spaceflight Dynamics, Pitman, London, 1979.
[17] Gurfil, P., and Seidelmann, P. K., Celestial Mechanics and Astrodynamics Theory and Practice, Springer, 2016.
[18] Avendaño, M., Martín-Molina, V., Martín-Morales, J., and Ortigas-Galindo, J., “Algebraic Approach to the Minimum-Cost
Multi-Impulse Orbit-Transfer Problem,” Journal of Guidance, Control, and Dynamics, Vol. 39, No. 8, 2016, pp. 1734–1743.
[19] Curtis, H. D., Orbital Mechanics for Engineering Students, 3rd ed., Elsevier Ltd., 2014.
[20] James M. Longuski, J. J. G., and Prussing, J. E., Optimal Control with Aerospace Applications, Springer, 2014.
[21] Avriel, M., Nonlinear Programming: Analysis and Methods, an unabridged republication of the edition published by
Prentice-Hall in 1976 ed., Dover Publications Inc., 2003.
[22] Bertsekas, D., Nonlinear Programming, 2nd ed., Athena Scientific, 1999.
[23] Bryson Jr., A. E., and Ho, Y.-C., Applied Optimal Control, Hemisphere Publishing Corp., London, 1975.
[24] Zefran, M., Desai, J. P., and Kumar, V., “Continuous Motion Plans for Robotic Systems with Changing Dynamic Behavior,”
The second Int. Workshop on Algorithmic Fundations of Robotics, Toulouse, France, 1996.
[25] Shampine, L., Gladwell, I., and Thompson, S., Solving ODEs with Matlab, Cambridge University Press, 2003.
[26] Kierzenka, J., “Studies in the numerical solution of ordinary differential equations,” Ph.D. thesis, Department of Mathematics,
Southern Methodist University, Dallas, TX., 1998.

20

