813 lines
45 KiB
TeX
813 lines
45 KiB
TeX
\documentclass[defaultstyle,11pt]{thesis}
|
||
|
||
\usepackage{graphicx}
|
||
\usepackage{amssymb}
|
||
\usepackage{hyperref}
|
||
\usepackage{amsmath}
|
||
|
||
\title{Designing Optimal Low-Thrust Interplanetary Trajectories Utilizing Monotonic Basin Hopping}
|
||
\author{Richard C.}{Johnstone}
|
||
\otherdegrees{B.S., Unviersity of Kentucky, Mechanical Engineering, 2016 \\
|
||
B.S., University of Kentucky, Physics, 2016}
|
||
\degree{Master of Science}{M.S., Aerospace Engineering}
|
||
\dept{Department of}{Aerospace Engineering}
|
||
\advisor{Prof.}{Natasha Bosanac}
|
||
\reader{TBD: Kathryn Davis}
|
||
\readerThree{TBD: Daniel Scheeres}
|
||
|
||
\abstract{ \OnePageChapter
|
||
There are a variety of approaches to finding and optimizing low-thrust trajectories in
|
||
interplanetary space. This thesis analyzes one such approach, Sims-Flanagan transcriptions, and
|
||
its applications in a multiple-shooting non-linear solver for the purpose of finding valid
|
||
low-thrust trajectory arcs between planets given poor initial conditions. These valid arcs are
|
||
then fed into a Monotonic Basin Hopping (MBH) algorithm, which combines these arcs in order to
|
||
find and optimize interplanetary trajectories, given a set of flyby planets. This allows for a
|
||
fairly rapid searching of a very large solution space of low-thrust profiles via a medium
|
||
fidelity inner-loop solver and a well-suited optimization routine. The trajectories found by
|
||
this method can then be optimized further by feeding the solutions back, once again, into the
|
||
non-linear solver, this time allowing the solver to perform optimization.
|
||
}
|
||
|
||
\dedication[Dedication]{
|
||
Dedicated to some people.
|
||
}
|
||
|
||
\acknowledgements{ \OnePageChapter
|
||
This will be an acknowledgement.
|
||
}
|
||
|
||
\LoFisShort
|
||
\emptyLoT
|
||
|
||
\begin{document}
|
||
|
||
\input macros.tex
|
||
|
||
\chapter{Introduction}
|
||
|
||
Continuous low-thrust arcs utilizing technologies such as Ion propulsion, Hall thrusters, and
|
||
others can be a powerful tool in the design of interplanetary space missions. They tend to be
|
||
particularly suited to missions which require very high total change in velocity or $\Delta V$
|
||
values and take place over a particularly long duration. Traditional impulsive thrusting
|
||
techniques can achieve these changes in velocity, but they typically have a far lower specific
|
||
impulse and, as such, are much less efficient and use more fuel, costing the mission valuable
|
||
financial resources that could instead be used for science. Because of their inherently high
|
||
specific impulse (and thus efficiency), low-thrust fuels are well-suited to interplanetary
|
||
missions.
|
||
|
||
For instance, low thrust ion propulsion was used on the Bepi-Colombo, Dawn, and Deep
|
||
Space 1 missions. In general, anytime an interplanetary trajectory is posed, it is advisable to
|
||
first explore the possibility of low-thrust technologies. In an interplanetary mission, the
|
||
primary downside to low-thrust orbits (that they require significant time to achieve large
|
||
$\Delta V$ changes) is made irrelevant by the fact that interplanetary trajectories take such a
|
||
long time as a matter of course.
|
||
|
||
Another technique often leveraged by interplanetary trajectory designers is the gravity assist.
|
||
Gravity assists cleverly utilize the inertia of a large planetary body to ''slingshot`` a
|
||
spacecraft, modifying the direction of its velocity with respect to the central body, the Sun.
|
||
This technique lends itself very well to impulsive trajectories. The gravity assist maneuver
|
||
itself can be modeled very effectively by an impulsive maneuver with certain constraints, placed
|
||
right at the moment of closest approach to the (flyby) target body. Because of this,
|
||
optimization with impulsive trajectories and gravity assists are common.
|
||
|
||
% TODO: Might need to remove the HOCP stuff
|
||
However, there is no physical reason why low-thrust trajectories can't also incorporate gravity
|
||
assists. The optimization problem becomes much more complicated. The separate problems of
|
||
optimizing flyby parameters (planet, flyby date, etc.) and optimizing the low-thrust control
|
||
arcs don't combine very easily. In this paper, a technique is explored by setting the
|
||
dual-problem up as a Hybrid Optimal Control Problem (HOCP).
|
||
|
||
This thesis will explore these concepts in a number of different sections. Section
|
||
\ref{traj_opt} will explore the basic principles of trajectory optimization in a manner agnostic
|
||
to the differences between continuous low-thrust and impulsive high-thrust techniques. Section
|
||
\ref{low_thrust} will then delve into the different aspects to consider when optimizing a low
|
||
thrust mission profile over an impulsive one. Section \ref{interplanetary} provides more detail
|
||
on the interplanetary considerations, including force models and gravity assists. Section
|
||
\ref{algorithm} will cover the implementation details of the HOCP optimization algorithm
|
||
developed for this paper. Finally, section \ref{results} will explore the results of some
|
||
hypothetical missions to Saturn.
|
||
|
||
\chapter{Trajectory Optimization} \label{traj_opt}
|
||
|
||
Trajectory optimization is concerned with a narrow problem (namely, optimizing a spaceflight
|
||
trajectory to an end state) with a wide range of possible techniques, approaches, and even
|
||
solutions. In this section, the foundations for direct optimization of these sorts of problems
|
||
will be explored by first introducing the Two-Body Problem, then an algorithm for directly
|
||
solving for states in that system, then exploring approaches to Non-Linear Problem (NLP) solving
|
||
in general and how they apply to spaceflight trajectories.
|
||
|
||
\section{The Two-Body Problem}
|
||
The motion of a spacecraft in space is governed by a large number of forces. When planning and
|
||
designing a spacecraft trajectory, we often want to use the most complete (and often complex)
|
||
model of these forces that is available. However, in the process of designing these
|
||
trajectories, we often have to compute the path of the spacecraft many hundreds, thousands, or
|
||
even millions of times. Utilizing very high-fidelity force models that account for aerodynamic
|
||
pressures, solar radiation pressures, multi-body effects, and many others may be infeasible
|
||
for the method being used if the computations take too long.
|
||
|
||
Therefore, a common approach (and the one utilized in this implementation) is to first look
|
||
simply at the single largest force governing the spacecraft in motion, the gravitational force
|
||
due to the primary body around which it is orbiting. This can provide an excellent
|
||
low-to-medium fidelity model that can be extremely useful in categorizing the optimization
|
||
space as quickly as possible. In many cases, including the algorithm used in this paper, it is
|
||
unlikely that local cost-function minima would be missed due to the lack of fidelity of the
|
||
Two Body Problem.
|
||
|
||
In order to explore the Two Body Problem, we must first examine the full set of assumptions
|
||
associated with the force model. Firstly, we are only concerned with the nominative two
|
||
bodies: the spacecraft and the planetary body around which it is orbiting. Secondly, both of
|
||
these bodies are modeled as simple point masses. This removes the need to account for
|
||
non-uniform densities and asymmetry. The third assumption is that the mass of the spacecraft
|
||
($m_2$) is much much smaller than the mass of the planetary body ($m_1$) and enough so as to be
|
||
considered negligible. The only force acting on this system is then the force of gravity that
|
||
the primary body enacts upon the secondary. Lastly, we'll assume a fixed inertial frame. This
|
||
isn't necessary for the formulation of a solution, but will simplify the derivation.
|
||
|
||
Reducing the system to two point masses with a single gravitational force acting between them
|
||
(and only in one direction) we can model the force on the secondary body as:
|
||
|
||
\begin{equation}
|
||
\ddot{\vec{r}} = - \frac{G \left( m_1 + m_2 \right)}{r^2} \frac{\vec{r}}{\left| r \right|}
|
||
\end{equation}
|
||
|
||
Where $\vec{r}$ is the position of the spacecraft, $G$ is the universal gravitational
|
||
parameter, $m_1$ is the mass of the planetary body, and $m_2$ is the mass of the spacecraft.
|
||
Due to our assumption that the mass of the spacecraft is significantly smaller than the mass
|
||
of the primary body ($m_1 >> m_2$) we can reduce that formulation to simply:
|
||
|
||
\begin{equation}
|
||
\ddot{\vec{r}} = - \frac{\mu}{r^2} \hat{r}
|
||
\end{equation}
|
||
|
||
Where $\mu = G m_1$ is the specific gravitational parameter for our primary body of interest.
|
||
|
||
\subsection{Kepler's Laws and Equations}
|
||
|
||
% TODO: Can I segue better from 2BP to Keplerian geometry?
|
||
|
||
Now that we've fully qualified the forces acting within the Two Body Problem, we can concern
|
||
ourselves with more practical applications of this as a force model. It should be noted,
|
||
firstly, that the spacecraft's position and velocity (given an initial position and velocity
|
||
and of course the $\mu$ value of the primary body) is actually analytically solvable for all
|
||
future points in time. This can be easily observed by noting that there are three
|
||
one-dimensional equations (one for each component of the three-dimensional position) and
|
||
three unknowns (the three components of the second derivative of the position).
|
||
|
||
In the early 1600s, Johannes Kepler produced just such a solution. By taking advantages of
|
||
what is also known as ``Kepler's Laws'' which are:
|
||
|
||
\begin{enumerate}
|
||
\item Each planet's orbit is an ellipse with the Sun at one of the foci. This can be
|
||
expanded to any orbit by re-wording as ``all orbital paths follow a conic section
|
||
(circle, ellipse, parabola, or hyperbola) with a primary body at one of the foci''.
|
||
\item The area swept out by the imaginary line connecting the primary and secondary
|
||
bodies increases linearly with respect to time. This implies that the magnitude of the
|
||
orbital speed is not constant.
|
||
\item The square of the orbital period is proportional to the cube of the semi-major
|
||
axis of the orbit, regardless of eccentricity. Specifically, the relationship is: $T = 2
|
||
\pi \sqrt{\frac{a^3}{\mu}}$ where $T$ is the period and $a$ is the semi-major axis.
|
||
\end{enumerate}
|
||
|
||
\section{Analytical Solutions to Kepler's Equations}
|
||
|
||
Kepler was able to produce an equation to represent the angular displacement of an orbiting
|
||
body around a primary body as a function of time, which we'll derive now for the elliptical
|
||
case. Since the total area of an ellipse is the product of $\pi$, the semi-major axis, and
|
||
the semi-minor axis ($\pi a b$), we can relate (by Kepler's second law) the area swept out
|
||
by an orbit as a function of time:
|
||
|
||
\begin{equation}\label{swept}
|
||
\frac{\Delta t}{T} = \frac{k}{\pi a b}
|
||
\end{equation}
|
||
|
||
This leaves just one unknown variable $k$, which we can determine through use of the
|
||
geometric auxiliary circle, which is a circle with radius equal to the ellipse's semi-major
|
||
axis and center directly between the two foci, as in Figure~\ref{aux_circ}.
|
||
|
||
\begin{figure}
|
||
\centering
|
||
\includegraphics[width=0.8\textwidth]{fig/kepler}
|
||
\caption{Geometric Representation of Auxiliary Circle}\label{aux_circ}
|
||
\end{figure}
|
||
|
||
In order to find the area swept by the spacecraft, $k$, we can take advantage of the fact
|
||
that that area is the triangle $k_1$ subtracted from the elliptical segment $PCB$:
|
||
|
||
\begin{equation}\label{areas_eq}
|
||
k = area(seg_{PCB}) - area(k_1)
|
||
\end{equation}
|
||
|
||
Where the area of the triangle $k_1$ can be found easily using geometric formulae:
|
||
|
||
\begin{align}
|
||
area(k_1) &= \frac{1}{2} \left( ae - a \cos E \right) \left( \frac{b}{a} a \sin E \right) \\
|
||
&= \frac{ab}{2} \left(e \sin E - \cos E \sin E \right)
|
||
\end{align}
|
||
|
||
Now we can find the area for the elliptical segment $PCB$ by first finding the circular
|
||
segment $POB'$, subtracting the triangle $OB'C$, then applying the fact that an ellipse is
|
||
merely a vertical scaling of a circle by the amount $\frac{b}{a}$.
|
||
|
||
\begin{align}
|
||
area(PCB) &= \frac{b}{a} \left( \frac{a^2 E}{2} - \frac{1}{2} \left( a \cos E \right)
|
||
\left( a \sin E \right) \right) \\
|
||
&= \frac{abE}{2} - \frac{ab}{2} \left( \cos E \sin E \right) \\
|
||
&= \frac{ab}{2} \left( E - \cos E \sin E \right)
|
||
\end{align}
|
||
|
||
By substituting the two areas back into Equation~\ref{areas_eq} we can get the $k$ area
|
||
swept out by the spacecraft:
|
||
|
||
\begin{equation}
|
||
k = \frac{ab}{2} \left( E - e \sin E \right)
|
||
\end{equation}
|
||
|
||
Which we can then substitute back into the equation for the swept area as a function of
|
||
time (Equation~\ref{swept}):
|
||
|
||
\begin{equation}
|
||
\frac{\Delta t}{T} = \frac{E - e \sin E}{2 \pi}
|
||
\end{equation}
|
||
|
||
Which is, effectively, Kepler's equation. It is commonly known by a different form:
|
||
|
||
\begin{align}
|
||
M &= E - e \sin E \\
|
||
&= \sqrt{\frac{\mu}{a^3}} \Delta t
|
||
\end{align}
|
||
|
||
Where we've defined the mean anomaly as $M$ and used the fact that $T =
|
||
\sqrt{\frac{a^3}{\mu}}$. This provides us a useful relationship between Eccentric Anomaly
|
||
($E$) which can be related to spacecraft position, and time, but we still need a useful
|
||
algorithm for solving this equation.
|
||
|
||
\subsection{LaGuerre-Conway Algorithm}\label{laguerre}
|
||
For this application, I used an algorithm known as the LaGuerre-Conway algorithm, which was
|
||
presented in 1986 as a faster algorithm for directly solving Kepler's equation and has been
|
||
in use in many applications since. This algorithm is known for its convergence robustness
|
||
and also its speed of convergence when compared to higher order Newton methods.
|
||
|
||
This thesis will omit a step-through of the algorithm itself, but the code will be present
|
||
in the Appendix.
|
||
|
||
\section{Non-Linear Problem Optimization}
|
||
|
||
Now we can consider the formulation of the problem in a more useful way. For instance, given a
|
||
desired final state in position and velocity we can relatively easily determine the initial
|
||
state necessary to end up at that desired state over a pre-defined period of time by solving
|
||
Kepler's equation. In fact, this is often how impulsive trajectories are calculated since,
|
||
other than the impulsive thrusting event itself, the trajectory is entirely natural.
|
||
|
||
However, often in trajectory design we want to consider a number of other inputs. For
|
||
instance, a low thrust profile, a planetary flyby, the effects of rotating a solar panel on
|
||
solar radiation pressure, etc. Once these inputs have been accepted as part of the model, the
|
||
system is generally no longer analytically solvable, or, if it is, is too complex to calculate
|
||
directly.
|
||
|
||
Therefore an approach is needed, in trajectory optimization and many other fields, to optimize
|
||
highly non-linear, unpredictable systems such as this. The field that developed to approach
|
||
this problem is known as Non-Linear Problem (NLP) Optimization.
|
||
|
||
There are, however, two categories of approaches to solving an NLP. The first category,
|
||
indirect methods, involve declaring a set of necessary and/or sufficient conditions for declaring
|
||
the problem optimal. These conditions then allow the non-linear problem (generally) to be
|
||
reformulated as a two point boundary value problem. Solving this boundary value problem can
|
||
provide a control law for the optimal path. Indirect approaches for spacecraft trajectory
|
||
optimization have given us the Primer Vector Theory.
|
||
|
||
The other category is the direct methods. In a direct optimization problem, the cost function
|
||
itself is calculated to provide the optimal solution. The problem is usually thought of as a
|
||
collection of dynamics and controls. Then these controls can be modified to minimize the cost
|
||
function. A number of tools have been developed to optimize NLPs via this direct method in the
|
||
general case. For this particular problem, direct approaches were used as the low-thrust
|
||
system dynamics adds too much complexity to quickly optimize indirectly and the individual
|
||
optimization routines needed to proceed as quickly as possible.
|
||
|
||
\subsection{Non-Linear Solvers}
|
||
For these types of non-linear, constrained problems, a number of tools have been developed
|
||
that act as frameworks for applying a large number of different algorithms. This allows for
|
||
simple testing of many different algorithms to find what works best for the nuances of the
|
||
problem in question.
|
||
|
||
One of the most common of these NLP optimizers is SNOPT, which is a proprietary package
|
||
written primarily using a number of Fortran libraries by the Systems Optimization Laboratory
|
||
at Stanford University. It uses a sparse sequential quadratic programming approach.
|
||
|
||
Another common NLP optimization packages (and the one used in this implementation) is the
|
||
Interior Point Optimizer or IPOPT. It can be used in much the same way as SNOPT and uses an
|
||
Interior Point Linesearch Filter Method and was developed as an open-source project by the
|
||
organization COIN-OR under the Eclipse Public License.
|
||
|
||
Both of these methods utilize similar approaches to solve general constrained non-linear
|
||
problems iteratively. Both of them can make heavy use of derivative Jacobians and Hessians
|
||
to improve the convergence speed and both have been ported for use in a number of
|
||
programming languages, including in Julia, which was used for this project.
|
||
|
||
This is by no means an exhaustive list, as there are a number of other optimization
|
||
libraries that utilize a massive number of different algorithms. For the most part, the
|
||
libraries that port these are quite modular in the sense that multiple algorithms can be
|
||
tested without changing much source code.
|
||
|
||
\subsection{Linesearch Method}
|
||
As mentioned above, this project utilized IPOPT which leveraged an Interior Point Linesearch
|
||
method. A linesearch algorithm is one which attempts to find the optimum of a non-linear
|
||
problem by first taking an initial guess $x_k$. The algorithm then determines a step
|
||
direction (in this case through the use of automatic differentiation to calculate the
|
||
derivatives of the non-linear problem) and a step length. The linesearch algorithm then
|
||
continues to step the initial guess, now labeled $x_{k+1}$ after the addition of the
|
||
``step'' vector and iterates this process until predefined termination conditions are met.
|
||
|
||
In this case, the IPOPT algorithm was used, not as an optimizer, but as a solver. For
|
||
reasons that will be explained in the algorithm description in Section~\ref{algorithm} it
|
||
was sufficient merely that the non-linear constraints were met, therefore optimization (in
|
||
the particular step in which IPOPT was used) was unnecessary.
|
||
|
||
\chapter{Low-Thrust Considerations} \label{low_thrust}
|
||
|
||
Thus far, the techniques that have been discussed can be equally useful for both impulsive and
|
||
continuous thrust mission profiles. In this section, we'll discuss the intricacies of continuous
|
||
low-thrust trajectories in particular. There are many methods for optimizing such profiles and
|
||
we'll briefly discuss the difference between a direct and indirect optimization of a low-thrust
|
||
trajectory as well as introduce the concept of a control law and the notation used in this
|
||
thesis for modelling low-thrust trajectories more simply.
|
||
|
||
\section{Low-Thrust Control Laws}
|
||
|
||
In determining a low-thrust arc, a number of variables must be accounted for and, ideally,
|
||
optimized.
|
||
|
||
Firstly, we must determine the presence or absence of thrust. Often, this is a question of
|
||
preference in the arsenal of the mission designer. Generally speaking, there are points along
|
||
an orbit at which thrusting in order to achieve the final orbit are more or less efficient.
|
||
For instance, in a classic orbit raising, if increasing the semi-major axis is the only goal,
|
||
then thrusting nearer to the periapsis is far more efficient than thrusting near the apoapsis.
|
||
For this reason, a mission designer may choose to reduce the thrust or turn it off altogether
|
||
during certain segments of the trajectory.
|
||
|
||
Secondly, the direction of thrust must also be determined. The methods for determining this
|
||
direction varies greatly depending on the particular control law chosen for that mission.
|
||
Generally speaking, a control law determines these two parameters: thrust presence and thrust
|
||
direction, at each point along the arc.
|
||
|
||
This is, of course, also true for impulsive trajectories. However, since the thrust presence
|
||
for those trajectories are generally taken to be impulse functions, the control laws can
|
||
afford to be much less complicated for a given mission goal, by simply thrusting only at the
|
||
moment on the orbit when the transition will be most efficient. For a low-thrust mission,
|
||
however, the control law must be continuous rather than discrete and therefore the control law
|
||
inherently gains a lot of complexity.
|
||
|
||
\section{Sims-Flanagan Transcription}
|
||
|
||
The major problem with optimizing low thrust paths is that the control law must necessarily be
|
||
continuous. Also, since indirect optimization approaches are quite difficult, the problem must
|
||
necessarily be reformulated as a discrete one in order to apply a direct approach. Therefore,
|
||
this thesis chose to use a model well suited for discretizing low-thrust paths: the
|
||
Sims-Flanagan transcription (SFT).
|
||
|
||
The SFT is actually quite a simple method for discretizing low-thrust arcs. First the
|
||
continuous arc is subdivided into a number ($N$) of individual consistent timesteps of length
|
||
$\frac{tof}{N}$. The control thrust is then applied at the center of each of these time steps.
|
||
|
||
Using the SFT, it is relatively straightforward to propagate a state (in the context of the
|
||
Two-Body Problem) that utilizes a continuous low-thrust control, without the need for
|
||
computationally expensive numeric integration algorithms, by simply solving Kepler's equation
|
||
(using the LaGuerre-Conway algorithm introduced in Section~\ref{laguerre}) $N$ times. This
|
||
greatly reduces the computation complexity, which is particularly useful for cases in which
|
||
low-thrust trajectories need to be calculated many millions of times, as is the case in this
|
||
thesis. The fidelity of the model can also be easily fine-tuned. By simply increasing the
|
||
number of sub-arcs, one can rapidly approach a fidelity equal to a continuous low-thrust
|
||
trajectory within the Two-Body Problem, with only linearly-increasing computation time.
|
||
|
||
\chapter{Interplanetary Trajectory Considerations} \label{interplanetary}
|
||
|
||
The question of interplanetary travel opens up a host of additional new complexities. While
|
||
optimizations for simple single-body trajectories are far from simple, it can at least be
|
||
said that the assumptions of the Two Body Problem remain fairly valid. In interplanetary
|
||
travel, the primary body most responsible for gravitational forces might be a number of
|
||
different bodies, dependent on the phase of the mission. In the ideal case, every relevant
|
||
body would be considered as an ``n-body'' perturbation during the entire trajectory. For
|
||
some approaches, this method is sufficient and preferred. However, for other uses, a more
|
||
efficient model is necessary. The method of patched conics can be applied in this case to
|
||
simplify the model.
|
||
|
||
Interplanetary travel does not simply negatively impact trajectory optimization. The
|
||
increased complexity of the search space also opens up new opportunities for orbit
|
||
strategies. The primary strategy investigated by this thesis will be the gravity assist, a
|
||
technique for utilizing the gravitational energy of a planet to modify the direction of
|
||
solar velocity.
|
||
|
||
\section{Patched Conics}
|
||
|
||
The first hurdle to deal with is the problem of reconciling the Two-Body problem with
|
||
the presence of multiple and varying planetary bodies. The most common method for
|
||
approaching this is the method of patched conics. In this model, we break the
|
||
interplanetary trajectory up into a series of smaller sub-trajectories. During each of
|
||
these sub-trajectories, a single primary is considered to be responsible for the
|
||
trajectory of the orbit, via the Two-Body problem.
|
||
|
||
The transition point can be calculated a variety of ways. The most typical method is to
|
||
calculate the gravitational force due to the two bodies separately, via the Two-Body
|
||
models. Whichever primary is a larger influence on the motion of the spacecraft is
|
||
considered to be the primary at that moment. This effectively breaks the trajectory into
|
||
a series of orbits defined by the Two-Body problem (conics), patched together by
|
||
distinct transition points.
|
||
|
||
\section{Gravity Assist Maneuvers}
|
||
|
||
As previously mentioned, there are methods for utilizing the orbital energy of the other
|
||
planets in the Solar System. This is achieved via a technique known as a Gravity Assist,
|
||
or a Gravity Flyby. During a gravity assist, the spacecraft enters into the
|
||
gravitational sphere of influence of the planet and, because of its excess velocity,
|
||
proceeds to exit the sphere of influence. Relative to the planet, the speed of the
|
||
spacecraft increases as it approaches, then decreases as it departs. From the
|
||
perspective of the planet, the velocity of the spacecraft is unchanged. However, the
|
||
planet is also orbiting the Sun.
|
||
|
||
From the perspective of a Sun-centered frame, though, this is effectively an elastic
|
||
collision. The overall momentum remains the same, with the spacecraft either gaining or
|
||
losing some in the process (dependent on the directions of travel). The planet also
|
||
loses or gains momentum enough to maintain the overall system momentum, but this amount
|
||
is negligible compared to the total momentum of the planet. The overall effect is that
|
||
the spacecraft arrives at the planet from one direction and, because of the influence of
|
||
the planet, leaves in a different direction.
|
||
|
||
This effect can be used strategically. The ``bend'' due to the flyby is actually
|
||
tunable via the exact placement of the fly-by in the b-frame, or the frame centered at
|
||
the planet, from the perspective of the spacecraft at $v_\infty$. By modifying the
|
||
turning angle of this bend. In doing so, one can effectively achieve a (restricted) free
|
||
impulsive thrust event.
|
||
|
||
\section{Multiple Gravity Assist Techniques}
|
||
|
||
Naturally, therefore, one would want to utilize these gravity flybys to reduce the fuel
|
||
cost to arrive at their destination target state. However, these flyby maneuvers are
|
||
quite restricted. The incoming hyperbolic velocity must be equal in magnitude to the
|
||
outgoing hyperbolic velocity. Also, the turning angle $\delta$, in the following
|
||
equation, correlates with the radius of periapsis of the hyperbolic trajectory crossing
|
||
the planet:
|
||
|
||
\begin{equation}
|
||
r_p = \frac{\mu}{v_\infty^2} \left[ \frac{1}{\sin\left(\frac{\delta}{2}\right)} - 1 \right]
|
||
\end{equation}
|
||
|
||
Where $v_\infty$ is the magnitude of hyperbolic velocity. Naturally, the radius of
|
||
periapsis must not fall below some safe value, in order to avoid the risk of the
|
||
spacecraft crashing into the planet or its atmosphere.
|
||
|
||
In order to visualize which trajectories are possible within these constraints, porkchop
|
||
plots are often employed, such as the plot in Figure~\ref{porkchop}. These plots outline
|
||
various incoming and outgoing qualities of the trajectory arc between two planetary
|
||
bodies. For instance, during an arc from launch at Earth to a flyby one might plot the
|
||
launch C3 against the Mars arrival $v_\infty$ for a variety of launch and arrival dates.
|
||
|
||
\begin{figure}
|
||
\centering
|
||
\includegraphics[width=\textwidth]{fig/porkchop}
|
||
\caption{A sample porkchop plot of an Earth-Mars transfer}
|
||
\label{porkchop}
|
||
\end{figure}
|
||
|
||
This is made possible by solving Lambert's problem for the planetary ephemeris at the
|
||
epochs plotted. Lambert's problem is concerned with determining the orbit between two
|
||
positions at two different times in space. There are a number of different Lambert's
|
||
problem algorithms that allow a mission designer to determine the velocity needed (and
|
||
thus the $\Delta V$) required to achieve a position at a later time. From this, the
|
||
designer can algorithmically determine trajectory properties in the porkchop plot for
|
||
easy visualization.
|
||
|
||
However, this is an impulsive thrust-centered approach. The solution to Lambert's
|
||
problem assumes a natural trajectory. However, to the low-thrust designer, this is
|
||
needlessly limiting. A natural trajectory is unnecessary when the trajectory can be
|
||
modified by a continuous thrust profile along the arc. Therefore, for the hybrid problem
|
||
of optimizing both flyby selection and thrust profiles, porkchop plots are less helpful,
|
||
and an algorithmic approach is preferred.
|
||
|
||
% \chapter{Genetic Algorithms}
|
||
% I will probably give only a brief overview of genetic algorithms here. I don't personally know
|
||
% that much about them. Then in the following subsections I can discuss the parts that are
|
||
% relevant to the specific algorithm that I'm using.
|
||
|
||
% \section{Decision Vectors}
|
||
% Discuss what a decision vector is in the context of an optimization problem.
|
||
|
||
% \section{Selection and Fitness Evaluation}
|
||
% Discuss the costing being used as well as the different types of fitness evaluation that are
|
||
% common. Also discuss the concept of generations and ``survival''.
|
||
|
||
% \subsection{Tournament Selection}
|
||
% Dive deeper into the specific selection algorithm being used here.
|
||
|
||
% \section{Crossover}
|
||
% Discuss the concept of crossover and procreation in a genetic algorithm.
|
||
|
||
% \subsection{Binary Crossover}
|
||
% Discuss specific crossover algorithm used here.
|
||
|
||
% \subsection{Mutation}
|
||
% Discuss both the necessity for mutation and the mutation algorithm being used.
|
||
|
||
\chapter{Algorithm Overview} \label{algorithm}
|
||
|
||
In this section, we will review the actual execution of the algorithm developed. As an
|
||
overview, the routine was developed to enable the determination of an optimized spacecraft
|
||
trajectory from the selection of some very basic mission parameters. Those parameters
|
||
include:
|
||
|
||
\begin{itemize}
|
||
\item Spacecraft dry mass
|
||
\item Thruster Specific Impulse
|
||
\item Thruster Maximum Thrusting Force
|
||
\item Thruster Duty Cycle Percentage
|
||
\item Number of Thruster on Spacecraft
|
||
\item Total Starting Weight of the Spacecraft
|
||
\item A Maximum Acceptable $V_\infty$ at arrival and $C_3$ at launch
|
||
\item The Launch Window Timing and the Latest Arrival
|
||
\item A cost function relating the mass usage, $v_\infty$ at arrival, and $C_3$ at
|
||
launch to a cost
|
||
\item A list of flyby planets starting with Earth and ending with the destination
|
||
\end{itemize}
|
||
|
||
Which allows for extremely automated optimization of the trajectory, while still providing
|
||
the mission designer with the flexibility to choose the particular flyby planets to
|
||
investigate.
|
||
|
||
This is achieved via an optimal control problem in which the ``inner loop'' is a
|
||
non-linear programming problem to determine the optimal low-thrust control law and flyby
|
||
parameters given a suitable initial guess. Then an ``outer loop'' monotonic basin hopping
|
||
algorithm is used to traverse the search space and more carefully optimize the solutions
|
||
found by the inner loop.
|
||
|
||
\section{Trajectory Composition}
|
||
|
||
In this thesis, a specific nomenclature will be adopted to define the stages of an
|
||
interplanetary mission in order to standardize the discussion about which aspects of the
|
||
software affect which phases of the mission.
|
||
|
||
Overall, a mission is considered to be the entire overall trajectory. In the context of
|
||
this software procedure, a mission is taken to always begin at the Earth, with some
|
||
initial launch C3 intended to be provided by an external launch vehicle. This C3 is not
|
||
fully specified by the mission designer, but instead is optimized as a part of the
|
||
overall cost function (and normalized by a designer-specified maximum allowable value).
|
||
|
||
This overall mission can then be broken down into a variable number of ``phases''
|
||
defined as beginning at one planetary body with some excess hyperbolic velocity and
|
||
ending at another. The first phase of the mission is from the Earth to the first flyby
|
||
planet. The final phase is from the last flyby planet to the planet of interest.
|
||
|
||
Each of these phases are then connected by a flyby event at the boundary. Each flyby
|
||
event must satisfy the following conditions:
|
||
|
||
\begin{enumerate}
|
||
\item The planet at the end of one phase must match the planet at the beginning of
|
||
the next phase.
|
||
\item The magnitude of the excess hyperbolic velocity coming into the planet (at the
|
||
end of the previous phase) must equal the magnitude of the excess hyperbolic
|
||
velocity leaving the planet (at the beginning of the next phase).
|
||
\item The flyby ``turning angle'' must be such that the craft maintains a safe
|
||
minimum altitude above the surface or atmosphere of the flyby planet.
|
||
\end{enumerate}
|
||
|
||
These conditions then effectively stitch the separate mission phases into a single
|
||
coherent mission, allowing for the optimization of both individual phases and the entire
|
||
mission as a whole. This nomenclature is similar to the nomenclature adopted by Jacob
|
||
Englander in his Hybrid Optimal Control Problem paper, but does not allow for missions
|
||
with multiple targets, simplifying the optimization.
|
||
|
||
\section{Inner Loop Implementation}
|
||
|
||
The optimization routine can be reasonable separated into two separate ``loops'' wherein
|
||
the first loop is used, given an initial guess, to find valid trajectories within the
|
||
region of the initial guess and submit the best. The outer loop is then used to traverse
|
||
the search space and supply the initial loop with a number of well chosen initial
|
||
guesses.
|
||
|
||
Figure~\ref{nlp} provides an overview of the process of breaking a mission guess down
|
||
into an NLP, but there are essentially three primary routines involved in the inner
|
||
loop. A given state is propagated forward using the LaGuerre-Conway Kepler solution
|
||
algorithm, which itself is used to generate powered trajectory arcs via the
|
||
Sims-Flanagan transcribed propagator. Finally, these powered arcs are connected via a
|
||
multiple-shooting non-linear optimization problem. The trajectories describing each
|
||
phase complete one ``Mission Guess'' which is fed to the non-linear solver to generate
|
||
one valid trajectory within the vicinity of the original Mission Guess.
|
||
|
||
\begin{figure}
|
||
\centering
|
||
\includegraphics[width=\textwidth]{flowcharts/nlp}
|
||
\caption{A flowchart of the Non-Linear Problem Solving Formulation}
|
||
\label{nlp}
|
||
\end{figure}
|
||
|
||
\subsection{LaGuerre-Conway Kepler Solver}
|
||
|
||
The most basic building block of any trajectory is a physical model for simulating
|
||
natural trajectories from one point forward in time. The approach taken by this
|
||
paper uses the solution to Kepler's equation put forward by
|
||
Conway\cite{laguerre_conway} in 1986 in order to provide simple and very
|
||
processor-efficient propagation without the use of integration. The code logic
|
||
itself is actually quite simple, providing an approach similar to the Newton-Raphson
|
||
approach for finding the roots of the Battin form of Kepler's equation.
|
||
|
||
The following pseudo-code outlines the approach taken for the elliptical case. The
|
||
approach is quite similar when $a<0$:
|
||
|
||
% TODO: Some symbols here aren't recognized by the font
|
||
\begin{singlespacing}
|
||
\begin{verbatim}
|
||
i = 0
|
||
# First declare some useful variables from the state
|
||
σ0 = (position ⋅ velocity) / √(μ)
|
||
a = 1 / ( 2/norm(position) - norm(velocity)^2/μ )
|
||
coeff = 1 - norm(position)/a
|
||
|
||
# This loop is essentially a second-order Newton solver for ΔE
|
||
ΔM = ΔE_new = √(μ/a^3) * time
|
||
ΔE = 1000
|
||
while abs(ΔE - ΔE_new) > 1e-10
|
||
ΔE = ΔE_new
|
||
F = ΔE - ΔM + σ0 / √(a) * (1-cos(ΔE)) - coeff * sin(ΔE)
|
||
dF = 1 + σ0 / √(a) * sin(ΔE) - coeff * cos(ΔE)
|
||
d2F = σ0 / √(a) * cos(ΔE) + coeff * sin(ΔE)
|
||
ΔE_new = ΔE - n*F / ( dF + sign(dF) * √(abs((n-1)^2*dF^2 - n*(n-1)*F*d2F )))
|
||
i += 1
|
||
end
|
||
|
||
# ΔE can then be used to determine the F/Ft and G/Gt coefficients
|
||
F = 1 - a/norm(position) * (1-cos(ΔE))
|
||
G = a * σ0/ √(μ) * (1-cos(ΔE)) + norm(position) * √(a) / √(μ) * sin(ΔE)
|
||
r = a + (norm(position) - a) * cos(ΔE) + σ0 * √(a) * sin(ΔE)
|
||
Ft = -√(a)*√(μ) / (r*norm(position)) * sin(ΔE)
|
||
Gt = 1 - a/r * (1-cos(ΔE))
|
||
|
||
# Which provide transformations from the original position and velocity to the
|
||
# final
|
||
final_position = F*position + G*velocity
|
||
final_velocity = Ft*position + Gt*velocity
|
||
\end{verbatim}
|
||
\end{singlespacing}
|
||
|
||
This approach was validated by generating known good orbits in the 2 Body Problem.
|
||
For example, from the orbital parameters of a certain state, the orbital period can
|
||
be determined. If the system is then propagated for an integer multiple of the orbit
|
||
period, the state should remain exactly the same as it began. In
|
||
Figure~\ref{laguerre_plot} an example of such an orbit is provided.
|
||
|
||
\begin{figure}
|
||
\centering
|
||
\includegraphics[width=\textwidth]{fig/laguerre_plot}
|
||
\caption{Example of a natural trajectory propagated via the Laguerre-Conway
|
||
approach to solving Kepler's Problem}
|
||
\label{laguerre_plot}
|
||
\end{figure}
|
||
% TODO: Consider adding a paragraph about the improvements in processor time
|
||
|
||
\subsection{Sims-Flanagan Propagator}
|
||
|
||
Until this point, we've not yet discussed how best to model the low-thrust
|
||
trajectory arcs themselves. The Laguerre-Conway algorithm efficiently determines
|
||
natural trajectories given an initial state, but it still remains, given a control
|
||
law, that we'd like to determine the trajectory of a system with continuous input
|
||
thrust.
|
||
|
||
For this, we leverage the Sims-Flanagan transcription mentioned earlier. This allows
|
||
us to break a single phase into a number of ($n$) different arcs. At the center of
|
||
each of these arcs we can place a small impulsive burn, scaled appropriately for the
|
||
thruster configured on the spacecraft of interest. Therefore, for any given phase,
|
||
we actually split the trajectory into $2n$ sub-trajectories, with $n$ scaled
|
||
impulsive thrust events. As $n$ is increased, the trajectory becomes increasingly
|
||
accurate as a model of low-thrust propulsion in the 2BP. This allows the mission
|
||
designer to trade-off speed of propagation and the fidelity of the results quite
|
||
effectively.
|
||
|
||
\begin{figure}
|
||
\centering
|
||
\includegraphics[width=\textwidth]{fig/spiral_plot}
|
||
\caption{An example trajectory showing that classic continuous-thrust orbit
|
||
shapes, such as this orbit spiral, are easily achievable using a Sims-Flanagan
|
||
model}
|
||
\label{sft_plot}
|
||
\end{figure}
|
||
|
||
Figure~\ref{sft_plot} shows that the Sims-Flanagan transcription model can be used
|
||
to effectively model these types of orbit trajectories. In fact, the Sims-Flanagan
|
||
model is capable of modeling nearly any low-thrust trajectory with a sufficiently
|
||
high number of $n$ samples.
|
||
|
||
\subsection{Non-Linear Problem Solver}
|
||
|
||
Now that we have the basic building blocks of a continuous-thrust trajectory, we can
|
||
leverage one of the many non-linear optimization packages to find solutions near to
|
||
a (proposed) trajectory. This trajectory need not be valid.
|
||
|
||
For the purposes of discussion in this Section, we will assume that the inner-loop
|
||
algorithm starts with just such a ''Mission Guess``, which represents the proposed
|
||
trajectory. However, we'll briefly mention what quantities are needed for this
|
||
input:
|
||
|
||
A Mission Guess object contains:
|
||
\begin{singlespacing}
|
||
\begin{itemize}
|
||
\item The spacecraft and thruster parameters for the mission
|
||
\item A launch date
|
||
\item A launch $v_\infty$ vector representing excess Earth velocity
|
||
\item For each phase of the mission:
|
||
\begin{itemize}
|
||
\item The planet that the spacecraft will encounter (either flyby or
|
||
complete the mission) at the end of the phase
|
||
\item The $v_{\infty,out}$ vector representing excess velocity at the
|
||
planetary flyby (or launch if phase 1) at the beginning of the phase
|
||
\item The $v_{\infty,in}$ vector representing excess velocity at the
|
||
planetary flyby (or completion of mission) at the end of the phase
|
||
\item The time of flight for the phase
|
||
\item The unit-thrust profile in a sun-fixed frame represented by a
|
||
series of vectors with each element ranging from 0 to 1.
|
||
\end{itemize}
|
||
\end{itemize}
|
||
\end{singlespacing}
|
||
|
||
From this information, as can be seen in Figure~\ref{nlp}, we can formulate the
|
||
mission in terms of a non-linear problem. Specifically, the Mission Guess object can
|
||
be represented as a vector, $x$, the propagation function as a function $F$, and the
|
||
constraints as another function $G$ such that $G(x) = \vec{0}$.
|
||
|
||
This is a format that we can apply directly to the IPOPT solver, which Julia (the
|
||
programming language used) can utilize via bindings supplied by the SNOW.jl
|
||
package\cite{snow}.
|
||
|
||
IPOPT also requires the derivatives of both the $F$ and $G$ functions in the
|
||
formulation above. Generally speaking, a project designer has two options for
|
||
determining derivatives. The first option is to analytically determine the
|
||
derivatives, guaranteeing accuracy, but requiring processor time if determined
|
||
algorithmically and sometimes simply impossible or mathematically very rigorous to
|
||
determine manually. The second option is to numerically derive the derivatives,
|
||
using a technique such as finite differencing. This limits the accuracy, but can be
|
||
faster than algorithmic symbolic manipulation and doesn't require rigorous manual
|
||
derivations.
|
||
|
||
However, the Julia language has an excellent interface to a new technique, known as
|
||
automatic differentiation\cite{RevelsLubinPapamarkou2016}. Automatic differentiation
|
||
takes a slightly different approach to numerical derivation. It takes advantage of
|
||
the fact that any algorithmic function, no matter how complicated, can be broken
|
||
down into a series of smaller arithmetic functions, down to the level of simple
|
||
arithmetic. Since all of these simple arithmetic functions have a known derivative,
|
||
we can define a new datatype that carries through the function both the float and a
|
||
second number representing the derivative. Then, by applying (to the derivative) the
|
||
chain rule for every minute arithmetic function derivative as that arithmetic
|
||
function is applied to the main float value, the derivative can be determined,
|
||
accurate to the machine precision of the float type being used, with a processing
|
||
equivalent of two function calls (this of course depends on the simplicity of the
|
||
chained derivatives compared to the function pieces themselves). Generally speaking
|
||
this is much faster than the three or more function calls necessary for accurate
|
||
finite differencing and removes the need for the designer to tweak the epsilon value
|
||
in order to achieve maximum precision.
|
||
|
||
\section{Outer Loop Implementation}
|
||
Overview the outer loop. This may require a final flowchart, but might potentially be too
|
||
simple to lend itself to one.
|
||
|
||
\subsection{Inner Loop Calling Function}
|
||
The primary reason for including this section is to discuss the error handling.
|
||
|
||
\subsection{Monotonic Basin Hopping}
|
||
Outline the MBH algorithm, going into detail at each step. Mention the long-tailed PDF being
|
||
used and go into quite a bit of detail. Englander's paper on the MBH algorithm specifically
|
||
should be a good guide. Mention validation.
|
||
|
||
\chapter{Results Analysis} \label{results}
|
||
Simply highlight that the algorithm was tested on a sample trajectory to Saturn.
|
||
|
||
\section{Sample Trajectory to Saturn}
|
||
Give an overview of the trajectory that was ultimately chosen.
|
||
|
||
\subsection{Comparison to Less Optimal Solutions}
|
||
I should have a number of elite but less-optimal solutions. Honestly, I may write the
|
||
algorithm to keep all of the solutions to provide many points of comparison here.
|
||
|
||
\subsection{Cost Function Analysis}
|
||
Give some real-world context for the mass-use, time-of-flight, etc.
|
||
|
||
\subsection{Comparison to Impulsive Trajectories}
|
||
I may also remove this section. I could do a quick comparison (using porkchop plots) to
|
||
similar impulsive trajectories. Honestly, this is a lot of work for very little gain,
|
||
though, so probably the first place to chop if needed.
|
||
|
||
\chapter{Conclusion} \label{conclusion}
|
||
\section{Overview of Results}
|
||
Quick re-wording of the previous section in a paragraph or two for reader's convenience.
|
||
|
||
\section{Applications of Algorithm}
|
||
Talk a bit about why this work is valuable. Missions that could have benefited, missions that
|
||
this enables, etc.
|
||
|
||
\section{Recommendations for Future Work}
|
||
Recommend future work, obviously. There are a \emph{ton} of opportunities for improvement
|
||
including parallelization, cluster computing, etc.
|
||
|
||
\bibliographystyle{plain}
|
||
\nocite{*}
|
||
\bibliography{thesis}
|
||
|
||
\appendix
|
||
% \input appendixA.tex
|
||
|
||
\end{document}
|