University of Southern California Dissertations and Theses

Spacecraft trajectory optimization: multiple-impulse to time-optimal finite-burn conversion

(USC Thesis Other)

Spacecraft trajectory optimization: multiple-impulse to time-optimal finite-burn conversion

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content SPACECRAFT TRAJECTORY OPTIMIZATION: MULTIPLE-IMPULSE TO TIME-OPTIMAL FINITE-BURN CONVERSION
A Dissertation
by
JOSHUA FOGEL
Presented to the FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In partial fulfillment of the requirements for the degree
DOCTOR OF PHILOSOPHY
Department of Astronautical Engineering
Viterbi School of Engineering
Advisory Committee:
Mike Gruntman (Chair)
Joseph Wang
Henryk Flashner
December 2019 Copyright 2019 Joshua Fogel
ABSTRACT
A novel algorithm is presented that provides an improvement over traditional parameter
optimization methods when solving time-optimal, finite-burn pseudo-rendezvous spacecraft
trajectory problems. A hybrid optimization procedure is described that converts a set of multiple-
impulses, representing high- or low-thrust maneuvers, to an exact time-optimal finite-burn
trajectory for a thrust limited, constant exhaust velocity spacecraft. The Hybrid Method applies a
control law derived from the Euler-Lagrange system of equations within the classical Indirect
Method to a modern Direct Method. An iterative adjoint-control transformation and an evolving
constraint vector are introduced to solve the optimal control two-point boundary value problem.
This method requires no prior knowledge of the solution, which adds simplicity to the trajectory
design process and aids automation. Examples are shown for low-thrust apogee raise maneuvers,
non-coplanar Earth orbit transfers, a two-burn Earth-Mars transfer, and a modified Deep Space 1
low-thrust trajectory. A numerically significant improvement to objective cost is shown across all
application problems compared to a traditional solution method, as well as a significant
improvement to convergence speed for a select class of problems.
!2
DEDICATION
To my parents —
Whose love and support launched me on the right trajectory to achieve my dreams.
To my friends —
Whose bright light helped guide me through the darkest of times.
And to Ivonne —
Who believed in me when I could not,
Supported me when I fell,
and Challenged me to be better.
I would not have made it here without you.
!3
ACKNOWLEDGEMENTS
Thank you to my academic advisor Mike Gruntman for shepherding me through two
degrees and finally out the other side. Your guidance has helped correct my trajectory many
times over the years, and for that I will always be grateful.
My sincerest appreciation to everyone in the ASTE department for always being there to
help me out of a jam. Thank you for being amazing!
And finally, a huge thank you to the NASA-JSC Pathways Office, USRA Office, and the
people of the EG5 Branch. Frank, Jeff, Jacob, Tim, Max, Amelia, Jerry, Bharat— you have been
some of the best colleagues I have ever had the honor of working with. Thank you for helping
me grow as a budding astrodynamicist, and showing me that there is indeed a bright light at the
end of the tunnel. Thank you for giving me a chance. Ricky— thank you for going above and
beyond as a mentor, and for helping me navigate through tough seas. I hope to one day be as
good of a role model as you. Cesar— thank you for your willingness to take me under your wing
and guide me to my topic. You helped open my eyes to a fascinating realm of optimization that I
was excited to explore. !4
TABLE OF CONTENTS
ABSTRACT 2 ....................................................................................................................................
DEDICATION 3 .................................................................................................................................
ACKNOWLEDGEMENTS 4 .............................................................................................................
LIST OF FIGURES 7 .........................................................................................................................
LIST OF TABLES 9 ...........................................................................................................................
ABBREVIATIONS 10 .......................................................................................................................
1. INTRODUCTION 11 .....................................................................................................................
1.1 Introduction to Low-Thrust Trajectory Optimization 11 ...............................................................
1.2 Indirect vs. Direct Methods 11 ......................................................................................................
1.3 Motivation: Bridging the Gap 13 ...................................................................................................
1.4 Research Goals 15 .........................................................................................................................
1.5 Advancing the State of the Art 15 ..................................................................................................
2. PROBLEM STATEMENT 16 .......................................................................................................
2.1 Overview of the Low-Thrust Trajectory Optimization Problem 16 ..............................................
2.2 The Impulsive Burn Model 18 .......................................................................................................
2.3 The Finite-Burn Optimal Control Problem 22 ...............................................................................
2.4 The Multi-Impulse Finite Burn Conversion Problem 27 ..............................................................
3. PRIOR WORK 30 ..........................................................................................................................
3.1 Survey of Existing Optimization Methods 30 ...............................................................................
3.2 Single Impulse-to-Finite Burn Conversion 42 ...............................................................................
4. EXISTING GAPS 45 .....................................................................................................................
4.1 General Gaps in Finite-Burn & Low-Thrust Optimization 45 ......................................................
4.2 Specific Gaps within the Impulse-to-Finite Burn Conversion Problem 52 ...................................
5. HYPOTHESES 55 ..........................................................................................................................
6. INITIAL STUDIES 57 ...................................................................................................................
6.1 The MID-MPCM Method for Solving Perturbed Two-Body Problems 57 ...................................
7. MULTI-IMPULSE TO FINITE-BURN TRAJECTORY CONVERSION 73 ...............................
7.1 Assumptions 73 ..............................................................................................................................
!5
7.2 The Multiple Finite-Burn Conversion Algorithm 76 .....................................................................
7.3 Application Problems 86 ...............................................................................................................
7.4 Proof of Numerically Significant Improvement 109 .....................................................................
8. SUMMARY AND CONCLUSION 113 .........................................................................................
BIBLIOGRAPHY 116 .......................................................................................................................
!6
LIST OF FIGURES
Figure 1: Trajectory Optimization Schema....................................................................................17
Figure 2: Single Impulse-to-Finite-Burn Conversion via Pseudo-Rendezvous Targeting.............19
Figure 3: Multi-Impulse to Multi-Finite-Burn Conversion...........................................................27
Figure 4: Multi-Impulse to Multi-Finite-Burn Conversion, with Finite-Burn Arc Overlap..........29
Figure 5: The Two-Body Multi-Impulse Trajectory Structure Used in MALTO..........................31
Figure 6: Example Exponential-Sinusoid Shapes..........................................................................33
Figure 7: Intersections of an Exponential-Sinusoid Arc in the Ecliptic.........................................34
Figure 8: Three- and Four-Body Problems....................................................................................35
Figure 9: Sample Low-Thrust Capture Trajectory at Mars via Invariant Manifold......................36
Figure 10: Possible Maneuver Combinations for a Segment in Copernicus.................................37
Figure 11: An Example Trajectory in Copernicus, Comprised of Multiple Segments..................38
Figure 12: Example Edelbaum LEO-to-GEO Continuous Low-Thrust Transfer..........................42
Figure 13: Converged GEO-to-GEO Non-Coplanar Transfer using MID-MPCM.......................62
Figure 14: Comparison of Thrust-Acceleration Unit Vectors for a LEO-to-GEO Transfer...........64
Figure 15: Structure of the Initial Guess Trajectory Generation Algorithm via Edelbuam...........66
Figure 16: Example LEO-to-LLO Patched-Edelbaum Initial Guess Trajectory............................68
Figure 17: Parametric Rotation of Initial Mean Anomalies of Edelbaum Spirals.........................70
Figure 18: Initial Guess SOI-Crossing Trajectories.......................................................................71
Figure 19: Conceptual Overview of Multi-Impulse to Multi-Finite-Burn Conversion.................78
Figure 20: Conversion Algorithm Burn Model Parameterizations................................................80
Figure 21: Implementation of Single Impulse Conversion, with added Pseudo-Rendezvous.......81
Figure 22: Multi-Finite-Burn Conversion Algorithm....................................................................82
Figure 23: Artist depiction of the Deep Space Transport...............................................................86
Figure 24: Solution Comparison of Two-Burn, Apogee Raise Conversion...................................88
Figure 25: Change in Radius and Velocity vs. Time for Apogee Raise Conversion.....................89
Figure 26: Two-Impulse Inclination Change Conversion, Optimization Sub-Problem 1..............90
Figure 27: Finite-Burn Arc Overlap Problem for Two-Impulse Inclination Change.....................92
!7
Figure 28: Mitigated Finite-Burn Arc Overlap (Optimization Sub-Problem 3 Solution)..............92
Figure 29: Feasibility and Optimality Metrics vs Iteration for Hybrid OCT: SP-1 Solution.......96
Figure 30: Feasibility and Optimality Metrics vs Iteration for Hybrid OCT: SP-2 Solution.......97
Figure 31: Comparison of Objective Function Value vs Iteration.................................................98
Figure 32: Feasibility Metric vs Iteration for Direct SOC: SP-2 Solution....................................98
Figure 33: Feasibility & Optimality Metrics vs Iterations for Direct SOC: SP-2 Solution..........99
Figure 34: Two-Burn, Non-Coplanar, Earth-Mars Low-Thrust Conversion, Sub-Problem 1.....100
Figure 35: Two-Burn, Non-Coplanar, Earth-Mars Low-Thrust Conversion, Sub-Problem 3.....102
Figure 36: Switching Function for the Earth-Mars Low-Thrust Conversion, Sub-Problem 3....103
Figure 37: Impulsive Approximation of Modified DS1 Trajectory, with increased ∆V ..............106
Figure 38: Modified DS1 Trajectory, SP-1 Hybrid OCT Solution..............................................107
Figure 39: Modified DS1 Trajectory, SP-2 Hybrid OCT and Direct SOC Solutions..................108
Figure 40: Flown Deep Space 1 Trajectory to Asteroid 1992 KD (now called 9969 Braille).....109
!8
LIST OF TABLES
Table 1: Survey of Indirect & Direct Methods for Low-Thrust Trajectory Optimization.............12
Table 2: Comparison of Final Orbital Elements, Total ∆V & Flight Time for GEO-GEO............62
Table 3: Comparison of Final Orbital Elements, Total ∆V & Flight Time for LEO-GEO............64
Table 4: Optimization Sub-Problem 1...........................................................................................83
Table 5: Optimization Sub-Problem 2...........................................................................................84
Table 6: Optimization Sub-Problem 3...........................................................................................85
Table 7: Two-Burn, Low-Thrust Apogee Raise Solution Burn Parameters...................................88
Table 8: Two-Burn, Low-Thrust Inclination Change Sub-Problem 1 Solution Parameters..........91
Table 9: Two-Burn, Low-Thrust Inclination Change Sub-Problem 3 Solution Parameters..........93
Table 10: Comparison of Optimization Parameters Between Hybrid OCT and Direct SOC........94
Table 11: Two-Burn, Non-Coplanar, Earth-Mars Low-Thrust Sub-Problem 1 Solution.............100
Table 12: Two-Burn, Non-Coplanar, Earth-Mars Low-Thrust Sub-Problem 1 Solution.............102
Table 13: Modified DS1 Trajectory, SP-2 Hybrid OCT and Direct SOC Burn Parameters........108
Table 14: Comparison of Optimization Convergence Violations................................................112
!9
ABBREVIATIONS
FB Finite-Burn
LEP Low-thrust Electric Propulsion
SEP Solar Electric Propulsion
LEO Low Earth Orbit
GEO Geosynchronous Orbit
LTGA Low-Thrust Gravity Assist
TPBVP Two-Point Boundary Value Problem
IVP Initial Value Problem
NLP Non-Linear Programming
∆V Delta-Velocity
V∞ Hyperbolic Excess Velocity
SNOPT Sparse Nonlinear Optimizer
[26]

CRTBP Circularly Restricted Three-Body Problem
MID-MPCM Multi-Impulse Discretization via the Modified Picard-Chebyshev Method
ECI Earth Centered Inertial
E/M-SOI Earth/Moon Sphere of Influence
FPA Flight Path Angle
SOC Sub-Optimal Control
OCT Optimal Control Theory
P-R Pseudo-Rendezvous
SQP Sequential Quadratic Programming
SP Sub-Problem
DS1 Deep Space 1
NEA Near Earth Asteroid !10
1. INTRODUCTION
1.1 Introduction to Low-Thrust Trajectory Optimization
Recent technological advancements in the realm of low-thrust electric propulsion (LEP)
have allowed such systems to offer significant performance benefits over other in-space
propulsion systems in certain applications. Shorter flight times, increased payload mass fraction,
and smaller launch vehicles are all potential benefits to LEP missions
[3]
. The growing dominance
of LEP systems requires the development of robust, accurate, and practical methods of
computing optimal low-thrust trajectories, and more generally, finite-burn trajectories. The
performance index typically minimized in these problems is transfer time or propellant
expenditure.
Unfortunately, finding global optimal solutions for these problems is a formidable
challenge, as no analytic, closed-form solutions exist for the full perturbed N-body problem.
While many indirect & direct numerical optimization methods have been developed
[40]
, they are
often tailored for specific transfer problems, with sensitive convergence properties as the
objective functions are often highly non-linear. This makes many existing methods difficult to
generalize on their own, which restricts the total change in inclination and eccentricity that can
be computed for a given set of maneuvers. However several of these methods can be readily
applied as initial guess solvers for a larger, more generalized optimization routine.
1.2 Indirect vs. Direct Methods
The methods commonly used to solve general, finite-burn trajectory optimization
problems, including low-thrust problems, fall into two main categories: indirect and direct
methods
[1, 2, 10]
. Table 1 showcases a selection of Indirect and Direct optimization methods
applicable to the finite-burn problem, which will be explored further in Section 2.1. The
selection of appropriate solution methods depends on the problem formulation, how much time is
willing to be taken to implement, and access to relevant software tools.
!11
Indirect methods are based on calculus of variations, and convert the trajectory into a
parameter optimization two-point boundary value problem (TPBVP) with N unknowns and N
constraints, solved using a non-linear programming (NLP) method
[10]
. The controls are
eliminated from the problem, and indirectly solved as feedbacks of adjoint Lagrange multipliers
known as “co-states”. This is where the term Indirect Method is derived. The solution to this
TPBVP is exactly the optimal solution for the specified initial conditions, objective function, and
final conditions (i.e. targeting constraints). However, finding such solutions in general is
difficult, as the trajectory requires six time-varying co-state variables (for a three-dimensional
problem expressed in Cartesian form), with unintuitive initial values, to be propagated along
with the spacecraft’s state. The convergence behavior of this problem is extremely sensitive to
!12
Indirect Methods Direct Methods
Kechichian
[4]
analyzed the low-thrust minimum time
rendezvous problem for constant acceleration spacecraft
- Used UNCMIN (Unconstrained Minimization Package) to calculate 3D LEO-GEO transfer trajectories
- Method involves propagating six costate parameters,
whose initial values are tuned to the specific problem
MALTO/GALLOP
[6,17,18]
: Interplanetary optimization
algorithms that numerically simulate continuous thrust
trajectories with sets of impulses and match points
- 2-body, instantaneous gravity assists, medium fidelity
Thorne
[5]
derived a minimum time, constant thrust orbit
transfer algorithm with non-circular boundary conditions
- Planar, 2-body assumption yields analytic expressions of Lagrange costates
- Solution is still sensitive to the initial guess values of the
costates
Mystic
[12]
: Non-linear optimal control method for
optimizing static & dynamic variables simultaneously
- Automatically finds and exploits gravity assists
- Trajectories are fully integrated using finite burns, multi- body propagation, and select perturbations; high fidelity
- Produces continuous trajectories with no “match” points
Nah and Vadali
[16]
investigated 3D, fuel optimal Earth-
Mars trajectories
- Used multiple-shooting approach and Newton’s method to optimize individual trajectory segments
STOUR-LTGA
[6,18]
: Shape-based algorithm for
automated design of initial guess LTGA trajectories
- Planar 2-body (extended to 3D); models gravity assists as instantaneous perturbations to an exponential sinusoid
SEPTOP
[2]
: Uses calculus of variations to form the state
& costate equations, integrated numerically to solve the
two-point BVP
GPOPS-II
[11]
: Direct collocation method for solving low
thrust orbit transfer optimization problems
- Uses NLP solver IPOPT
TfMin
[38]
: 3D minimum time orbit transfer, with free
final longitude
- Numerical single & multiple-shooting capability
- Automatic differentiation & BVP generation
MID-MPCM
[1]
: Multi-Impulse Discretization via the
Modified Picard-Chebyshev Method
- Finds minimum time orbit transfer trajectory using a series of impulses at intermediate waypoints
- Uses NLP solver SNOPT
[26]
and GPOPS-II
[11]

- Includes numerical perturbation and N-body force models
ITOP/HYTOP
[39]
: Indirect Trajectory Optimization
Program / Hybrid Trajectory Optimization Program
- Low-thrust orbit transfer optimization
- Developed by Aerospace Corp. (restricted use)
COPERNICUS
[2, 45]
: N-body, multiple shooting and
direct integration for targeting & state propagation
- Multi-spacecraft modeling, capable of optimization for both constant and variable specific impulse trajectories
Table 1: Survey of indirect & direct methods for low-thrust trajectory optimization
[40]
.
the initial values of these parameters, which is made worse by the fact that the co-states hold no
intuitive mapping to the physical dynamics. Adding gravity assists and other perturbations to the
trajectory compounds this sensitivity
[17]
. While successfully crafted indirect methods yield high
fidelity solutions, they are not often used during the conceptual design phase as they are difficult
to implement within automated, global optimization tools due to long execution times, small
regions of convergence, and their highly specified nature
[2]
.
Conversely, Direct Methods parameterize the optimization problem and use NLP solvers
to directly optimize a performance index (i.e. objective function) via the actual control. Direct
methods can be applied to a wide variety of optimal control problems
[2]
, and can be based on
embedded collocation, direct transcription, and differential inclusion schemes
[13-15]
. Born out of
the computer age, they rely on the power of modern processors to solve optimization problems
with “brute force” computations. Advantages to using direct methods include coding efficiency
and more robust convergence, as the solutions are less sensitive to the initial guesses compared
to indirect methods
[2]
. Direct methods are also more easily adaptable to complex force models,
including N-body gravity models and J2 perturbations, due to the problem formulation keeping
the controls (typically thrust directions and magnitude) independent of the dynamics
[40]
.
The main disadvantage with using direct methods with continuous low-thrust problems is
that they discretize the problem into many unknowns, e.g. impulsive nodes, or are controlled by
shape-based steering control laws which impose a pre-defined structure to the solution, meaning
they are inherently sub-optimal. However, since the accuracies of direct methods are typically
acceptable, and generally have more physically intuitive initial guess trajectories, they are most
often applied in automated global search algorithms
[2]
, and are the basis for most modern
trajectory design.
1.3 Motivation: Bridging the Gap
This dissertation is motivated by the current challenge in astrodynamics of bridging the
gap between Indirect and Direct optimization methods. Most industry astrodynamicists rely on
Direct Methods for designing spacecraft trajectories due to their adaptability, ease of use, and
!13
relative high fidelity solutions. There is a desire, however, to find improved solutions which
yield a savings in flight time, propellant, or other mission resources, or can be computed faster of
more efficiently. In order to leverage the highest fidelity solutions available from Indirect
Methods, a new methodology must be derived.
The research presented herein provides an improvement over traditional parameter
optimization methods when solving time-optimal, finite-burn spacecraft trajectory problems
based on existing solutions computed by readily available techniques. The new conversion
algorithm takes an approximate finite-burn trajectory comprised of n-number of discretized,
impulsive nodes, or sub-optimal parametrically-steered finite burns, and through automated
procedure, replaces the approximate trajectory with a set of time-optimal finite-burns separated
by ballistic arcs, whose pre- and post- maneuver states match that of the original approximation.
The methodology of this conversion algorithm can be broken down into two cases: the
single impulse, and the multi-impulse. Inspired by existing state-of-the-art work
[51]
, the single
impulse case formulates an optimal control problem that is treated as a bounded, pseudo-
rendezvous, targeting initial and final states upstream and downstream of a single ∆V on the
initial impulsive trajectory. This becomes a TPBVP, which can be solved by applying the control
law from the classical indirect method, utilizing calculus of variations and primer vector theory
to minimize the Hamiltonian with respect to the controls.
The multi-impulse case is an extension of the first to a more complex conversion
problem. This new approach chains together single-impulse conversions, separated by ballistic
arcs, thereby forming a larger parameter optimization problem around the sequence of
individually solved burn arcs. To complete the conversion, focus must be given to mitigating the
finite-burn arc overlap problem, which occurs when interior ballistic arcs between optimized
burns result in negative times of flight. This is mitigated by evolving the boundary constraints
across multiple sub-problems as part of an iterative adjoint-control transformation. The hybrid
conversion algorithm requires no prior knowledge of the solution beyond the supply of an
approximate finite-burn trajectory, which adds simplicity to the design process and aids
automation.
!14
1.4 Research Goals
• Improve upon the traditional Direct Method by applying an Indirect Method’s control law
• Formulate the optimal control problem for multi-impulse conversion, using pseudo-
rendezvous targeting on the impulsive initial guess
• Mitigate the finite-burn overlap issue by introducing additional constraints on the ballistic
arcs, and iteratively introducing the full constraint vector
1.5 Advancing the State of the Art
The presented work advances the state of the art with regards to enabling a robust
capability that is otherwise unavailable to industry mission designers, promising a moderate
savings in objective cost (flight time or propellant) when compared to conventional Direct
Methods for a wide variety of trajectory problems. Moreover, the new Hybrid Method is shown
to yield superior convergence behavior in terms of iterations and convergence time for
conversion problems where both ! , and the burn is changing inclination. This
specific class of problem will therefore reap larger benefits than other problems when the Hybrid
Method is employed in lieu of the classical Direct.
The presented algorithm also advances the state of the art by solving the finite-burn arc
overlap problem via re-targeting and re-defining the optimization problem through an evolving
constraint vector. This process improves all solutions, regardless of the existence of any overlaps,
providing substantial objective cost savings in the output solution compared to earlier sub-
problem solutions. The method requires no prior knowledge of the solution, merely a multi-
impulse or sub-optimal finite-burn initial guess, which aids automation.
The following text discusses these advancements in detail. The results presented herein
were completed in preparation for presentation at the 2020 AIAA SciTech Forum. Parts of this
thesis were presented at the 2015 International Astronautical Congress
[40]
, and presented with
special publication at the 2015 International Conference on Computational & Experimental
Engineering and Sciences
[1]
. ΔT
burn
>
T
Period
2
!15
2. PROBLEM STATEMENT
One of the fundamental questions in trajectory design is: “What is the best way to steer
my rocket to get me where I want to go using the least amount of fuel?” Or, in more formal
parlance, “What is the required time history of thrust magnitudes and directions necessary to
minimize a propellant or time cost?” The problem sees us taking a spacecraft with some
propulsion model (e.g. impulsive, finite-burn, power limited), subject it to some force field
(dynamics), and we wish to fly it from one place to another (a specified state, Keplerian orbit, or
pre-defined ephemeris), in some optimal sense (parameter or functional optimization, minimize
∆V or time of flight, maximize final spacecraft mass), and under some conditions (equality and
inequality constraints, boundary constraints, path constraints, control constraints).
While many types of orbital maneuvers are able to be represented by either simple
impulsive ∆Vs or parameterized sub-optimal thrust control laws, these modeling assumptions are
often only sufficient for burns whose durations are significantly less than the spacecraft’s
instantaneous orbital period about the primary central gravitating body. For finite-burns with
durations that span a significant fraction of their orbital period (e.g. low-thrust spiral
trajectories), or for the case where an exact time history of optimal thrust directions and
magnitudes is desired, we must turn to Optimal Control Theory (OCT).
2.1 Overview of the Low-Thrust Trajectory Optimization Problem
The low-thrust problem, and more generally the low-thrust N-body problem, is a
fundamental challenge in spacecraft trajectory design. The goal is to determine the time- or fuel-
optimal trajectory necessary for a spacecraft to achieve a target orbit or state. A procedure for
solving this kind of optimization problem is illustrated in Figure 1. Typical performance indexes
include maximizing the final spacecraft mass, minimizing time of flight or total fuel usage (or
alternatively total ∆V), or minimizing some other user defined infeasibility, i.e. the magnitude of
constraint violation. In the case of continuous thrust-acceleration problems, minimizing time of
flight yields the same solution as minimizing ∆V .
!16
Optimization is used to determine the required thrust profile (thrust directions and
magnitudes), i.e. the control, in order to satisfy a constraint vector and minimize an objective
function for a spacecraft under the influence of some gravity model. It should be noted that this
variety of NLP problem is highly non-convex, meaning a trajectory optimized for a set of
solution invariant parameters does not guarantee a globally optimal solution over the complete
trajectory design space
[2]
. While many optimization methods require a feasible start (i.e. an
initial trajectory which satisfies the constraint vector) to converge on a solution, the resulting
locally optimal solutions must then be checked for global optimality.
For problems involving multiple planetary encounters (i.e. gravity assists or planetary
captures), the sequence of fly-bys must either be provided or known a-priori, or determined in
situ by the optimization tool. In some solution methods, the sequence itself may be included as
design variables to be optimized, however the sequence is often defined prior to optimization
based on previously known feasible configurations, as well as trial and error. As with the single
Low-Thrust Gravity Assist (LTGA) problem, automated broad-search algorithms have been
developed for multi-LTGAs
[2, 6, 7, 8]
. However, a full search of the design space during
preliminary optimization can be computationally prohibitive, depending on the complexity of the
problem and chosen force models
[2]
. This means the trajectory analyst must bring prior
knowledge to the design process to inform the search if such methods are to be used.
The following sections present the single impulse-to-finite burn conversion theory
developed by Jesick and Ocampo
[51]
. This serves as the foundational theory upon which the
presented multi-impulse conversion algorithm is derived. Sections 2.2 and 2.3 detail the
established single impulse to finite-burn conversion problem, while Section 2.4 describes the
new problem of multi-impulse to finite-burn conversion, including the associated finite-burn arc
overlap problem. Note that in this text, a bold letter with an arrow (e.g. ! ) denotes a vector.
⃗ X
!17
Formulate
Problem
Generate Feasible Ini2al Guess
Trajectory
Solve Op2mal Control Problem:
find thrust direc-ons that
minimize a performance index
Figure 1: Trajectory Optimization Schema: A common process for conducting trajectory optimization.
2.2 The Impulsive Burn Model
Before expressing the finite-burn optimal control problem, we must first define the
impulsive maneuver that we wish to convert into a finite-burn (as shown in Figure 2). Consider
the state vector:
! (2.1)
with associated state equations:
! (2.2)
where ! is the gravitational acceleration term, ! is the thrust acceleration term, and ! is
the parameterized thrust direction unit vector defined by spherical angles ! and ! :
! (2.3)
Additional perturbing accelerations may be included in the state dynamics as well, but for the
simplicity of this work, only two-body gravitational and thrust acceleration terms are included.
The initial conditions at t0 are , and m0. The spacecraft’s thrust magnitude is bounded:
! (2.4)
This model assumes two arcs of pure ballistic motion (T = 0), separated by discontinuities in m
and ! at some interior point ! , representing an impulsive change in velocity and mass. This
impulse occurs at position ! . The states immediately before and after the impulse are
! , and ! , respectively.
⃗ X =
[
⃗ r
⊤
⃗ v
⊤
m]
⊤
·
⃗ X =
⃗ v(t)
⃗ g(t, ⃗ r )+
T
m
̂ u(t)
−
T
c
⃗ g(t, ⃗ r)
T
m
̂ u(t)
α β
̂ u(t) =
cosα cosβ
sinα cosβ
sinβ
⃗ r
0
, ⃗ v
0
0≤ T≤ T
max
⃗ v t
i
≥ t
0
⃗ r
i
[
⃗ r
i
⃗ v
i
−
m
−
i ]
⊤
[
⃗ r
i
⃗ v
i
+
m
+
i
]
⊤
!18
The impulsive ! is then defined as:
! (2.5)
! (2.6)
! (2.7)
with the rocket’s effective exhaust velocity c. Since ! , the instantaneous change in mass
! , such that:
! (2.8)
The final time ! is appropriately chosen such that the yet-to-be-found finite-burn with T > 0
that will replace the impulse will terminate at some ! .
⃗ ΔV
⃗ v
i
+
= ⃗ v
i
−
+
⃗ ΔV
i
m
+
i
= m
−
i
+e
Δv
i
c
m
−
i
= m
0
ΔV
i
> 0
Δm
i
< 0
Δm
i
= m
+
i
−m
−
i
t
f
> t
i
t
i
f
≤ t
f
!19
Figure 2: Single Impulse to Finite-Burn Conversion via Pseudo-Rendezvous targeting.
The pre-impulse ballistic trajectory is defined as a continuous arc from t0 to ti, and the
post-impulse ballistic trajectory is defined as a continuous arc from ti to . This is important, as
these pre- and post-impulse ballistic arcs (expressed as propagated functions of time) will
become the boundary targeting constraints once the finite-burn TPBVP is formulated. Herein the
point-mass spacecraft on the pre-impulse and post-impulse trajectories will be referred to as the
chase and target spacecraft, respectively.
The goal is to replace the impulse with a finite burn with bounded thrust magnitude ! . In
general, ! is typically taken to be a control variable. However, it is well known that for any
minimum time finite-burn problem,
[41]
.
We seek an estimate for the finite burn duration:
! (2.9)
Estimates of the finite burn start and stop times can be similarly computed, assuming the time of
impulse ti is also at the midpoint of the finite burn:
! (2.10)
! (2.11)
The state on the post-impulse trajectory at ! is ! , where ! . The total
transfer time is ! , and the total finite burn time is ! . The optimization variables are
then simply ! and ! , since the target state is a function of time.
The optimization problem is to find the minimum-time finite burn maneuver that begins
on the pre-impulse trajectory at ! and matches exactly the position and velocity of the target
spacecraft at ! . Since the trajectory of the target spacecraft is defined by a pair of ballistic arcs
separated by a pre-defined discontinuity in velocity and mass, its state is a known function of
t
f
T
T
T = T
max
Δt
i
=
Δm
i
T
max
c
t
i
0
= t
i
−
Δt
i
2
t
i
f
= t
i
+
Δt
i
2
t
i
f [
⃗ r*
i
f
⃗ v*
i
f
m*
i
f]
m*
i
f
= m
+
i
t
i
f
−t
0
t
i
f
−t
i
0
t
i
0
t
i
f
t
i
0
t
i
f
!20
time: ! and ! . As such, the finite-burn optimization problem is more accurately described
as a Pseudo-Rendezvous (P-R) problem, as the boundary conditions are known functions of time.
The unknowns of our problem are the start and stop times of the finite burn, ! and ! ,
respectively, and the thrust direction as a function of time ! for ! . The target
conditions for the finite-burn (denoted with *) are:
! (2.12)
! (2.13)
The cost functional J which we seek to maximize is chosen to be a negative quantity:
! (2.14)
for a given ! , and c.
While this method can be derived using any general, time and position dependent gravity
field ! , a central body inverse square force model is employed for simplicity. The two-
body gravitational acceleration is defined as:
! (2.15)
with its associated gravity gradient matrix:
! (2.16)
This minimum time, pseudo-rendezvous problem, with the optimal control function ! and the
finite-burn start and stop times ( ! and ! ), can be solved using direct or indirect methods. The
application of the direct method is a parameter finite-burn model which uses a simple set of
parameters to define the thrust directions ! , whose sub-optimal solution can be found using
standard NLP algorithms via multi-shooting numerical integration. The indirect method leads to
a TPBVP associated with the optimal control problem. Either a nonlinear equation algorithm or
an NLP algorithm can be used to solve for the time-optimal controls.
⃗ r(t) ⃗ v(t)
t
i
0
t
i
f
̂ u(t) t
i
0
≤ t ≤ t
i
f
⃗ r(t
i
f
) =
⃗ r*(t
i
f
)
⃗ v(t
i
f
) =
⃗ v*(t
i
f
)
J = max k(t
i
0
−t
i
f
)
k > 0, ⃗ g(t, ⃗ r), t
0
,
⃗ X
0
, t
i
,
⃗ ΔV
i
, T
max
⃗ g(t, ⃗ r )
⃗ g(t, ⃗ r ) =−
μ
r
3
⃗ r
G( ⃗ r ) =
∂ ⃗ g( ⃗ r )
∂ ⃗ r
=
μ
r
3
(
3 ⃗ r ⃗ r
⊤
r
2
−I)
̂ u(t)
t
i
0
t
i
f
̂ u(t)
!21
2.3 The Finite-Burn Optimal Control Problem
Using established theory, one of the fundamental optimal control problems in orbital
mechanics can be formulated
[41-43]
. Expressing the cost functional in the Mayer form:
! (2.17)
Assuming there exists some u*(t) that maximizes J, then the Euler-Lagrange theorem states that
there must exist a time-varying Lagrange multiplier (i.e. “co-state”) vector
! and a constant Lagrange multiplier vector ! such
that a Hamiltonian function, H, can be formed for which first order necessary conditions for
optimality are valid. For our pseudo-rendezvous problem, the Hamiltonian can be expressed as:
! (2.18)
where the time-varying Lagrange multiplier (co-state) vector, ! , and state dynamics,
! , are defined as:
! (2.19)
! (2.20)
Note that the position and velocity vectors, and mass scalar, each have associated co-states.
Hence the Hamiltonian becomes:
! (2.21)
The co-states are required to satisfy the following necessary conditions:
! (2.22)
! (2.23)
! (2.24)
J = max ϕ(t
i
0
, t
i
f
)
⃗ λ
⊤
(t) =
[
λ
1
,λ
2
, … λ
n]
⃗ ν
⊤
=
[
ν
1
,ν
2
, … ν
q
]
H =
⃗ λ
⊤
⃗ f
⃗ λ
⊤
(t)
⃗ f =
·
⃗ X
⃗ λ =
[
⃗ λ
r
⊤
⃗ λ
v
⊤
λ
m
]
⃗ f =
·
⃗ X =
⃗ v
⃗ g( ⃗ r )+
T
m
̂ u
−
T
c
H =
⃗ λ
r
⊤
⃗ v +
⃗ λ
v
⊤
[ ⃗ g( ⃗ r )+
T
m
̂ u]+λ
m
T
c
·
⃗ λ
r
=−(
∂H
∂ ⃗ r
)
⊤
=−G
⃗ λ
v
·
⃗ λ
v
=−(
∂H
∂ ⃗ v
)
⊤
=−
⃗ λ
r
·
λ
m
=−(
∂H
∂m
)
⊤
=−
T
m
2
⃗ λ
v
⊤
̂ u
!22
Equations (2.22—2.24) form the co-state equations of motion, or “co-state dynamics”, which are
propagated along with the state dynamics’ equations of motion. The necessary condition for
optimality of the controls, ! , that maximizes the Hamiltonian is the well known primer vector
unit direction
[44]
:
! (2.25)
The Hamiltonian can be rewritten in terms of the well known switching function, ! (t):
! (2.26)
with
! (2.27)
and the thrust throttling function T determined from the following bang-bang control law:
! (2.28)
The augmented state vector to be numerically integrated is then:
! (2.29)
governed by the corresponding first order vector differential equation:
! (2.30)
Note that the explicit thrust direction ! has been successfully removed from the problem, and
replaced with the primer vector, governed by the co-state dynamics.
̂ u(t)
̂ u(t) = +
⃗ λ
v
(t)
λ
v
(t)
$
H =
⃗ λ
r
⊤
⃗ v +
⃗ λ
v
⊤
⃗ g( ⃗ r )+$(t)T
$(t) =
λ
v
m
−
λ
m
c
T =
0 if $ < 0
T
max
if $ > 0
0≤ T≤ T
max
if $ = 0 (i.e. singular arc, undetermined)
⃗ X =
[
⃗ r
⊤
⃗ v
⊤
m
⃗ λ
r
⊤
⃗ λ
v
⊤
λ
m
]
⊤
14x1
·
⃗ X =
⃗ v
⃗ g( ⃗ r )+
T
m
⃗ λ
v
(t)
λ
v
(t)
−
T
c
−G
⃗ λ
v
−
⃗ λ
r
−
T
m
2
λ
v
14x1
̂ u(t)
!23
The last step in defining the optimal control problem is to derive the transversality
(boundary) conditions. The initial kinematic boundary conditions at ! are
! (2.31)
where ! . The terminal kinematic boundary conditions at ! are:
! (2.32) The scalar Bolza function can be written as:
! (2.33) ! (2.34)
The required necessary conditions for the initial co-state vector and the Hamiltonian at ! are:
! (2.35)
! (2.36)
which reduces to:
! (2.37)
t
i
0
⃗ ϑ =
⃗ r
i
0
− ⃗ r(t
i
0
)
⃗ v
i
0
− ⃗ v(t
i
0
)
m
i
0
−m(t
i
0
)
7x1
=
⃗ 0
m(t
i
0
) = m
0
t
i
f
⃗ ψ =
⃗ r
i
f
−
⃗ r*(t
i
f
)
⃗ v
i
f
−
⃗ v*(t
i
f
)
6x1
=
⃗ 0
B = k(t
i
0
−t
i
f
)+
⃗ ξ
⊤ ⃗ ϑ + ⃗ ν
⊤
⃗ ψ
B = k(t
i
0
−t
i
f
)+
[
⃗ ξ
r
⊤
⃗ ξ
v
⊤
ξ
m
]
⃗ r
i
0
− ⃗ r(t
i
0
)
⃗ v
i
0
− ⃗ v(t
i
0
)
m
i
0
−m(t
i
0
)
+
[
⃗ ν
r
⊤
⃗ ν
v
⊤
]
⃗ r
i
f
⊤
− ⃗ r(t
i
f
)
⃗ v
i
f
⊤
− ⃗ v(t
i
f
)
t
i
0
⃗ λ
i
0
=−
∂B
∂
⃗ X
0
=
−
⃗ ξ
r
−
⃗ ξ
v
−ξ
m
7x1
H*
i
0
= +
∂B
∂t
i
0
= k−
⃗ ξ
r
⊤
⃗ v(t
i
0
)+
⃗ ξ
v
⊤
⃗ g(t
i
0
)
H*
i
0
= k +
⃗ λ
r
⊤
i
0
⃗ v(t
i
0
)+
⃗ λ
v
⊤
i
0
⃗ g(t
i
0
)
!24
Similarly, the necessary conditions for the final co-state vector and Hamiltonian at ! are:
! (2.38)
! (2.39)
which reduces to:
! (2.40)
It is now apparent that the only co-state known explicitly at the terminal state is ! .
Since ! , the necessary target values for H at ! and ! are ! and ! ,
respectively. Recall that H at ! is
!
with Switching Function
!
Because ! and ! , then:
! (2.41)
! (2.42)
leading to the conditions:
! and ! (2.43)
t
i
f
⃗ λ
i
f
= +
∂B
∂
⃗ X
f
=
− ⃗ ν
r
− ⃗ ν
v
0
7x1
H*
i
f
=−
∂B
∂t
i
f
= k + ⃗ ν
r
⊤
⃗ v*(t
i
f
)+ ⃗ ν
v
⊤
⃗ g*(t
i
f
)
H*
i
f
= k +
⃗ λ
r
⊤
i
f
⃗ v*(t
i
f
)+
⃗ λ
v
⊤
i
f
⃗ g*(t
i
f
)
λ
m
i
f
= 0
⃗ ϑ = ⃗ ψ =
⃗ 0 t
i
0
t
i
f
H*
i
0
H*
i
f
t
i
f
H
i
f
= k +
⃗ λ
r
⊤
⃗ v(t
i
f
)+
⃗ λ
v
⊤
⃗ g(t
i
f
)+$
i
f
(t)T
i
f
$
i
f
=
λ
v
i
f
m
i
f
−
λ
m
i
f
c
λ
v
i
f
> 0 λ
m
i
f
= 0
$
i
f
=
λ
v
i
f
m
i
f
> 0
T
i
f
= T
max
H*
i
f
= H
i
f
H*
i
0
= H
i
0
!25
Through continued reduction and expansion, it can be shown that ! will always be positive
and ! in the interval ! .
Finally, the conditions
! and ! (2.44)
require that
! (2.45)
! (2.46)
The described TPBVP can be solved either via forward time propagation ( ! ) or
backward time propagation ( ! ). A forward time propagation will require specifying
! and solving for ! to satisfy ! , and ! . To solve the TPBVP
for the forward propagation (forward shooting) case, the optimization variable vector becomes:
! (2.47)
with the inequality constraint ! , and equality constraint vector:
! (2.48)
which can be used with an NLP solver to determine the time-minimum optimal finite burn
trajectory.
$(t)
T = T
max
t
i
0
≤ t ≤ t
i
f
$
i
0
T
max
= k $
i
f
T
max
= k
$
i
0
= $
i
f
λ
v
i
0
m
i
0
−
λ
m
i
0
c
=
λ
v
i
f
m
i
f
t
i
0
→ t
i
f
t
i
0
← t
i
f
m
i
0
= m
0
, λ
v
i
0
= 1.0, λ
m
i
0
λ
m
i
f
= 0 $
i
0
= $
i
f
⃗ X
opt
=
[
t
i
0
t
i
f
α
i
0
·
α
i
0
β
i
0
·
β
i
0
·
λ
v
i
0
λ
m
i
0
]
⊤
8x1
λ
v
i
f
> 0
⃗ c(
⃗ X
opt
) =
⃗ r
i
f
−
⃗ r*(t
i
f
)
⃗ v
i
f
−
⃗ v*(t
i
f
)
λ
m
f
$
i
0
−$
i
f
8x1
=
⃗ 0
!26
2.4 The Multi-Impulse Finite Burn Conversion Problem
The main thrust of the presented research is an extension of the current cutting-edge
method for single impulse-to-finite burn conversion by Jesick & Ocampo
[51]
. As described in
Sections 2.2 and 2.3, a single impulsive ∆V can be converted to an time-optimal, continuous
finite-burn arc within the bounds of ! and ! . Now we wish to extend the single impulse case to
the multiple-impulse case by chaining together single impulse conversions, as illustrated in
Figure 3.
Let there be a set of n-number of ∆V impulses separated by ballistic arcs. We can say ∆Vn
occurs at ! at the state! . Though this trajectory is more complex than the single
impulse case, we can treat each impulse as independent from the rest of the sequence. This is
t
0
t
f
t
n
[
⃗ r
n
⊤
⃗ v
n
⊤
]
⊤
!27
Figure 3: Multi-Impulse to multi-finite burn conversion. The ballistic arcs between finite-burn arcs match
time, position, and velocity of the original impulsive trajectories.
because, as described, the single impulse conversion problem assumes bounded ballistic arcs on
either side of the finite-burn ! , and ! . In others words, there exists some
coast time before the start of the finite burn:
!
and some coast time after the finite burn:
!
If we consider the multiple-impulse case to be composed of a chain of single impulse
optimization problems (still minimizing finite-burn time), we can set the boundary constraints on
each impulse conversion to allow for the longest possible coast-time around either end of the
finite burns. This is the same as saying that for some impulse ∆Vk at tk within the multi-impulse
chain, the time bounds (t0, tf) become:
!
!
with corresponding states ! , ! , and ! , ! . This essentially constrains each
individual impulse conversion to have finite burn arcs which occur between subsequent
impulses. For a flyable multi-impulse converted trajectory to exist, the coast times between
finite-burns must be ! , where:
!
However, in general, it is possible within this formulation that overlaps in time-state
space may occur between subsequent finite burn arcs. In the event this occurs, a finite-burn arc
overlap between finite-burn k and k+1 means that the coast time between the two finite-burn arcs
is negative, i.e.:
!
(t
0
, ⃗ r
0
, ⃗ v
0
) (t
f
, ⃗ r
f
, ⃗ v
f
)
Δt
coast
0
= t
0
−t
i
0
≥ 0
Δt
coast
f
= t
f
−t
i
f
≥ 0
t
0
(k)
= t
f
(k−1)
t
f
(k)
= t
0
(k+1)
⃗ r(t
0
(k)
) ⃗ v(t
0
(k)
) ⃗ r(t
f
(k)
) ⃗ v(t
f
(k)
)
Δt
c
(k→k+1)
≥ 0
Δt
c
(k→k+1)
= t
(k+1)
0
−t
(k)
f
t
(k+1)
0
−t
(k)
0
=Δt
c
(k→k+1)
< 0
!28
This is because each impulse in the multi-impulse chain is treated as an independent conversion
problem, only bounded on the ballistic arcs in time and state by neighboring impulses. This
crucial issue within the multiple-impulse conversion problem is illustrated in Figure 4. The
solution to this issue is illustrated in Section 7. !29
∆t
c
1-2
< 0
Figure 4: Multi-Impulse to multi-finite burn conversion, with finite-burn arc overlap between burns 1 & 2.
3. PRIOR WORK
3.1 Survey of Existing Optimization Methods
An overview of Indirect & Direct Methods for solving continuous low-thrust orbit
transfers, and more generally finite-burn problems, is illustrated in Table 1. A key selection of
relevant optimization methods and tools that are representative of published state-of-the-art work
is discussed in greater detail in the following subsections. These methods provide context for the
presented conversion algorithm.
3.1.1 MALTO (NASA-JPL)
The Mission Analysis Low Thrust Optimization (MALTO) software package, developed
by Sims et al.
[17]
at NASA’s Jet Propulsion Laboratory (JPL), is a robust medium-fidelity tool
used for preliminary design of low-thrust trajectories. This direct method discretizes continuous
two-body trajectories into multiple-impulse segments, with flybys of secondary bodies modeled
as instantaneous changes in the direction of the hyperbolic excess velocity, V∞, vectors.
The underlying trajectory structure used in MALTO is shown in Figure 5. MALTO
employs a forward/backward multiple shooting technique with intermediate control nodes, which
can represent planets, small bodies, or free points in space. Each leg begins and ends at a control
node, and has a single match point where the two halves of a segment are patched together. The
trajectory is propagated forward in time from the leg’s previous control node to the match point,
and backward in time from the leg’s next control node to the match point. The match points are
constrained so that through subsequent numerical optimization, the discontinuities in state at
these points are eliminated. MALTO’s developers found that this technique significantly reduces
the propagated trajectory’s sensitivity to intermediate flybys compared to strict forward
propagation
[17]
. Each control node has associated optimization variables of the spacecraft’s
velocity relative to the node body, the spacecraft’s mass, and corresponding epoch.
For optimization, MALTO uses the constrained, NLP software tool Sparse Nonlinear
Optimizer (SNOPT) developed by Gill et al.
[26]
at the University of California San Diego and
!30
Stanford University. The partial derivatives of the constraint and cost functions are computed
analytically, i.e. not by finite differencing. The developers note that this represents the bulk of the
method’s computational effort and code, as well as a key component leading to fast execution
and robust convergence
[17]
. The locally optimal solutions found by MALTO depend on the initial
guess values of the node variables. Initial guess trajectories are generated using a fast-converging
patched conic trajectory optimization program. Vector theory or a spiral guess control law are
employed that assume thrust directions along or against the velocity vector on each segment. The
spiral guess was found to be useful for generating feasible trajectories with multiple revolutions
about the Sun
[17]
. MALTO has been used in the development of mission proposals and advanced
concept studies, most notably the Jupiter Icy Moons Orbiter (JIMO) mission
[27-28]
.
3.1.2 Mystic (NASA-JPL)
The Mystic software, developed by G. Whiffen
[12]
at JPL, is a high-fidelity direct
optimization method designed to compute, analyze, and visualize optimal low-thrust trajectories.
!31
Figure 5: The two-body multi-impulse trajectory structure used in MALTO
[17]
.
Mystic employs a non-linear optimization method that can simultaneously optimize both “static
variables” (constant parameters, e.g. time of flight) and “dynamic variables” (functions of time,
e.g. thrust vectors). The tool seeks to maximize final spacecraft mass or minimize a user defined
infeasibility, i.e. the magnitude of constraint violation. Trajectories are fully integrated using
finite-burns, multi-body propagation, gravitational harmonics, and solar radiation pressure. The
trajectories computed by Mystic are continuous from beginning to end, rendering solutions
closer to the true LEP dynamics than multiple-impulse methods. Mystic is best suited for multi-
body or non-spherical body gravity problems, and automatically seeks out and exploits gravity
assists during optimization
[2]
.
Mystic uses three classes of optimization variables: state variables, dynamic control
variables, and static control variables. State variables include 3D position, velocity, and
spacecraft mass as functions of time. Dynamic control variables include the 3D thrust vectors,
and the static control variables include the transfer start time, total flight time or arrival date, and
initial state. Mystic’s optimization constraint settings allow for general orbital element targeting,
as well as geometry constraints on the trajectory, propellant mass, and thrust beam cone angle.
Trajectory geometry constraints include minimum fly-by distances, forced fly-bys, launch
asymptote & energy, and intermediate targeting and rendezvous constraints. A separate tool
generates a sub-optimal initial guess trajectory comprised of patched two-body problems using
simple feedback laws or Q-Law guidance
[12]
.
Mystic has been successfully used to design NASA’s JIMO mission’s reference trajectory

[29-30]
. The limitations of the types of continuous trajectories Mystic can optimize correspond
directly to the stability and accuracy limitations for multi-body propagation. There is a limit of
approximately three close fly-bys of massive interplanetary bodies, or six close flybys of massive
moons around a single planet, within a single trajectory. While Mystic’s continuous trajectory
modeling is an attractive feature, initial guesses can take time to set up. Spiraling problems are
limited to 250 orbit revolutions for practical run times, however sub-optimal trajectories with
300 or more revolutions are possible using sophisticated control laws
[12]
.
!32
3.1.3 STOUR-LTGA with GALLOP & SEPTOP (Purdue)
The Satellite Tour Design Program for Low-Thrust Gravity Gravity Assists (STOUR-
LTGA) is a program developed by Petropoulos and Longuski
[6]
from Purdue University. The
analytic shape-based method employs exponential-sinusoid thrust arcs and automatically
searches for LTGA trajectories. However, a sequence of gravity-assisting bodies, launch dates,
and launch V∞’s must be provided
[6]
. Two-body motion is assumed between fly-bys, with
gravity-assists modeled as discontinuities in velocity from instantaneous turning of the
spacecraft’s V∞ vector with respect to the assisting body. The motions of the assisting bodies are
modeled by polynomial representations or higher fidelity ephemeris data. Coast arcs are
permitted between constant specific impulse thrusting arcs on the transfers between bodies. This
approach allows for rapid & broad searches over the design space of trajectories, whose solutions
can be used as initial guesses in optimization routines.
The exponential-sinusoid shape used to represent the spacecraft's trajectory is given in
2D polar coordinates of position and true anomaly (r, θ):
! (3.1)
where k0, k1, k2, and are constants. The dynamic range parameter, k1, controls the ratio of the
apoapsis distance to the periapsis distance. Parameter k2 is the called the winding parameter, k0 is
a scaling factor, and the phase angle ! controls the in-plane orientation of the exponential-
sinusoid (Figure 6). The tuning of these parameters, namely k1 k2, is done using a search-
r = k
0
e
k
1
sin(k
2
θ+ϕ)
ϕ
ϕ
!33
Figure 6: Example exponential-sinusoid shapes, from periapsis to apoapsis, for k1 = 0.5
[6]
.
direction-controlled Newton’s method, while the flyby geometries (i.e. miss angles) are
computed by quadrature. This allows intersections between the exponential-sinusoid arc in the
ecliptic and planar projections of the assisting bodies to be found (Figure 7). Such intersections,
computed via root-finding, can be used to identify inbound and outbound gravitational body
encounters. The out-of-plane motion is determined based on whether the spacecraft encounters a
target body on a thrust or coast arc. The target’s out-of-plane position at the time of the in-plane
encounter is matched by using additional thrust acceleration acting along or against the
spacecraft’s angular momentum vector for some final portion of that particular thrust arc
[6]
.
STOUR-LTGA trajectories have been used to seed Purdue’s in-house Gravity-Assist,
Low-Thrust, Local Optimization Program (GALLOP)
[18]
. GALLOP uses the same multi-impulse
discretization optimization method as in MALTO. STOUR-LTGA trajectories have also been
validated using the indirect optimization tool Solar Electric Propulsion Trajectory Optimization
Program (SEPTOP)
[2]
. The initial guesses for the co-state values used in SEPTOP were found by
trial and error, as neither STOUR-LTGA nor GALLOP are able to estimate initial co-state values.
!34
Figure 7: Intersections of an exponential-sinusoid arc in the ecliptic with the projection of Ceres’s orbit (Points B, C, G are outbound, D is inbound)
[6]
.
3.1.4 Exploiting Invariant Manifolds via the Circularly Restricted Three-Body Problem
An alternate means of computing feasible initial guess trajectories is to exploit the natural
dynamics of the planar Circularly Restricted Three-Body Problem (CRTBP)
[19, 20, 21]
. Figure 8
depicts the classic 3-body problem as well as the 4-body problem. The CRTBP model provides a
framework for well understood low energy formulations such as invariant manifolds and periodic
orbits to help guide low-thrust trajectory design. Anderson and Lo
[19]
explored the role invariant
manifolds play in the dynamics of low-thrust trajectories moving through unstable regions in
three-body systems. They showed that an optimization algorithm incorporating no knowledge of
invariant manifolds converges on low-thrust trajectories that use the invariant manifolds of
unstable resonant orbits to traverse resonances. Their results indicate that resonant orbits and
their invariant manifolds may be fundamentally important for understanding the optimal
pathways that may exploited by low-thrust trajectories.
Similarly, Mingotti et. al.
[20]
investigated using patched-planar three- and four-body
systems to design Earth-to-Mars transfers. Low-thrust propulsion was incorporated into standard
invariant manifold techniques by defining special attainable sets to fill in the gaps between two
patched-CRTBP systems, i.e. connect nonintersecting manifolds with low-thrust propulsion. An
application of their method is comprised of three phases: ballistic escape, low-thrust rendezvous,
and low-thrust capture. The models include a Moon-perturbed Sun-Earth problem (CR4BP),
!35
Figure 8: Three- (left) and four- (right) body problems. For example, an Earth-Moon problem and a Moon-
perturbed Sun-Earth problem, respectively
[20]
.
used for ballistic escape via lunar gravity assist, and a Mars-Sun CRTBP used for low-thrust
capture (Figure 9). The patched-CRNBP trajectory was then optimized numerically using a
direct approach, developed by Hargraves and Paris
[22]
, using the full N-body problem. While
both Mongotti & Anderson’s approaches demonstrate significant advantages to exploiting the
CRTBP and low energy formulations in feasible trajectory generation, a non-trivial effort is
required by experienced analysts to navigate the design process to find suitable manifolds.
3.1.5 Copernicus (NASA-JSC)
Developed at NASA’s Johnson Space Center (JSC) primarily for use in trajectory design
for Human exploration missions, Copernicus is a general trajectory design and optimization
system that can be used to solve a wide range of spacecraft trajectory problems in a robust and
efficient manner
[45]
. The software tool can be applied to the following three main problems:
!36
Figure 9: Sample low-thrust capture trajectory at Mars via an invariant manifold in the Sun-Mars CRTBP
rotating frame
[20]
.
Y
X
1. Modeling and Open Loop Simulation: Copernicus can simulate open loop
trajectories in any forcefield (N-body, spherical gravity, etc.), inclusive of multiple
spacecraft, and impulsive and/or finite-burn maneuvers.
2. General Targeting via Non-linear Root Finding: Copernicus can search for the
values of the independent variables required to satisfy a set of constraint functions, or,
more generally, formulate and solve general systems of nonlinear equations that
depend on user-defined independent variables.
3. Optimization: Copernicus can extremize the value of any computable function,
typically for maneuver design and optimization, which requires determining the
maneuver parameters that describe the burns (impulsive or finite-burn).
The trajectory model employed by Copernicus uses a fundamental building block called a
segment, as shown in Figure 10. A segment is a trajectory arc that can have velocity impulses
and/or finite-burn maneuvers. Each segment is defined by two terminal node points, which are
tagged with epochs relative to a reference epoch. The segment duration is the difference in time
between the initial and final nodes, which is allowed to be negative or positive. Specific
!37
Figure 10: Possible maneuver combinations for a segment in Copernicus
[63]
.
parameters within the states at the initial and final nodes can either be user defined, or chosen to
be optimization variables to be determined by the internal solver. Both the orbital state (position
and velocity) and spacecraft vehicle state (mass and power) can be explicitly modeled. The
segment’s trajectory is propagated real-time between the nodes, or for a combination of an initial
or final time and a user given time of flight (dt), which can also be positive or negative.
Chaining together multiple segments constructs a complete trajectory for one or more
spacecraft, an example of which is illustrated in Figure 11. Discontinuities in velocity and mass
between subsequent nodes can represent impulsive burns or module separation. To create a
continuous trajectory across multiple segments, terminal continuity constraints may be assigned
to the node’s orbital and spacecraft states. This will patch two segments together, matching the
constrained states on the final and initial nodes, respectively. Forward and backward propagation
is a commonly employed tactic when designing missions in Copernicus, as it makes the
convergence behavior of the optimizer more favorable. Often one segment will be assigned a
positive dt (propagated forward in time), with the next chronological segment defined with a
negative dt (propagated backwards in time), with continuity constraints on the final and initial
nodes to match states (orbital position, velocity, and spacecraft mass).
Copernicus allows for manually crafted initial guess trajectories via its intuitive graphical
user interface (or autonomously via scripting through Python plug-ins). As such, more consistent
!38
Figure 11: An example trajectory in Copernicus, comprised of multiple segments, maneuvers, and plug-ins
[63]
.
and expedited convergence times have been found by placing match points between segments in
relative “deep-space”, roughly at the midpoint of a coast arc. This is due to the sensitivity that
small changes in optimization variables have on the propagated state when closer to a gravitating
body. For instance, a small change in final node orbital state at periapsis on a planetary flyby will
have a dramatic affect on the initial node orbital state when propagated backwards in time.
Finite-burns (including low-thrust) are modeled during the propagation between nodes on
a segment. In Copernicus, finite-burn arcs are modeled in two ways. The first is a parameter
model that uses a finite set of parameters to describe the time evolution of the control variables.
The thrust direction unit vector is commonly parameterized by two spherical angles, which can
be fixed to any desired frame (typically the spacecraft body frame, an osculating frame, inertial
frame, etc.). Copernicus is capable of modeling these angles with constant, linear, quadratic, and
sinusoidal terms as functions of time. In most applications, analysts use only the constant and
linear terms. However in longer finite-burns, such as those associated with low-thrust problems,
both the quadratic and sinusoidal terms can be useful.
The second method used to model finite-burns is via optimal control. As described in
Section 2, this involves augmenting the state vector to include the co-state vector of variable
Lagrange multipliers adjoined to the state, and iterating on the initial co-state vector. In all
maneuver problems (impulsive and finite-burn), the optimization architecture is represented in
standard form currently used by state of the art nonlinear constrained optimization codes. As
such, Copernicus is designed to be modular, and can employ a variety of NLP solvers, most
notably SNOPT. Equality and inequality constraints can be defined for every segment node and
design parameter, but only thorough experience reveals how best to set these constraints in order
to not over constrain the problem, which can lead to poor or failed convergence.
!39
3.1.6 Edelbaum’s Low-Thrust Transfer Solution
Edelbaum derived an analytic method for computing continuous low-thrust transfers
between non-coplanar circular orbits in the 2-body problem
[35]
. It should be noted that this
solution is only valid for relatively low values of constant thrust-acceleration, such that the
spacecraft is considered to always be in a circular orbit. This means that the Edelbaum solution is
not applicable to high-thrust, finite-burn problems, but is useful for modeling low-thrust
transfers, namely 3D spiral orbits with inclination change (but not eccentricity change, due to the
circular assumption). For this reason it is an excellent choice for producing initial guess solutions
for many continuous, low-thrust problems. It should be noted that the algorithm is only valid for
inclination changes of less than 114.6˚
For numerical solutions based on Edelbaum’s method
[36-37]
, the initial thrust vector yaw
angle is computed by:
! (3.2) where: ! (3.3)
The total velocity change required for a low-thrust transfer is found from:
! (3.4)
With the total transfer time described by:
! (3.5)
where at is the magnitude of the constant thrust-acceleration vector. The yaw angle as a function
of time is described by:
! (3.6)
tan(β
0
) =
sin(
π
2
Δi)
V
0
V
f
−cos(
π
2
Δi)
V
0
=
μ
r
0
V
f
=
μ
r
f
ΔV = V
2
0
−2V
0
V
f
cos(
π
2
Δi)+V
2
f
t
f
=
ΔV
a
t
β(t) = tan
−1
(
V
0
sin(β
0
)
V
0
cos(β
0
)−a
t
t
)
!40
The change in inclination with respect to time is:
! (3.7)
The instantaneous velocity is given by:
! (3.8)
Using the orbital-energy invariance law, the semi-major axis with respect to time is:
! (3.9)
where µ is the standard gravitational parameter of the central body (i.e. Earth, Moon, Mars, etc.).
The thrust-acceleration unit vector is calculated from:
! (3.10)
where [r, t, n] is the spacecraft-centered radial-tangential-normal coordinate frame defined by:
! (3.11)
The out-of-plane pitch angle with respect to time, ! , and in-plane yaw angle, ! , can then
computed by:
! (3.12)
! (3.13)
An example LEO-to-GEO transfer computed using this method is shown in Figure 12.
The example spacecraft has a constant thrust-acceleration of 150 µg, an initial circular orbit
altitude of 350 km inclined at 28.5° (corresponding to a launch from Cape Canaveral), and a final
circular equatorial altitude of 35,786 km (GEO). The Edelbaum solution finds the total ∆V to be
5.924 km/s, with a total transfer time of 46.67 days.
Δi(t) =
2
π
[tan
−1
(
a
t
t−V
0
cos(β
0
)
V
0
cos(β
0
)
)+
π
2
−β
0
]
V(t) = V
0
sin(β
0
)
1+tan
2
(β(t))
tan(β(t))
a(t) =
μ
V(t)
2
̂ u
T
(t) =
sin(β(t))
cos(β(t))cos(ψ(t))
cos(β(t))sin(ψ(t))
=
U
T
r
U
T
t
U
T
n
rtn
R
eci
rtn
=
[
⃗ r
| ⃗ r |
( ⃗ r × ⃗ v )× ⃗ r
| ⃗ r × ⃗ v || ⃗ r |
⃗ r × ⃗ v
| ⃗ r × ⃗ v |]
ψ(t) β(t)
ψ(t) = tan
−1
(U
T
n
,U
T
t
)
β(t) = sin
−1
(U
T
r
)
!41
3.2 Single Impulse-to-Finite Burn Conversion
One of the earliest investigations into converting an impulsive trajectory into an optimal
control finite burn trajectory was made by Pines in 1964
[47]
. In that study, four constants of
motion were derived for optimal thrust trajectories in a central inverse squared force field. These
constants were applied to the impulsive thrust case in order to estimate the minimum propellant
optimal finite burn thrust arc via an indirect method. It was already known that solving the finite
burn TPBVP via shooting method is made difficult by having to guess an initial set of co-states,
and rely on an iterative procedure to converge on a solution that satisfies the defined boundary
conditions. The issue for Pines, as it still is for us today, was the sensitivity of the co-state
dynamics to initial guess values, making convergence difficult-to-near-impossible to obtain if
insufficient guesses were provided. To avoid this difficulty, Pines proposed iteration starting
from ! from an already determined optimal impulsive trajectory, similar in ∆V to that of the
desired finite-burn. In other words, to use the impulsive adjoint vector as an initial guess
[49]
.
However, numerical validation of this approach was not provided by Pines
[48]
.
⃗ λ
i
(0)
!42
Figure 12: Example Edelbaum LEO-to-GEO continuous low-thrust transfer.
A few years later, Handelsman
[48]
sought to validate Pines’ conjecture by using a sweep
method to approximate the optimal control solution of a two-dimensional, two maneuver
trajectory with an initial impulsive solution. Handelsman calculated the set of ! from the
orbital elements of the ballistic 2-body transfer between two terminal states (i.e. an elliptical arc).
The initial thrust-acceleration was assumed to be a relatively large value, which initiates the
iteration with a non-impulsive case where the co-states differ slightly from the impulsive case.
This ensures a converged trajectory. The solution is then iterated upon with decreasing steps of
thrust-acceleration (seeded by the previously converged iterate) until the desired value of thrust-
acceleration is obtained. Handelsman’s approach was shown to be successful in solving for 2-
dimensional, 2-body trajectories with “bang-bang” optimal steering control.
Next in 1970, Hazelrigg and Lion
[49]
developed a method for generating more accurate
estimates of the initial adjoint vector for use in the widely used iterative shooting solution
strategy. Building upon the work done by both Pines and Handelsman, their approach was based
on the idea that the impulsive trajectory represents the limit of the finite-burn trajectory as the
thrust level is increased to infinity. Instead of just using the impulsive adjoint vector the initial
guess for the finite-burn adjoint vector, a solution for the finite-burn adjoint vector is sought in
the form of a series expansion in 1/a, where a is the maximum thrust-acceleration
[49]
. This
technique was shown to be most accurate for high thrust problems, and through numerical
experiments Hazelrigg et. al. showed their approximation can work for low-thrust ranges as well.
Soon thereafter, Kornhauser, Lion, and Hazelrigg
[50]
documented a comprehensive study
of a minimum propellant transfer using two and three finite-burn arcs for an idealized Earth-
Mars transfer. Their solution method is obtained in series form as an expansion about the
optimal, unbounded thrust impulsive trajectory. This expansion is done in terms of two
independent spacecraft parameters: thrust-acceleration, and the exhaust velocity. Since this
expansion is once again about the impulsive case, this technique is also most accurate for high
thrust problems.
These early studies are indeed trailblazing in the topic of impulse-to-finite burn
conversion, but their work does not address the core questions at the center of the presented
⃗ λ
i
(0)
!43
research. The main distinction is that here, the formulation and algorithms are extended to the
multi-impulse conversion problem, in a general position dependent gravity field. The presented
research is an extension of this cutting-edge work, taking direct inspiration from the efforts of
Jesick and Ocampo
[51]
. They present an automated procedure for converting a single impulsive
maneuver into a minimum fuel finite-burn trajectory.
The process involves first converting the impulsive maneuver into sub-optimal finite-
burn solution using parameter optimization, and then using this solution to seed a solution to the
full optimal control problem. Their hybrid optimization method uses a parameter optimization
algorithm to directly optimize a performance index, namely minimize propellant, and then the
control is found using optimal control theory. An advantage they cite of this approach is that the
state and control constraints can be easily applied in a parameter optimization framework
[51]
.
They employ an adjoint control transformation to determine the optimal thrust vectors along the
finite-burn arc (see Section 2 for more details). The general algorithm of this hybrid method is as
follows:
1. Solve the impulsive optimization problem. 1a. Convert each impulse into a linearly steered finite burn arc.
2. Solve linearly steered finite burn optimization problem. 2a. Estimate co-states on each thrust arc.
3. Solve optimal control finite burn optimization problem. !44
4. EXISTING GAPS
4.1 General Gaps in Finite-Burn & Low-Thrust Optimization
As detailed in Section 3, the topic of low-thrust trajectory optimization, and in general
finite-burn optimization, has been studied for decades. While the underlying optimal control
theory has be known for far longer, due to the highly non-convex, nonlinear nature of the finite-
burn optimal control problem
[40]
, its methodology has not been widely employed in modern
applications. Beyond the sheer numerical difficulties in converging on indirect solutions, the
industry has in recent years shifted towards direct optimization methods, largely due to the
overwhelming computational power that has become available in the digital age. Though direct
methods do not guarantee globally optimal solutions, the ease and speed to which direct
solutions may be found using modern processors have made the search for complete indirect
solutions no longer a priority for most industry players. This is because accuracies attainable
using direct methods have proven sufficient for most mission design applications. While many
highly lauded optimization tools such as MALTO, MYSTIC, and Copernicus rely on direct
parameter optimization, many aspects of the finite-burn optimal control problem still remain
unsolved or difficult to implement.
In summation of the reviewed literature, the following represent 8 key knowledge gaps,
theoretical deficiencies, or difficulties in applying existing methods and theory used to solve
low-thrust, and in general finite-burn, trajectory optimization problems. Seven general gaps are
presented in this section, with the eighth (representing the chosen topic for this research) detailed
in Section 4.2. These knowledge gaps represent distinct areas of research within the field, and
could each individually serve as the backbone of a unique PhD dissertation.
Gap 1: There are no generalized methods for solving finite-burn optimal control
problems that are easily adaptable to problems with different state dynamics. Existing methods for solving finite-burn optimal control problems from scratch are
difficult to adapt to similar problems with different force models
[2, 4, 5, 39, 42-45, 47-51]
.
!45
While the underlying mathematical theory used to describe the optimal control
TPBVB problem is well known, in orbit trajectory applications the TPBVP must be
re-derived for each unique dynamical context. This requires the re-calculation of the
partial derivatives of the objective and constraints, which can be a formidable
challenge with complex state dynamics. It is not uncommon to see orbit transfer
problems that involve anywhere from 2- to N-gravitating bodies, that may or may
not be represented by central-body force fields, spherical harmonics, or other
complex models. Additional perturbations can include extra gravitational effects
(e.g. J2), aerodynamic drag, Solar radiation pressure, collisions with micro-meteors
and orbital debris, spacecraft outgassing and leaks, among others. Within the
spacecraft model itself, specific vehicle parameters may be desired to be part of the
optimization objective function (e.g. power limited, propulsion system thermal
constraints, vehicle slew time to orient the engines, etc.). For each added
acceleration term, constraint, and objective, the optimal control problem must be
reformulated and the partial derivatives of the constraints and objective found by
hand or symbolic computation, since the co-state dynamics are intimately linked to
the unique state dynamics. This is in contrast to direct methods, where the controls
(e.g. thrust vectors) are directly used as optimization variables, with derivatives
estimated internally by standard optimization algorithms; each perturbation,
objective, or constraint can then easily be added to the corresponding function,
without explicitly computing any gradients. This is why direct methods have been
historically more adaptable to different orbit problems than indirect methods. Gap 2: There is no straightforward method for confidently determining the
appropriate initial co-state values to ensure convergence for general problems. In existing indirect methods, determining the appropriate initial values for the co-
states is difficult, leading to indirect methods not being widely used in existing tools
[1-2, 7, 10-16, 23, 39-40, 42-45, 47-51, 57-59]
.
!46
It is well known that in the finite-burn optimal control problem, convergence on an
optimal solution is extremely sensitive to the initial values of the co-state vector.
This difficulty is amplified by the fact that estimating the initial co-states can be
challenging, as they do not hold any intuitive, physical representation. The common
approach used to estimate the initial co-state vector is to employ an adjoint control
transformation by solving a simpler problem with either analytic or numerical results
which can then be used to approximate the optimal initial co-states in the full
problem
[58]
. This idea was first introduced by Dixon and Biggs
[59]
, who suggested
estimating physical control variables and their derivatives instead of the initial co-
states in order to reduce problem sensitivity and bake in more physical significance
to the co-states. Multiple mappings can connect the initial co-states to the initial
evolution of the thrust vector or initial parameter solution, and a clever selection of
reference frame can also reduce the problem’s sensitivity to the initial co-state
vector. While these kinds of methods have been proven to successfully estimate the
initial co-states, such transformations must be considered and derived for each
unique dynamical context. For instance, a unique transformation for the problem of
continuous-thrust, coplanar, circular to circular orbit transfer, compared with a 3-
body, Lagrange point or halo orbit transfer problem. Gap 3: It is difficult to determine optimal trajectories for low-thrust problems with
large inclination change. Existing low-thrust methods assume planar trajectories
[5, 16, 18]
, or have limited
ability to optimize large changes in inclination
[4, 16, 38, 18]
.
There are many classes of finite-burn transfer problems that can be realistically
modeled using planar, or near-planar assumptions. However for high-revolution,
low-thrust spiral orbits with high inclination change, the switching function within
the optimal control problem runs into problems due to the oscillatory change in
thrust direction as the spacecraft repeatedly passes through the ascending and
!47
descending nodes. Some methods instead employ shape-based estimations of the
spirals, such as Edelbaum’s method, which can accommodate large inclination
changes. However in the case of Edelbaum’s method, this leads to the implicit
constraint that during the finite-burn, the spacecraft is assumed to always be in an
instantaneous circular orbit. This essentially requires the thrust-acceleration to
remain well within the realm of ultra-low-thrust (proportional to the instantaneous
gravitational acceleration), and is not well suited for general finite-burn modeling.
3D thrust vector modeling using parameter sub-optimal control, such as in
Copernicus, easily gets around this, but is then an inherently locally-optimal, direct
optimization approach. Gap 4: How to handle bang-bang thrust within the switching function. The switching function in the finite-burn optimal control problem leads to bang-bang
thrust with the possibility for singular arcs, causing issues with optimizer and
integrator performance
[41-45, 51, 56-57]
.
When the switching function described in Section 2.3 is exactly zero (which can
occur during a sign change, signaling the thrust to turn on or off), the problem has
reached a non-smooth singular arc and the thrust is indeterminate. While in practice
this condition is rare in finite-burn trajectory problems, it can exist. Through each
shooting attempt, it is necessary for the optimizer to integrate the system of ODEs
defining the state and co-state dynamics, which due to the switching function is a
discontinuous function of time. As such, numerical integration schemes like Runge-
Kutta with adaptive step size develop issues when encountering the discontinuity
caused by a singular thrust arc. Additionally, minimization theorems for Newton’s
method require the equations to be twice continuously differentiable, and require
their Jacobian matrices to be non-singular in the vicinity of the solution. This means
Newton’s method often fails to converge in the case of non-smooth functions, as is
the case with the piece-wise switching function
[56]
.
!48
There are a variety of existing methods to handle bang-bang thrust, summarized
below. They each approach the singular arc problem from different perspectives, and
therefore leave much room for further study and investigation into novel methods. Use quadratic cost: This removes bang-bang thrust structure, making it easier to
handle numerically, but is not mass-optimal for the constant specific impulse case.
Quadratic cost is only mass-optimal for the variable specific impulse case. Pre-define thrust structure: Divide finite-burn arc into N sub-legs, alternating forced
thrust-coast-thrust arcs (still optimal per Potryagin principal, considering admissible
controls). The times where the switching function changes sign must be determined,
and targeted at junctions. The problem is then kept smooth, and variable step
integrators can be used on each sub-leg. The disadvantages of this method include:
thrust level is no longer a control, and the solution is no longer optimal for the
arbitrary thrust level case or those with thrust arcs in the middle of a coast arc (i.e. it
is only optimal for the fixed thrust-coast-thrust structure). Approximate the step function with Heaviside function approximations: Define the Heaviside function (H) approximation:
! (4.1) ! (4.2) ! (4.3) The approximate step function switch from ! to ! then becomes:
! (4.4)
where ! is the switching function defined in Section 2.3.
H(x) =
{
1 if x≥ 0
0 if x < 0
H = lim
k→∞
(H*)
H*(x) =
1
2
+
1
2
tanh(kx) =
1
1+e
−2kx
T
min
T
max
T = T
min
+H$(t)(T
max
−T
min
)
$(t)
!49
As the value of the constant k increases, the Heaviside function H* becomes closer
to a pure step function. The main issue associated with employing the Heaviside
approximation is that the transversality conditions and Hamiltonian of the optimal
control problem must be redefined so that this new control solves the problem.
Another issue is that gradient based optimizers run into problems with exponential
functions since they are so non-linear. The cliff on H*(x) can be very steep
depending on the choice(s) of k, which can cause high sensitivities in the indirect
problem where initial co-states control multiple cliffs. Homotopy / Continuation Method “Smoothing”: Iterate on the Heaviside function
approximation for higher values of k , using the previously converged solution as a
guess for the next sub-problem. This iteratively guides the optimizer closer to the
true bang-bang problem. Alternatively, the problem could be first solved using
quadratic cost (for q = 1 in the following objective function involving thrust-
acceleration Γ), then used to iteratively solve new problems that approach the bang-
bang problem (q = 0): ! (4.5) Employ an integrator with an event function to stop exactly at $(t) = 0: The simplest
approach, this allows for an arbitrary number of sub-legs without growing the
number of problem unknowns or constraints. This also avoids having to compute
analytic derivatives across boundaries. Gap 5: More robust and accurate methods for computing finite-burn initial guess
trajectories. Existing direct, parameter optimization methods commonly employ sub-optimal,
shape-based methods to approximate finite-burn thrust arcs
[2, 6, 18, 35, 57]
.
Shape-based and curve-fitting approaches improve the ease of initial guess
generation and encourage convergence behavior, but their solutions are inherently
sub-optimal (or more specifically, the solutions are local minima, with no guarantee
for global optimality). Further study may reveal more diverse and accurate methods
for approximating finite-burn arcs in varying force field models. J =
∫
Γ−qΓ(1−Γ)dt
!50
Gap 6: Low-thrust, multi-revolution spiral problems often require prohibitively long
simulation times and high computational power to solve. Optimization run-times of high revolution spiral orbits or multi-gravity assist low-
thrust transfers become computationally prohibitive, requiring sophisticated global
solvers or heavy simplifying assumptions
[6-8, 13, 31-34, 51-55]
.
Due to the computational intensity involved with simulating high revolution count
spiral trajectories, as well as the complexity of multiple-low-thrust gravity assist
(LTGA) problems, low-thrust optimization methods are often slow to converge
simply due to the large propagation times and difficulty in estimating the state and
co-state dynamics of the problems. Existing methods rely on one of two techniques
for expedited optimization. The first, employing sophisticated global solvers that
utilize genetic algorithms, evolutionary neurocontrollers, or robust parallel
computing processors to intelligently “brute-force” through the computationally
taxing problem. The second technique is to use shape-based approximations to
represent a portion or the totality of the spiral transfer. Either option represents a
different approach to grappling with the size of the spiral and multiple-LTGA
problems, and leaves room for further work to improve solution strategies for these
classes of problems. This could involve the development and utilization of highly
parallelizeable subroutines, such as the Picard-Chebyshev iteration method
[31-34]
, or
reformulating the optimization problem with an eye for efficiency. Gap 7: There is a prohibitively steep learning curve to implementing finite-burn
optimal control methods. Existing indirect methods require such complexity and rigor to derive and implement
in generic solvers, that they are often overlooked by industry orbit designers.
In practice, trajectory optimization has become “an art”, meaning the process of
designing a mission requires significant experience and intuition with each relevant
unique problem type and tool. Design and optimization efforts could be vastly
improved by developing solution methods that are more adaptable to different types
!51
of finite-burn problems and force models, which make use of more intuitive initial
starts. On top of this, many existing high-fidelity tools are either restricted use, or are
based upon underlying theory that is not well documented
[12, 39]
. In order for the
optimization power of indirect methods to be more widely taken advantage of, the
methodology could be reformulated to maximize Human-in-the-loop usability.
4.2 Specific Gaps within the Impulse-to-Finite Burn Conversion Problem
One final knowledge gap in finite-burn trajectory optimization is the lack of a robust
method for converting a multiple-impulse trajectory into a time-optimal, multi-finite-burn
trajectory, and addressing the potential finite-burn arc overlap problem. Such an impulsive
trajectory may represent an approximation of a finite-burn or low-thrust transfer, which can be
easier and faster to derive than computing a more accurate finite-burn initial guess.
Existing tools for modeling and optimizing finite-burn problems are powerful, but limited
in scope to either work with direct methods (e.g. Copernicus
[2, 45]
, MALTO/GALLOP
[6,17,18]
), or
indirect methods (e.g. Kechichian’s method
[4]
), with no clear conversion strategy from one to the
other. This is important because many tools, such as Copernicus, make it relatively easy to model
and optimize impulsive trajectories (due to the reduced number and sensitivity of optimization
variables), as well as sub-optimal controlled, parameter-steered finite-burns (due to the
optimizer’s direct control of the thrust steering law). This means that in real-world mission
design problems, these tools make the effort of constructing initial trajectories far easier and
more intuitive by staying in the direct optimization realm.
The issue arises once the initial guess has been optimized. Since the finite-burns have
been modeled by either multiple impulses or parameter steering laws, the converged solutions
are inherently sub-optimal. In many real-world design scenarios, the attained accuracy via direct
methods is sufficient, i.e. the inefficiency in flight time or propellant (or some other cost
functional to be minimized or maximized) is acceptable. However if a larger savings in flight
!52
time or propellant is desired, the true solution to the finite-burn optimal control problem must be
sought. It is clear that within practical design problems, it would be beneficial to create an initial
approximation of a desired finite-burn trajectory using a simple set of impulsive ∆Vs, and then
have the methodology to then generate a solution to the true finite-burn optimal control problem.
Due to the aforementioned difficulty in determining an estimate of the initial co-state
vector (a result of the solution’s extreme sensitivity to the initial co-states), the problem of taking
an existing sub-optimal finite-burn or multiple-impulse trajectory and converting it to an optimal
finite-burn trajectory is non-trivial. As described in Section 3.2, researchers have been working
for many decades to solve the single impulse-to-finite-burn optimal control problem. Though
they have shown great success, there still remain key gaps in the theory and methodology
regarding the application to the multiple-impulse case, as well as generalizing to the conditions
of practical, industry problems.
Handelsman
[48]
sought to validate the idea of iterating from an already determined
optimal impulsive trajectory with similar in ∆V to that of the desired finite burn. He used a
sweep method to approximate the optimal control solution of a two-dimensional, two maneuver
trajectory with an impulsive adjoint vector as an initial guess
[49]
. Handelsman’s approach was
shown to be successful in solving for two-dimensional, two-body trajectories with “bang-bang”
optimal steering control, however his method does not account for three-dimensional problems,
which is of course a key requirement for a general real-world optimization tool.
Kornhauser, Lion, and Hazelrigg
[50]
applied series expansions about an impulsive
approximation. Their method was found to be most accurate for high-thrust problems in the
single-impulse conversion case. This leaves a gap in the low-thrust regime, and does not cover in
any way the multiple-impulse conversion problem, or when the impulsive ∆Vs are relatively low
(i.e. when a set of multiple-impulses represents a continuous low-thrust arc).
Finally, Jesick and Ocampo
[51]
contributed heavily to the foundational theory for single
impulse conversion in full 3D, assuming two-body gravity. Their hybrid method employs an
interim parameter solver which takes an impulsive trajectory and converts each impulse into a
linearly steered finite-thrust arc (via sub-optimal control, direct optimization method). This
!53
parameter solution is then used in an adjoint control transformation to estimate the initial co-
states of the thrust arc in the full optimal control problem. Their work suggests extension to
multiple-impulse conversion, but does not address this new problem at all, or the resulting finite-
burn arc overlap issue that can arise. !54
5. HYPOTHESES
To reiterate, the goals of the presented research are threefold:
• Improve upon the traditional Direct Method by applying an Indirect
Method’s control law
• Formulate the optimal control problem for multi-impulse conversion, using
pseudo-rendezvous targeting on the impulsive initial guess
• Mitigate the finite-burn overlap issue by introducing additional constraints
on the ballistic arcs, and iteratively introducing the full constraint vector
To accomplish these goals, the presented research tested the following hypotheses:
1. Hypothesis: The new hybrid conversion algorithm will yield solutions to the
finite-burn optimal control problem with lower cost functionals than those found
using a traditional direct method. The new conversion algorithm will result in a reduction in propellant or time cost.
This was likely to be true for most cases as direct method solutions are often local
minima. For the direct cases involving parameter-steered approximate finite-burns,
the difference in cost with the indirect hybrid solution will be dependent on the sub-
optimal control law employed. How it was tested: Compare hybrid conversion solutions to parameter steered,
Direct Method solutions. For each test problem, solutions found via the new hybrid conversion algorithm were
validated against those from a direct parameter optimization schema. The metrics of
comparison include cost functional value, feasibility and optimality violation, and
spacecraft performance. !55
2. Hypothesis: Evolving the optimization problem definition (ballistic arc and
targeting constraints) in an iterative fashion will eliminate the overlap problem
and result in improved convergence. The added constraints to the multi-finite-burn optimization problem may reduce the
optimality of the specific overlapping impulse conversions, but will result in a
converged solution for the entire sequence of converted finite-burns. As the
optimization of the full set of burns is the goal of the procedure, its minimization
supersedes the optimality of the individual finite-burn conversions. How it was tested: Solutions were compared between conversion steps. For the presented application problems, each conversion step yields different
solutions as the optimization problem evolves. These intermediate solutions were
compared to show that the iterative constraint vector progressively reduces objective
cost throughout the conversion algorithm.
!56
6. INITIAL STUDIES
In order for the presented conversion algorithm to work, an initial trajectory comprised of
multiple-impulses or sub-optimal finite-burn arcs must be provided. To investigate such initial
trajectories, Section 6.1 describes preliminary studies that were conducted on multiple-impulse
direct optimization, where an impulsive trajectory approximates a low-thrust transfer
[1, 40]
. These
studies were presented at the International Conference on Computational & Experimental
Engineering and Sciences (ICCES ‘15), as well as the International Astronautical Congress (IAC ‘15). The proven architecture from these works served as the basis for deriving the initial
guess impulsive trajectories to be converted by the conversion algorithm.
6.1 The MID-MPCM Method for Solving Perturbed Two-Body Problems
6.1.1 Overview of the MID-MPCM
The Multi-Impulse Discretization via the Modified Picard-Chebyshev Method (MID-
MPCM) is a direct numerical optimization algorithm developed by Koblick, Xu, Fogel, and
Shankar
[1]
. The MID-MPCM finds the minimum time orbit transfer trajectory of a satellite using
a series of multiple impulses at intermediate waypoints, which may represent a continuous low-
thrust transfer or a sequence of high-thrust ∆Vs. These waypoint nodes correspond to
instantaneous impulses that are computed using a nonlinear constrained optimization routine,
such as SNOPT, MATLAB’s fmincon, or IPOPT. This is coupled with the Modified Picard-
Chebyshev Method (MPCM), a numerical scheme used to solve initial value problems (IVP) and
BVPs
[31-34]
. Previous studies
[1, 40]
explored two-body LEO-to-GEO and GEO-to-GEO transfers
with J2 perturbations, computed using the MID-MPCM, with initial guess trajectories generated
by both Edelbaum’s analytic solution
[23]
and numerical solutions provided by Kechichian
[4]
.
As described in Section 3.1.6, Edelbaum’s analytic solution solves the circular-to-circular
two-body, non-coplanar, continuous thrust transfer problem. The algorithm is valid for
inclination changes of less than 114.6˚, and assumes that the spacecraft’s thrust-acceleration
magnitude is small enough such that the instantaneous orbits during the transfer remain circular.
!57
Kechichian’s indirect numerical two-body method adheres to the standard optimal control
problem formulation discussed in Section 2.3, which seeks to minimize the Hamiltonian of the
problem’s cost functional. The Hamiltonian is a function of adjoint variables λ(t), the time-
varying co-states, and a vector function v, constant co-states, in which the first order necessary
conditions for optimality are valid
[24]
:
! (6.1)
It can be shown that to minimize the performance index of total flight time, and indirectly ∆V ,
H(t) = 1 ∈ [t0, tf] such that the transversality condition is met and the solution is time optimal
[1]
.
Kechichian derived the initial and final Lagrange multiplier (co-state) values corresponding to a
specific LEO-GEO transfer, however these co-states are difficult to generalize to other transfer
problems. The full derivations of these methods can be found in
[4, 23, 24]
. The procedure for
conducting continuous low-thrust optimization using the MID-MPCM is as follows:
Step 1: Generate Approximate Continuous, Low-Thrust Orbit Transfer: Edelbaum’s solution for two-body, non-coplanar circular orbit transfers is used as
an initial approximation for continuous low-thrust trajectories. Other feasible
initial guess solutions may be used as well, such as those found by running a
shape-based method using parameter steered sub-optimal control laws.
Step 2: Discretization of the Orbit Transfer: The low-thrust trajectory computed in Step 1 is discretized into N-number of node
points, equally spaced in time. The motivation for this is to reduce the number of
optimization variables for a full, continuous burn problem to a smaller, discrete
set of impulses, each with only thrust magnitude and direction as the controls.
Either a two-body BVP solver such as Lambert’s solution or an N-body numerical
BVP solver such as the MPCM is used to compute the instantaneous impulses
required to travel from point rN-1 to point rN. Shooting methods using IVPs can be
used to perform this optimization as well. If the MPCM is used, the two-body
Lambert solver serves as an initial guess for the MPCM.
H =
⃗ λ(t)
⊤
⃗ f + ⃗ ν
⊤
⃗ g
!58
Step 3: Direct, Non-Linear Optimization: Given a set of desired constraints, a non-linear optimization routine such as
SNOPT, fmincon, or IPOPT is used to determine the ∆V necessary for each
impulse, as well as the total time of flight, in order to minimize the objective
function (defined as either total time of flight or total ∆V).
Since this approach discretizes a continuous trajectory into multiple impulsive segments, it is a
perfect application for the MPCM, which performs best on partial orbital segments as opposed to
many revolutions
[25]
. Previous efforts demonstrated the ability to improve runtime performance
of SNOPT by over a factor of three when using the MPCM over MATLAB’s ODE113 for
numerical propagation
[1]
.
6.1.2 Trajectory Optimization using the MID-MPCM
The objective function to be minimized for the perturbed two-body problem is either total
∆V or total time of flight, tf, but not both. This is because the maximal ∆V corresponds to the
minimal tf for the multiple impulse approximation. The time of flight is discretized evenly into N
intervals. The impulse vector for each interval is represented as ! , where ! and !
are spherical azimuth and elevation angles with respect to the Earth Centered Inertial (ECI)
coordinate frame, and , where at is the magnitude of the constant low-thrust
acceleration. Using a spherical coordinate system to express the ! vectors allows the search
space to be within the range of (0—π radians) in elevation and (0—2π radians) in azimuth.
Upper and lower bounds on ηi can be provided to ensure that very large impulses are not in the
feasible region. The optimization variables (control) are written as:
! (6.2)
⃗ ΔV
i
= [η
i
, ϕ
i
,δ
i
] ϕ δ
|
⃗ ΔV
i
| =
η
i
t
f
a
t
N
⃗ ΔV
i
⃗ Z =
[
ΔV
⊤
1
ΔV
⊤
2
ΔV
⊤
3
...ΔV
⊤
N
,t
f
]
⊤
!59
The performance index for the minimum ∆V objective is then:
! (6.3)
with the gradient given by:
! (6.4) Similarly, the minimum time of flight performance index and its gradient are written as:
! (6.5)
! (6.6)
For the optimization performed in the following examples, the minimum time of flight was
chosen as the performance index, with a penalty given for when ∆Vi is larger than the maximum
low-thrust equivalent ∆Vi for the segment, i.e., when ηi > 1.
The constraints associated with orbital maneuvers depend on which mission parameters
are important to the designer. The initial conditions will remain the same regardless of which
optimization objectives one chooses. For full orbital rendezvous optimization, it may be
necessary to fix the final state vector such that the position and velocity at the end of the transfer
match that of a target spacecraft or destination state. For orbit transfers similar to those described
by Kechichian, it may be desirable to fix only five of the six orbital elements, e.g. let true or
mean anomaly be free at the final state.
6.1.3 Results from Initial Numerical Experiments
The MID-MPCM was used to perform low-thrust trajectory optimization for two cases
[1]
:
Case 1: Non-coplanar Circular Transfer (GEO-to-GEO) Case 2: Non-coplanar LEO-to-GEO Transfer
J =
∫
t
f
0
|
⃗ ΔV |dt =
N
∑
i=1
|
⃗ ΔV
i
|
∇J =
[
ΔV
i
x
|
⃗ ΔV
i
|
ΔV
i
y
|
⃗ ΔV
i
|
ΔV
i
z
|
⃗ ΔV
i
|
]
J =
∫
t
f
0
dt = t
f
∇t
f
= 1
!60
For Case 1, a constant thrust-acceleration transfer from an inclined circular GEO to equatorial
circular GEO was optimized. The initial orbital conditions are: a0 = 42000 km, e0 = 0.0, i0 =
10.0˚, Ω0 = 0.0˚, ω0 = 0.0˚, M0 = -220˚, with constant thrust-acceleration at = 1x10
-5
km/s
2
. The
final conditions of the target orbit are: af = 42,000 km, ef = 0.0, if = 0.0˚; Ωf, ωf, and Mf are
determined by initial approximation; tf is to be minimized.
Edelbaum’s analytic solution to this problem corresponds to a ∆V of 0.8419 km/s and a
time of flight of 84193.833 seconds. It was found that there is a clear local optimal setting to
achieve the highest optimality in SNOPT and the most realistic low-thrust segmentation (all ηi < 1.0), which occurs around eight impulses per revolution. The corresponding optimized
total ∆V = 0.84183 km/s, and optimized tf = 84,183.1071 sec. A 3D view of the converged
transfer is shown in Figure 13 with the node impulsive thrust vectors.
These results are compared with Edelbaum’s analytic solution, as well as the solution
found using GPOPS-II, in Table 2. It should be noted that for these comparisons, the
optimization routines minimize the total time of flight required for the transfer instead of the
minimum ∆V . Theoretically, the minimum transfer time should correspond to the minimum ∆V
for the continuous, constant-thrust-acceleration, low-thrust case. Multiplying the total transfer
time by the thrust-acceleration should produce the ∆V associated with the transfer. As the total
number of impulses increases, the solution should approach that of the constant thrust-
acceleration case (i.e. all ηi should approach unity).
Case 2 optimizes a non-coplanar LEO-to-GEO transfer described by Kechichian
[4]
. The
initial conditions for this problem are: a0 = 7000 km, e0 = 0.0, i0 = 28.5°, Ω0 = 0.0°, ω0 = 0.0°, M0 = -220°, with constant thrust-acceleration at = 9.8x10
-5
km/s
2
. The final conditions of the
target GEO are: af = 42000 km, ef = 0.001, if = 1.0°; Ωf = 0.0°, ωf = 0.0°, Mf = -43.779715°, and
tf = 58624.094 sec. This equates to ∆VLow-Thrust = (at tf )= 5.745161 km/s.
With the Lagrange multipliers specified by Kechichian as input to the low-thrust orbit
propagation IVP, the trajectory was discretized into 30 nodes which were then optimized using
SNOPT. The converged solutions showed sufficient optimality and feasibility, with a total ∆V = 5.7346 km/s, tf = 58516.19 sec, and all ηi < 1.0. These results can be compared to
!61
!62
Figure 13: Converged GEO-to-GEO non-coplanar transfer using MID-MPCM, with impulsive thrust vectors
[1]
.
Final Conditions Edelbaum’s Solution MID-MPCM GPOPS-II
af 42014.065 km 42000.000 km 42000.000 km
ef 0.0260 0.00 0.00
if 0.667º 0.000º 0.000º
Ωf 0.000º 0.000º 0.000º
ωf 0.000º 0.000º 0.000º
Mf 124.390º 124.390º 356.036º
∆Vtotal 841.94 m/s 841.83 m/s 848.25 m/s
tf 84193.833 sec 84183.107 sec 84824.731 sec
Table 2: Comparison of final orbital elements, total ∆V & flight time for the GEO-GEO transfer
[1]
.
Kechichian’s indirect solution and that of GPOPS-II (Table 3). All methods converged to the
target orbital state, with MID-MPCM and GPOPS-II yielding smaller differences in final a, e, i. As expected, the MID-MPCM total ∆V is slightly lower than using the other methods (by 0.0154
km/s) due to the method’s impulsive discretization assumption.
The unit vectors of the ∆V impulses expressed in the ECI reference frame are compared
in Figure 14. Note that the black arrows represent the solution found via Kechichian’s indirect
method, while the blue arrows are the thrust vectors at each node point via the MID-MPCM at
the same time values. These results show that the MID-MPCM’s discretized solution
successfully approximates Kechichian’s continuous indirect solution. Further analysis on runtime
performance and parallelization, as well as experiments inclusive of J2 perturbation force models,
can be found in previously published work by Koblick et al.
[1]
.
!63
!64
Final Conditions Kechichian MID-MPCM GPOPS-II
af 42000.007 km 42000.000 km 42000.000 km
ef 1.00022 x 10
-3
1.00000 x 10
-3
1.00000 x 10
-3
if 1.000012º 1.000000º 1.000000º
Ωf 359.999963º 359.999994º 359.999863º
ωf 1.966524º x 10
-2
1.67º x 10
-4
4.445109º
Mf 43.779715º 43.858910º 39.5072109º
∆Vtotal 5.75 km/s 5.7346 km/s 5.75 km/s
tf 58624.09 sec 58624.09 sec 58622.90 sec
Table 3: Comparison of final orbital elements, total ∆V , and total flight time for the converged LEO-to-GEO
transfer
[1]
.
X Y Z
Figure 14: Comparison of thrust-acceleration unit vectors (bottom) for a LEO-to-GEO transfer (top) using
Kechichian’s indirect method (blue), and the MID-MPCM seeded by Kechichian’s solution (black)
[1]
.
6.1.4 Extending the MID-MPCM to N-Body Problems: Patched Conic Edelbaum Solution
The MID-MPCM has demonstrated its robust ability to discretize a feasible initial guess
trajectory into a set of impulsive waypoints which are then optimized under the influence of
additional perturbations and force models. The natural extension of this method is to the
interplanetary, N-body regime. While the MID-MPCM optimization architecture may stay
largely unchanged, the methods for generating feasible start trajectories must be revised to suit
the types of N-body problems under consideration.
Edelbaum’s two-body, non-coplanar circular-to-circular orbit transfer solution was
adapted to create a patched-conic representation of an interplanetary low-thrust capture problem.
For the purposes of this preliminary study, a LEO-to-Low Lunar Orbit (LLO), planar-circularly
restricted three-body transfer was explored. The structure of the initial guess method is shown in
Figure 15.
First, Edelbaum’s two-body solution is computed for two separate phases: one solution
for the outward spiral from the Earth to the Earth/Moon sphere of influence (E/M-SOI) radius, with the Earth as the central body (Phase I), and one for a downward spiral
from the E/M-SOI to the target LLO with the Moon as the central body (Phase II).
When µ1 >> µ2, the SOI radius is calculated by
[23]
:
! (6.7)
The two phases are then patched together in the inertial coordinate frame of the primary body (in
this case, ECI).
Next, the Moon is forward and backward propagated in its planar-circular orbit about the
Earth from the two trajectories’ match, or “junction”, point. Due the underlying circular-to-
circular orbit assumption in Edelbaum’s solution, there is an inherent discontinuity in velocity
(magnitude and direction) between the two phases at the junction point. This is because the
circular velocity about any central body at a given position is not necessarily equal to the circular
r
SOI
= a
2
[
μ
2
μ
1
]
2
5
!65
velocity about a secondary body orbiting the primary when transformed into the second’s frame.
For example, assuming a circular aMoon = 384,399 km, the circular orbit velocity about the Earth
at the E/M-SOI is roughly 1.1192 km/s. Comparably, the circular velocity about the Moon at the
E/M-SOI is 0.2722 km/s, which when added to the Moon’s velocity yields an ECI velocity of
1.291 km/s. Clearly there is a velocity discontinuity at this match point due to the Edelbaum’s
circular assumption.
!66
2. Trajectory Patching
3. Discre2za2on &
Tuning of the Ini2al
Guess Trajectory
1A. Earth 3D Spiral (Phase I)
1B. Moon 3D Spiral (Phase II)
Figure 15: Structure of the initial guess trajectory generation algorithm via Edelbaum patched conics for an
Earth-Moon transfer (with steps 2 & 3 in ECI)
This means that the flight path angle (FPA) and position difference at the junction point
can be minimized together, but they cannot both be eliminated simultaneously. If the position
discontinuity is eliminated, the difference in FPA at the match point will result in a large
additional ∆V at the nearby discretized nodes. If the FPA discontinuity is removed, there will be
a non-trivial gap in position at the match point which dramatically increases node ∆Vs. Instead
of minimizing either discontinuity, rotational tuning of the Phase I & II trajectories must be
conducted to minimize the total ∆V after discretization and N-body propagation via the MPCM.
This is described in the next subsection. Once tuned, the patched, discretized trajectory is
optimized using the same architecture as in previous MID-MPCM cases.
6.1.5 Extending the MID-MPCM to N-Body Problems: Tuning the Patched Conic Edelbaum Solution
Without appropriately tuning the patched trajectory, the two Phases may be matched to
minimize both the flight path angle and position discontinuity simultaneously (Figure 16).
However, as described in the previous section, this geometry may not yield the lowest total ∆V
initial guess solution after discretization and propagation, making it difficult for the optimizer to
converge. Instead, the Phase I & II trajectories are parametrically rotated about their respective
central bodies, discretized, then propagated to find geometries that minimize the total node ∆V .
For the LEO-LLO problem presented below, it was observed that the MID-MPCM
optimization method was able to rapidly converge to a solution for segments of the transfer well
within each of the body’s SOIs. Similar to the previous GEO-to-GEO and LEO-to-GEO
experiments, the third-body perturbation effects in these domains are easily accommodated by
the MID-MPCM, as Edelbam’s circular transfer approximation is close enough to the final
optimized trajectory. Conversely, it was found that the patched-Edelbaum segments approaching
and crossing the SOI radius diverge significantly from the optimized solution. As such, attention
was focused on this transition region, with the majority of the outward spiral from LEO and
downward spiral to LLO post-capture assumed to be optimized separately.
!67
!68
Figure 16: Example LEO-to-LLO patched-Edelbaum initial guess trajectory in ECI (top) and MCI (bottom),
tuned to minimize both ∆flight path angle and ∆position at the junction point. Alt0 = 300 km (ECI), Altf =
500 km (MCI), constant thrust-acceleration = 50 µg, total ∆V = 8.1 km/s, transfer time = 6.6 months.
When selecting the trajectory matching geometry, two cases are be considered: prograde
capture about the Moon, and retrograde capture. In this context, a prograde capture corresponds
to a transfer whose Phase II two-body angular momentum vector is in the same direction as in
Phase I, i.e. the spacecraft passes behind the Moon during the encounter. Conversely, a
retrograde capture is one where the Phase II angular momentum vector is opposite to that of
Phase I, i.e. the spacecraft passes in-front of the Moon. Both cases produce feasible solutions that
result in different total ∆Vs. Finding both types of solutions is critical for LTGA trajectory
design, as this helps control the spacecraft’s velocity vector change over the course of the fly-by.
6.1.6 Extending the MID-MPCM to N-Body Problems: Results from Patched Conic Edelbaum Experiments
A parametric study was conducted for a 5 µg thrust-acceleration LEO-to-LLO patched-
Edelbaum trajectory, with an initial circular Earth orbit radius of 187,730.1 km, and a final
circular lunar orbit radius of 29,999.1 km. The two Edelbaum phases were calculated for a match
point 18,114.6 km into the Moon’s SOI. While the match point can be set to the E/M-SOI or any
other radius, initial optimization attempts indicated that a match point closer to the Moon
improved N-body convergence. This is likely due to the low thrust-acceleration magnitude used
in the transfer, which reduces the ∆V per node available to aid in lunar capture during the
relatively fast SOI approach & crossing (compared to the complete transfer time of flight).
The trajectory was discretized into nine node points equally spaced in time that span the
SOI crossing region, using the MPCM for N-body propagation between the nodes. The number
of node points was selected based on previous experience with the MID-MPCM and SNOPT. To
tune the patched-conic Edelbaum initial guess, the two trajectory phases were parametrically
rotated about their respective central body by indexing through a range of initial mean anomalies
(M0) for the Phase I spiral in ECI, and Phase II spiral in MCI (Moon Centered Inertial). The
resulting total ∆Vs of the discretized waypoints versus the initial M0’s are shown in Figure 17,
indicating the minimum ∆V configurations of the initial guesses for both prograde and retrograde
capture at the Moon. The corresponding feasible solutions for minimum total ∆V in Figure 18.
!69
!70
Figure 17: Parametric rotation of initial mean anomalies of Edelbaum spirals in ECI & MCI for prograde
capture (top) and retrograde capture (bottom) versus the resulting total ∆V upon discretization via MID-
MPCM.
Minimum ∆V = 0.483640 km/s
Retrograde Capture
Minimum ∆V = 0.531426 km/s
Prograde Capture
!71
1.
2.
3.
4.
5.
6.
7.
8.
9.
1.
2.
3.
4.
5.
6.
7.
8.
9.
--- ∆V for Constant Low-Thrust
--- ∆V for Constant Low-Thrust
1.
9.
1.
9.
2.
2.
3.
4.
5.
6.
7.
8.
3.
4.
5.
6.
7.
8.
−− Raw Edelbaum Solu2ons
-•- Discre2zed Nodes
-•- Moon Orbit Nodes
−− MPCM N-Body Propaga2on
→ ∆V Impulse Unit Vector
1.
9.
1.
9.
2.
2.
3.
4.
5.
6.
7.
8.
3.
4.
5.
6.
7.
8.
ηmax = 20.17
ηmin = 0.38
ηmax = 24.78
ηmin = 0.36
tf = 370.98 hrs
∆V = 483.640 m/s
tf = 370.98 hrs
∆V = 531.426 m/s
Retrograde Capture
ECI Trajectory ECI Trajectory
MCI Trajectory MCI Trajectory
Prograde Capture
Figure 18: Initial guess SOI-crossing trajectories with minimum total ∆V for Prograde Lunar Capture (left) and
Retrograde Lunar Capture (right). The discretized trajectory, with N-body propagation between nodes using
the MPCM, is shown in ECI (top) and MCI (middle), along with the corresponding ∆V impulse magnitudes
for each node (bottom). The red dash-line represents the nodal ∆V required for continuous, low-thrust
approximation.
For the prograde case, M0 / ECI (Phase I) was varied between 221.2˚ and 251.2˚ in steps
of 0.5˚, and M0 / MCI (Phase II) was varied between 87.5˚ and 177.5˚ in steps of 0.5˚. The
minimum total ∆V for the prograde case was found to be 0.531426 km/s, which corresponded to
an M0 / ECI of 237.7˚ and M0 / MCI of 153.5˚. For the retrograde case, M0 / ECI (Phase I) was
varied between 197.6˚ and 277.6˚ in steps of 0.5˚, and M0 / MCI (Phase II) was varied between
-14.0˚ and 136˚ in steps of 0.5˚. The minimum total ∆V for the retrograde case was found to be
0.483640 km/s, which corresponded to an M0 / ECI of 237.6˚ and M0 / MCI of 61˚. A ∆V
savings of 0.047784 km/s over the SOI-crossing segment can be achieved by selecting the
retrograde capture. This is as expected, as passing in front of a perturbing body during a fly-by
encounter reduces the spacecraft’s velocity, which in this case aids in lunar capture.
In both capture cases, there is a distinct gap in ECI position between the two Edelbaum
phases of 46,049.6 km in the prograde case, and 6,467.7 km in the retrograde case. Upon discretization, the MPCM BVP solver and propagator smooths out this discontinuity by
distributing additional ∆V across multiple nodes. While the added ∆V decreases the initial
trajectory’s optimality, as the node ∆Vs are generally larger than the discretized low-thrust
equivalent ! , it is now without discontinuities and can serve as a starting point for
the optimization routine. The MID-MPCM objective function can be set to include a factor that
drives down the ratio of the propagated node ∆Vs to the discretized low-thrust ∆V to unity (i.e.,
minimize ∑|ηi - 1|). Similarly, a constraint can be added to enforce |ηi = 1|.
(|Δ
⃗ V
i
|) =
t
f
a
t
N
!72
7. MULTI-IMPULSE TO FINITE-BURN TRAJECTORY CONVERSION
7.1 Assumptions
The presented conversion algorithm was developed under the following assumptions:
Conversion Solution Structure: The conversion algorithm is derived for a pseudo-
rendezvous optimal control problem composed of multiple burns. The output solution to
this algorithm is constructed with a coast-thrust-coast structure that repeats for the
number of burns input to the algorithm. The coast arcs between the computed finite-
burns are allowed to have flight times greater than or equal to zero.
Input Trajectory: The trajectory input to the conversion algorithm is assumed to be
comprised of a set of impulsive and/or sub-optimal, parametrically steered finite-burns,
separated by ballistic arcs. The spacecraft is assumed to have a constant specific impulse,
variable thrust engine model.
Force Model: A two-body, point-mass gravity field is used to derive the dynamics and
boundary constraints associated with the pseudo-rendezvous optimal control problem.
Primer Vector Alignment: The initial thrust directions provided by the input trajectory
are assumed to be closely aligned with the optimal Primer Vector direction, as dictated
by the co-state dynamics control law. This enables a series of adjoint-control
transformations to be used to estimate the initial co-states from the input trajectory. This
represents the fundamental assumption of this work, which then requires the input
trajectory to be a reasonable approximation of the true time-optimal finite-burn solution.
This does not require a priori knowledge of the solution beyond a simple calculation to
ensure that the engine model has enough time to provide the input ∆V before the next
impulse. Even when this is not the case, the algorithm is designed to mitigate the
resulting overlap problem its omission creates.
!73
Optimization Settings: The software program Copernicus 4.6
[45]
was used to simulate
all trajectories. A variable order, variable step-size Adams-Basforth-Moulton integration
method
[61]
was utilized for propagation, with absolute and relative error tolerances of
1.0E-12. The constrained optimization system SNOPT
[26]
was employed with the
following key settings which were found to improve convergence:
Major Optimality Tolerance: 1.0E-03
Major Feasibility Tolerance: 1.0E-06
Major Step Limit: 0.01
Gradient Perturbation Step Size, Difference Interval: 1.0E-07
Gradient Perturbation Step Size, Central Difference Interval: 1.0E-06
Gradient Sparsity Mode: Assume Dense
Lastly, a few comments on the optimization solution scheme, namely second-order
necessary conditions for optimality. Neither the work of Jesick & Ocampo
[51]
nor the work in
this dissertation compute the second-order necessary conditions for optimality. This is omitted
purposely, as such second derivatives (i.e. the Hessian matrix) are often difficult to derive, or
numerically taxing to compute or store for each re-definition of the optimization problem at
hand. This decision is acceptable, and in no way detracts from the merits of this work, as
justified by the following points.
First, to reiterate, the prime goal of this dissertation is to provide an improvement over
traditional Sub-Optimal Control (SOC) Direct Methods specific to the multi-impulse conversion
problem, not to replicate a full Indirect Method, or make any comparative statement regarding
the two. The well-known
[40]
disadvantage of Indirect Methods is their inability to be easily
adapted to different trajectory problems, as a re-derivation of the second-order necessary
conditions for optimality is required, which can become prohibitively complex. Conversely,
Direct Methods often do not compute second-order information, and are therefore both
inherently sub-optimal, and easier to apply to diverse problems. It is this second fact that
!74
motivated this dissertation. The objective of this dissertation was to replace the SOC control law
in a parameter optimization-based solver (i.e. Direct Method) with an improved control law
derived from the Euler-Lagrange equations from the finite-burn optimal control TPBVP (i.e. an
Indirect Method). The intent was to apply a key element of an Indirect Method to a traditional
Direct Method, in order to improve upon the Direct solution.
This leads to the second justification for the omission of second-order information in this
work. Direct Methods, as well the Hybrid Method developed in this dissertation, often utilize
sequential quadratic programming (SQP) algorithms to solve constrained optimization problems,
which do not always require second derivative information. Such optimization tools approximate
the Hessian matrix at each iterate, and check to confirm that the solution has converged to at
least a local minima. The SQP solver SNOPT
[24]
was used for just this purpose in this
dissertation. Selections from the SNOPT user’s guide
[64]
are presented below to provide a brief
description of the algorithm’s solution approach.
“SNOPT finds solutions that are locally optimal, and ideally any nonlinear
functions should be smooth and users should provide gradients. It is often more
widely useful. For example, local optima are often global solutions, and
discontinuities in the function gradients can often be tolerated if they are not too
close to an optimum. Unknown gradients are estimated by finite differences.
SNOPT uses a sequential quadratic programming (SQP) algorithm. Search
directions are obtained from QP subproblems that minimize a quadratic model of
the Lagrangian function subject to linearized constraints. An augmented
Lagrangian merit function is reduced along each search direction to ensure
convergence from any starting point.
Sequential quadratic programming (SQP) methods have proved highly effective
for solving constrained optimization problems with smooth nonlinear functions in
the objective and constraints. Here we consider problems with general inequality
constraints (linear and nonlinear)… Second derivatives are assumed to be
unavailable or too expensive to calculate.
!75
…A point is regarded as a satisfactory solution if it satisfies the first-order
optimality conditions to within certain tolerances. We emphasize that this point
may not be a constrained local minimizer for the problem. Verifying second-order
conditions requires second derivatives.
…The basic structure of the SQP algorithm involves major and minor iterations.
The major iterations generate a sequence of iterates{xk} that satisfy the linear
constraints and converge to a point that satisfies the nonlinear constraints and the
first-order conditions for optimality. At each xk a QP subproblem is used to
generate a search direction toward what will be the next iterate xk+1. The
constraints of the subproblem are formed from the linear constraints ALx − sL = 0
and the linearized constraint f (xk) + f0 (xk)(x − xk ) − sN = 0, where f0(xk)
denotes the Jacobian matrix, whose elements are the first derivatives of f (x)
evaluated at xk.
SNOPT uses a quasi-Newton approximation to the Hessian of the Lagrangian. A
BFGS (Broyden-Fletcher-Goldfarb-Shanno) update is applied after each major
iteration.”
As detailed above, SNOPT handles the local minimization internally, meaning the
provision of derivative information is not required for successful convergence to at least a local
minima. While this would not be sufficient to solve the full optimal control TPBVP (via Indirect
Method), the goal of this dissertation is simply to improve upon the Direct Method, as
demonstrated later sections.
7.2 The Multiple Finite-Burn Conversion Algorithm
A novel conversion algorithm is presented that extends the hybrid optimization method
by Jesick & Ocampo
[51]
to a new problem domain. When attempting to solve general optimal
control finite-burn trajectory optimization problems (e.g. low-thrust orbit transfers) using a
!76
control law derived from the Euler-Lagrange equations, difficulty arises in autonomously
deriving a sufficient initial co-state vector at ignition to enable solution convergence. Due to the
extreme sensitivity of convergence to these initial co-states (i.e. the control), such optimal
control methods are often difficult to generalize and automate, as prior knowledge of the problem
domain is required.
The presented method mitigates the difficulty in generating the approximate initial co-
state vector by evolving the optimal control Two-Point Boundary Value Problem (TPBVP) across
multiple solution iterations, culminating in the desired problem definition. This is done by
iteratively applying adjoint control transformations and slowly introducing the full constraint
vector. This provides an increasingly more accurate approximation of the initial co-states, which
by the end of the process become sufficiently close to the optimal values to enable solution
convergence for the desired problem. The method requires no prior knowledge of the solution,
merely a multi-impulse or suboptimal finite-burn initial guess, which aids automation.
An algorithm is herein described that converts a set of multiple-impulses or parameter
steered sub-optimal finite-burns, representing the entirety or a portion of a high- or low-thrust
maneuver, to an exact time optimal finite-burn trajectory for a thrust limited, constant exhaust
velocity spacecraft. This can be called a Pseudo-Rendezvous (P-R) problem, as it yields a finite-
burn solution whose final time, position and velocity state is equal to that of the original post-
multi-impulse trajectory.
A conceptual overview of the conversion procedure is shown in Figure 19. First, an
approximate multi-impulse (or sub-optimal finite burn) trajectory is used as input to the
algorithm. Each burn is then replaced with a time-optimal finite-burn computed using the single
impulse conversion method by Jesick and Ocampo
[51]
, targeting the ballistic arcs between
neighboring burns. In some circumstances, this may yield a Finite-Burn Arc Overlap, which is
said to exist when the coast time between two sequential finite-burn segments is negative. This
problem can arise when patching together a series of single-impulse conversions with large ΔVs
too close together in time (relative to the engine performance), meaning the bounded time
between impulses is not sufficient for the engine to produce the required finite-burn ∆V .
!77
The Finite-Burn Arc Overlap problem is mitigated through a process of introducing
additional boundary constraints, and by relaxing the pseudo-rendezvous state targeting on all
interior conversions, leaving the state targeting for only the final impulse conversion. These
redefinitions of the optimal control TPBVP allow the finite-burns to be re-optimized, while
ensuring the final pseudo-rendezvous targeting is upheld, which is the stated goal of the
procedure. This process also has the benefit of improving optimality for the entire sequence of
!78
Figure 19: Conceptual Overview of Multi-Impulse to Multi-Finite-Burn Conversion Algorithm. 1) Multi-
impulse trajectory input. 2a) Set of converted impulses, targeting interior ballistic arcs. 2b) Existence of
finite-burn arc overlap, with negative coast time on interior ballistic arc. 3) Re-targeted solution, relaxing
targeting constraints and adding positive interior coast time constraints.
burns, as by definition each finite-burn was optimized individually, not as a member of the larger
set of conversions. By evolving the targeting constraints for the entire set, regardless if any
overlaps exist, the optimizer is afforded greater freedom when minimizing the full multi-finite-
burn problem, as the burn out states of each finite-burn are no longer constrained to be on the
inherently sub-optimal impulsive input trajectory.
The adjoint-control transformation central to the solution method represents a mapping
across three fundamental burn models: Impulsive, Sub-Optimal Control (SOC), and Optimal
Control Theory (OCT) control laws [Figure 20]. The impulsive model has the thrust vector
comprised of ΔV magnitudes and spherical angular directions, with instantaneous mass change.
The SOC model represents a sub-optimal finite-burn, with thrust directions described as a
function of polynomial-based spherical steering angles. The OCT model represents an optimal
control finite-burn, which requires the numerical integration of the Euler-Lagrange system of
equations from the standard formulation of the optimal control TPBVP
[51]
(Section 2.3).
Selecting the proper parameterizations within these models allows for the transfer of
solution information from one to another, which can be used to estimate the initial co-state vector
for each finite-burn conversion. Instead of attempting to estimate the full vector (7x1), selecting
the parameterization shown in Figure 20 replaces either 2 or 4 of these values with thrust angles
which can be inferred from the control parameterizations of the input trajectory. This leaves only
3 additional values that are needed to complete the estimate of the initial co-state vector, one of
which is known from Jesick & Ocampo’s single impulse conversion strategy
[51]
: .
Through experimentation, it was found that setting ! provided suitable
convergence for the problems explored. Monotonic Basin Hopping (MBH), i.e. a bounded
randomization of the control variables
[62]
, was utilized to encourage convergence in the event
that the optimizer exited prematurely.
The last piece of theory needed to construct the multi-finite-burn conversion algorithm is
how to express the pseudo-rendezvous targeting within the TPBVP. As discussed in Section 2.2,
the post-impulse or SOC ballistic arc is assumed to be defined by the input trajectory. As such,
|λ
v
i
0
| = 1.0
|
·
λ
v
i
0
| = λ
m
i
0
= 0
!79
!80
Figure 20: Conversion algorithm burn model parameterizations. Mapping from Impulsive & SOC thrust
control to OCT initial co-state vector.
this ballistic arc can be treated as a known function of time, enabling a simplification of the
boundary constraints. However, to implement this within the conversion algorithm, a phantom
segment was added to the trajectory structure to provide a variable state for the finite-burn to
target. This Pseudo-Rendezvous Segment is shown in Figure 21. The P-R segment is defined as
inheriting the time and state immediately after the ∆V impulse, and has a final time TF equal to
the TF of the OCT finite-burn segment. This means the OCT burn-out time will always be
identical to the TF of the P-R segment. With identical TFs, the OCT FB segment can now target
its final state to be equal to that of the P-R segment. Satisfying this constraint ensures that the
post-OCT FB state and time are identical to that of the original input trajectory, satisfying the
requirement that the conversion yield a pseudo-rendezvous trajectory.
Now the full multi-finite-burn conversion algorithm can be constructed, as shown in
Figure 22. The algorithm is comprised of 3 key evolutions of the optimal control problem
(Tables 4—6), which each seed the next Sub-Problem (SP). This structure was found to aid in
solution convergence, as each interim solution guides the optimizer closer to the optimal initial
co-state vector. Attempts to skip immediately to the final optimization sub-problem were met
with peril due to the sensitivity of the initial co-states, and significant difference between the
input trajectories and the true solutions to the TPBVP.
!81
Figure 21: Implementation of Single Impulse Conversion, with added Pseudo-Rendezvous Segment used for
burn-out targeting.
!82
Figure 22: Multi-Finite-Burn Conversion Algorithm, with example low-thrust inclination change problem
(detailed in Section 7.3.2). Note the addition of the final equality constraint, switching function $(t), in
Step 5, to avoid over-constraining the earlier sub-optimization problems. The 3 Optimization Sub-
Problems are identified (right).
Optimization 1
Optimization 2
Optimization 3
Algorithm Step
∆Vs
∆Vs
∆V
∆T
∆T
!83
Objective Function Control Vector
Variable Initial Guess
∆t
0.000
0.000
Size of Control Vector: (8 x # of burns) x 1
Constraint Vector Size
(6 x # of burns) x 1
(1 x # of burns) x 1
(1 x # of burns) x 1
Size of Constraint Vector: (8 x # of burns) x 1
0 OR ! from SOC input
·
α
i
o
SOC
!|
⃗ λ
v
i
f
|≥ 0.0
!α
i
o
!J = min
∑
(t
i
f
−t
i
0
)
n
!λ
M
i
o
! estimate from input (t
i
f
−t
i
0
)
! ⃗ X
i
f
=
⃗ X
PR
f
!t
i
0
! ·
β
i
o
! from input β
i
o
! ·
α
i
o
! from input α
i
o
!|
·
⃗ λ
v
i
o
|
!λ
M
i
f
= 0.0
! estimate from input t
i
0
0 OR ! from SOC input
·
β
i
o
SOC
!β
i
o
Table 4: Optimization Sub-Problem 1, with state & co-state dynamics as defined in Section 2, for a set of n-
number of burn conversions. The selection of constraints represents a majority of the full TPBVP
constraint vector, reduced to not over-constrain the problem. Note an equal number of controls and
constraints (i.e., problem is square). EXAMPLE: For a 3-burn conversion, the control vector would be 24x1, with the constraint vector 24x1.
!84
Objective Function Control Vector
Variable Initial Guess
from Sub-Problem 1 solution
∆t
from Sub-Problem 1 solution
from Sub-Problem 1 solution
from Sub-Problem 1 solution
from Sub-Problem 1 solution
from Sub-Problem 1 solution
from Sub-Problem 1 solution
Size of Control Vector: (8 x # of burns) x 1
Constraint Vector Size
6 x 1 (final burn only)
(1 x # of burns - 1) x 1
(1 x # of burns) x 1
(1 x # of burns) x 1
Size of Constraint Vector: (3 x # of burns + 5) x 1
!λ
M
i
f
= 0.0
!α
i
o
!J = min
∑
(t
i
f
−t
i
0
)
n
!λ
M
i
o
! from Sub-Problem 1 solution (t
i
f
−t
i
0
)
! ⃗ X
i
f
Final−Burn
=
⃗ X
PR
f
Final−Burn
!t
i
0
! ·
β
i
o
!|
⃗ λ
v
i
f
|≥ 0.0
! ·
α
i
o
!ΔT
coast
≥ 0
!|
·
⃗ λ
v
i
o
|
!β
i
o
Table 5: Optimization Sub-Problem 2, with state & co-state dynamics as defined in Section 2, for a set of n-
number of burn conversions. The P-R targeting constraints are removed for all but the final burn, and a
constraint is added to ensure interior coast times are positive. This frees the optimizer to find cheaper
neighboring solutions to those found in Sub-Problem 1, and reduces the number of constraints to always be
less than the number of controls (for multiple-finite-burn conversions), in preparation for Sub-Problem 3. EXAMPLE: For a 3-burn conversion, the control vector would be 24x1, with the constraint vector 14x1.
!85
Objective Function Control Vector
Variable Initial Guess
from Sub-Problem 2 solution
∆t
from Sub-Problem 2 solution
from Sub-Problem 2 solution
from Sub-Problem 2 solution
from Sub-Problem 2 solution
from Sub-Problem 2 solution
from Sub-Problem 2 solution
Size of Control Vector: (8 x # of burns) x 1
Constraint Vector Size
6 x 1 (final burn only)
(1 x # of burns - 1) x 1
(1 x # of burns) x 1
(1 x # of burns) x 1
(1 x # of burns) x 1
Size of Constraint Vector: (4 x # of burns + 5) x 1
!β
i
o
!λ
M
i
f
= 0.0
!α
i
o
!J = min
∑
(t
i
f
−t
i
0
)
n
!λ
M
i
o
! from Sub-Problem 2 solution (t
i
f
−t
i
0
)
! ⃗ X
i
f
Final−Burn
=
⃗ X
PR
f
Final−Burn
!t
i
0
! ·
β
i
o
!|
⃗ λ
v
i
f
|≥ 0.0
! ·
α
i
o
!ΔT
coast
≥ 0
!|
·
⃗ λ
v
i
o
|
!$(t
i
f
) = $(t
i
o
)
Table 6: Optimization Sub-Problem 3, with state & co-state dynamics as defined in Section 2, for a set of n-
number of burn conversions. Note the addition of the switching function constraint, completing the
multiple-finite-burn optimal control TPBVP. The reduced P-R constraints are maintained, ensuring the
number of constraints will still be always less than the number of controls (for multi-finite-burn
conversions), and the interior coast time constraint is maintained. EXAMPLE: For a 3-burn conversion, the control vector would be 24x1, with the constraint vector 17x1.
7.3 Application Problems
For the subsequent application problems, the conversion algorithm was modeled in
NASA JSC’s in-house orbit trajectory optimization software Copernicus 4.6
[45]
, however the
procedure may be implemented within any trajectory optimization schema. Examples are shown
for low-thrust apogee raise maneuvers, non-coplanar Earth orbit transfers, a two-burn Earth-
Mars transfer, and a modified Deep Space 1 low-thrust trajectory. Results are compared to
conventional Direct SOC solutions. The Hybrid OCT solutions are shown to offer a numerically
significant improvement to objective cost savings across all application problems, as well as a
notable improvement to convergence speed for a select class of problems.
Due to the current industry push to provide rapid advancement in the realm of high-
powered Solar Electric Propulsion (SEP), one of NASA’s state-of-the-art SEP spacecraft designs
was selected as the low-thrust test article for several of the application problems. NASA’s Deep
Space Transport (DST)
[65]
shown in Figure 23 is a hybrid SEP/chemical crewed spacecraft
designed to ferry astronauts to and from destinations such as Mars on missions lasting in excess
of 1000 days. The single-launch reusable system will be aggregated in cis-lunar space (namely
the Gateway space station), where it will depart and arrive from its deep space missions.
!86
Figure 23: Artist depiction of the Deep Space Transport
[65]
. Note the two arrays of hall effect thrusters on
angled trusses near the center of the image, and the two deployable solar arrays.
The DST power generation system is comprised of two 250 kW (BOL, @1 AU) solar
arrays with IMM4J solar cells (33% efficient). This provides ~470 kW solar array power at the
start of mission, or ~300 kW of electric propulsion (EP) power. The hybrid propulsion system is
split between 4 main LOx/MCH4 (methane) thrusters, providing a total thrust of 17.8 kN at a
specific impulse of 351 sec, and two gimbaled arrays of 24 Hall Effect Thrusters that provide a
total thrust of 19.2 N at 2000 — 3000 sec of specific impulse using Xenon propellant.
The DST has ~24t of EP propellant, and ~16t of chemical propellant. The dry mass of the
vehicle is estimated at ~21.5t. For the purposes of this dissertation, the chemical propellant was
omitted as the trajectories computed were purely low-thrust, as the mission can be assumed to be
a test flight exercising the SEP systems. The total wet mass of the modeled spacecraft was set to
40t, with 19.2 N of max thrust at an average specific impulse of 2500 sec. The propellant mass
fraction was not set a priori for the problems explored, but verified after optimization
convergence. All computed trajectories were found to be well within the SEP propellant mass
budget laid out in the DST design.
7.3.1 Low-Thrust: Apogee Raise
The first introductory test problem is a simple two-burn, low-thrust apogee raise
maneuver sequence. An input trajectory was generated with two ∆V impulses occurring at
perigee on consecutive revolutions, starting from an initial orbit with elements:
! , with a final orbit with
elements: ! ! . The ∆V vectors are aligned in the velocity
vector direction, ensuring pure apogee raise. The identical ∆V magnitudes provide a total change
in apogee altitude of 2,954.159 km. A SOC solution and OCT solution were generated using the
multi-finite-burn conversion algorithm, whose solutions are detailed in Figure 24, Table 7 and
Figure 25. For this simple problem, the thrust profiles are nigh identical, however the OCT
solution manages to save 7.2 sec of burn time, or an equivalent 5 grams of propellant mass.
While not a substantial savings in cost, this does validate the conversion algorithm as being able
a
0
= 16,378.14 km, e
0
= 0.0, i
0
= 28.474 deg,Ω
0
= ω
0
= 0.0 deg
a
f
= 17,855.220 km, R
p
f
= 16,378.14 km, R
a
f
= 19,332.299 km, e
f
= 0.082725,
i
f
= 28.474 deg,Ω
f
= 0.0 deg, ω
f
= 170.856 deg
!87
!88
Figure 24: Solution comparison of two-burn, apogee raise conversion from multi-impulse input.
Impulsive Parameters
2,954.159 km
0.2 km/s
0.0 deg
0.0 deg
SOC Parameters OCT Parameters
Thrust 19.2 N 19.2 N
2,500 sec 2,500 sec
Total Burn ∆t 3 hrs 11 min 51.5 sec 3 hrs 11 min 44.2 sec
1,000 kg 1,000 kg
Propellant Used 9.015 kg 9.009 kg
Effective ∆V 221. 1 m/s 222.0 m/s
!I
sp
!ΔR
a
β
i
o
!α
i
o
!M
i
o
Total !|
⃗ ΔV |
Table 7: Two-burn, low-thrust apogee raise solution burn parameters.
to provide a solution at least as good as a SOC solution, and in this case marginally improved.
Moreover, it was found that converted finite-burns whose burn ∆t are less than about a quarter of
a revolution are more likely to yield comparable SOC and OCT solutions. Only when longer
duration burns (relative to orbit period) or more complex maneuvers are demanded does the
hybrid conversion algorithm yield appreciable benefits over SOC-type solutions.
!89
Figure 25: Change in radius (top) and velocity (bottom) vs. time for apogee raise conversion (OCT solution).
Initial orbit (purple), OCT finite-burns (red), intermediate coast orbit (blue), and final orbit (green).
7.3.2 Low-Thrust: Inclination Change
The next application problem explored was a two-impulse to low-thrust inclination
change conversion from an initial circular orbit with elements: ! ! . Two ∆V impulses (100 m/s each) were place half a
revolution apart, with thrust directions -/+ 90 deg out-of-plane to yield a change in inclination of
! per impulse, or 2.292 deg total. Conversion solutions were generated, with
the optimization SP-1 solution shown in Figure 26 and Table 8. As with the first application
problem, the OCT and SOC solutions are comparable, with a slight benefit to OCT in terms of
iteration time and optimality. Note the OCT solution’s savings in cost of 43.4 sec of burn time, or
34 grams of propellant, over the SOC solution. This represents a larger cost savings compared to
the apogee raise problem, since the converted finite-burns for this problem have longer durations
for the same spacecraft model.
a
0
= 16,378.14 km, e
0
= 0.0,
i
0
= 28.474 deg,Ω
0
= ω
0
= 0.0 deg
Δi
imp
= +1.146 deg
!90
Figure 26: Two-impulse inclination change conversion, Optimization Sub-Problem 1 solution. OCT thrust
vectors (red) and SOC thrust vectors (purple) are shown, with initial orbit (blue dashed) and final orbit
(magenta dashed). A planar projection (right) shows the two finite-burn arcs (yellow).
To fully exercise the conversion algorithm, it was desired to induce a finite-burn arc
overlap. This was done by increasing the first ∆V impulse to 150 m/s, resulting in a longer finite-
burn and larger change in inclination (! ). The overlap is shown in Figure 27
(SP-1 solution). After solving SP-2 and 3 (evolving the constraint vector by relaxing the P-R
targeting, and then adding the switching function constraints), the overlap is eliminated and the
full multi-finite-burn TPBVB solved (Figure 28, Table 9). Note that in addition to mitigating the
finite-burn arc overlap problem, the Hybrid OCT Sub-Problem 3 solution also yields a lower cost
compared to that of Sub-Problem 1, illustrating the benefits of solving all 3 sub-problems,
regardless of the existence of an overlap. Note that the Direct SOC solution was unable to
converge to the first sub-problem, but was able to converge to SP-2. This is discussed in greater
detail in subsequent sections.
Δi = +1.718 deg
!91
Impulsive Parameters (per impulse)
1.146 deg
100 m/s
0.0 deg
SOC Parameters OCT Parameters
Thrust 19.2 N 19.2 N
2,500 sec 2,500 sec
Total Burn ∆t 3 hrs 18 min 57.2 sec 3 hrs 18 min 13.8 sec
1,000 kg 1,000 kg
Propellant Used 9.349 kg 9.315 kg
Effective ∆V 229.8 m/s 229.51 m/s
!Δi
!|
⃗ ΔV |
!α
i
o
β
i
o
!I
sp
! deg −/+90.0
!M
i
o
Table 8: Two-burn, low-thrust inclination change Sub-Problem 1 solution burn parameters.
!92
Figure 28: Mitigated Finite-Burn Arc Overlap (Optimization Sub-Problem 3 solution) for two-impulse
inclination change conversion.
Figure 27: Finite-Burn Arc Overlap Problem for two-impulse inclination change conversion. Optimization
Sub-Problem 1 solution. Note the inefficiency in the first burn thrust direction slew (red).
To further provide proof of advancement, optimization data was collected that clearly
reinforces the improvement to the state of the art offered by the new Hybrid Method. Table 10
compares optimization solution data between the Hybrid and Direct Methods for the central
conversion problem: 2-burn LEO ∆i with an overlap. The first observation is that the Hybrid
Method converges for all three of the sub-problems, while the Direct Method only converges for
the second sub-problem (and is incapable of solving the third, as it involves adding the switching
function constraint, which is not a part of the SOC control law). As such, the remainder of the
table only compares solutions after the first and second sub-problems. Despite solving two
problems instead of one, the Hybrid Method still converges to a final SP-2 solution in a total of
128 iterations, which is 25 iterations faster than the Direct Method. This also equates to a 38.5%
reduction in convergence time, a substantial benefit when it comes to automation. This
improvement can be attributed to the overlap inherent to the SP-1 solution, which also explains
the inability of the Direct Method to converge on SP-1. The overlap in this problem yields an
OCT SP-1 solution with a burn spanning greater than half a revolution:
! ΔT
burn
>
T
Period
2
!93
Impulsive Parameters (per impulse)
1.146 deg, 1.718 deg
150 m/s, 100 m/s
0.0 deg
SOC Parameters OCT Parameters
Thrust 19.2 N 19.2 N
2,500 sec 2,500 sec
2 hrs 27 min 55.7 sec 2 hrs 27 min 55.1 sec
2 hrs 27 min 46.1 sec 2 hrs 27 min 44.2 sec
1,000 kg 1,000 kg
Propellant Used 14.013 kg 13.892 kg
!Δi
!|
⃗ ΔV |
!α
i
o
β
i
o
!I
sp
! deg −/+90.0
Finite-Burn !Δt
1
Finite-Burn !Δt
2
!M
i
o
Table 9: Two-burn, low-thrust inclination change Sub-Problem 3 solution burn parameters.
This means that due to the targeted plane change with each burn roughly centered around
the line of nodes, the converged OCT thrust vector slews nearly 180º during the burn arc as the
spacecraft nears the opposite node to maintain consistent change in inclination (which can be
seen in Figure 27). This rapid change in thrust direction is easily handled by the Hybrid OCT
control law (derived from the Euler-Lagrange equations), but the linearly-steered SOC control
law is incapable of such a slew in such a short fraction of the orbit period. This demonstrates that
the Hybrid OCT Method will yield superior convergence behavior for conversion problems
where both ! , and the burn is changing inclination. This specific class of
problem will therefore reap larger benefits than other problems when the Hybrid OCT Method is
employed in lieu of Direct SOC.
Figures 29– 33 compare the convergence behavior for the Hybrid OCT and Direct SOC
Methods for this overlap problem. In SNOPT, the feasibility metric is a measure of how
infeasible the constraint vector is at each iterate, or in other words, the maximum constraint
ΔT
burn
>
T
Period
2
!94
Application Problem: 2-burn LEO ∆i, w/ Overlap
Hybrid OCT Direct SOC
SP-1 Solution Convergence YES NO
SP-2 Solution Convergence YES YES
SP-3 Solution Convergence YES N/A
SP-1 # of Iterations 57 N/A
SP-2 # of Iterations 71 153
Total # of Iterations after SP-2 128 153
SP-1 Convergence Time 4.82 sec N/A
SP-2 Convergence Time 4.50 sec 15.16 sec
Total Convergence Time after SP-2 9.32 sec 15.16 sec
Min Prop Mass after SP-2 13.926 kg 14.013 kg
Table 10: Comparison of optimization parameters between Hybrid OCT and Direct SOC for 2-burn LEO ∆i
application problem w/ FB overlap. Green = superior solution, red = lesser solution. Note the faster
convergence of the Hybrid OCT solution in terms of iteration count and computer sim time, as well its
cheaper propellant mass solution.
residual. This specifies how accurately the nonlinear constraints are satisfied. The maximum
nonlinear constraint violation is normalized by the size of the solution, and is required to satisfy:
!
where is the violation of the i
th
nonlinear constraint, and is the major feasibility tolerance.
The default value of 1.0E-6 is appropriate when the linear and nonlinear constraints contain data
to about that accuracy
[64]
. This maximum violation is recorded as the feasibility metric.
The optimality metric measures the violation of the optimality tolerance (1.0E-02). This
metric specifies the final accuracy of the dual variables (Lagrange multipliers for the general
constraints) based on an estimate of the maximum complementarity slackness variable
(maxComp). This maxComp is computed from the final QP solution using the reduced gradients
derived from the objective gradient, the associated column of the constraint matrix, and the set of
the QP dual variables. Both of these convergence metrics approach zero near a solution (i.e. local
minima).
It can be seen that the infeasibility metric decreases nearly monotonically for both
methods, while the optimality metric periodically increases and decreases, with amplitude
trending towards zero (optimality). Keep in mind that the optimization algorithm SNOPT
prioritizes the feasibility of an iterate over the optimality when computing the next step. This is
done by selecting search directions formed by the linear constraints and the linearized constraint
computed from the Jacobian matrix of first derivatives of the objective function. Over the course
of major and minor iterations, the solution trends towards a point that satisfies the nonlinear
constraints and the first-order conditions for optimality. The effect results in a more chaotic
optimality violation than the max constraint residual. The final type of figure plots the value of
the objective function for each iterate, which for all problems considered is the minimum total
finite-burn time.
max
i
viol
i
/∥x∥≤ ϵ
r
viol
i
ϵ
r
!95
!96
0.00E+00
1.00E-03
2.00E-03
3.00E-03
4.00E-03
5.00E-03
0 5 10 15 20 25 30 35 40 45 50 55 60
Max Constraint Residual
Iteration
Hybrid OCT: SP-1 Feasibility
0.00E+00
1.00E-05
2.00E-05
3.00E-05
4.00E-05
5.00E-05
6.00E-05
7.00E-05
8.00E-05
0 5 10 15 20 25 30 35 40 45 50 55 60
Max Constraint Residual
Iteration
Feasibility Tolerance = 1.0E-06
0.00E+00
2.00E-01
4.00E-01
6.00E-01
8.00E-01
1.00E+00
1.20E+00
1.40E+00
1.60E+00
0 5 10 15 20 25 30 35 40 45 50 55 60
Optimality Tolerance Violation
Iteration
Hybrid OCT: SP-1 Optimality
Optimality Tolerance = 1.0E-02
Figure 29: Feasibility (top) and Optimality (bottom) metrics vs. iteration for Hybrid OCT: SP-1 solution.
!97
0.00E+00
2.00E-03
4.00E-03
6.00E-03
8.00E-03
1.00E-02
1.20E-02
0 10 20 30 40 50 60 70 80
Max Constraint Residual
Iteration
Hybrid OCT: SP-2 Feasibility
Feasibility Tolerance = 1.0E-06
0.00
0.50
1.00
1.50
2.00
2.50
3.00
3.50
4.00
4.50
5.00
0 10 20 30 40 50 60 70 80
Optimality Tolerance Violation
Iteration
Hybrid OCT: SP-2 Optimality
Optimality Tolerance = 1.0E-02
Figure 30: Feasibility (top) and Optimality (bottom) metrics vs. iteration for Hybrid OCT: SP-2 solution.
!98
0.14
0.15
0.16
0.17
0.18
0.19
0.20
0.21
0.22
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160
Total Burn Time (days)
Iteration
SP-2 Objective Function Comparison
SOC Obective
Hybrid OCT
0.00
0.01
0.02
0.03
0.04
0.05
0 20 40 60 80 100 120 140 160
Max Constraint Residual
Iteration
SOC: SP-2 Feasibility
Feasibility
Feasibility Tolerance = 1.0E-06
Figure 31: Comparison of objective function value vs. iteration for Hybrid OCT & Direct SOC: SP-2 solution.
Figure 32: Feasibility metric vs. iteration for Direct SOC: SP-2 solution.
7.3.3 Low-Thrust: Earth-Mars Hohmann-Type Transfer
The next application problem explored was a two-burn, non-coplanar, Earth—Mars low-
thrust transfer. Using a Lambert solver, two ∆V impulses were computed to travel from Earth’s
orbit to Mars’ orbit (ignoring planetary gravity). These ∆V’s may be thought of as ! vectors
departing and arriving at Earth and Mars, respectively, but in this problem they represent the
spacecraft departing with a ! of zero. The Optimization Sub-Problem 1 OCT solution is shown
in Figure 34 and Table 11, with comparisons to the SOC solution. It should be noted that the
initial velocity co-state magnitude was set experimentally to be ! , as this scaling
was found to greatly improve convergence behavior compared to the previously value of 1.0.
This is likely due to the larger geometric scale of this specific interplanetary problem, however
further study is required to better understand the relationship between the initial co-states.
V
∞
C
3
|λ
v
i
0
| = 2,000
!99
0.00
0.05
0.10
0.15
0.20
0.25
0 20 40 60 80 100 120 140 160
Tolerance Violation
Iteration
SOC: SP-2 Feasibility & Optimality
Feasibility
Optimality
Feasibility Tolerance = 1.0E-06
Optimality Tolerance = 1.0E-02
Figure 33: Feasibility & Optimality metrics vs. iteration for Direct SOC: SP-2 solution.
!100
Figure 34: Two-burn, non-coplanar, Earth—Mars low-thrust conversion, Optimization Sub-Problem 1.
Impulsive solution (orange), OCT solution (red), SOC solution (purple).
Impulsive Parameters
4,338.0 m/s
3,151.0 m/s
7,489.0 m/s
SOC Parameters OCT Parameters
Thrust 19.2 N 19.2 N
2,500 sec 2,500 sec
104.2 days 103.9 days
58.5 days 58.5 days
40,000 kg 40,000 kg
Propellant Used 11,006.7 kg 10,992.8 kg
Effective ∆V 7,889.8 m/s 7,878.0 m/s
!|
⃗ ΔV
1
|
!|
⃗ ΔV
2
|
!|
⃗ ΔV
Total
|
!I
sp
Finite-Burn !Δt
1
Finite-Burn !Δt
2
!M
i
o
Table 11: Two-burn, non-coplanar, Earth—Mars low-thrust Sub-Problem 1 solution burn parameters.
The Hybrid OCT solution to the final Optimization SP-3 is shown in Figure 35 and
Table 12, with comparisons to the SOC solution to SP-1. This was due to the surprisingly poor
convergence behavior of SOC in SP-2, yielding no converged solution despite OCT
demonstrating successful minimization. This further illustrates the benefit of the Hybrid
conversion method over Direct SOC, as the Hybrid had no greater difficulty in solving the larger,
more complex interplanetary problem than any of the simpler problems explored. Comparing the
Hybrid OCT SP-3 solution to the Direct SOC SP-1 solution reveals a significant cost savings of
2,320 kg of propellant, a 21.1% savings, when using Hybrid OCT over Direct SOC.
It is clear that on this large interplanetary problem scale, the savings in cost over the
course of the conversion algorithm is substantial. Comparing across a SP-1 OCT solution and an
SP-3 OCT solution yields a difference of 10,992.8 kg vs. 8,686.7 kg of used propellant,
respectively. Proof that the SP-3 solution satisfies the previously stated boundary constrains can
be seen by directly comparing the switching function across Sub-Problems. Figure 36 shows that
after the second round of optimization, the switching function for both burns does not have equal
end point values. After applying the final constraint, the third optimization does indeed yield
equal switching function values at the ignition and burn-out times for each converted finite-burn.
Note that the Sub-Problem 1 solution yields a switching function that is structurally close to that
of the third solution. This provides reassurance that the strategy of multi-impulse to finite-burn
conversion via pseudo-rendezvous is a suitable approach that aids in convergence of the optimal
control TPBVP.
!101
!102
Figure 35: Two-burn, non-coplanar, Earth—Mars low-thrust conversion, Optimization Sub-Problem 3.
Impulsive solution (orange), OCT solution (red), with SOC SP-1 solution (purple) for comparison.
Impulsive Parameters
4,338.0 m/s
3,151.0 m/s
7,489.0 m/s
SOC Parameters (SP-1) OCT Parameters (SP-3)
Thrust 19.2 N 19.2 N
2,500 sec 2,500 sec
104.2 days 87.367 days
58.5 days 41.015 days
40,000 kg 40,000 kg
Propellant Used 11,006.7 kg 8,686.7 kg
!|
⃗ ΔV
1
|
!|
⃗ ΔV
2
|
!|
⃗ ΔV
Total
|
!I
sp
Finite-Burn !Δt
1
Finite-Burn !Δt
2
!M
i
o
Table 12: Two-burn, non-coplanar, Earth—Mars low-thrust Sub-Problem 2 (SOC) and 3 (OCT) solution
burn parameters.
7.3.4 Low-Thrust: Modified Deep Space 1 Mission to Near-Earth Asteroid 9969 Braille
In order to further validate the capabilities of Hybrid OCT Method, an interplanetary
conversion problem was derived based closely on the flown trajectory of NASA’s Deep Space 1
(DS1), launched October 24th, 1998
[66-68]
. The technology demonstration mission was designed
to test high risk systems, namely an electrostatic ion thruster. During the mission, DS1’s
propulsion system was shown to provide 90 mN of thrust in vacuum, with a demonstrated
specific impulse of 3,100 seconds
[68]
. The launch mass of the vehicle was 486 kg, of which there
was 82 kg of Xenon propellant. DS1 was launched using a Boeing Delta II 7326 rocket, stacked
with a Star-37FM third-stage apogee kick motor (AKM). This third-stage provided the ∆V
necessary to escape from its LEO parking orbit and onto a heliocentric trajectory. The Star-37FM
stage provided 63 seconds of burn time at a thrust of 47.26 kN, a specific impulse of 289.8
seconds, a launch mass of 1,147.4 kg, and a burn-out mass of 81.5 kg for the stage.
!103
Figure 36: Switching function for the two-burn, non-coplanar, Earth—Mars low-thrust conversion;
Optimization Sub-Problem 2 Solution (top), Optimization Sub-Problem 3 Solution (bottom). The first
OCT burn (left), and second OCT burn (right).
The first destination of DS1 was a fly-by of the near-Earth asteroid (NEA) 9969 Braille.
During the first outbound interplanetary transfer leg, DS1 fired its ion thruster twice, first for
only 4.5 min on November 10th, 1998, and then for 14 days starting November 24th, 1998. DS1
made its closest approach to 9969 Braille on July 29, 1999 at a relative distance and velocity of
26 km and 15.5 km/s, respectively. For the purposes of this dissertation, two larger burns were
designed for this outbound leg to better exercise the hybrid conversion method, while
maintaining a comparable fly-by geometry of 9969 Braille. Note that the goal was not to
replicate the DS1 trajectory, but to develop a test problem comparable to a flown mission.
Taking inspiration from the DS1 trajectory, a set of Earth-escape orbital parameters was
computed using the spacecraft design of DS1, and making use of the Star-37FM AKM. Two
impulsive ∆Vs of equal magnitude (based on achievable propellant expenditure during the
outbound time interval) were then set for November 10, 1998, and March 22, 1999. The
spherical angles , dictating the ∆V and V∞ directions in heliocentric space, were then
optimized to target a fly-by geometry close to that of the flown DS1 trajectory. The parameters of
this computed impulsive trajectory are shown below.
Earth Escape !
!
!
where ! are expressed in the Sun J2000 VUW frame.
!
!
!
!
!
where ! are expressed in the Sun J2000 IJK frame.
α & β
V
∞
= 2.9489 km/s
α
∞
=−14.8047 deg
β
∞
=−62.9735 deg
α
∞
& β
∞
ΔV
1
=ΔV
2
= 1.650806 km/s
α
1
= 103.6443 deg
β
1
= 3.1642 deg
α
2
= 97.5617 deg
β
2
= 102.7469 deg
α
1,2
& β
1,2
!104
This two-burn impulsive trajectory can be seen in Figure 37, along with the orbit of
asteroid 9969 Braille. The approximation of the modified DS1 trajectory was then input to the
Hybrid Conversion Algorithm, and separately solved using the traditional Direct SOC Method.
The SP-1 Hybrid OCT solution is shown in Figure 38, while the SP-2 solution compared to the
Direct SOC solution is shown in Figure 39 and Table 13.
It should be noted that the modified DS1 impulsive trajectory was designed to purposely
induce a FB overlap to further exercise the conversion method, as seen in Figure 38. The overlap
between burns 1 and 2 is found to be 39.069 days, significantly longer than the overlap in
application problem 2. For comparison, Figure 40 shows the first leg of the flown DS1 trajectory
(note the former name of asteroid 9669 Braille, asteroid 1992 KD). After executing two burns
with its electrostatic thruster during this outbound segment (4.5 min and 14 days, respectively),
this equates to roughly 3.581 kg of propellant used, compared to over 50 kg in the computed
trajectories. Recall that the purpose of this work was not to replicate the DS1 trajectory, but to
use it as a foundation to validate the conversion algorithm with a realistic, flown spacecraft. This
required increasing the burn time to make the conversion process more significant. That being
said, the modified trajectory is still comparable in terms of the mission objectives considered for
this study (targeting the NEA fly-by), and does not exceed the flown propellant mass budget.
After SP-2 conversion, the difference in used propellant shows a 0.12% improvement when
using the Hybrid OCT Method over Direct SOC. The numerical validity of this will be
investigated in the next section.
Before that critical discussion, however, it must first be emphasized that the Direct SOC
solution failed to converge for SP-1 and SP-2, just as Direct SOC failed within the SP-2 Earth-
Mars application problem. In both cases, this is due to the inability of the SOC control law
(linearly steered) to model the thrust direction flip shown to be associated with low-thrust FB
overlaps for burns spanning a significant portion of a revolution. The only way the modified DS1
trajectory conversion was able to be solved using Direct SOC (in SP-2) was to use the converged
Hybrid OCT SP-2 solution in a reverse-adjoint control transformation (i.e. use the Hybrid OCT
solution as an initial guess for the Direct SOC solution). This re-seeded, “assisted” Direct SOC
!105
!106
Figure 37: Impulsive Approximation of Modified DS1 Trajectory, with increased ∆V . Earth, Mars, and
asteroid 9969 Braille orbits are shown in blue, red, and white, respectively.
∆V1
∆V2
∆V∞
∆V1
∆V2
∆V∞
Asteroid
9969 Braille
Asteroid
9969 Braille
DS1 Flyby
!107
Figure 38: Modified DS1 Trajectory, SP-1 Hybrid OCT Solution. Note the FB overlap of 39.069 days.
OCT FB 1
∆V∞
Asteroid
9969 Braille
DS1 Flyby
OCT FB 2
FB Overlap
]
FB Overlap
]
OCT FB 1
Asteroid
9969 Braille
DS1 Flyby
OCT FB 2
!108
Figure 39: Modified DS1 Trajectory, SP-2 Hybrid OCT Solution (red & green), and assisted SP-2 Direct SOC
Solution (purple).
Impulsive Parameters
1,650.806 m/s
1,650.806 m/s
3,301.612 m/s
SOC Parameters (SP-2) OCT Parameters (SP-2)
Thrust 90 mN 90 mN
3,100 sec 3,100 sec
61.113 days 60.574 days
138.052 days 138.357 days
486 kg 486 kg
Propellant Used 50.943 kg 50.883 kg
!|
⃗ ΔV
1
|
!|
⃗ ΔV
2
|
!|
⃗ ΔV
Total
|
!I
sp
Finite-Burn !Δt
1
Finite-Burn !Δt
2
!M
i
o
Table 13: Modified DS1 Trajectory, SP-2 Hybrid OCT and Direct SOC Solution Burn Parameters.
solution is presented in Table 13, but it should be noted that such a solution would not be
available were it not for the successful performance of the Hybrid OCT Method. This means that
for all practical considerations, the Direct SOC Method completely failed to solve the modified
DS1 conversion problem, while the Hybrid OCT Method showed no issue, once again
demonstrating its superiority for the classes of problems explored.
7.4 Proof of Numerically Significant Improvement
Through the four presented application problems, it can be seen that using the Hybrid
OCT Method instead of the traditional Direct SOC Method yields a roughly consistent objective
cost savings of ~0.1% for all problems explored. The one exception to this is in the Earth-Mars
transfer problem, where the Direct SOC Method failed to converge on a solution to SP-2,
resulting in its SP-1 solution netting a ~21% cost increase compared to the significantly cheaper
Hybrid OCT SP-3 solution. In order for the generally small improvement of the Hybrid OCT
Method to be considered numerically significant, i.e. not a result of numerical error or noise in
the simulation, focus must be given to the integration error tolerances, and the optimization
convergence tolerances.
!109
Figure 40: Flown Deep Space 1 Trajectory to Asteroid 1992 KD (now called 9969 Braille)
[66]
.
First, as described in the Section 7.1, a variable order, variable step-size Adams-Basforth-
Moulton integration method
[61]
was utilized for propagation, with absolute and relative error
tolerances of 1.0E-12. The Simulink user’s guide
[69]
provides a clear description of what these
standard integration error tolerances represent for such solvers:
“Variable-step solvers use standard control techniques to monitor the local
error at each time step. During each time step, the solvers compute the state
values at the end of the step and determine the local error —the estimated error
of these state values. They then compare the local error to the acceptable error,
which is a function of both the relative tolerance (rtol) and the absolute
tolerance (atol). If the local error is greater than the acceptable error for any
one state, the solver reduces the step size and tries again.
…Relative tolerance measures the error relative to the size of each state. The
relative tolerance represents a percentage of the state value. The default, 1e-3,
means that the computed state is accurate to within 0.1%.
Absolute tolerance is a threshold error value. This tolerance represents the
acceptable error as the value of the measured state approaches zero.
The solvers require the error for the i
th
state, ei, to satisfy:
ei ≤ max(rtol × |xi|, atoli)”
This tells us that by setting these integration error tolerances to 1E-12, the computed state
is accurate to within 0.0000000001%. As this error tolerance is nine orders of magnitude smaller
than ~0.1%, the difference in objective cost between the Hybrid OCT Method and Direct SOC
Method cannot be attributed to integration error or noise in the propagated trajectory.
The next pieces needed to prove numerical significance are the optimization convergence
tolerances, namely, we wish to verify that the converged solutions have acceptable accuracy
relative to the observed difference in cost between solution types. Recall that the optimization
package SNOPT utilizes an SQP algorithm with two key tolerances. First, a feasibility metric:
the measure of how infeasible the constraint vector is at each iterate, or in other words, the
maximum constraint residual. This specifies how accurately the nonlinear constraints are
!110
satisfied. The maximum nonlinear constraint violation is normalized by the size of the solution,
and is required to satisfy:
!
where is the violation of the i
th
nonlinear constraint, and is the major feasibility tolerance.
The default value of 1.0E-06 is appropriate when the linear and nonlinear constraints contain
data to about that accuracy
[64]
. This violation is recorded as the feasibility metric for each iterate.
The second metric is the optimality tolerance (set to 1.0E-03). This metric specifies the
final accuracy of the dual variables (Lagrange multipliers for the general constraints) based on an
estimate of the maximum complementarity slackness variable (maxComp). This maxComp is
computed from the final QP solution using the reduced gradients derived from the objective
gradient, the associated column of the constraint matrix, and the set of the QP dual variables.
Both of these convergence metrics approach zero near a solution (i.e. local minima).
What this means, in essence, is that we must compare the feasibility and optimality
metrics of the converged solutions with the size of the observed difference in objective cost.
Table 14 presents the feasibility and optimality metrics for the final iterate of each of the four
application problems (i.e. the converged solutions for each). It can be clearly seen that without
exception, all converged solutions boast feasibility and optimality violations that are several
orders of magnitude smaller than the observed difference in objective cost (0.1%), indicating
sufficient convergence. This fact combined with the demonstrated numerical integration
accuracy puts a nail in the coffin for the notion that the difference in cost between solutions is a
result of numerical error or noise. It can now be confidently stated that the Hybrid OCT Method
yields a numerically significant improvement in objective cost over the Direct SOC Method for
all conversion problems explored. This benefit is in addition to the aforementioned 38.4%
improvement in convergence time for conversion problems where both ! , and
the burn is changing inclination.
max
i
viol
i
/∥x∥≤ ϵ
r
viol
i
ϵ
r
ΔT
burn
>
T
Period
2
!111
!112
Optimization Convergence Metrics for all Application Problems
Apogee Raise Feasibility Optimality
Hybrid OCT Method 9.2E-07 3.3E-05
Direct SOC Method 6.6E-07 1.5E-05
LEO ∆i, w/ FB Overlap Feasibility Optimality
Hybrid OCT Method 2.1E-07 7.3E-04
Direct SOC Method 1.9E-08 5.6E-05
Earth-Mars Transfer Feasibility Optimality
Hybrid OCT Method 3.1E-10 1.8E-06
Direct SOC Method 1.5E-08 4.0E-11
Modified Deep Space 1 Feasibility Optimality
Hybrid OCT Method 1.0E-06 8.8E-07
Direct SOC Method 1.9E-08 1.1E-06
Table 14: Comparison of optimization convergence violations via Hybrid OCT and Direct SOC Methods for
all four application problems: apogee raise, LEO ∆i w/ FB overlap, Earth-Mars Transfer, and Modified
DS1. Note that for all problems, feasibility tolerance = 1.0E-06, and optimality tolerance = 1.0E-03.
8. SUMMARY AND CONCLUSION
A novel finite-burn trajectory optimization algorithm has been presented that offers
objective cost savings, and optimization convergence benefits for certain class of problems, over
conventional methods. The procedure combines the high fidelity control law of an indirect
optimization method with the generality of direct methods to improve upon solving time-
optimal, finite-burn (FB) Pseudo-Rendezvous (P-R) spacecraft trajectory problems. The
traditional difficulty in approximating the initial co-state vector within the optimal control Two-
Point Boundary Value Problem (TPBVP) is mitigated by a unique series of adjoint-control
transformations. These mappings evolve the constraint vector across multiple sub-problems (SP),
culminating in the desired problem definition. This provides an increasingly more accurate
approximation of the initial co-state vector, which by the end of the process become sufficiently
close to the optimal values, enabling solution convergence for the desired P-R problem.
Focus was given to solving the FB Arc Overlap problem, which can arise when patching
together a series of single-impulse conversions with insufficiently spaced ΔV magnitudes,
relative to engine performance. The FB overlap problem results in the ballistic arc connecting
two converged finite-burns having a negative time of flight, yielding an infeasible trajectory. The
conversion algorithm mitigates this through a process of introducing additional boundary
constraints, and by relaxing the P-R state targeting on all interior conversions, leaving the state
constraints on only the final impulse conversion. This process also has the benefit of improving
optimality for the entire sequence of burns, as by definition each FB was first optimized
individually, not as a member of the larger set of conversions, which is the true problem at hand.
By evolving the targeting constraints and growing the problem to encompass the entire set,
regardless if any overlaps exist, the optimizer is afforded greater freedom when minimizing the
full multi-FB problem. This is because the burn-out states of each FB are no longer constrained
to be on the inherently sub-optimal impulsive input trajectory.
The conversion method requires no prior knowledge of the solution, merely a multi-
impulse or sub-optimal finite-burn initial guess, which aids automation. Numerical examples
were shown for a variety of 3D orbit transfers, illustrating the adaptability of the method. These
!113
include a co-planar apogee raise maneuver, low-thrust non-coplanar Earth orbit transfers with
and without FB arc overlaps, a low-thrust Earth-Mars transfer, and a modified low-thrust
interplanetary trajectory based upon the flown Deep Space 1 mission. Hybrid conversion results
are shown with comparisons to standard sub-optimal control (SOC) solutions. The Hybrid OCT
solutions offer faster convergence under certain circumstances, and in all cases provide
moderately cheaper propellant cost when compared to SOC solutions. Simulations were run
using common industry tools, demonstrating the potential new capability that can be added to
conventional mission design practice by applying the presented method.
The prime goal of this dissertation was to advance the state-of-the-art by improving upon
the traditional Direct Method used to solve the multi-impulse to time-optimal FB conversion
problem. The presented results demonstrate that the new Hybrid Method yields cheaper solutions
over the traditional Direct Method for every problem examined, a roughly 0.1% improvement.
The numerical significance of this was demonstrated by examining the numerical integration
error and optimization convergence tolerances. First, the numerical integration error tolerance
was shown to be nine orders of magnitude smaller than the observed difference in objective cost
between solutions, ruling out the notion that it can be attributed to integration error or noise in
the propagated trajectory. Next, the feasibility and optimality violations of all converged
solutions were shown to be several orders of magnitude smaller than the observed difference in
cost, indicating sufficient convergence and finally confirming numerical significance. It can now
be confidently stated that the Hybrid OCT Method yields a numerically significant improvement
in objective cost over the Direct SOC Method for all conversion problems explored. In addition
to this wide-spread benefit, explorations into the FB-arc overlap problem indicate that the Hybrid
OCT Method will yield superior convergence behavior (i.e. number of iterations and
convergence time) for conversion problems where both ! , and inclination is
changed during the burn. This specific class of problem will therefore reap larger benefits than
other problems when the Hybrid OCT Method is employed in lieu of Direct SOC.
Having advanced the state of the art with regards to hybrid finite-burn trajectory
optimization, the work described in this dissertation can be continued in a multitude of
ΔT
burn
>
T
Period
2
!114
directions. The first step would be to apply the hybrid conversion algorithm to more complex
transfer problems with a larger number of impulses than those presented in Section 7.3. A perfect
example of this is a multi-impulse representation of a low-thrust spiral with many revolutions,
such as those described in Section 6. This represents a key transfer problem found in the realm
of low-thrust and solar electric propulsion mission design, and one that would benefit
substantially from the robust convergence properties of the hybrid conversion algorithm. As low-
thrust spirals are often employed to place commercial satellites in their operational orbits, a
savings in objective cost (transfer time and propellant) could represent a significant savings in
financial cost as well. This would be thanks to the satellite’s reduction in launch mass, and more
importantly reaching its operational orbit sooner than if a sub-optimal spiral trajectory was
employed.
Another way this research can be advanced by future work is to focus on automation and
implementation. First, effort could be given to derive and compute the analytic gradients of the
objective and constraints, as opposed to leaving them to the optimization routine to estimate.
This can be done by applying linear perturbation theory, namely the computation of the State
Transition Matrix. This has the potential to save considerable computer processing time
[51]
.
Second, due to the inherently parallelizable nature of the first conversion sub-problem
(individually optimizing a set of smaller P-R problems), significant improvements to
computational processing and convergence times may be achieved by leveraging subroutines
such as the Modified-Picard Chebyshev Method (MPCM) for propagation. Since the hybrid
conversion algorithm computes trajectories comprised of discretized, sub-revolution sized
impulsive and finite-burns, it is a perfect application for the MPCM, which performs best on
partial orbital segments as opposed to many revolutions
[25]
. Previous efforts demonstrated the
ability to improve runtime performance of SNOPT by over a factor of three when using the
MPCM over MATLAB’s ODE113 for numerical propagation
[1]
. Such a subroutine combined
with the cost savings and case-dependent convergence improvements afforded by the hybrid
conversion algorithm will likely yield substantial overall benefits compared to using standard
industry optimization tools for FB problems. !115
BIBLIOGRAPHY
[1] Koblick, D., Xu, S., Fogel, J., Shankar, P., “Low Thrust Minimum Time Orbit Transfer
Nonlinear Optimization Using Multi-Impulse Discretization via the Modified Picard-
Chebyshev Method”, International Conference on Computational & Experimental
Engineering and Sciences, 2015.
[2] Alemany, K., Braun, R., “Survey of Global Optimization Methods for Low-Thrust, Multiple
Asteroid Tour Missions”, AAS 07-211.
[3] Williams, S.N., Coverstone-Carrol, V ., “Benefits of Solar Electric Propulsion for the Next
Generation of Planetary Exploration Missions,” The Journal of the Astronautical Sciences,
V ol. 45, No. 2, April-June 1997, pp. 143-159.
[4] Kechichian, J., “Optimal Low-Earth-Orbit-Geostationary-Earth-Orbit Intermediate
Acceleration Orbit Transfer,” AIAA Journal of Guidance, Control, and Dynamics, 1997, pp.
803-811.
[5] Thorne, J., “Minimum-Time, Constant-Thrust Orbit Transfers with Noncircular Boundary
Conditions,” Institute for Defense Analyses, Alexandria, V A.
[6] Petropoulos, A.E., Longuski, J.M., “Shape-Based Algorithm for Automated Design of Low-
Thrust, Gravity-Assist Trajectory,” Journal of Spacecraft and Rockets, V ol. 41, No. 5,
September-October 2004, pp 787-796.
[7] Patel, P.R., Gallimore, A.D., Zurbuchen, T.H., Scheeres, D.J., “An Algorithm for Generating
Feasible Low Thrust Interplanetary Trajectories,” IEPC-2009-219, September 2009.
[8] Carnelli, I., Dachwald, B., Vasile, M., “Evolutionary Neurocontrol: A Novel Method for
Low-Thrust Gravity-Assist Trajectory Optimization,” Journal of Guidance, Control, and
Dynamics, V ol. 32, No. 2, March-April 2009, pp. 615-624.
[9] Vallado, D., McClain, W., “Fundamentals of Astrodynamics and Applications”, Managinf
Forest Ecosystems, Springer, 2001.
[10] V on Stryk, O., Bulirsch, R., “Direct and Indirect Methods for Trajectory Optimization”,
Annals of Operations Research, 1992, vol. 37, no. 1, pp. 357-373.
!116
[11] Patterson, M.A., Rao, A.V ., “GPOPS-II: A MATLAB Software for Solving Multiple-Phase
Optimal Control Problems Using hp-Adaptive Gaussian Quadrature Collocation Methods
and Sparse Nonlinear Programming,” ACM Transactions on Mathematical Software, V ol.
41, No. 1, Article 1, October 2014.
[12] Whiffen, G.J., and Sims, J., “Mystic: Implementation of the Static Dynamic Optimal
Control Algorithm for High-Fidelity, Low-Thrust Trajectory Design,” AIAA Paper
2006-6741, August 2006.
[13] Betts, J.T., “Survey of Numerical Methods for Trajectory Optimization,” Journal of
Guidance, Control, and Dynamics, V ol. 21, No. 2, March-April 1998, pp. 193-207.
[14] Hargraves, C.R., Paris, S.W., “Direct Trajectory Optimization Using Nonlinear
Programming and Collocation,” Journal of Guidance, Control, and Dynamics, V ol. 10, No.
4, 1987, 338-342.
[15] Tang, S., Conway, B.A., “Optimization of Low-Thrust Interplanetary Trajectories Using
Collocation and Nonlinear Programming,” Journal of Guidance, Control, and Dynamics,
V ol. 18, No. 3, 1995, pp. 599-604.
[16] Nah, R.S., Vadali, S.R., “Fuel-Optimal, Low-Thrust, Three-Dimensional Earth-Mars
Trajectories,” Journal of Guidance, Control, and Dynamics, V ol. 24, No. 6, Nov-Dec 2001.
[17] Sims, J.A., Finlayson, P.A., Rinderle, E.A., Vavrina, M.A., Kowalkowski, T.D.,
“Implementation of a Low-Thrust Trajectory Optimization Algorithm for Preliminary
Design,” AIAA/AAS Astrodynamics Specialist Conference and Exhibit, August 2006.
[18] McConaghy, T.T., Debban, T.J., Petropoulos, A.E., Longuski, J.M., “Design and
Optimization of Low-Thrust Trajectories with Gravity Assists,” Journal of Spacecraft and
Rockets, V ol. 40, No. 3, May-June 2003, pp. 380-387.
[19] Anderson, R.L., Lo, M.W., “Role of Invariant Manifolds in Low-Thrust Trajectory Design,”
Journal of Guidance, Control, and Dynamics, V ol. 32, No. 6, November-December 2009.
[20] Mingotti, G., F. Topputo, Bernelli-Zazzera, F., “Earth-Mars transfers with ballistic escape
and low-thrust capture,” Celest. Mech. Dyn. Astr., May 2011, pp 169-188.
!117
[21] Zhang, C., Topputo, F., Bernelli-Zazzera, F., Zhao, Y .S., “Low-Thrust Minimum-Fuel
Optimization in the Circular Restricted Three-Body Problem,” Journal of Guidance,
Control, and Dynamics, V ol. 38, No. 8, August 2015.
[22] Hargraves, C.R., Paris, S.W., “Direct Trajectory Optimization Using Nonlinear
Programming and Collocation,” J. Guidance., V ol. 10, No. 4, 1987.
[23] Chobotov, V ., “Orbital Mechanics,” AIAA Education Series, American Institute of
Aeronautics and Astronautics, 1996.
[24] V on Stryk, O., Bulirsch, R., “Direct and Indirect Methods for Trajectory Optimization,”
Annals of Operations Research, V ol. 37, No. 1, pp. 357-373, 1992.
[25] Koblick, D., Shankar, P., “Evaluation of the Modified Picard-Chebyshev Method for High-
Precision Orbit Propagation,” Journal of Aerospace Engineering, 2014.
[26] Gill, P.E., Murray, W., Saunders, M.A., “User’s Guide for SNOPT Version 7: Software for
Large-Scale Nonlinear Programming,” June 16, 2008.
[27] Parcher, D.W., Sims, J.A., “Gravity-Assist Trajectories to Jupiter Using Nuclear Electric
Propulsion,” AAS/AIAA Astrodynamics Specialist Conference, AAS 05-398, Lake Tahoe,
California, August 7-11, 2005.
[28] Kowalkowski, T.D., Kangas, J.A., Parcher, D.W., “Jupiter Icy Moons Orbiter Interplanetary
Injection Period Analysis,” AAS/AIAA Space Flight Mechanics Meeting, AAS 06-187,
Tampa, Florida, January 22-26, 2006.
[29] Whiffen, G.J., “An Investigation of a Jupiter Galilean Moon Orbiter Trajectory,” AAS/
AIAA Astrodynamics Specialists Conference, AAS 03-544, Big Sky, Montana, August 3-7,
2003.
[30] Whiffen, G.J., Lam, T., “Jupiter Icy moons Orbiter Reference Trajectory,” AAS/AIAA
Space Flight Mechanics Conference, AAS 06-186, Tampa, Florida, January 22-26, 2006.
[31] Bai, X., Junkins, J.L., “Modified Chebyshev-Picard Iteration Methods for Orbit
Propagation,” The Journal of the Astronautical Sciences, 2011, V ol. 58, No. 4, pp. 583-613.
!118
[32] Bai, X., Junkins, J.L., “Modified Chebyshev-Picard Iteration Methods for Solution of
Boundary Value Problems,” The Journal of the Astronautical Sciences, 2011, V ol. 58, No. 4,
pp. 615-644.
[33] Bai, X., Junkins, J.L., “Modified Chebyshev-Picard Iteration Methods for Solution of Initial
Value Problems,” The Journal of the Astronautical Sciences, 2012, V ol. 59, No. 1-2, pp.
327-351.
[34] Bai, X., Junkins, J.L., “Modified Chebyshev-Picard Iteration Methods for Solution for
Station-Keeping of Translunar Halo Orbits,” Mathematical Problems in Engineering, V ol.
2012.
[35] Edelbaum, T.N., “Propulsion Requirements for Controllable Satellites,” American Rocket
Society Journal, 1961, V ol. 31, pp. 1079-1089.
[36] Kechichian, J., “Reformulation of Edelbaum’s Low-Thrust Transfer Problem using Optimal
Control Theory,” Journal of Guidance, Control, and Dynamics, 1997.
[37] Eagle, D., “Low Thrust Transfer Between Non-Coplanar Circular Orbits,” Orbital
mechanics with MATLAB.
[38] Caillau, J.B., Gergaud, J., Noailles, J., “TfMin: Short Reference Manual,” ENSEEIHT-IRIT
Technical Report RT/APO/01/3 July 2001.
[39] Garza, D., “Hybrid Trajectory Optimization Program,” “Indirect Trajectory Optimization
Program,” Aerospace Corp.
[40] Fogel, J., Koblick, D., Shah, N., “Survey of Low-Thrust Gravity Assist Trajectory
Optimization Methods, with Comparisons to a Novel, Multi-Impulse Discretization
Approach”, International Astronautical Congress, Oct 2015.
[41] Longuski, M., Guzmán, J., Prussing, J., Optimal Control with Aerospace Applications,
Space Technology Library, Microcosm Press and Springer Science, 2014.
[42] Bryson, A., Ho, Y ., Applied Optimal Control, Hemisphere Publishing 1975.
[43] Hull, D., Optimal Control Theory for Applications, Springer Verlag, 2003.
!119
[44] Jezewski, D., “Primer Vector Theory and Applications,” NASA TRR-454, Lyndon B.
Johnson Space Center, Houston, TX, Nov. 1975
[45] Ocampo, C., Senent, J., Williams, J., “Theoretical Foundation of Copernicus: A Unified
System for Trajectory Design and Optimization”, 4th International Conference on
Astrodynamics Tools and Techniques, 3-6 May 2010, Madrid, Spain.
[46] “Apollo 13 Technical Air-to-Ground V oice Transcription”, Test Division, Apollo Spacecraft
Program Office, National Aeronautics and Space Administration, Manned Spacecraft
Center, Houston, Texas, April 1970.
[47] Pines, S., “Constant of Motion for Optimum Thrust Trajectories in a Central Force Field”,
AIAA Journal, V ol. 2, No. 11, 1964.
[48] Handelsman, M., “Optimal Free-Space Fixed-Thrust Trajectories Using Impulsive
Trajectories as Starting Iterates”, AIAA Journal, V ol. 4, No. 6, 1966.
[49] Hazelrigg, G., Lion, P., “Analytical Determination of the Adjoint Vector for Optimum Space
Trajectories”, Journal of Spacecraft and Rockets, V ol. 7, No. 10, 1970.
[50] Komhauser, A., Lion, P., Hazelrigg, G., “An Analytic Solution for Constant-Thrust,
Optimal-Coast, Minimum-Propellant Space Trajectories”, AIAA Journal, V ol. 9, No. 7,
1971.
[51] Jesick, M., Ocampo, C., “Automated Design of Optimal Finite Thrust Orbit Insertion with
Ballistic Flyby Constraints”, Journal of Spacecraft and Rockets, V ol. 51, No. 6, 2014.
[52] Bai, X., Junkins, J., “Solving initial value problems by the Picard-Chebyshev method with
nvidia GPUs”, Proceedings of the 20th Spaceflight Mechanics Meeting, San Diego, CA,
AAS 10-197, 2010.
[53] Koblick, D., Poole, M., Shankar, P., “Parallel high-precision orbit propagation using the
modified picard-chebyshev method”, ASME International Mechanical Engineering
Congress and Exposition, Houston, TX, November 9-15, 2012.
[54] Macomber, B., Probe, A., Woollands, R., Junkins, J., “Parallel modified-Chebyshev Picard
iteration for orbit catalog propagation and monte-carlo analysis”, Proceedings of the 38th
Annual AAS/AIAA Guidance and Control Conference, Breckenridge, CO, AAS 15-009,
2015.
!120
[55] Read, J., “Modified Chebyshev Picard Iteration: Integration of Perturbed Motion Using
Modified Equinoctial Elements”, PhD Dissertation, Texas A&M University, 2016.
[56] Bertrand, R., Epenor, R., “New smoothing techniques for solving bang-bang optimal control
problems—numerical results and statistical interpretation”, Optimal Control Applied
Methods, vol. 23, 171-197, 2002.
[57] Russell, R., “Primer Vector Theory Applied to Global Low-Thrust Trade Studies”, Journal
of Guidance, Control, and Dynamics, V ol. 30, No. 2, 460-472, 2007.
[58] Thorne, J., Hall, C., “Approximate Initial Lagrange Costates for Continuous-Thrust
Spacecraft”, Journal of Guidance, Control, and Dynamics, V ol. 19, No. 2, 1996, pp.
283-288.
[59] Dixon, L., Biggs, M., “The Advantages of Adjoint-Control Transformations when
Determining Optimal Trajectories by Pontryagin’s Maximum Principle”, Aeronautical
Journal, V ol. 76, No. 735, 1972, pp. 169-174.
[60] VanderVeen, A., “Triple-Planet Ballistic Flybys of Mars and Venus”, Journal of Spacecraft
and Rockets, V ol. 6, No. 4, 1969, pp. 383-389.
[61] Shampine, L.F., Watts, H.A., “DEPAC - Design of a User Oriented Package of ODE
Solvers”, Technical Report SAND-79-2374, Sandia National Labs, September 1980.
[62] Englander, J.A., Englander, A.C., “Tuning Monotonic Basin Hopping: Improving the
Efficiency of Stochastic Search as Applied to Low-Thrust Trajectory Optimization”,
Technical Report GSF-E-DAA-TN14154, NASA Technical Report Server, July 30, 2014.
[63] Williams, J, “Copernicus 4.6 User Guide”, NASA Johnson Space Center, Aeroscience and
Flight Mechanics Division, Flight Mechanics and Trajectory Design Branch (EG5), April
30, 2018.
[64] Gill, P.E., Murray, W., Saunders, M.A., “SNOPT: An SQP Algorithm for Large-Scale
Constrained Optimization,” Society for Industrial and Applied Mathematics, V ol. 47, No. 1,
2005.
[65] Connolly, J., “Deep Space Transport (DST) and Mars Mission Architecture,” NASA Mars
Study Capability Team, October 17, 2017.
!121
[66] “Deep Space 1 Launch Press Kit October 1998,” National Aeronautics and Space
Administration, October 1998.
[67] Rayman, M., Chadbourne, P., Culwell, J., Williams, S., “Mission Design for Deep Space 1:
A Low-Thrust Technology Validation Mission,” Acta Astronautica, V ol. 45, Nos. 4-9, 1999.
[68] Siddiqi, A., “Beyond Earth: A Chronicle of Deep Space Exploration, 1958-2016,” United
States, NASA History Division, NASA SP-2018-4041, V ol. 45, Nos. 4-9, 1999.
[69] “Simulink User’s Guide” MathWorks R2019b, The Math Works, Inc., 2019.
!122

Abstract (if available)

Abstract A novel algorithm is presented that provides an improvement over traditional parameter optimization methods when solving time-optimal, finite-burn pseudo-rendezvous spacecraft trajectory problems. A hybrid optimization procedure is described that converts a set of multiple-impulses, representing high- or low-thrust maneuvers, to an exact time-optimal finite-burn trajectory for a thrust limited, constant exhaust velocity spacecraft. The Hybrid Method applies a control law derived from the Euler-Lagrange system of equations within the classical Indirect Method to a modern Direct Method. An iterative adjoint-control transformation and an evolving constraint vector are introduced to solve the optimal control two-point boundary value problem. This method requires no prior knowledge of the solution, which adds simplicity to the trajectory design process and aids automation. Examples are shown for low-thrust apogee raise maneuvers, non-coplanar Earth orbit transfers, a two-burn Earth-Mars transfer, and a modified Deep Space 1 low-thrust trajectory. A numerically significant improvement to objective cost is shown across all application problems compared to a traditional solution method, as well as a significant improvement to convergence speed for a select class of problems.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Mission design study of an RTG powered, ion engine equipped interstellar spacecraft

PDF

Relative-motion trajectory generation and maintenance for multi-spacecraft swarms

PDF

Trajectory mission design and navigation for a space weather forecast

PDF

Optimal guidance trajectories for proximity maneuvering and close approach with a tumbling resident space object under high fidelity J₂ and quadratic drag perturbation model

PDF

The space environment near the Sun and its effects on Parker Solar Probe spacecraft and FIELDS instrumentation

PDF

Optimization of relative motion rendezvous via modified particle swarm and primer vector theory

PDF

Periodic motion analysis in spacecraft attitude dynamics

PDF

Techniques for analysis and design of temporary capture and resonant motion in astrodynamics

PDF

On practical network optimization: convergence, finite buffers, and load balancing

PDF

Exploitation of sparse and low-rank structures for tracking and channel estimation

PDF

Iterative path integral stochastic optimal control: theory and applications to motor control

PDF

Application of the fundamental equation to celestial mechanics and astrodynamics

PDF

Kinetic studies of collisionless mesothermal plasma flow dynamics

PDF

Cooperative localization of a compact spacecraft group using computer vision

PDF

Imposing classical symmetries on quantum operators with applications to optimization

PDF

Trajectory planning for manipulators performing complex tasks

PDF

Landscape analysis and algorithms for large scale non-convex optimization

PDF

Aerodynamics at low Re: separation, reattachment, and control

PDF

Hydrogen peroxide vapor for small satellite propulsion

PDF

Learning and control for wireless networks via graph signal processing

Asset Metadata

Creator Fogel, Joshua (author)

Core Title Spacecraft trajectory optimization: multiple-impulse to time-optimal finite-burn conversion

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Astronautical Engineering

Publication Date 11/13/2019

Defense Date 10/10/2019

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag low-thrust,OAI-PMH Harvest,optimal control theory,optimization,spacecraft,trajectory

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Gruntman, Mike (committee chair), Flashner, Henryk (committee member), Wang, Joseph (committee member)

Creator Email fogelfever@gmail.com,joshuafo@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-234177

Unique identifier UC11675602

Identifier etd-FogelJoshu-7919.pdf (filename),usctheses-c89-234177 (legacy record id)

Legacy Identifier etd-FogelJoshu-7919.pdf

Dmrecord 234177

Document Type Dissertation

Rights Fogel, Joshua

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA