Robert D. Falck

Glenn Research Center, Cleveland, Ohio

John W. Dankanich

AeroDank, Inc., Cleveland, Ohio

Optimization of Low-Thrust Spiral Trajectories

by Collocation

NASA/TM—2012-217699

September 2012

AIAA–2012–4423
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

 443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
Robert D. Falck

Glenn Research Center, Cleveland, Ohio

John W. Dankanich

AeroDank, Inc., Cleveland, Ohio

Optimization of Low-Thrust Spiral Trajectories

by Collocation

NASA/TM—2012-217699

September 2012

AIAA–2012–4423

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

2012 Astrodynamics Specialists Conference

sponsored by the American Institute of Aeronautics and Astronautics (AIAA) and 

American Astronautical Society (AAS)

Minneapolis, Minnesota, August 13–16, 2012
Acknowledgments

The work described in this paper was funded in whole or in part by the In-Space Propulsion Technology Program, which 

is managed by NASA’s Science Mission Directorate in Washington, D.C., and implemented by the In-Space Propulsion 

Technology Project at the NASA Glenn Research Center in Cleveland, Ohio.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://www.sti.nasa.gov

Trade names and trademarks are used in this report for identifi cation 

only. Their usage does not constitute an offi cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

Level of Review: This material has been technically reviewed by technical management. 

This report is a formal draft or working 

paper, intended to solicit comments and 

ideas from a technical peer group.

This report contains preliminary fi ndings, 

subject to revision as analysis proceeds.
Optimization of Low-Thrust Spiral Trajectories by

Collocation

As NASA examines potential missions in the post space shuttle era, there has been

a renewed interest in low-thrust electric propulsion for both crewed and uncrewed missions. While much progress has been made in the field of software for the optimization

of low-thrust trajectories, many of the tools utilize higher-fidelity methods which, while

excellent, result in extremely high run-times and/or poor convergence when dealing with

planetocentric spiraling trajectories deep within a gravity well. Conversely, faster tools

like SEPSPOT provide a reasonable solution but typically fail to account for other forces

such as third-body gravitation, aerodynamic drag, solar radiation pressure. SEPSPOT is

further constrained by its solution method, which may require a very good guess to yield

a converged optimal solution. Here the authors have developed an approach using collocation intended to provide solution times comparable to those given by SEPSPOT while

allowing for greater robustness and extensible force models.

Nomenclature

t Time, days

x State variable vector

u Time-varying control variable vector

φg Global design variable vector

φp Phase design variable vector

ψ0 Initial boundary constraint

ψf Final boundary constraint

ψp Path constraint

τ Non-dimensional time within a segment polynomial

∆c Continuity defect

∆s Differential defect

Li Lagrange interpolation matrix for interior nodes

Di Lagrange derivative matrix for interior nodes

Dc Lagrange derivative matrix for cardinal nodes

n Number of state variables

m Number of time-varying control variables

ncnp Number of cardinal nodes per phase

ninp Number of interior nodes per phase

z Independent variable vector in the nonlinear optimization problem

g Boundary constraint vector for the nonlinear optimization problem

a Semi-major axis, km

e Eccentricity

i Inclination, deg

Robert D. Falck 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

John W. Dankanich 

AeroDank, Inc. 

Cleveland, Ohio 44135

NASA/TM—2012-217699 1
Ω Right ascension of ascending node, deg

ω Argument of periapsis, deg

θ True anomaly, deg

h Specific angular momentum, km2

sec

p Semilatus rectum, km

f X-component of eccentricity vector in the equinoctial frame

g Y-component of eccentricity vector in the equinoctial frame

j X-component of ascending node vector in the equinoctial frame

k Y-component of ascending node vector in the equinoctial frame

L True longitude, deg

r Spacecraft radius from central body, km

v Spacecraft velocity in the EME2000 frame, km

sec

m Spacecraft mass, kg

rev Number of orbital revolutions

Tp Orbital period, days

P QW Perifocal coordinate frame

NTW Velocity-vector aligned spacecraft coordinate frame

RSW Local-vertical, local-horizontal spacecraft coordinate frame

EME2000 Earth Mean Equator of J2000 inertial coordinate frame

rcb Central body equatorial radius, km

μ Gravitational parameter of the central body, km3

s2

P0 Spacecraft propulsion system input power, kW

η Spacecraft propulsion system efficiency

Isp Spacecraft propulsion system specific impulse, sec

Subscript

L Lower bound

U Upper bound

i Interior node

c Cardinal node

0 Initial condition

f Final condition

I. Introduction

For several decades, SEPSPOT has been NASA’s primary tool for the optimization of planetocentric

low-thrust trajectories.1 SEPSPOT takes advantage of orbital-averaging techniques to greatly increase

the speed of computation at the expense of fidelity. Generally, orbital averaging techniques provide an

accurate answer regarding the transfer of one orbit to another, but lack the precision to put a spacecraft at

a given position at a specific time. Although relatively fast, SEPSPOT suffers from some weaknesses that

result in difficulty when using it to analyze state-of-the-art low-thrust transfers. SEPSPOT utilizes indirect

optimization techniques based on Hamilton-Lagrange theory. In addition to providing an initial guess at

values of the orbital state variables, one must also provide a fairly accurate guess for initial and final values

of the costates, which can be challenging. More critically, extending the equations of motion to include

additional force models (atmospheric drag, solar pressure, perturbations by the Moon and Sun, etc.) would

require the equations of motion of the costate variables to be rederived from their current form.

To address these issues, the authors have embarked on developing an approach which uses a collocation

technique to transform an initial-value-problem (IVP) or two-point boundary value problem (TPBVP) into a

nonlinear programming (NLP) problem, which may be solved utilizing an off-the-shelf optimization software

package. Since this approach utilizes direct optimization method, it has the advantage that new environmental and vehicle related force models may be added with relative ease. Collocation and pseudospectral

optimization techniques were demonstrated by Dickmanns and Well2 and Hargraves and Paris3 and have

been used with success in software such as OTIS4,4 DIDO,5 and SOCS.6 Collocation and pseudospectral

methods have previously been applied to the low-thrust spiraling problem by Betts.7 However, in his work

he directly transcribed the equinoctial equations of motion with respect to time, which led to a mesh grid

NASA/TM—2012-217699 2
of over 16000 points and a run time of several hours.7 Rather than take this approach, which sees highly

oscillatory behavior in the time-histories of the states and controls, the authors utilize orbital averaging

techniques in conjunction with hybrid control formulations which largely remove the oscillatory behavior

from the problem. This substantially reduces the number and order of the collocating polynomials, which

reduces the size of the NLP problem and results in run times on the order of a few seconds to minutes.

Kluever demonstrated an orbital averaging problem which utilizes the Gauss form of the Lagrange planetary equations in conjunction with a hybrid control technique.8 Using only a handful of mesh points for the

controls, he demonstrated the ability to optimize a LEO to GEO transfer in good agreement with SEPSPOT.

For this approach the authors elected to use collocation of the states and controls since the method tends

to show an ability to find optimal solutions given a relatively poor initial guess. The authors use equations

of motion based on a modified set of equinoctial elements, which don’t suffer from as many singularities as

the Lagrange planetary equations.

Like Kluever, the authors utilize a hybrid set of control variables. These control variables are the costates

of the Hamilton-Lagrange formulation. They show a relatively smooth time-history throughout an orbital

transfer, which is beneficial for collocation, but internally are transformed into a set of pitch and yaw angles

which may oscillate greatly over the course of any single orbit. Furthermore, due to the extensible nature

of the collocated form of this problem, alternative guidance strategies can be easily implemented so long as

they can be parameterized such that their control variables do not exhibit ”excessively” oscillatory behavior

throughout the orbital transfer. As such, this method allows users to determine both the optimal (minimum

time, minimum propellant, etc.) transfer guidance history, but also how it compares to a custom guidance

strategy. The use of collocation to solve the optimal control problem should, in general, result in a greater

radius of convergence such that the initial guesses to the costates need not be as accurate as those supplied

to SEPSPOT.

This approach does sacrifice fidelity for speed, and thus is intended for the preliminary analysis of lowthrust trajectories and guidance algorithm development. For example, orbital averaging is incapable of

analyzing a lunar insertion maneuver as demonstrated by Betts. However, the run-time is so small that

using a collocated averaging technique to ”dig” out of the deepest part of the gravity well would greatly

reduce analysis time. Once at a sufficiently high intermediate orbit, the oscillatory behavior of the state

and control time histories is low enough to permit collocation of non-averaged equations of motion from the

intermediate orbit to the final target.

II. Formulation

A. The Collocation Problem

Collocation is used to simultaneously simulate and optimize the trajectory of a dynamical system. As

a fallback, the implementation also supports a more traditional explicitly integrated shooting method for

trajectory optimization, though experience has shown that implicit integration is not only faster, but often

capable of finding more optimal solutions.

In general terms, the optimization problem is of the following form:

Minimize fobj (x(t0|f), t0|f , u(t0|f), φg, φp) (1)

s.t. xL ≤ x ≤ xU (2)

tL ≤ t ≤ tU (3)

uL ≤ u ≤ uU (4)

φg,L ≤ φg ≤ φg,U (5)

φp,L ≤ φp ≤ φp,U (6)

ψ0,L ≤ fψ0 (x(t0), t0, u(t0), φg, φp) ≤ ψ0,U (7)

ψf,L ≤ fψ0 (x(tf), tf , u(tf), φg, φp) ≤ ψf,U (8)

ψp,L ≤ fpath(x(t), t, u(t), φg, φp) ≤ ψp,U (9)

Where x and u are the time-varying states and controls, respectively. The vectors φg and φp contain

static optimization variables (design variables) for the entire problem, and a specific phase, respectively.

NASA/TM—2012-217699 3
Variable t is the time at the point of the objective function evaluation; either the beginning or end of a

phase. All states, controls, design variables, and time may have simple bounds. Furthermore, we may

impose boundary constraints at the start or end of a phase (ψ0 and ψf ), or path constraints to be assessed

throughout a phase (ψp).

Figure 1. The state time-history of a single state in a phase consisting of three 3rd-order segments, including

the differential defects (∆si,j ) and continuity defects (∆ci,j )

The collocation problem consists of a series of phases in which the forces acting upon the spacecraft

are consistent such that the states are C1 continuous within a phase. Each phase is defined on an interval

[t0, t0 + tp], and the equations of motion dictate some set of states x and controls u. Each phase is further

subdivided into a series of segments along which the states and controls are represented by a polynomial of

at least order 3, and thus C2 continuous. For each state and control along each segment, this polynomial is

defined at a series of nodes called the cardinal nodes. Internally, each segment of order o is defined on the

interval [−1, 1] (so called τ -space) and has cardinal nodes at the (o+ 1) Legendre-Gauss-Lobatto points. For

example, a 3rd-order polynomial segment has four cardinal nodes. The collocation engine uses a nonlinear

optimization routine to vary the cardinal values of the states and controls, the phase start time and duration,

and design variable values such that state time histories accurately reflect the equations of motion and the

objective function (1) is satisfied.

The accuracy of the collocated polynomials in representing the dynamics of the system is measured by

comparing the slope of the polynomial for a given state to the derivative of that state as given by the

equations of motion. These differences are called differential defects. A well-posed initial-value-problem has

a unique solution and requires that the number of independent variables is equal to the number of equality

constraints. Figure 1 illustrates the setup of a collocation problem for one state variable in a phase of three

3rd-order polynomial segments. With just one state this problem has 14 variables: the initial time and phase

duration (2), and the state values at four cardinal nodes in each of three segments (12). We require that

NASA/TM—2012-217699 4
there be no discontinuity in the value of a state variable at the segment boundaries.

∆ci,j = xc,j − xc,i = 0 (10)

Furthermore, the initial time, duration, and initial state value in this phase are fixed via bounds, giving

a total of five equality constraints. The remaining nine equality constraints are achieved by requiring the

differential defects at three points in each segment to be zero. The points in each polynomial where the

defects are assessed are referred to as the interior nodes. Within each segment we may construct a Lagrange

interpolation matrix Li and a Lagrange derivative matrix Di such that the values and derivatives of the

states and controls at the interior nodes may be obtained by simple matrix-vector products:

xi = Lixc (11)

dxi

dt = Dixc

dτ

dt (12)

(13)

Finally, the differential defects of the states at the interior node of the phase are:

∆si = feom (xi, ti, ui, φg, φp) − dxi

dt = 0 (14)

Some preliminary tests suggest that a good choice of interior nodes is to determine the LGL points for

an (o + 1) polynomial and remove the endpoints. Reusing the first (or last) o cardinal nodes as the interior

nodes is also possible, but generally exhibits poorer convergence. On the other hand, if the convergence

issues can be fixed, reusing the cardinal nodes as interior nodes would reduce equation of motion evaluations

by nearly a factor of two, and no matrix multiplications (11) would be required to determine the state and

control values at the interior nodes.

Experience has also shown that convergence is aided by requiring both value and rate continuity to

be imposed on the control values at the segment boundaries. Control rates at the segment boundaries

are determined with a Lagrange derivative matrix constructed such that it returns the derivatives of the

polynomial at the cardinal nodes:

duc

dt = Dcuc

dτ

dt (15)

B. The NLP Interface

The optimization problem above is transcribed by the collocation engine into the following form:

Minimize fobj (z) (16)

s.t. gL ≤ g(z) ≤ gU (17)

zL ≤ z ≤ zU (18)

This transformed optimization problem may be solved by a variety of off-the-shelf nonlinear optimization

software packages. Currently the authors are using the IPOPT9 which is able to capitalize on the sparsity

of the Jacobian matrix that is characteristic of the collocaiton problem.

The times, states, controls, and design parameters are scaled and packed into the independent variable

array for the NLP (z). Currently a basic scaling scheme is used whereby each independent variable is scaled

by the inverse of its expected magnitude such that the corresponding component of z is roughly on the order

of one. The independent parameter array is packed in the following order:

1. The scaled global design parameter values comprise the first elements of z

2. The scaled design parameter values of the first phase comprise the next elements of z

NASA/TM—2012-217699 5
3. The scaled initial time and duration of the first phase are the next two elements of z

4. The cardinal values of the n states and m controls of the first phase are packed as the next (n+m)·ncnp

elements, where ncnp is the number of cardinal nodes in the phase. These values are packed such that

the states and columns at any given cardinal node are contiguous. Conceptually, this is the same as

constructing the following ncnp by n + m matrix XcUc and unraveling its values in row-major order.

XcUc =













x0,0 x1,0 ··· xn−1,0 u0,0 ··· um−1,0

x0,1 x1,1 ··· xn−1,1 u0,1 ··· um−1,1

.

.

. .

.

. ... .

.

. .

.

. .

.

. .

.

.

x0,ncnp−1 x1,ncnp−1 ··· xn−1,ncnp−1 u0,ncnp−1 ··· um−1,ncnp−1













(19)

5. Items 2 through 4 are repeated for each phase in the problem.

z =

'

[φg]

T [φp]T t0 tp [XcUc]

0 ··· [XcUc]ncnp−1

(

(20)

Similarly, the constraints of the collocation problem are packed into the NLP constraint array (g) array

in the following order:

1. The phase linkage constraints comprise the first elements of g

2. The ninp · n state defects for the first phase are stored node-by-node as the next elements of g.

3. The state and control continuity conditions at each of the num seg −1 segment boundaries in the first

phase make up the next (num seg − 1) · (n + m) elements of g.

4. The control continuity rate conditions at each of the num seg−1 segment boundaries in the first phase

make up the next (num seg − 1) · m elements of g.

5. The initial and final boundary constraints of the first phase make up the next elements of g

6. Finally, the path constraints are evaluated at each cardinal node and comprise the final ncnp · npc

elements of g, where ncp is the number of path constraints in the phase.

7. Items 2 through 6 are repeated for each phase in the problem.

g =

'

[∆s]

T

0:ninp−1 [∆c]

T

0:num seg−1 [∆ ̇c]

T

0:num seg−1 [ψ0]

T [ψf ]T [ψp]T

0 ··· [ψp]

T

ncnp−1

(

(21)

By ordering the z and g arrays in such a way the sparsity pattern of the Jacobian matrix is largely

block-diagonal. The use of a nonlinear programming routine which is able to capitalize on the sparsity of

the Jacobian matrix is a critical feature which enables the solution to the collocation problem to be obtained

rapidly.

C. Equations of Motion

In his work, Kluever showed that, by applying orbital averaging to the Gauss form of the Lagrange planetary

equations,10 he was able to get very good agreement with SEPSPOT results by parameterizing the controls

at some set of control nodes and explicitly integrating the solution.8 However, the Lagrange planetary

equations are problematic in the context of a collocation formulation since they exhibit singularities at zero

eccentricity and zero inclination, a point which is of interest for the purpose of examining trajectories to or

from geostationary orbit. Instead, the authors utilized a set of modified equinoctial elements as given by

Walker.11 Betts also utilized this set of equations of motion in his low-thrust orbit transfer example, though

he did not employ an averaging technique.12

NASA/TM—2012-217699 6
x ̇ =







































0 2p

w

)p

μ

0

)p

μsin (L)

)p

μ

[(w+1)cos(L)+f]

w −

)p

μ

g(jsin(L)−kcos(L))

w

−

)p

μcos (L)

)p

μ

[(w+1)sin(L)+g]

w

)p

μ

f(jsin(L)−kcos(L))

w

0 0 )p

μ

s2

2w cos (L)

0 0 )p

μ

s2

2w sin (L)

0 0 )p

μ

jsin(L)−kcos(L)

w

00 0

00 0













































ar

as

aw





 +































0

0

0

0

0

√μp *q

p

+2

−m ̇

rev ̇































(22)

where the state vector consists of the modified equinoctial elements and spacecraft mass

x =

'

p f g j k L m rev(T

(23)

the accelerations

a =







ar

as

aw





 (24)

are the radial, local-horizontal, and orbit normal perturbing accelerations where

s2 =1+ j2 + k2 (25)

w = p

r

=1+ fcos (L) + gsin (L) (26)

Note we refer to the equinoctial element h as j to avoid confusion with the specific angular momentum

of the orbit. The equations of motion (22) account for central-body gravitation in the absence of perturbing

accelerations. Adding different perturbing accelerations is as simple as converting them to the RSW frame

and including them in (24). For the purposes of this paper, the only non-two-body accelerations are due to

thrust and J2 perturbations.

The mass flow rate (  ̇m) is based on the propulsion system model (29).

The number of orbital revolutions performed by the spacecraft can be approximated by treating it as an

integrated state variable with a derivative function equal to the frequency of the spacecraft orbit.

rev ̇ = 1

Tp

=

√μ

2πa 3

2

(27)

D. Solar Electric Propulsion Model

The solar electric propulsion system generates power based on its distance to the Sun and the angle between

the array plane and the Sun vector. In the analysis below the Sun is assumed to be 1 A.U. from the spacecraft

at all times and the arrays are always perfectly pointed at the Sun. If the nominal array power at 1 A.U. is

given by P0, then the spacecraft acceleration due to thrust is8

aT = 2ηP0

mgIsp

(28)

and the rate of propellant expenditure due to engine firing is derived from the definition of thrust:

FT =  ̇mgIsp −→ m ̇ = 2ηP0

g2I2

sp

(29)

NASA/TM—2012-217699 7
E. Gravitational Perturbations due to Oblateness of the Central Body

Kechichian quantified J2 perturbations in the RSW frame as:13

aJ2 =









−3μJ2r2

cb

2r4

*

1 − 12(jsin(L)−kcos(L))2

s4

+

−12μJ2r2

cb

r4

(jsin(L)−kcos(L))(jcos(L)+ksin(L))

s4

−6μJ2r2

cb

r4

(jsin(L)−kcos(L))(1−k2−j2)

s4









(30)

F. Orbital Averaging

The equations of motion (22) are subject to oscillatory behavior which can greatly increase the number of

collocation segments required for a converged solution. To enable better performance from the collocation

routine, the equations of motion are averaged. The averaged equations of motion are obtained by converting

 ̇x to dx

dθ , integrating for an entire orbital revolution, and dividing by the orbital period to obtain the average

rates of change in the modified equinoctial elements for a given orbital state.8 The anomaly term (true

longitude for the modified equinoctial elements) is not included among the averaged orbital elements.

x ̄ ̇ = 1

Tp

, π

−π

dx

dt

dt

dθ

dθ (31)

The transformation from time to true anomaly (θ) is approximated using the perturbation-independent term

of the Gauss form of the Lagrange planetary equation for true anomaly.10

dθ

dt ≈ hr2 (32)

G. Eclipse Arcs

While accelerations due to solar electric propulsion only need to compute the integral in (31) through the

illuminated orbital arcs, forces due to gravitational harmonics, solar radiation pressure, drag, third-body

effects, and trapped radiation impingement require that the entire orbit be integrated. To determine those

points at which the spacecraft enters and exits the shadow of the central body, the authors employed Escobal’s

quartic shadow function:1415

S = α0cos (θ)

4 + α1cos (θ)3 + α2cos (θ)2 + α3cos (θ) + α4 (33)

where















α1

α2

α3

α4

α5















=















e4 −2(β2

2 − β21 )e2 (β22 + β21 )2

4e3 −4e(β2

2 − β21 ) 0

6e2 −2((β2

2 − β21 ) + e2 ∗ (1 − β22 )) 2(β22 − β21 )(1 − β22 ) − 4(β21 β22 )

4e −4e(1 − β2

2 ) 0

1 −2(1 − β2

2 ) (1 − β22 )2





















α4

α2

1





 (34)

α = rcb

p (35)

and β is the unitized vector from the central (shadowing) body to the Sun in the perifocal frame (PQW):15

ˆrs = β1ˆp + β2ˆq + β3 ˆw (36)

NASA/TM—2012-217699 8
The solution of the quartic shadow funtion (33) will yield false positives, which must be filtered out.

Furthermore, since we seek θ and the solution is in the form cos (θ), an even function, we must account for

all possible combinations. This yields as many as eight potential solutions for the roots, expressed in terms

of sin (θ) and cos (θ) :

'

cos (θ1) cos (θ2) cos (θ3) cos (θ4) cos (θ1) cos (θ2) cos (θ3) cos (θ4)

(

(37)

'

sin (θ1) sin (θ2) sin (θ3) sin (θ4) −sin (θ1) −sin (θ2) −sin (θ3) −sin (θ4)

(

(38)

The true roots in cos (θ) are obtained by passing (37) through the shadow function as given by Vallado

(39).15 Here the shadow function has been normalized by the central body radius to help with numerical

errors:

S(θ) = (1 + ecos (θ))2 +

- p2

r2

cb

.

(β1cos (θ) + β2sin (θ))2 −

- p2

r2

cb

.

(39)

Both Escobal’s quartic shadow function (33) and the form given in Vallado (39) will also have valid roots

on the illuminated side of the central body. These are filtered by assuring that we only take those roots for

values of true anomaly at which central body is between the spacecraft and the Sun. Using equation 5-3

from Vallado,15 the shadow entry and exits must be such that:

β1cos (θ) + β2sin (θ) ≤ 0 (40)

Now the roots of the shadow function in true anomaly are known, and the derivative of the shadow

function at the roots indicates whether the root is a shadow entry ( dS

dθ > 0) or shadow exit ( dSdθ < 0).

dS

dθ = 2 - p2r2

cb

.

(β1cos (θ) + β2sin (θ))(−β1sin (θ) + β2cos (θ)) − 2esin (θ) (1 + ecos (θ)) (41)

We can introduce breakpoints in the integral from (31) and evaluate the integral along each illuminated

or shadowed segment, applying thrust due to a solar electric propulsion system accordingly. For each interval

in a given orbit, the integrand in (31) is evaluated using a ninth-order Legendre-Gauss quadrature.

H. Control Parameterization

The forces due to thrust are determined by the guidance scheme of the spacecraft. In his non-averaged

scheme, Betts used a unitized thrust vector in the RSW frame as the control.12 This is not possible using

an orbital-averaging approach, since a single control vector is needed that will provide the appropriate

oscillatory behavior of the thrust direction within a single averaged orbit. Kluever showed that the costates

of the classical orbital elements can be used as the time-varying guidance parameters for an averaging

approach.8 Kluever’s approach (42) assumes we only have costates governing the ”slow” orbital elements,

and provides a thrust vector in the velocity vector-aligned frame (NTW):

aTMcoe ˆu =



























0 2 a2v

μ

0

rsin(θ)

a·v

2(e+cos(θ))

v 0

0 0 r cos(ω+θ)

h

0 0 r sin(ω+θ)

hsin(inc)

−2e+(r cos(θ)

a ) e·v

2sin(θ)

e·v

−rsin(ω+θ)cos(inc)

hsin(inc)







































an

at

aw













(42)

NASA/TM—2012-217699 9
Since our equations of motion are formulated using modified equinoctial elements, the classical elements

must first be computed (see reference 12). The guidance parameters are a given as a vector of the classical

orbital element costates, excluding true anomaly.

λcoe =

'

λa λe λi λΩ λω

(T

(43)

and the thrust unit vector in the NTW frame is:

ˆuNTW = −MT

coeλcoe

'MT

coeλcoe' (44)

In his work, Kluever divides λa by the semi-major axis and parameterizes the other costates as a function

of the semi-major axis.8 In this approach λa is normalized by the semi-major axis before being used in (44),

but all costates are functions of time. Notably, since the guidance parameters here are not functions of

semi-major axis, the posibility exists to use these controls for a single-phase ”round-trip” trajectory or other

trajectories which involve both increasing and decreasing the semi-major axis of the orbit.

We may derive a similar guidance scheme where the equinoctial equations of motion and costates are

used in place of the classical orbital element equations and costates.

aTMmee ˆu =































0 2p

w

)p

μ

0

)p

μsin (L)

)p

μ

[(w+1)cos(L)+f]

w −

)p

μ

g(jsin(L)−kcos(L))

w

−

)p

μcos (L)

)p

μ

[(w+1)sin(L)+g]

w

)p

μ

f(jsin(L)−kcos(L))

w

0 0 )p

μ

s2

2w cos (L)

0 0 )p

μ

s2

2w sin (L)











































ar

as

aw













(45)

The guidance parameters are a given as a vector of the modified equinoctial element costates, excluding true

longitude.

λmee =

'

λp λf λg λj λk

(T

(46)

and the thrust unit vector in the RSW frame is:

ˆuRSW = −MT

meeλmee

'MT

meeλmee' (47)

In the modified equinoctial element control formulation, λp is divided by the semi-latus rectum before

being passed to (47).

The modified equinoctial element-based guidance scheme is similar to that used in SEPSPOT1 and is

less prone to singularities, except when an orbit is nearly retrograde. However, if one only cares to control

semi-major axis, eccentricity, and inclination of the spacecraft, this requires only three time-varying control

variables (λa, λe, λi) in the case of the classical element formulation, but five time-varying control variables

(λp, λf , λg, λj , λk) in the case of of the modified equinoctial element formulation. The additional control

variables increase the size of the collocation problem, and may adversely affect runtime, though this effect

has not yet been quantified.

III. Results

The collocation formulation of a planetocentric low-thrust trajectory optimization problem shown here

is demonstrated with comparisons to the LEO to GEO and GTO to GEO cases in Kluever’s demonstration

of his solution method.8

NASA/TM—2012-217699 10
A. Minimum Time LEO to GEO

The minimum time LEO to GEO spiral transfer uses the same assumptions as those given by Kluever.8 The

spacecraft initially has a mass of 1200 kg and the propulsion system has a nominal power level of 10 kW,

a specific impulse of 3300 sec, and a propulsive efficiency of 65%. The initial time of the propagation is

January 1, 2000 and the initial state vector is given in Table 1.

Variable Initial Value

p 6927 km

f 1.0E-6

g 0

j 0.2539676

k 0

m 1200 kg

rev 0

Table 1. LEO to GEO initial variable values (fixed)

The solution shown here uses the classical orbital element costates as guidance parameters. The semimajor axis costate (λaa) is fixed at -1, giving a tangential thrust vector when λe = 0 and λi = 0. Costates

for right ascension of ascending node (λΩ) and argument of periapsis (λω) are fixed at 0.

The costates corresponding to eccentricity and inclination are time-varying optimal controls. Initially

the value of these controls is a linear fit of the guessed initial and final values.

Parameter Initial Value (guess) Final Value (guess) Lower Bound Upper Bound Scale Factor

λe 0 0.5 -1 1 100

λi 0 3 -10 10 100

Table 2. LEO to GEO time varying guidance parameters

Initially, the trajectory is simulated starting at the initial state for the guessed elapsed time using linear

fits to the time-varying controls. This explicit simulation provides values for the state variables at the

cardinal nodes in the phase. At this point, the solution is physically accurate (the defects are approximately

zero) but the constraints are not yet satisfied.

Variable Lower Bound Upper Bound Scale Factor

a 42164 (km) 42164 (km) 1

42164

e 0.0001 0.001 100

i 0 (deg) 0.01 (deg) 100

Table 3. LEO to GEO final boundary constraints

In addition to the boundary constraints above, a path constraint is used to ensure the value of periapsis

altitude is at least 300 km. Without this path constraint, the optimizer sometimes attempts to push periapsis

of the orbit below the surface of the Earth.

With the initial trajectory of the vehicle reasonably defined from the explicit simulation and the boundary

and path constraints in place, the solution is solved using IPOPT.9 Reasonable convergence was achieved

with the phase broken into 15 equal 3rd-order polynomial segments. The solution converged in 91 iterations.

The resulting minimum transfer time is given in Table 4 including comparisons to results obtained by Kluever

and SEPSPOT.

The time histories of semi-major axis, eccentricity, and inclination are shown in figure 2. They show good

agreement with the results achieved by Kluever.8 Figure 3 shows the time history for the eccentricity and

inclination costates used as optimal control variables. Both exhibit some ”wagging” at the beginning and end

of the phase. This indicates that the scaling in the collocation problem needs improvement, but the very good

agreement with other results suggests that the solution is insensitive to the values of the control variables at

NASA/TM—2012-217699 11
those times. Better scaling may reduce or eliminate the ”wagging” behavior and improve convergence but

not have a significant impact on the solution.

Figure 2. Semi-major axis, eccentricity, and inclination histories for the minimum-time LEO to GEO transfer.

Markers indicate cardinal nodes of the collocation problem.

Figure 3. Eccentricity and inclination costates used as the time-varying optimal controls for the minimum-time

LEO to GEO transfer.

Solution Technique Transfer Time (days)

Collocation 198.6

Kluever Control Parameterization8 199.0

SEPSPOT Result8 198.8

Table 4. LEO to GEO results

NASA/TM—2012-217699 12
Figure 4 shows the evolution of the orbit from LEO to GEO for the minimum time solution, with orbits

plotted at the collocation segment boundaries, roughly once every 13.2 days.

Figure 4. Evolution from the initial orbit (blue) to the final orbit (red) for the minimum time LEO to GEO

transfer.

B. Minimum Time GTO to GEO

As in the previous case, the minimum time GTO to GEO spiral transfer uses the same assumptions as those

given by Kluever.8 The spacecraft initially has a mass of 1200 kg and the propulsion system has a nominal

power level of 5 kW, a specific impulse of 1800 sec, and a propulsion system efficiency of 55%. The initial

time of the propagation is March 22, 2000 and the initial state vector is fixed to the values given in Table 5.

The costates corresponding to eccentricity and inclination are time-varying optimal controls with their

values linearly fit to guessed initial and final values (Table 6), and the initial state values at the cardinal

nodes are obtained through an explicit simulation. The final boundary constraints are the same as those

given in Table 3 and the path constraint on periapsis radius has the minimum value constrained to 185 km.

In this case the phase was broken into 10 equal 3rd-order polynomial segments. Despite the reduction

in the number of variables in the problem due to the reduced number of cardinal nodes, the GTO to GEO

transfer required 299 iterations to achieve convergence. The resulting minimum transfer time is given in

Table 7 with comparisons to results obtained by Kluever and results generated by SEPSPOT.

NASA/TM—2012-217699 13
Variable Initial Value

p 11359.07 km

f 0.7306

g 0

j 0.2539676

k 0

m 1200 kg

rev 0

Table 5. GTO to GEO initial time and state values (fixed)

Parameter Initial Value (guess) Final Value (guess) Lower Bound Upper Bound Scale Factor

λe 0.5 6 0 15 100

λi 1 10 -15 15 100

Table 6. GTO to GEO time-varying guidance parameters

Solution Technique Transfer Time (days)

Collocation 118.29

Kluever Control Parameterization8 118.36

SEPSPOT 118.29

Table 7. GTO to GEO results

Figure 5 shows the time histories of semi-major axis, eccentricity, and inclination for the minimum time

GTO to GEO solution, which again compare favorably with the results achieved by Kluever.8

Figure 5. Semi-major axis, eccentricity, and inclination history for the minimum-time GTO to GEO transfer.

Figure 6 shows the time history for the eccentricity and inclination costates used as optimal control

variables. Again, ”wagging” behavior is present, especially in the inclination costate. The insensitivity of

the result to the terminal value of the costates is a likely source of some of the convergence issues experienced

in this case. One possible reason for this is scaling. However, the fact that the behavior seems to be exhibited

where the orbit is nearly circular or near zero inclination suggest the singularities in the classical orbit element

control formulation may be responsible.

NASA/TM—2012-217699 14
Figure 6. Eccentricity and inclination costates used as the time-varying optimal controls for the minimum-time

GTO to GEO transfer.

Figures 7 and 8 show the evolution of the orbit from GTO to GEO for the minimum time solution. Orbits

are plotted at the collocation segment boundaries, or approximately once every 12 days.

Figure 7. Evolution from the initial orbit (blue) to the final orbit (red) for the minimum time GTO to GEO

transfer.

NASA/TM—2012-217699 15
Figure 8. Evolution from the initial orbit (blue) to the final orbit (red) for the minimum time GTO to GEO

transfer.

IV. Conclusions and Future Work

The results shown here indicate that a collocation technique, combined with orbital averaging of the

modified equinoctial elements, is capable of generating optimal solutions for planetocentric low-thrust orbit

transfer problems relatively quickly. While the results given here were set up for comparisons with Kluever’s

results,8 the implementation can be extended to include other forces on the spacecraft such as atmospheric

drag and perturbation by the Moon and Sun.

Despite the fact that this approach uses costates of the classical or modified equinoctial orbital elements,

the use of collocation instead of a shooting method generally allows for greater convergence. Rough guesses

at the guidance parameters can yield an optimal solution.

The issues of high nonlinearity in the terminal control conditions clearly need to be addressed. Better

scaling of the NLP will likely help, but experience with OTIS4 has shown that other measures may be

necessary. In the use of OTIS4, a common approach to such issues is to add control rate constraints, which

often forces the control time history to be more smooth. A grid refinement algorithm should also be put

into place to ensure that the converged solution adequately matches an explicitly integrated trajectory.

The ease of which this method can be extended to include other control formulations lends itself to

experimenting with more control laws in the future. Solving these problems using the modified equinoctial

element control formulations may also prove to be more robust due to the lack of singularities in those

equations. Other guidance laws which do not require high-frequency time-varying controls, such as Q-Law,16

will also be explored.

Modeling the spacecraft subsystem in greater fidelity will also be a priority. Degradation of the solar

arrays due to the impingement of trapped particles can be implemented as it was in SEPSPOT. Reference 1

includes a great amount of detail as to that implementation. The authors intend to implement more detailed

radiation modeling in the future using more recent and higher-fidelity models.

So far this method has been demonstrated for a single phase trajectory. Using a multiple phase approach

in which later phases employ non-averaged equations of motion may allow for fast, robust solutions to

planetocentric spiral trajectories which target lunar orbit or lagrange points in a terminal phase.

NASA/TM—2012-217699 16
References

1Edelbaum, T. N., Sackett, L. L., and Malchow, H. L., “Solar Electric Geocentric Transfer With Attitude Constraints:

Analysis,” NASA CR 134927, National Aeronautics and Space Administration, August 1975. 2Dickmanns, E. and Well, K., “Approximate solution of optimal control problems using third order hermite polynomial

functions,” Optimization Techniques IFIP Technical Conference Novosibirsk, July 17, 1974, edited by G. Marchuk, Vol. 27 of

Lecture Notes in Computer Science, Springer Berlin / Heidelberg, 1975, pp. 158–166. 3Hargraves, C. R. and Paris, S. W., “Direct trajectory optimization using nonlinear programming and collocation,” Journal

of Guidance, Control, and Dynamics, Vol. 10, Aug. 1987, pp. 338–342. 4Paris, S. W., Riehl, J. P., and Sjauw, W. K., “Enhanced Procedures for Direct Trajectory Optimization Using Nonlinear

Programming and Implicit Integration,” Proceedings of the AIAA/AAS Astrodynamics Specialists Conference and Exhibit, No.

2006–6309, 2006. 5Ross, I., “User’s Manual for DIDO: A MATLAB Application Package for Solving Optimal control Problems,” Tech. rep.,

TOMLAB Optimization, Sweden, 2004. 6Betts, J. T. and Huffman, W., “Sparse Optimal Control Software SOCS,” Mathematics and Engineering Analysis Technical Document MEA-LR-085, The Boeing Company, Seattle, WA, 1997. 7Betts, J. T. and Erb, S. O., “Optimal Low Thrust Trajectories to the Moon,” SIAM Journal on Applied Dynamical

Systems, Vol. 2, January 2003, pp. 144–170. 8Conway, B. A. and Kluever, C. A., Spacecraft Trajectory Optimization, chap. 5, Cambridge University Press, New York,

New York, 2010, pp. 113–137. 9Wachter, A. and Biegler, L. T., “On the implementation of an interior-point filter line-search algorithm for large-scale

nonlinear programming,” Mathematical Programming, Vol. 106, 2006, pp. 25–57. 10Roy, A., Orbital Motion, A. Hilger, 1988, p. 192. 11Walker, M. J. H., Ireland, B., and Owens, J., “A Set of Modified Equinoctial Elements,” Celestia Mechanics, Vol. 36,

July 1985, pp. 409–419. 12Betts, J., Practical Methods for Optimal Control and Estimation Using Nonlinear Programming, Advances in Design

and Control, Society for Industrial and Applied Mathematics, 2010. 13Kechichian, J. A., “The treatment of the earth oblateness effect in trajectory optimization in equinoctial coordinates,”

Acta Astronautica, Vol. 40, No. 1, 1997, pp. 69 – 82. 14Escobal, P., Methods of Orbit Determination, J. Wiley and Sons, New York, New York, 1985, pp. 155–159. 15Vallado, D. A. and McClain, W. D., Fundamentals of Astrodynamics and Applications, chap. 5, Springer, New York,

New York, 2007, pp. 304–306. 16Petropoulos, A. E., “Low-Thrust Orbit Transfers Using Candidate Lyapunov Functions with a Mechanism for Coasting,”

Proceedings of the AIAA/AAS Astrodynamics Specialists Conference and Exhibit, No. 2004–5089, 2004.

NASA/TM—2012-217699 17
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-09-2012 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Optimization of Low-Thrust Spiral Trajectories by Collocation 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Falck, Robert, D.; Dankanich, John, W. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 346620.02.01.11.01.01 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-18386 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2012-217699 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 13 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

As NASA examines potential missions in the post space shuttle era, there has been a renewed interest in low-thrust electric propulsion for 

both crewed and uncrewed missions. While much progress has been made in the field of software for the optimization of low-thrust 

trajectories, many of the tools utilize higher-fidelity methods which, while excellent, result in extremely high run-times and poor 

convergence when dealing with planetocentric spiraling trajectories deep within a gravity well. Conversely, faster tools like SEPSPOT 

provide a reasonable solution but typically fail to account for other forces such as third-body gravitation, aerodynamic drag, solar radiation 

pressure. SEPSPOT is further constrained by its solution method, which may require a very good guess to yield a converged optimal 

solution. Here the authors have developed an approach using collocation intended to provide solution times comparable to those given by 

SEPSPOT while allowing for greater robustness and extensible force models.

15. SUBJECT TERMS

Trajectory optimization; Low-thrust; Collocation; Orbital mechanics 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

24 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18