connor/thesis
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'main' of https://gitlab.rconnorjohnstone.com/school/thesis

Browse Source
This commit is contained in:
Connor Johnstone
2022-03-19 13:58:21 -06:00
parent aed8f2d78b 5329d15067
commit e731ed16f5
18 changed files with 332 additions and 664 deletions
Download Patch File Download Diff File Expand all files Collapse all files

76
LaTeX/approach.tex
View File

@@ -1,11 +1,10 @@
\chapter{Algorithm Overview} \label{algorithm} \chapter{Algorithm Overview} \label{algorithm}
This thesis will attempt to develop an algorithm for the preliminary analysis of feasibility in This thesis focuses on designing a low-thrust interplanetary mission to an outer planet by
designing a low-thrust interplanetary mission to an outer planet by leveraging a monotonic basin leveraging a monotonic basin hopping algorithm. This section will review the actual execution of
hopping algorithm. In this section, we will review the actual execution of the algorithm the algorithm developed. As an overview, the routine is designed to enable the determination of
developed. As an overview, the routine was designed to enable the determination of an optimized an optimal spacecraft trajectory that minimizes propellant usage and $C_3$ from the selection of
spacecraft trajectory from the selection of some very basic mission parameters. Those parameters some very basic parameters. Those parameters include:
include:
\begin{itemize} \begin{itemize}
\setlength\itemsep{-0.5em} \setlength\itemsep{-0.5em}
@@ -25,7 +24,7 @@
\end{itemize} \end{itemize}
Which allows for an automated approach to optimization of the trajectory, while still providing Which allows for an automated approach to optimization of the trajectory, while still providing
the mission designer with the flexibility to choose the particular flyby planets to investigate. the 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'' involves solving a This is achieved via an optimal control problem in which the ``inner loop'' involves solving a
TPBVP to find the optimal solution given a suitable initial guess. Then an ``outer loop'' TPBVP to find the optimal solution given a suitable initial guess. Then an ``outer loop''
@@ -115,7 +114,7 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{flowcharts/nlp} \includegraphics[width=\textwidth]{LaTeX/flowcharts/nlp}
\caption{A flowchart of the TPBVP Solution Approach} \caption{A flowchart of the TPBVP Solution Approach}
\label{nlp} \label{nlp}
\end{figure} \end{figure}
@@ -133,12 +132,11 @@
The following pseudo-code outlines the approach taken for the elliptical case. The The following pseudo-code outlines the approach taken for the elliptical case. The
approach is quite similar when $a<0$: approach is quite similar when $a<0$:
% TODO: Some symbols here aren't recognized by the font
\begin{singlespacing} \begin{singlespacing}
\begin{verbatim} \begin{verbatim}
i = 0 i = 0
# First declare some useful variables from the state # First declare some useful variables from the state
sig0 = (position ⋅ velocity) / √(mu) sig0 = dot(position, velocity) / √(mu)
a = 1 / ( 2/norm(position) - norm(velocity)^2/mu ) a = 1 / ( 2/norm(position) - norm(velocity)^2/mu )
coeff = 1 - norm(position)/a coeff = 1 - norm(position)/a
@@ -178,13 +176,13 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{fig/laguerre_plot} \includegraphics[width=\textwidth]{LaTeX/fig/laguerre_plot}
\caption{Example of a natural trajectory propagated via the Laguerre-Conway \caption{Example of a natural trajectory propagated via the Laguerre-Conway
approach to solving Kepler's Problem} approach to solving Kepler's Problem}
\label{laguerre_plot} \label{laguerre_plot}
\end{figure} \end{figure}
\subsection{Sims-Flanagan Propagator} \subsection{Propagating with Sims-Flanagan Transcription}
Until this point, we've not yet discussed how best to model the low-thrust Until this point, we've not yet discussed how best to model the low-thrust
trajectory arcs themselves. The Laguerre-Conway algorithm efficiently determines trajectory arcs themselves. The Laguerre-Conway algorithm efficiently determines
@@ -204,7 +202,7 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{fig/spiral_plot} \includegraphics[width=\textwidth]{LaTeX/fig/spiral_plot}
\caption{An example trajectory showing that classic continuous-thrust orbit \caption{An example trajectory showing that classic continuous-thrust orbit
shapes, such as this orbit spiral, are easily achievable using a Sims-Flanagan shapes, such as this orbit spiral, are easily achievable using a Sims-Flanagan
model} model}
@@ -228,11 +226,9 @@
and the mass flow rate (a function of the duty cycle percentage ($d$), thrust ($f$), and the mass flow rate (a function of the duty cycle percentage ($d$), thrust ($f$),
and the specific impulse of the thruster ($I_{sp}$), commonly used to measure and the specific impulse of the thruster ($I_{sp}$), commonly used to measure
efficiency)\cite{sutton2016rocket}: efficiency)\cite{sutton2016rocket}:
\begin{equation} \begin{equation}
\Delta m = \Delta t \frac{f d}{I_{sp} g_0} \Delta m = \Delta t \frac{f d}{I_{sp} g_0}
\end{equation} \end{equation}
Where $\Delta m$ is the fuel used in the sub-trajectory, $\Delta t$ is the time of Where $\Delta m$ is the fuel used in the sub-trajectory, $\Delta t$ is the time of
flight of the sub-trajectory, and $g_0$ is the standard gravity at the surface of flight of the sub-trajectory, and $g_0$ is the standard gravity at the surface of
Earth. From knowledge of the mass flow rate, we can then decrement the mass Earth. From knowledge of the mass flow rate, we can then decrement the mass
@@ -264,7 +260,7 @@
\item The $v_{\infty,in}$ vector representing excess velocity at the \item The $v_{\infty,in}$ vector representing excess velocity at the
planetary flyby (or completion of mission) at the end of the phase planetary flyby (or completion of mission) at the end of the phase
\item The time of flight for the phase \item The time of flight for the phase
\item The unit-thrust profile in a sun-fixed frame represented by a \item The unit-thrust profile in a sun-centered frame represented by a
series of vectors with each element ranging from 0 to 1. series of vectors with each element ranging from 0 to 1.
\end{itemize} \end{itemize}
\end{itemize} \end{itemize}
@@ -272,10 +268,9 @@
From this information, as can be seen in Figure~\ref{nlp}, we can formulate the mission From this information, as can be seen in Figure~\ref{nlp}, we can formulate the mission
in terms of a non-linear programming problem. Specifically, the variables describing the in terms of a non-linear programming problem. Specifically, the variables describing the
trajectory contained within the Guess object can be represented as an input vector, trajectory from the free variable, $\vec{x}$, the cost function produced by an entire
$\vec{x}$, the cost function produced by an entire trajectory propagation as $F$, and trajectory propagation, $F$, and the constraints that the trajectory must satisfy as
the constraints that the trajectory must satisfy as another function $\vec{G}$ such that another function $\vec{G}$ such that $\vec{G}(\vec{x}) = \vec{0}$.
$\vec{G}(\vec{x}) = \vec{0}$.
This is a format that we can apply directly to the IPOPT solver, which Julia (the 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 programming language used) can utilize via bindings supplied by the SNOW.jl
@@ -325,7 +320,7 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{flowcharts/mbh} \includegraphics[width=\textwidth]{LaTeX/flowcharts/mbh}
\caption{A flowchart visualizing the steps in the monotonic basin hopping \caption{A flowchart visualizing the steps in the monotonic basin hopping
algorithm} algorithm}
\label{mbh_flow} \label{mbh_flow}
@@ -333,12 +328,12 @@
\subsection{Random Trajectory Generation}\label{random_gen_section} \subsection{Random Trajectory Generation}\label{random_gen_section}
At a basic level, the algorithm needs to produce a guess (represented by all of the At a basic level, the algorithm needs to produce a guess for the free variable vector
values described in Section~\ref{inner_loop_section}) that contains random values within (represented by all of the values described in Section~\ref{inner_loop_section}) that
reasonable bounds in the space. However, that still leaves the determination of which contains random values within reasonable bounds in the space. However, that still leaves
distribution function to use for the random values over each of those variables, which the determination of which distribution function to use for the random values over each
bounds to use, as well as the possibilities for any improvements to a purely random of those variables, which bounds to use, as well as the possibilities for any
search. improvements to a purely random search.
Currently, the first value set for the mission guess is that of $n$, which is the Currently, the first value set for the mission guess is that of $n$, which is the
number of sub-trajectories that each arc will be broken into for the Sims-Flanagan number of sub-trajectories that each arc will be broken into for the Sims-Flanagan
@@ -372,18 +367,18 @@
missions with more flybys. missions with more flybys.
Then, the internal components for each phase are generated. It is at this step, that Then, the internal components for each phase are generated. It is at this step, that
the mission guess generator splits the outputs into two separate outputs. The first the trajectory guess generator splits the outputs into two separate outputs. The first
is meant to be truly random, as is generally used as input for a monotonic basin is meant to be truly random, as is generally used as input for a monotonic basin
hopping algorithm. The second utilizes a Lambert's solver to determine the hopping algorithm. The second utilizes a Lambert's solver to determine the
appropriate hyperbolic velocities (both in and out) at each flyby to generate a appropriate hyperbolic velocities (both in and out) at each flyby to generate a
natural trajectory arc. For this Lambert's case, the mission guess is simply seeded natural trajectory arc. For this Lambert's case, the trajectory guess is simply seeded
with zero thrust controls and outputted to the monotonic basin hopper. The intention with zero thrust controls and outputted to the monotonic basin hopper. The intention
here is that if the time of flights are randomly chosen so as to produce a here is that if the time of flights are randomly chosen so as to produce a
trajectory that is possible with a control in the vicinity of a natural trajectory, trajectory that is possible with a control in the vicinity of a natural trajectory,
we want to be sure to find that trajectory. More detail on how this is handled is we want to be sure to find that trajectory. More detail on how this is handled is
available in Section~\ref{mbh_subsection}. available in Section~\ref{mbh_subsection}.
However, for the truly random mission guess, there are still the $v_\infty$ values However, for the truly random trajectory guess, there are still the $v_\infty$ values
and the initial thrust guesses to generate. For each of the phases, the incoming and the initial thrust guesses to generate. For each of the phases, the incoming
excess hyperbolic velocity is calculated in much the same way that the launch excess hyperbolic velocity is calculated in much the same way that the launch
velocity was calculated. However, instead of multiplying the randomly generate unit velocity was calculated. However, instead of multiplying the randomly generate unit
@@ -397,18 +392,12 @@
non-powered flyby. non-powered flyby.
From these two velocity vectors the turning angle, and thus the periapsis of the flyby, From these two velocity vectors the turning angle, and thus the periapsis of the flyby,
can then be calculated by Equation~\ref{turning_angle_eq} and the following equation: can then be calculated by Equation~\ref{turning_angle_eq} and
Equation~\ref{periapsis_eq}. If this radius of periapse is then found to be less than
\begin{equation} the minimum safe radius (currently set to the radius of the planet plus 100 kilometers),
r_p = \frac{\mu}{\vec{v}_{\infty,in} \cdot \vec{v}_{\infty,out}} \cdot \left( then the process is repeated with new random flyby velocities until a valid seed flyby
\frac{1}{\sin(\delta/2)} - 1 \right) is found. These checks are also performed each time a mission is perturbed or generated
\end{equation} by the NLP solver.
If this radius of periapse is then found to be less than the minimum safe radius
(currently set to the radius of the planet plus 100 kilometers), then the process is
repeated with new random flyby velocities until a valid seed flyby is found. These
checks are also performed each time a mission is perturbed or generated by the NLP
solver.
The final requirement then, is the thrust controls, which are actually quite simple. The final requirement then, is the thrust controls, which are actually quite simple.
Since the thrust is defined as a 3-vector of values between -1 and 1 representing some Since the thrust is defined as a 3-vector of values between -1 and 1 representing some
@@ -477,14 +466,11 @@
Because of this, the perturbation used in this implementation follows a Because of this, the perturbation used in this implementation follows a
bi-directional, long-tailed Pareto distribution generated by the following bi-directional, long-tailed Pareto distribution generated by the following
probability density function\cite{englander2014tuning}: probability density function\cite{englander2014tuning}:
\begin{equation} \begin{equation}
1 + 1 +
\left[ \frac{s}{\epsilon} \right] \cdot \left[ \frac{s}{\epsilon} \right] \cdot
\left[ \frac{\alpha - 1}{\frac{\epsilon}{\epsilon + r}^{-\alpha}} \right] \left[ \frac{\alpha - 1}{\frac{\epsilon}{\epsilon + r}^{-\alpha}} \right]
\end{equation} \end{equation}
\noindent
Where $s$ is a random array of signs (either plus one or minus one) with dimension Where $s$ is a random array of signs (either plus one or minus one) with dimension
equal to the perturbed variable and bounds of -1 and 1, $r$ is a uniformly equal to the perturbed variable and bounds of -1 and 1, $r$ is a uniformly
distributed random array with dimension equal to the perturbed variable and bounds distributed random array with dimension equal to the perturbed variable and bounds

7
LaTeX/conclusion.tex
View File

@@ -10,11 +10,8 @@
In performing this examination, two results were selected for further analysis. These In performing this examination, two results were selected for further analysis. These
results are outlined in Table~\ref{results_table}. As can be seen in the table, both results are outlined in Table~\ref{results_table}. As can be seen in the table, both
resulting trajectories have trade-offs in mission length, launch energy, fuel usage, and resulting trajectories have trade-offs in mission length, launch energy, fuel usage, and
more. However, both results show very interesting trajectories that could indicate some more. Each of these trajectories appear to be within the capabilities of existing launch
favorable possibilities for such a mission profile. Each of these trajectories should be vehicles in terms of $C_3$.
within the capabilities of existing launch vehicles in terms of $C_3$.
\section{Recommendations for Future Work}\label{improvement_section}
In the course of producing this algorithm, a large number of improvement possibilities were In the course of producing this algorithm, a large number of improvement possibilities were
noted. This work was based, in large part, on the work of Jacob Englander in a number of noted. This work was based, in large part, on the work of Jacob Englander in a number of

BIN
LaTeX/fig/EMJS_thrust_components.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 40 KiB

After

Width:  |  Height:  |  Size: 43 KiB

BIN
LaTeX/fig/EMS_thrust_components.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 38 KiB

After

Width:  |  Height:  |  Size: 42 KiB

75
LaTeX/introduction.tex
View File

@@ -1,16 +1,16 @@
\chapter{Introduction} \chapter{Introduction}
Continuous low-thrust engines utilizing technologies such as Ion propulsion, Hall thrusters, and Continuous low-thrust engines utilizing technologies such as Ion propulsion, Hall thrusters, and
others can be a powerful system in the enabling of long-range interplanetary missions with fuel others enable long-range interplanetary missions with fuel efficiencies unrivaled by those that
efficiencies unrivaled by those that employ only impulsive thrust systems. The challenge in employ only impulsive thrust systems. The challenge in utilizing these systems, then, is the
utilizing these systems, then, is the design of trajectories that effectively utilize this design of trajectories that effectively utilize this technology. Continuous thrust propulsive
technology. Continuous thrust propulsive systems tend to be particularly suited to missions systems tend to be particularly suited to missions which require very high total change in
which require very high total change in velocity ($\Delta V$) values and take place over a velocity ($\Delta V$) values and take place over a particularly long duration. Traditional
particularly long duration. Traditional impulsive thrusting techniques can achieve these changes impulsive thrusting techniques can achieve these changes in velocity, but typically have a far
in velocity, but typically have a far lower specific impulse and, as such, are much less fuel lower specific impulse and, as such, are much less fuel efficient, costing the mission valuable
efficient, costing the mission valuable financial resources that could instead be used for financial resources that could instead be used for science. Because of their inherently high
science. Because of their inherently high specific impulse (and thus efficiency), low-thrust specific impulse (and thus efficiency), low-thrust propulsion systems are well-suited to
propagation systems are well-suited to interplanetary missions. interplanetary missions.
The first attempt by NASA to use an electric ion-thruster for an interplanetary mission was the The first attempt by NASA to use an electric ion-thruster for an interplanetary mission was the
Deep Space 1 mission\cite{brophy2002}. This mission was designed to test the ``new'' technology, Deep Space 1 mission\cite{brophy2002}. This mission was designed to test the ``new'' technology,
@@ -29,16 +29,15 @@
in October 2018 and is projected to perform a flyby of Earth, two of Venus, and six of in October 2018 and is projected to perform a flyby of Earth, two of Venus, and six of
Mercury before inserting into an orbit around that planet. Mercury before inserting into an orbit around that planet.
A common theme in mission design is that there always exists a trade-off between efficiency A common theme in mission design is that there is a trade-off between efficiency (particularly
(particularly in terms of fuel use) and the time required to achieve the mission objective. Low in terms of fuel use) and the time required to achieve the mission objective. Low thrust systems
thrust systems in particular tend to produce mission profiles that sacrifice the rate of in particular tend to produce mission profiles that sacrifice the rate of convergence on the
convergence on the target state in order to achieve large increases in fuel efficiency. Often a target state in order to achieve large increases in fuel efficiency. Often a low-thrust transfer
low-thrust mission profile in Earth orbit will require multiple orbital periods to achieve the in Earth orbit will require multiple orbital periods to achieve the desired change in spacecraft
desired change in spacecraft state. Interplanetary missions, though, provide a particularly state. Interplanetary missions, though, provide a particularly useful case for continuous thrust
useful case for continuous thrust technology. The trajectory arcs in interplanetary space are technology. The trajectory arcs in interplanetary space are generally much, much longer than
generally much, much longer than orbital missions around the Earth. Because of this increase, orbital missions around the Earth. Because of this increase, even a small continuous thrust is
even a small continuous thrust is capable of producing large $\Delta V$ values over the course capable of producing large $\Delta V$ values over the course of a single trajectory arc.
of a single trajectory arc.
Another technique often leveraged by interplanetary trajectory designers is the gravity assist. Another technique often leveraged by interplanetary trajectory designers is the gravity assist.
Gravity assists utilize the inertia of a large planetary body to ``slingshot'' a spacecraft, Gravity assists utilize the inertia of a large planetary body to ``slingshot'' a spacecraft,
@@ -58,24 +57,22 @@
routine for producing unconstrained, globally optimal trajectories for realistic interplanetary routine for producing unconstrained, globally optimal trajectories for realistic interplanetary
mission development that utilizes both planetary flybys and efficient low-thrust electric mission development that utilizes both planetary flybys and efficient low-thrust electric
propulsion techniques. Similar studies have also been performed by a number of researchers propulsion techniques. Similar studies have also been performed by a number of researchers
including a team from JPL\cite{sims2006} as well as a Spanish team\cite{morante}, among several including a team from JPL\cite{sims2006}, among several others\cite{morante}.
others.
This thesis will attempt to develop an algorithm for the optimization of low-thrust enabled This thesis focuses on optimization of low-thrust enabled trajectories that use gravity assists.
trajectories for initial feasibility analysis in mission design. The algorithm will utilize The approach uses a non-linear programming solver to directly optimize a set of control thrusts
a non-linear programming solver to directly optimize a set of control thrusts for the for the user-provided flyby planets, for any provided cost function. A monotonic basin hopping
user-provided flyby planets, for any provided cost function. A monotonic basin hopping algorithm algorithm (MBH) is then employed to traverse the search space in an effort to find additional
(MBH) will then be employed to traverse the search space in an effort to find additional local local optima. This approach differs from the work produced earlier by Englander and the other
optima. This approach differs from the work produced earlier by Englander and the other teams, teams, but is largely meant to explore the feasibility of such techniques and propose a few
but is largely meant to explore the feasibility of such techniques and propose a few enhancements. The approach defined in this thesis is then used to design low thrust trajectories
enhancements. The approach defined in this thesis will then be used to investigate an example with gravity assits from the Earth to Saturn.
mission to Saturn.
This thesis will explore these concepts in a number of different sections. Section This thesis is organized as follows: Section \ref{traj_dyn} will explore the basic dynamical
\ref{traj_dyn} will explore the basic dynamical principles of trajectory design, beginning the principles of trajectory design, beginning the with fundamental system dynamics, then exploring
with fundamental system dynamics, then exploring interplanetary system dynamics and gravity interplanetary system dynamics and gravity flybys, and finally the dynamics that are specific to
flybys, and finally the dynamics that are specific to low-thrust enabled trajectories. Section low-thrust enabled trajectories. Section \ref{traj_optimization} will then discuss process of
\ref{traj_optimization} will then discuss process of optimizing spacecraft trajectories in optimizing spacecraft trajectories in general and the tool available for that. Section
general and the tool available for that. Section \ref{algorithm} will cover the implementation \ref{algorithm} will cover the implementation details of the optimization algorithm developed
details of the optimization algorithm developed for this paper. Finally, section \ref{results} for this paper. Finally, section \ref{results} will explore the results of some hypothetical
will explore the results of some hypothetical missions to Saturn. missions to Saturn.

50
LaTeX/presentation.tex Normal file
View File

@@ -0,0 +1,50 @@
\documentclass{beamer}
\usetheme{Luebeck}
\definecolor{color1}{HTML}{3A4040}
\definecolor{color2}{HTML}{F5F2F8}
\definecolor{color3}{HTML}{B65D4E}
\definecolor{color4}{HTML}{B6AD96}
\definecolor{color5}{HTML}{A96041}
\setbeamercolor*{structure}{bg=color3,fg=color3}
\setbeamercolor*{palette primary}{fg=color1,bg=color4}
\setbeamercolor*{palette secondary}{fg=color1,bg=color2}
\setbeamercolor*{palette tertiary}{fg=color1,bg=color2}
\setbeamercolor*{palette quaternary}{fg=color1,bg=color3}
\setbeamercolor{alerted text}{fg=color1,bg=color3}
\setbeamercolor{titlelike}{bg=color2,fg=color1}
\setbeamercolor*{titlelike}{bg=color2,fg=color1}
\setbeamercolor{frametitle}{bg=color1,fg=color2}
\setbeamercolor{background canvas}{bg=color2,fg=color2}
\title{Designing Optimal Low-Thrust Interplanetary Trajectories}
\subtitle{Utilizing Monotonic Basin Hopping}
\author{Richard Connor Johnstone}
\institute{University of Colorado -- Boulder}
\date{\today}
\begin{document}
\begin{frame}
\titlepage
\end{frame}
\section{Introduction}
\begin{frame} \frametitle{First Frame}
\begin{itemize}
\item Item 1
\item Item 2
\end{itemize}
\end{frame}
\section{Introduction}
\section{Introduction}
\section{Introduction}
\section{Introduction}
\end{document}

62
LaTeX/results.tex
View File

@@ -1,31 +1,29 @@
\chapter{Sample Saturn Trajectory Analysis} \label{results} \chapter{Application: Designing a Trajectory To Saturn} \label{results}
The algorithm described in this thesis is quite flexible in its design and could be used as To consider a relatively simple but representative mission design objective, a sample mission to
a tool for a mission designer on a variety of different mission types. However, to consider Saturn was investigated.
a relatively simple but representative mission design objective, a sample mission to Saturn
was investigated.
\section{Mission Constraints} \section{Mission Scenario}
The sample mission was defined to represent a general case for a near-future low-thrust The sample mission is defined to represent a general case for a near-future low-thrust
trajectory to Saturn. No constraints were placed on the flyby planets, but a number of trajectory to Saturn. No constraints are placed on the flyby planets, but a number of
constraints were placed on the algorithm to represent a realistic mission scenario. constraints were placed on the algorithm to represent a realistic mission scenario.
The first choice required by the application is one not necessarily designable to the The first choice required by the application is one not necessarily designable to the
initial mission designer (though not necessarily fixed in the design either) and is that initial mission designer (though not necessarily fixed in the design either) and is that of
of the spacecraft parameters. The application accepts as input a spacecraft object the spacecraft parameters. The application accepts as input a spacecraft object containing:
containing: the dry mass of the craft, the fuel mass at launch, the number of onboard the dry mass of the spacecraft, the fuel mass at launch, the number of onboard thrusters,
thrusters, and the specific impulse, maximum thrust and duty cycle of each thruster. and the specific impulse, maximum thrust and duty cycle of each thruster.
For this mission, the spacecraft was chosen to have a dry mass of only 200 kilograms for For this mission, the spacecraft was chosen to have a dry mass of only 200 kilograms for a
a fuel mass of 3300 kilograms. This was chosen in order to have an overall mass roughly fuel mass of 3300 kilograms. This was chosen in order to have an overall mass roughly in the
in the same zone as that of the Cassini spacecraft, which launched with 5712 kilograms same zone as that of the Cassini spacecraft, which launched with 5712 kilograms of total
of total mass, with the fuel accounting for 2978 of those kilograms\cite{cassini}. The mass, with the fuel accounting for 2978 of those kilograms\cite{cassini}. The dry mass of
dry mass of the craft was chosen to be extremely low in order to allow for a variety of the spacecraft was chosen to be extremely low in order to allow for a variety of
''successful`` missions in which the craft didn't run out of fuel. That way, the ''successful`` missions in which the spacecraft didn't run out of fuel. That way, the
delivered dry mass to Saturn could be thought of as a metric of success, without delivered dry mass to Saturn could be thought of as a metric of success, without discounting
discounting mission that may have delivered just under whatever more realistic dry mass mission that may have delivered just under whatever more realistic dry mass one might set,
one might set, in case those missions are in the vicinity of actually valid missions. in case those missions are in the vicinity of actually valid missions.
The thruster was chosen to have a specific impulse of 3200 seconds, a maximum thrust of The thruster was chosen to have a specific impulse of 3200 seconds, a maximum thrust of
250 millinewtons, and a 100\% duty cycle. This puts the thruster roughly in line with 250 millinewtons, and a 100\% duty cycle. This puts the thruster roughly in line with
@@ -121,7 +119,7 @@
The mission begins in late June of 2024 and proceeds first to an initial gravity assist with The mission begins in late June of 2024 and proceeds first to an initial gravity assist with
Mars after three and one half years to rendezvous in mid-December 2027. Unfortunately, the Mars after three and one half years to rendezvous in mid-December 2027. Unfortunately, the
launch energy required to effectively used the gravity assist with Mars at this time is launch energy required to effectively use the gravity assist with Mars at this time is
quite high. The $C_3$ value was found to be $60.4102 \frac{\text{km}^2}{\text{s}^2}$. However, quite high. The $C_3$ value was found to be $60.4102 \frac{\text{km}^2}{\text{s}^2}$. However,
for this phase, the thrust magnitudes are quite low, raising slowly only as the spacecraft for this phase, the thrust magnitudes are quite low, raising slowly only as the spacecraft
approaches Mars, allowing for a nearly-natural trajectory to Mars rendezvous. Note also that approaches Mars, allowing for a nearly-natural trajectory to Mars rendezvous. Note also that
@@ -130,7 +128,7 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.9\textwidth]{fig/EMS_plot} \includegraphics[width=0.9\textwidth]{LaTeX/fig/EMS_plot}
\caption{Depictions of the faster Earth-Mars-Saturn trajectory found by the \caption{Depictions of the faster Earth-Mars-Saturn trajectory found by the
algorithm to be most efficient; planetary ephemeris arcs are shown during the phase algorithm to be most efficient; planetary ephemeris arcs are shown during the phase
in which the spacecraft approached them} in which the spacecraft approached them}
@@ -139,7 +137,7 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.9\textwidth]{fig/EMS_plot_noplanets} \includegraphics[width=0.9\textwidth]{LaTeX/fig/EMS_plot_noplanets}
\caption{Another depiction of the EMS trajectory, without the planetary ephemeris \caption{Another depiction of the EMS trajectory, without the planetary ephemeris
arcs} arcs}
\label{ems_nop} \label{ems_nop}
@@ -158,14 +156,14 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.9\textwidth]{fig/EMS_thrust_mag} \includegraphics[width=0.9\textwidth]{LaTeX/fig/EMS_thrust_mag}
\caption{The magnitude of the unit thrust vector over time for the EMS trajectory} \caption{The magnitude of the unit thrust vector over time for the EMS trajectory}
\label{ems_mag} \label{ems_mag}
\end{figure} \end{figure}
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.9\textwidth]{fig/EMS_thrust_components} \includegraphics[width=0.9\textwidth]{LaTeX/fig/EMS_thrust_components}
\caption{The inertial x, y, and z components of the unit thrust vector over time for \caption{The inertial x, y, and z components of the unit thrust vector over time for
the EMS trajectory} the EMS trajectory}
\label{ems_components} \label{ems_components}
@@ -214,7 +212,7 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.9\textwidth]{fig/EMJS_plot} \includegraphics[width=0.9\textwidth]{LaTeX/fig/EMJS_plot}
\caption{Depictions of the slower Earth-Mars-Jupiter-Saturn trajectory found by the \caption{Depictions of the slower Earth-Mars-Jupiter-Saturn trajectory found by the
algorithm to be the second most efficient; planetary ephemeris arcs are shown during algorithm to be the second most efficient; planetary ephemeris arcs are shown during
the phase in which the spacecraft approached them} the phase in which the spacecraft approached them}
@@ -223,7 +221,7 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.9\textwidth]{fig/EMJS_plot_noplanets} \includegraphics[width=0.9\textwidth]{LaTeX/fig/EMJS_plot_noplanets}
\caption{Another depiction of the EMJS trajectory, without the planetary ephemeris \caption{Another depiction of the EMJS trajectory, without the planetary ephemeris
arcs} arcs}
\label{emjs_nop} \label{emjs_nop}
@@ -238,14 +236,14 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.9\textwidth]{fig/EMJS_thrust_mag} \includegraphics[width=0.9\textwidth]{LaTeX/fig/EMJS_thrust_mag}
\caption{The magnitude of the unit thrust vector over time for the EMJS trajectory} \caption{The magnitude of the unit thrust vector over time for the EMJS trajectory}
\label{emjs_mag} \label{emjs_mag}
\end{figure} \end{figure}
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.9\textwidth]{fig/EMJS_thrust_components} \includegraphics[width=0.9\textwidth]{LaTeX/fig/EMJS_thrust_components}
\caption{The inertial x, y, and z components of the unit thrust vector over time for \caption{The inertial x, y, and z components of the unit thrust vector over time for
the EMJS trajectory} the EMJS trajectory}
\label{emjs_components} \label{emjs_components}
@@ -306,8 +304,8 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{fig/c3} \includegraphics[width=\textwidth]{LaTeX/fig/c3}
\caption{Plot of Delta IV and Atlas V launch vehicle capabilities as they relate to \caption{Plot of Delta IV and Atlas V launch vehicle capabilities as they relate to
payload mass \cite{c3capabilities} from a source from 2007} payload mass \cite{c3capabilities} from Vardaxis, et al, 2007 }
\label{c3} \label{c3}
\end{figure} \end{figure}

25
LaTeX/thesis.tex
View File

@@ -1,4 +1,4 @@
\documentclass[defaultstyle,11pt]{thesis} \documentclass[defaultstyle,11pt]{LaTeX/thesis}
\usepackage{graphicx} \usepackage{graphicx}
\usepackage{amssymb} \usepackage{amssymb}
@@ -22,14 +22,15 @@
Much work has been performed recently to utilize the increasingly viable technology of Much work has been performed recently to utilize the increasingly viable technology of
low-thrust electric propulsion systems on missions of interplanetary scope. This thesis analyzes low-thrust electric propulsion systems on missions of interplanetary scope. This thesis analyzes
a technique for the initial analysis of feasibility of utilizing a combination of low-thrust a technique for designing trajectories for spacecraft with a low-thrust propulsion system that
propulsion systems and natural gravity flybys for missions to the outer planets. First, a method also use natural gravity flybys for missions to the outer planets. Often, the goal is to find
for finding local optima by utilizing an interior-point linesearch algorithm to directly feasible solutions that also minimize propellant mass requirements. First, locally optimal
optimize the entire trajectory as a Non-Linear Programming problem is presented. Then, a solutions are constructed by using an interior-point linesearch algorithm, along with multiple
Monotonic Basin Hopping algorithm is utilized to traverse the search space, improve the local shooting techniques for optimization. Then, Monotonic Basin Hopping is utilized to traverse the
optima determined by the internal optimizer, and determine the global optima. This allows for a search space, improve the local optima determined by the internal optimizer, and determine the
medium-fidelity, fully automated global optimization of the low thrust controls and flyby global optima. This approach allows for a medium-fidelity, fully automated global optimization
parameters for a given mission objective. of the low thrust controls and flyby parameters for a given target destination. As an
application of this method, two sample trajectories to Saturn are analyzed.
} }
@@ -60,7 +61,7 @@
\begin{document} \begin{document}
\input macros.tex \input LaTeX/macros.tex
\input LaTeX/introduction.tex \input LaTeX/introduction.tex
@@ -76,8 +77,6 @@
\bibliographystyle{plain} \bibliographystyle{plain}
\nocite{*} \nocite{*}
\bibliography{thesis} \bibliography{LaTeX/thesis}
\appendix
\end{document} \end{document}

310
LaTeX/trajectory_design.tex
View File

@@ -16,28 +16,28 @@
very high-fidelity force models that account for aerodynamic pressure, solar radiation very high-fidelity force models that account for aerodynamic pressure, solar radiation
pressure, multi-body effects, and other forces may be too time intensive for a pressure, multi-body effects, and other forces may be too time intensive for a
particular application. Initial surveys of the solution space often don't require such particular application. Initial surveys of the solution space often don't require such
complex models in order to gain valuable insight. complex models in order to gain valuable preliminary insight.
Therefore, a common approach (and the one utilized in this implementation) is to first A common approach (and the one utilized in this implementation) is to first use a
use a lower-fidelity dynamical model that captures only the gravitational force due to lower-fidelity dynamical model that captures only the gravitational force due to the
the primary body around which the spacecraft is orbiting. This approach can provide an primary body around which the spacecraft is orbiting. This approach can provide an
excellent low-to-medium fidelity model that is useful as an underlying model in an excellent low-to-medium fidelity model that is useful as an underlying model in an
algorithm for quickly categorizing a search space for initial mission feasibility algorithm for quickly categorizing a search space for initial mission feasibility
explorations. explorations.
In order to explore the Two Body Problem, we must first examine the full set of In order to explore the Two Body Problem, we must first examine the full set of
assumptions associated with the force model\cite{vallado2001fundamentals}. Firstly, we assumptions associated with the force model\cite{vallado2001fundamentals}. Firstly, we
are only concerned with the nominative two bodies: the spacecraft and the planetary body are only concerned with the gravitational influence between the nominative two bodies:
around which it is orbiting. Secondly, both of these bodies are modeled as point masses the spacecraft and the planetary body around which it is orbiting. Secondly, both of
with constant mass. This removes the need to account for non-uniform densities and these bodies are modeled as point masses with constant mass. This removes the need to
asymmetry. Finally, for convenience in notation at the end, we'll also assume that the account for non-uniform densities and asymmetry. Finally, for convenience in notation at
mass of the spacecraft ($m_2$) is much much smaller than the mass of the planetary body the end, we'll also assume that the mass of the spacecraft ($m_2$) is much much smaller
($m_1$) and enough so as to be considered negligible. The only force acting on this than the mass of the planetary body ($m_1$) and enough so as to be considered
system is then the force of gravity that the primary body enacts upon the secondary. negligible.
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.65\textwidth]{fig/2bp} \includegraphics[width=0.65\textwidth]{LaTeX/fig/2bp}
\caption{Figure representing the positions of the bodies relative to each other and \caption{Figure representing the positions of the bodies relative to each other and
the center of mass in the two body problem} the center of mass in the two body problem}
\label{2bp_fig} \label{2bp_fig}
@@ -45,7 +45,6 @@
Under these assumptions, the force acting on the body due to the law of universal Under these assumptions, the force acting on the body due to the law of universal
gravitation is: gravitation is:
\begin{align} \begin{align}
F_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\ F_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\
F_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|} F_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|}
@@ -53,7 +52,6 @@
And by Newton's second law (force is the product of mass and acceleration), we can And by Newton's second law (force is the product of mass and acceleration), we can
derive the following differential equations for $r_1$ and $r_2$: derive the following differential equations for $r_1$ and $r_2$:
\begin{align} \begin{align}
m_2 \ddot{\vec{r}}_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\ m_2 \ddot{\vec{r}}_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\
m_1 \ddot{\vec{r}}_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|} m_1 \ddot{\vec{r}}_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|}
@@ -65,7 +63,6 @@
inertial frame. $G$ is the universal gravitational parameter, $m_1$ is the mass of the inertial frame. $G$ is the universal gravitational parameter, $m_1$ is the mass of the
planetary body, and $m_2$ is the mass of the spacecraft. From these equations, we can planetary body, and $m_2$ is the mass of the spacecraft. From these equations, we can
then determine the acceleration of the spacecraft relative to the planet: then determine the acceleration of the spacecraft relative to the planet:
\begin{equation} \begin{equation}
\ddot{\vec{r}} = \ddot{\vec{r}}_2 - \ddot{\vec{r}}_1 = \ddot{\vec{r}} = \ddot{\vec{r}}_2 - \ddot{\vec{r}}_1 =
- \frac{G \left( m_1 + m_2 \right)}{r^2} \frac{\vec{r}}{\left| r \right|} - \frac{G \left( m_1 + m_2 \right)}{r^2} \frac{\vec{r}}{\left| r \right|}
@@ -76,27 +73,19 @@
negligible $m_2$ term. We can also introduce, for convenience, a gravitational parameter negligible $m_2$ term. We can also introduce, for convenience, a gravitational parameter
$\mu$ which represents the gravity constant for the system about the center of motion $\mu$ which represents the gravity constant for the system about the center of motion
($\mu = G (m_1 + m_2) \approx G m_1$). Doing so and simplifying produces: ($\mu = G (m_1 + m_2) \approx G m_1$). Doing so and simplifying produces:
\begin{equation} \begin{equation}
\ddot{\vec{r}} = - \frac{\mu}{r^2} \hat{r} \ddot{\vec{r}} = - \frac{\mu}{r^2} \hat{r}
\end{equation} \end{equation}
We may also wish to utilize the total orbital energy for a spacecraft within this model.
Since the spacecraft is acting only under the gravitational influence of the planet and Since the spacecraft is acting only under the gravitational influence of the planet and
no other forces, we can define the total specific mechanical energy as: no other forces, we can define the total specific mechanical energy as
\cite{vallado2001fundamentals}: \cite{vallado2001fundamentals}:
\begin{equation} \label{energy} \begin{equation} \label{energy}
\xi = \frac{v^2}{2} - \frac{\mu}{r} \xi = \frac{v^2}{2} - \frac{\mu}{r}
\end{equation} \end{equation}
\noindent
Where the first term represents the kinetic energy of the spacecraft and the second term Where the first term represents the kinetic energy of the spacecraft and the second term
represents the gravitational potential energy. represents the gravitational potential energy.
\subsection{Kepler's Laws}
Now that we've fully qualified the forces acting within the Two Body Problem, we can concern Now that we've fully qualified the forces acting within the Two Body Problem, we can concern
ourselves with more practical applications of it as a force model. It should be noted, ourselves with more practical applications of it as a force model. It should be noted,
firstly, that the spacecraft's position and velocity (given an initial position and velocity firstly, that the spacecraft's position and velocity (given an initial position and velocity
@@ -105,6 +94,8 @@
one-dimensional equations (one for each component of the three-dimensional space) and one-dimensional equations (one for each component of the three-dimensional space) and
three unknowns (the three components of the second derivative of the position). three unknowns (the three components of the second derivative of the position).
\subsection{Kepler's Laws}
In the early 1600s, Johannes Kepler produced just such a solution, by taking advantages of 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\cite{murray1999solar}: what is also known as ``Kepler's Laws'' which are\cite{murray1999solar}:
@@ -113,68 +104,61 @@
expanded to any orbit by re-wording as ``all orbital paths follow a conic section expanded to any orbit by re-wording as ``all orbital paths follow a conic section
(circle, ellipse, parabola, or hyperbola) with a primary mass at one of the foci''. (circle, ellipse, parabola, or hyperbola) with a primary mass at one of the foci''.
Specifically the path of the orbit follows the trajectory equation: The conic trajectory equation explains this observation and offers a description
of the path as:
\begin{equation} \begin{equation}
r = \frac{\sfrac{h^2}{\mu}}{1 + e \cos(\theta)} r = \frac{\sfrac{h^2}{\mu}}{1 + e \cos(\theta)}
\end{equation} \end{equation}
where $h$ is the angular momentum of the satellite, $e$ is the
Where $h$ is the angular momentum of the satellite, $e$ is the
eccentricity of the orbit, and $\theta$ is the true anomaly, or simply eccentricity of the orbit, and $\theta$ is the true anomaly, or simply
the angular distance the satellite has traversed along the orbit path. the angular distance the satellite has traversed along the orbit path from
periapsis.
\item The area swept out by the imaginary line connecting the primary and secondary \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 bodies increases linearly with respect to time. This implies that the magnitude of the
orbital speed is not constant. For the moment, we'll just take this orbital speed is not constant. For the moment, we'll just take this
value to be a constant: value to be a constant:
\begin{equation}\label{swept} \begin{equation}\label{swept}
\frac{\Delta t}{T} = \frac{k}{\pi a b} \frac{\Delta t}{T} = \frac{k}{\pi a b}
\end{equation} \end{equation}
where $k$ is the constant value, $a$ and $b$ are the semi-major and
Where $k$ is the constant value, $a$ and $b$ are the semi-major and
semi-minor axis of the conic section, and $T$ is the period. In the semi-minor axis of the conic section, and $T$ is the period. In the
following section, we'll derive the value for $k$. following section, we'll derive the value for $k$.
\item The square of the orbital period is proportional to the cube of the semi-major \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: axis of the orbit, regardless of eccentricity. For an elliptical orbit this
observation connects to the following known expression for the orbit period:
\begin{equation} \begin{equation}
T = 2 \pi \sqrt{\frac{a^3}{\mu}} T = 2 \pi \sqrt{\frac{a^3}{\mu}}
\end{equation} \end{equation}
where $T$ is the period and $a$ is the semi-major axis.
Where $T$ is the period and $a$ is the semi-major axis.
\end{enumerate} \end{enumerate}
\subsection{Kepler's Equation} \subsection{Kepler's Equation}
Kepler was able to produce an equation to represent the angular displacement of an 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 orbiting body around a primary body as a function of time, which we'll derive now for
the elliptical case\cite{vallado2001fundamentals}. Since the total area of an ellipse is the elliptical case\cite{vallado2001fundamentals}. Because the total area of an ellipse
the product of $\pi$, the semi-major axis, and the semi-minor axis ($\pi a b$), we can is the product of $\pi$, the semi-major axis, and the semi-minor axis ($\pi a b$), we
relate (by Kepler's second law) the area swept out by an orbit as a function of time, as can relate (by Kepler's second law) the area swept out by an orbit as a function of
we did in Equation~\ref{swept}. This leaves just one unknown variable $k$, which we can time, as we did in Equation~\ref{swept}. This leaves just one unknown variable $k$,
determine through use of the geometric auxiliary circle, which is a circle with radius which we can determine through use of the geometric auxiliary circle, which is a circle
equal to the ellipse's semi-major axis and center directly between the two foci, as in with radius equal to the ellipse's semi-major axis and center directly between the two
Figure~\ref{aux_circ}. foci, as in Figure~\ref{aux_circ}.
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.8\textwidth]{fig/kepler} \includegraphics[width=0.8\textwidth]{LaTeX/fig/kepler}
\caption{Geometric Representation of Auxiliary Circle}\label{aux_circ} \caption{Geometric representation of auxiliary circle}\label{aux_circ}
\end{figure} \end{figure}
In order to find the area swept by the spacecraft\cite{vallado2001fundamentals}, $k$, we In order to find the area swept by the spacecraft\cite{vallado2001fundamentals}, $k$, we
can take advantage of the fact that that area is the triangle $k_1$ subtracted from the can take advantage of the fact that that area is the triangle $k_1$ subtracted from the
elliptical segment $PCB$: elliptical segment $PCB$:
\begin{equation}\label{areas_eq} \begin{equation}\label{areas_eq}
k = area(seg_{PCB}) - area(k_1) k = area(seg_{PCB}) - area(k_1)
\end{equation} \end{equation}
\noindent
Where the area of the triangle $k_1$ can be found easily using geometric formulae: Where the area of the triangle $k_1$ can be found easily using geometric formulae:
\begin{align} \begin{align}
area(k_1) &= \frac{1}{2} \left( ae - a \cos E \right) \left( \frac{b}{a} a \sin E \right) \\ 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) &= \frac{ab}{2} \left(e \sin E - \cos E \sin E \right)
@@ -186,7 +170,6 @@
can find the area for the elliptical segment $PCB$ by first finding the circular segment can find the area for the elliptical segment $PCB$ by first finding the circular segment
$POB'$, subtracting the triangle $COB'$, then applying the fact that an ellipse is $POB'$, subtracting the triangle $COB'$, then applying the fact that an ellipse is
merely a vertical scaling of a circle by the amount $\frac{b}{a}$. merely a vertical scaling of a circle by the amount $\frac{b}{a}$.
\begin{align} \begin{align}
area(PCB) &= \frac{b}{a} \left( area(POB') - area(COB') \right) \\ area(PCB) &= \frac{b}{a} \left( area(POB') - area(COB') \right) \\
&= \frac{b}{a} \left( \frac{a^2 E}{2} - \frac{1}{2} \left( a \cos E \right) &= \frac{b}{a} \left( \frac{a^2 E}{2} - \frac{1}{2} \left( a \cos E \right)
@@ -197,26 +180,20 @@
By substituting the two areas back into Equation~\ref{areas_eq} we can get the $k$ area By substituting the two areas back into Equation~\ref{areas_eq} we can get the $k$ area
swept out by the spacecraft: swept out by the spacecraft:
\begin{equation} \begin{equation}
k = \frac{ab}{2} \left( E - e \sin E \right) k = \frac{ab}{2} \left( E - e \sin E \right)
\end{equation} \end{equation}
Which we can then substitute back into the equation for the swept area as a function of Which we can then substitute back into the equation for the swept area as a function of
time (Equation~\ref{swept}) for period of time since the spacecraft left periapsis: time (Equation~\ref{swept}) for period of time since the spacecraft left periapsis:
\begin{equation} \begin{equation}
\frac{\Delta t}{T} = \frac{t_2 - t_{peri}}{T} = \frac{E - e \sin E}{2 \pi} \frac{\Delta t}{T} = \frac{t_2 - t_{peri}}{T} = \frac{E - e \sin E}{2 \pi}
\end{equation} \end{equation}
Which is, effectively, Kepler's equation. It is commonly known by a different form: Which is, effectively, Kepler's equation. It is commonly known by a different form:
\begin{equation} \begin{equation}
M = \sqrt{\frac{\mu}{a^3}} \Delta t = E - e \sin E M = \sqrt{\frac{\mu}{a^3}} \Delta t = E - e \sin E
\end{equation} \end{equation}
where we've defined the mean anomaly as $M$ and used the fact that $T =
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
\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 ($E$) which can be related to spacecraft position, and time, but we still need a useful
algorithm for solving this equation in order to use this equation to propagate a algorithm for solving this equation in order to use this equation to propagate a
spacecraft. spacecraft.
@@ -224,34 +201,25 @@
\subsection{LaGuerre-Conway Algorithm}\label{laguerre} \subsection{LaGuerre-Conway Algorithm}\label{laguerre}
For this thesis, the algorithm used to solve Kepler's equation was the general numeric For this thesis, the algorithm used to solve Kepler's equation was the general numeric
root-finding scheme first developed by LaGuerre in the 1800s and first applied to root-finding scheme first developed by LaGuerre in the 1800s and first applied to Kepler's
Kepler's equation by Bruce Conway in 1985\cite{laguerre_conway}. In his paper, Conway equation by Bruce Conway in 1985\cite{laguerre_conway}. In his paper, Conway makes a
makes a compelling argument for utilizing the less common LaGuerre method over higher compelling argument for utilizing the less common LaGuerre method over higher order Newton
order Newton or Newton-Raphson methods. or Newton-Raphson methods. The Newton-Raphson methods, while found to generally have quite
impressive convergence rates (generally successfully solving Kepler's equation correctly
The Newton-Raphson methods, while found to generally have quite impressive convergence within 5 iterations), were prone to failures in convergence given certain specific initial
rates (generally successfully solving Kepler's equation correctly within 5 iterations), conditions. Therefore LaGuerre's algorithm is proposed as an alternative.
were prone to failures in convergence given certain specific initial conditions.
Therefore LaGuerre's algorithm is proposed as an alternative.
The algorithm can be relatively easily derived by examining the polynomial equation with
$m$ roots:
The algorithm can be derived by examining the polynomial equation with $m$ roots:
\begin{equation} \begin{equation}
g(x) = (x - x_1) (x - x_2) ... ( x - x_m) g(x) = (x - x_1) (x - x_2) ... ( x - x_m)
\end{equation} \end{equation}
\noindent
We can then generate some useful convenience functions as: We can then generate some useful convenience functions as:
\begin{align} \begin{align}
\ln|g(x)| &= \ln|(x - x_1)| + \ln|(x - x_2)| + ... + \ln|( x - x_m)| \\ \ln|g(x)| &= \ln|(x - x_1)| + \ln|(x - x_2)| + ... + \ln|( x - x_m)| \\
\frac{d\ln|g(x)|}{dx} &= \frac{1}{x - x_1} + \frac{1}{x - x_2} + ... + \frac{1}{x - \frac{d\ln|g(x)|}{dx} &= \frac{1}{x - x_1} + \frac{1}{x - x_2} + ... + \frac{1}{x -
x_m} = G_1(x) x_m} = G_1(x)
\end{align} \end{align}
and and
\begin{align} \begin{align}
\frac{-d^2\ln|g(x)|}{dx^2} &= \frac{1}{(x - x_1)^2} + \frac{1}{(x - x_2)^2} + ... + \frac{-d^2\ln|g(x)|}{dx^2} &= \frac{1}{(x - x_1)^2} + \frac{1}{(x - x_2)^2} + ... +
\frac{1}{(x - x_m)^2} = G_2(x) \frac{1}{(x - x_m)^2} = G_2(x)
@@ -259,42 +227,32 @@
Now we define the targeted root as $x_1$ and make the approximation that all of the Now we define the targeted root as $x_1$ and make the approximation that all of the
other roots are equidistant from the targeted root, which means: other roots are equidistant from the targeted root, which means:
\begin{equation} \begin{equation}
x - x_i = b, i=2,3,...,m x - x_i = b, i=2,3,...,m
\end{equation} \end{equation}
\noindent
We can then rewrite $G_1$ and $G_2$ as: We can then rewrite $G_1$ and $G_2$ as:
\begin{align} \begin{align}
G_1 &= \frac{1}{a} + \frac{n-1}{b} \\ G_1 &= \frac{1}{a} + \frac{n-1}{b} \\
G_2 &= \frac{1}{a^2} + \frac{n-1}{b^2} G_2 &= \frac{1}{a^2} + \frac{n-1}{b^2}
\end{align} \end{align}
\noindent
Which may be solved for $a$ in terms of $G_1$, $G_2$: Which may be solved for $a$ in terms of $G_1$, $G_2$:
\begin{equation} \begin{equation}
a = \frac{n}{G_1 \pm \sqrt{(n-1)(nG_2 - G_1^2)}} a = \frac{n}{G_1 \pm \sqrt{(n-1)(nG_2 - G_1^2)}}
\end{equation} \end{equation}
\noindent
With corresponding iteration function: With corresponding iteration function:
\begin{equation} \begin{equation}
x_{i+1} = x_i - \frac{n g(x_i)}{g'(x_i) \pm \sqrt{(n-1)^2 f'(x_i)^2 - n (n-1) f(x_i) x_{i+1} = x_i - \frac{n g(x_i)}{g'(x_i) \pm \sqrt{(n-1)^2 f'(x_i)^2 - n (n-1) f(x_i)
f''(x_i)}} f''(x_i)}}
\end{equation} \end{equation}
This iteration scheme can be shown to be globally convergent, regardless of the initial This iteration scheme can be shown to be globally convergent, regardless of the initial
guess. More relevantly, Conway also showed that the application of this method to guess. Conway also showed that the application of this method to Kepler's equation was shown
Kepler's equation was shown to converge with similar speed to many of the best common to converge with similar speed to many of the best common higher order Newton-Raphson
higher order Newton-Raphson solvers. However, LaGuerre's method was also found to be solvers. However, LaGuerre's method was also found to be incredibly robust, converging to
incredibly robust, converging to the correct value for every one of Conway's 500,000 the correct value for every one of Conway's 500,000 tests. Because of this robustness, it is
tests. Because of this robustness, it is very useful for propagating spacecraft states. useful for solving Kepler's equation.
\section{Interplanetary Considerations}\label{interplanetary} \section{Interplanetary Trajectories}\label{interplanetary}
In interplanetary travel, the primary body most responsible for gravitational forces might 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 fact, at some be a number of different bodies, dependent on the phase of the mission. In fact, at some
@@ -339,21 +297,22 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.8\textwidth]{fig/patched_conics} \includegraphics[width=0.8\textwidth]{LaTeX/fig/patched_conics}
\caption{Patched Conics Example Figure} \caption{Patched Conics Example Figure}
\label{patched_conics_fig} \label{patched_conics_fig}
\end{figure} \end{figure}
This effectively breaks the trajectory into a series of arcs each governed by a distinct This effectively breaks the trajectory into a series of arcs each governed by a distinct
Two-Body problem patched together by distinct transition points. These transition points Two-Body problem patched together by distinct transition points. These transition points
occur along the spheres of influence of the planets nearest to the spacecraft. occur along the spheres of influence of the planets nearest to the spacecraft. A
conceptual example of this process, labeled the method of patched conics, appears in
Figure~\ref{patched_conics_fig}.
Therefore, we must understand how to convert our spacecraft's state from the Sun frame Therefore, we must understand how to convert our spacecraft's state from the Sun frame
to the planetary frame as it crosses this boundary. An elliptical orbit about the sun to the planetary frame as it crosses this boundary. An elliptical orbit about the sun
will have enough orbital energy to represent a hyperbolic orbit around the planet. So we will have enough orbital energy to represent a hyperbolic orbit around the planet. So we
first need to determine the velocity of the spacecraft relative to the planet as it first need to determine the velocity of the spacecraft relative to the planet as it
crosses the SOI, which we can determine by subtraction \cite{vallado2001fundamentals}: crosses the SOI, which we can determine by subtraction \cite{vallado2001fundamentals}:
\begin{equation} \begin{equation}
\vec{v}_{sc/p} = \vec{v}_{sc/sun} - \vec{v}_{planet/sun} \vec{v}_{sc/p} = \vec{v}_{sc/sun} - \vec{v}_{planet/sun}
\end{equation} \end{equation}
@@ -361,24 +320,22 @@
Since the orbit around the planet is hyperbolic, in order to characterize the hyperbola Since the orbit around the planet is hyperbolic, in order to characterize the hyperbola
we must determine the velocity of the spacecraft when it has infinite distance relative we must determine the velocity of the spacecraft when it has infinite distance relative
to the planet. Since this never occurs, a further approximation is made that the to the planet. Since this never occurs, a further approximation is made that the
velocity that the spacecraft has (relative to the planet) as it crosses the SOI can be velocity of the spacecraft (relative to the planet) as it crosses the SOI can be modeled
modeled as the $\vec{v}_\infty$ of that hyperbolic arc. as the $\vec{v}_\infty$ of that hyperbolic arc.
As an example, we may wish to determine the velocity relative to the planet that the As an example, we may wish to determine the velocity relative to the planet that the
spacecraft has at the periapsis of its hyperbolic trajectory during the flyby. This spacecraft has at the periapsis of its hyperbolic trajectory during the flyby. This
could be useful, perhaps, for sizing the $\Delta V<$ required during the insertion stage could be useful, perhaps, for sizing the $\Delta V$ required during the insertion stage
of the mission if the spacecraft is intended to be captured into an elliptical orbit of the mission if the spacecraft is intended to be captured into an elliptical orbit
around its target planet. For a given incoming hyperbolic $\vec{v}_\infty$, we can first around its target planet. For a given incoming hyperbolic $\vec{v}_\infty$, we can first
determine the specific mechanical energy of the hyperbola at infinite distance by using determine the specific mechanical energy of the hyperbola at infinite distance by using
Equation~\ref{energy}: Equation~\ref{energy}:
\begin{equation} \begin{equation}
\xi = \frac{v^2}{2} - \frac{\mu}{r} = \frac{v_\infty^2}{2} \xi = \frac{v^2}{2} - \frac{\mu}{r} = \frac{v_\infty^2}{2}
\end{equation} \end{equation}
We can then leverage the conservation of energy to determine the velocity at a We can then leverage the conservation of energy to determine the velocity at a
particular point, $r_{ins}$: particular point, $r_{ins}$:
\begin{align} \begin{align}
\xi_{ins} &= \frac{v_{ins}^2}{2} - \frac{\mu}{r_{ins}} \\ \xi_{ins} &= \frac{v_{ins}^2}{2} - \frac{\mu}{r_{ins}} \\
\xi_{ins} &= \xi_\infty = \frac{v_\infty^2}{2} \\ \xi_{ins} &= \xi_\infty = \frac{v_\infty^2}{2} \\
@@ -387,31 +344,29 @@
\subsection{Launch Considerations} \subsection{Launch Considerations}
Generally speaking, an interplanetary mission begins with launch. For a satellite of For a satellite of given size, a certain amount of orbital energy can be imparted to the
given size, a certain amount of orbital energy can be imparted to the satellite by the satellite by the launch vehicle. In practice, this value, for a particular mission, is
launch vehicle. In practice, this value, for a particular mission, is actually actually determined as a parameter of the mission trajectory to be optimized. The excess
determined as a parameter of the mission trajectory to be optimized. The excess velocity velocity at infinity of the hyperbolic orbit of the spacecraft that leaves the Earth can
at infinity of the hyperbolic orbit of the spacecraft that leaves the Earth can be used be used to derive the launch energy. This is usually qualified as the quantity $C_3$,
to derive the launch energy. This is usually qualified as the quantity $C_3$, which is which is actually double the kinetic orbital energy with respect to the Sun, or simply
actually double the kinetic orbital energy with respect to the Sun, or simply the square the square of the excess hyperbolic velocity at infinity\cite{wie1998space}.
of the excess hyperbolic velocity at infinity\cite{wie1998space}.
This algorithm will assume that the initial trajectory at the beginning of the mission This algorithm will assume that the initial trajectory at the beginning of the mission
will be some hyperbolic orbit with velocity enough to leave the Earth. That initial will be some hyperbolic orbit with velocity enough to leave the Earth. That initial
$v_\infty$ will be used as a tunable parameter in the NLP solver. This allows the $v_\infty$ will be used as a tunable parameter in the optimization routine. This allows
mission designer to include the launch $C_3$ in the cost function and, hopefully, the mission designer to include the launch $C_3$ in the cost function and, hopefully,
determine the mission trajectory that includes the least initial launch energy. This can determine the mission trajectory that includes the least initial launch energy. This can
then be fed back into a mass-$C_3$ curve for prospective launch providers to determine then be fed back into a mass-$C_3$ curve for prospective launch providers to determine
what the maximum mass any launch provider is capable of imparting that specific $C_3$ what the maximum mass any launch provider is capable of imparting that specific $C_3$
to. to.
A similar approach is taken at the end of the mission. This algorithm doesn't attempt to A similar approach is taken at the end of the trajectory. This algorithm doesn't attempt
exactly match the velocity of the planet at the end of the mission. Instead, the excess to exactly match the velocity of the planet. Instead, the excess hyperbolic velocity is
hyperbolic velocity is also treated as a parameter that can be minimized by the cost also treated as a parameter that can be minimized by the cost function. If a trajectory
function. If a mission is to then end in insertion, a portion of the mass budget can is to then end in insertion, a portion of the mass budget can then be used for an
then be used for an impulsive thrust engine, which can provide a final insertion burn at impulsive thrust engine, which can provide a final insertion burn. This approach also
the end of the mission. This approach also allows flexibility for missions that might allows flexibility for missions that might end in a flyby rather than insertion.
end in a flyby rather than insertion.
\subsection{Gravity Assist Maneuvers} \subsection{Gravity Assist Maneuvers}
@@ -441,8 +396,8 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.8\textwidth]{fig/flyby} \includegraphics[width=0.8\textwidth]{LaTeX/fig/flyby}
\caption{Visualization of velocity changes during a gravity assist} \caption{Velocity changes during a gravity assist}
\label{grav_assist_fig} \label{grav_assist_fig}
\end{figure} \end{figure}
@@ -452,7 +407,7 @@
turning angle of this bend. In doing so, one can effectively achieve a (restricted) free turning angle of this bend. In doing so, one can effectively achieve a (restricted) free
impulsive thrust event. impulsive thrust event.
\subsection{Flyby Periapsis} \subsection{Flyby Periapsis Altitude}
Now that we understand gravity assists, the natural question is then how to leverage Now that we understand gravity assists, the natural question is then how to leverage
them for achieving certain velocity changes\cite{cho2017b}. But first, we must consider them for achieving certain velocity changes\cite{cho2017b}. But first, we must consider
@@ -461,7 +416,6 @@
mentioned in the previous section, given an excess hyperbolic velocity entering the mentioned in the previous section, given an excess hyperbolic velocity entering the
planet's sphere of influence ($\vec{v}_{\infty, in}$) and a target excess hyperbolic planet's sphere of influence ($\vec{v}_{\infty, in}$) and a target excess hyperbolic
velocity as the spacecraft leaves the sphere of influence ($\vec{v}_{\infty, out}$): velocity as the spacecraft leaves the sphere of influence ($\vec{v}_{\infty, out}$):
\begin{equation}\label{turning_angle_eq} \begin{equation}\label{turning_angle_eq}
\delta = \arccos \left( \frac{\vec{v}_{\infty,in} \cdot \delta = \arccos \left( \frac{\vec{v}_{\infty,in} \cdot
\vec{v}_{\infty,out}}{|\vec{v}_{\infty,in}| |\vec{v}_{\infty,out}|} \right) \vec{v}_{\infty,out}}{|\vec{v}_{\infty,in}| |\vec{v}_{\infty,out}|} \right)
@@ -471,12 +425,10 @@
that we must target in order to achieve the required turning angle. The periapsis of the that we must target in order to achieve the required turning angle. The periapsis of the
flyby, however, can provide a useful check on what turning angles are possible for a flyby, however, can provide a useful check on what turning angles are possible for a
given flyby, since the periapsis: given flyby, since the periapsis:
\begin{equation}\label{periapsis_eq}
\begin{equation}
r_p = \frac{\mu}{v_\infty^2} \left[ \frac{1}{\sin\left(\frac{\delta}{2}\right)} - 1 \right] r_p = \frac{\mu}{v_\infty^2} \left[ \frac{1}{\sin\left(\frac{\delta}{2}\right)} - 1 \right]
\end{equation} \end{equation}
cannot be lower than some safe value that accounts for the radius of the planet and
Cannot be lower than some safe value that accounts for the radius of the planet and
perhaps its atmosphere if applicable. perhaps its atmosphere if applicable.
\subsection{Multiple Gravity Assist Techniques} \subsection{Multiple Gravity Assist Techniques}
@@ -504,19 +456,21 @@
here for its robustness given any initial guess \cite{battin1984elegant}. here for its robustness given any initial guess \cite{battin1984elegant}.
Firstly, some geometric considerations must be accounted for. For any initial Firstly, some geometric considerations must be accounted for. For any initial
position, $\vec{r}_0$, and final position, $\vec{r}_f$, and time of flight $\Delta position, $\vec{r}_1$, and final position, $\vec{r}_2$, and time of flight $\Delta
t$, there are actually two separate transfer orbits that can connect the two points t$, there are actually two separate transfer orbits that can connect the two points
with paths that traverse less than one full orbit. For each of these, there are with paths that traverse less than one full orbit. Therefore, there are
actually then two trajectories that can connect the points actually then two trajectories that can connect the points
\cite{vallado2001fundamentals}. The first of the two will have a $\Delta \theta$ of \cite{vallado2001fundamentals}. The first of the two will have a $\Delta \theta$ of
less than 180 degrees, which we classify as a Type I trajectory, and the second will less than 180 degrees, which we classify as a Type I trajectory, and the second will
have a $\Delta \theta$ of greater than 180 degrees, which we call a Type II have a $\Delta \theta$ of greater than 180 degrees, which we call a Type II
trajectory. They will also differ in their direction of motion (clockwise or trajectory. They will also differ in their direction of motion (clockwise or
counter-clockwise about the focus). This can be seen in Figure~\ref{type1type2}. counter-clockwise about the focus). This can be seen in Figure~\ref{type1type2},
where both of the Lambert's solutions are presented for sample points in an orbit
around the Sun.
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.8\textwidth]{fig/lamberts} \includegraphics[width=0.8\textwidth]{LaTeX/fig/lamberts}
\caption{Visualization of the possible solutions to Lambert's Problem} \caption{Visualization of the possible solutions to Lambert's Problem}
\label{type1type2} \label{type1type2}
\end{figure} \end{figure}
@@ -524,7 +478,6 @@
The iteration used in this thesis will start by first calculating the change in true The iteration used in this thesis will start by first calculating the change in true
anomaly, $\Delta \theta$, as well as the cosine of this value, which can be found anomaly, $\Delta \theta$, as well as the cosine of this value, which can be found
by: by:
\begin{align} \begin{align}
\cos (\Delta \theta) &= \frac{\vec{r}_1 \cdot \vec{r}_2}{|\vec{r}_1| |\vec{r}_2|} \\ \cos (\Delta \theta) &= \frac{\vec{r}_1 \cdot \vec{r}_2}{|\vec{r}_1| |\vec{r}_2|} \\
\Delta \theta &= \arctan(y_2/x_2) - \arctan(y_1/x_1) \Delta \theta &= \arctan(y_2/x_2) - \arctan(y_1/x_1)
@@ -533,7 +486,6 @@
The direction of motion is then chosen such that counter-clockwise orbits are The direction of motion is then chosen such that counter-clockwise orbits are
considered, as travelling in the same direction as the planets is generally more considered, as travelling in the same direction as the planets is generally more
efficient. Next, the variable $A$ is defined: efficient. Next, the variable $A$ is defined:
\begin{equation} \begin{equation}
A = DM \sqrt{|r_1| |r_2| (1 - \cos(\Delta \theta))} A = DM \sqrt{|r_1| |r_2| (1 - \cos(\Delta \theta))}
\end{equation} \end{equation}
@@ -548,7 +500,6 @@
time of flight matches the expected value to within a provided tolerance. In order time of flight matches the expected value to within a provided tolerance. In order
to calculate the time of flight at each step, we must first calculate some useful to calculate the time of flight at each step, we must first calculate some useful
coefficients: coefficients:
\begin{equation}\label{loop_start} \begin{equation}\label{loop_start}
c_2 = \begin{cases} c_2 = \begin{cases}
\frac{1-\cos(\sqrt{\psi})}{\psi} \quad &\text{if} \, \psi > 10^{-6} \\ \frac{1-\cos(\sqrt{\psi})}{\psi} \quad &\text{if} \, \psi > 10^{-6} \\
@@ -556,31 +507,25 @@
1/2 \quad &\text{if} \, 10^{-6} > \psi > -10^{-6} 1/2 \quad &\text{if} \, 10^{-6} > \psi > -10^{-6}
\end{cases} \end{cases}
\end{equation} \end{equation}
\begin{equation} \begin{equation}
c_3 = \begin{cases} c_3 = \begin{cases}
\frac{\sqrt{\psi} - \sin sqrt{\psi}}{\psi^{3/2}} \quad &\text{if} \, \psi > 10^{-6} \\ \frac{\sqrt{\psi} - \sin \sqrt{\psi}}{\psi^{3/2}} \quad &\text{if} \, \psi > 10^{-6} \\
\frac{\sinh\sqrt{-\psi} - \sqrt{-\psi}}{(-\psi)^{3/2}} \quad &\text{if} \, \psi < -10^{-6} \\ \frac{\sinh\sqrt{-\psi} - \sqrt{-\psi}}{(-\psi)^{3/2}} \quad &\text{if} \, \psi < -10^{-6} \\
1/6 \quad &\text{if} \, 10^{-6} > \psi > -10^{-6} 1/6 \quad &\text{if} \, 10^{-6} > \psi > -10^{-6}
\end{cases} \end{cases}
\end{equation} \end{equation}
\noindent
Where the conditions of this piecewise function represent the elliptical, Where the conditions of this piecewise function represent the elliptical,
hyperbolic, and parabolic cases, respectively. Once we have these, we can calculate hyperbolic, and parabolic cases, respectively. Once we have these, we can calculate
another variable, $y$: another variable, $y$:
\begin{equation} \begin{equation}
y = |r_1| + |r_2| + \frac{A (c_3 \psi - 1)}{\sqrt{c_2}} y = |r_1| + |r_2| + \frac{A (c_3 \psi - 1)}{\sqrt{c_2}}
\end{equation} \end{equation}
We can then finally calculate the variable $\chi$, and from that, the time of We can then finally calculate the variable $\chi$, and from that, the time of
flight: flight:
\begin{equation} \begin{equation}
\chi = sqrt{\frac{y}{c_2}} \chi = \sqrt{\frac{y}{c_2}}
\end{equation} \end{equation}
\begin{equation} \begin{equation}
\Delta t = \frac{c_3 \chi^3 + A \sqrt{y}}{\sqrt{c_2}} \Delta t = \frac{c_3 \chi^3 + A \sqrt{y}}{\sqrt{c_2}}
\end{equation} \end{equation}
@@ -593,22 +538,17 @@
The resulting $f$ and $g$ functions (and the derivative of $g$, $\dot{g}$) can then The resulting $f$ and $g$ functions (and the derivative of $g$, $\dot{g}$) can then
be calculated: be calculated:
\begin{align} \begin{align}
f &= 1 - \frac{y}{|r_1|} \\ f &= 1 - \frac{y}{|r_1|} \\
g &= A \sqrt{\frac{y}{\mu}} \\ g &= A \sqrt{\frac{y}{\mu}} \\
\dot{g} &= 1 - \frac{y}{|r_2|} \dot{g} &= 1 - \frac{y}{|r_2|}
\end{align} \end{align}
And from these, we can calculate the velocities of the transfer points as: And from these, we can calculate the velocities of the transfer points as:
\begin{align} \begin{align}
\vec{v}_1 &= \frac{\vec{r}_1 - f \vec{r}_2}{g} \\ \vec{v}_1 &= \frac{\vec{r}_1 - f \vec{r}_2}{g} \\
\vec{v}_2 &= \frac{\dot{g} \vec{r}_2 - \vec{r}_1}{g} \vec{v}_2 &= \frac{\dot{g} \vec{r}_2 - \vec{r}_1}{g}
\end{align} \end{align}
Fully describing the connecting path with the specified flight time.
\noindent
Fully constraining the connecting orbit.
\subsubsection{Planetary Ephemeris} \subsubsection{Planetary Ephemeris}
@@ -621,8 +561,8 @@
The primary use of SPICE in this thesis, however, was to determine the planetary The primary use of SPICE in this thesis, however, was to determine the planetary
ephemeris at a known epoch. Using the NAIF0012 and DE430 kernels, ephemeris in the ephemeris at a known epoch. Using the NAIF0012 and DE430 kernels, ephemeris in the
ecliptic plane J2000 frame (ICRF) could be easily determined for a given epoch, provided as International Celestial Reference Frame could be easily determined for a given
a decimal Julian Day since the J2000 epoch. epoch, provided as a decimal Julian Day since the J2000 epoch.
\subsubsection{Porkchop Plots} \subsubsection{Porkchop Plots}
@@ -642,24 +582,17 @@
Using porkchop plots such as the one in Figure~\ref{porkchop}, mission designers can Using porkchop plots such as the one in Figure~\ref{porkchop}, mission designers can
quickly visualize which natural trajectories are possible between planets. Using the quickly visualize which natural trajectories are possible between planets. Using the
fact that incoming and outgoing $v_\infty$ magnitudes must be the same for a flyby, fact that incoming and outgoing $v_\infty$ magnitudes must be the same for a flyby,
a savvy mission designer can even begin to work out what combinations of flybys a mission designer can even begin to work out what combinations of flybys might be
might be possible for a given timeline, spacecraft state, and planet selection. possible for a given timeline, spacecraft state, and planet selection.
%TODO: Create my own porkchop plot
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{fig/porkchop} \includegraphics[width=\textwidth]{LaTeX/fig/porkchop}
\caption{A sample porkchop plot of an Earth-Mars transfer} \caption{A sample porkchop plot of an Earth-Mars transfer}
\label{porkchop} \label{porkchop}
\end{figure} \end{figure}
However, this is an impulsive thrust-centered approach. The solution to Lambert's \section{Modeling Low Thrust Control} \label{low_thrust}
problem assumes a natural trajectory. 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.
\section{Low Thrust Considerations} \label{low_thrust}
In this section, we'll discuss the intricacies of continuous low-thrust trajectories in 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 particular. There are many methods for optimizing such profiles and we'll briefly discuss
@@ -667,7 +600,7 @@
as introduce the concept of a control law and the notation used in this thesis for modelling as introduce the concept of a control law and the notation used in this thesis for modelling
low-thrust trajectories more simply. low-thrust trajectories more simply.
\subsection{Specific Impulse} \subsection{Engine Model}
The primary advantage of continuous thrust methods over their impulsive counterparts is The primary advantage of continuous thrust methods over their impulsive counterparts is
in their fuel-efficiency in generating changes in velocity. Put specifically, all in their fuel-efficiency in generating changes in velocity. Put specifically, all
@@ -679,45 +612,34 @@
This efficiency is often captured in a single variable called specific impulse, often This efficiency is often captured in a single variable called specific impulse, often
denoted as $I_{sp}$. We can derive the specific impulse by starting with the rocket denoted as $I_{sp}$. We can derive the specific impulse by starting with the rocket
thrust equation\cite{sutton2016rocket}: thrust equation\cite{sutton2016rocket}:
\begin{equation} \begin{equation}
F = \dot{m} v_e + \Delta p A_e F = \dot{m} v_e + \Delta p A_e
\end{equation} \end{equation}
\noindent
Where $F$ is the thrust imparted, $\dot{m}$ is the fuel mass rate, $v_e$ is the exhaust Where $F$ is the thrust imparted, $\dot{m}$ is the fuel mass rate, $v_e$ is the exhaust
velocity of the fuel, $\Delta p$ is the change in pressure across the exhaust opening, velocity of the fuel, $\Delta p$ is the change in pressure across the exhaust opening,
and $A_e$ is the area of the exhaust opening. We can then define a new variable and $A_e$ is the area of the exhaust opening. We can then define a new variable
$v_{eq}$, such that the thrust equation becomes: $v_{eq}$, such that the thrust equation becomes:
\begin{align} \begin{align}
v_{eq} &= v_e - \frac{\Delta p A_e}{\dot{m}} \\ v_{eq} &= v_e + \frac{\Delta p A_e}{\dot{m}} \\
F &= \dot{m} v_{eq} \label{isp_1} F &= \dot{m} v_{eq} \label{isp_1}
\end{align} \end{align}
\noindent
And we can then take the integral of this value with respect to time to find the total And we can then take the integral of this value with respect to time to find the total
impulse, dividing by the weight of the fuel to derive the specific impulse: impulse, dividing by the weight of the fuel to derive the specific impulse:
\begin{align} \begin{align}
I &= \int F dt = \int \dot{m} v_{eq} dt = m_e v_{eq} \\ I &= \int F dt = \int \dot{m} v_{eq} dt = m_e v_{eq} \\
I_{sp} &= \frac{I}{m_e g_0} = \frac{m_e v_{eq}}{m_e g_0} = \frac{v_{eq}}{g_0} I_{sp} &= \frac{I}{m_e g_0} = \frac{m_e v_{eq}}{m_e g_0} = \frac{v_{eq}}{g_0}
\end{align} \end{align}
Plugging Equation~\ref{isp_1} into the previous equation we can derive the following Plugging Equation~\ref{isp_1} into the previous equation we can derive the following
formula for $I_{sp}$: formula for $I_{sp}$:
\begin{equation} \label{isp_real} \begin{equation} \label{isp_real}
I_{sp} = \frac{F}{\dot{m} g_0} I_{sp} = \frac{F}{\dot{m} g_0}
\end{equation} \end{equation}
\noindent
Which is generally taken to be a value with units of seconds and effectively represents Which is generally taken to be a value with units of seconds and effectively represents
the efficiency with which a thruster converts mass to thrust. the efficiency with which a thruster converts mass to thrust.
\subsection{Sims-Flanagan Transcription} \subsection{Sims-Flanagan Transcription}
this thesis chose to use a model well suited for modeling low-thrust paths: the In this thesis the following approach is used for modeling low-thrust paths: the
Sims-Flanagan transcription (SFT)\cite{sims1999preliminary}. The SFT allows for Sims-Flanagan transcription (SFT)\cite{sims1999preliminary}. The SFT allows for
flexibility in the trade-off between fidelity and performance, which makes it very flexibility in the trade-off between fidelity and performance, which makes it very
useful for this sort of preliminary analysis. useful for this sort of preliminary analysis.
@@ -730,7 +652,7 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=0.6\textwidth]{fig/sft} \includegraphics[width=0.6\textwidth]{LaTeX/fig/sft}
\caption{Example of an orbit raising using the Sims-Flanagan Transcription with 7 \caption{Example of an orbit raising using the Sims-Flanagan Transcription with 7
Sub-Trajectories} Sub-Trajectories}
\label{sft_fig} \label{sft_fig}
@@ -753,7 +675,7 @@
continuous low-thrust trajectory within the Two-Body Problem, with only continuous low-thrust trajectory within the Two-Body Problem, with only
linearly-increasing computation time\cite{sims1999preliminary}. linearly-increasing computation time\cite{sims1999preliminary}.
\subsection{Low-Thrust Control Laws} \subsection{Low-Thrust Control Vector Description}
In determining a low-thrust arc, a number of variables must be accounted for and, In determining a low-thrust arc, a number of variables must be accounted for and,
ideally, optimized. Generally speaking, this means that a control law must be determined ideally, optimized. Generally speaking, this means that a control law must be determined
@@ -766,24 +688,23 @@
The methods for determining this direction varies greatly depending on the particular The methods for determining this direction varies greatly depending on the particular
control law chosen for that mission. Often, this process involves first determining a control law chosen for that mission. Often, this process involves first determining a
useful frame to think about the kinematics of the spacecraft. In this case, we'll use a useful frame to think about the kinematics of the spacecraft. In this case, we'll use a
frame often used in these low-thrust control laws: the spacecraft $\hat{R} \hat{\theta} frame often used in these low-thrust control laws: the spacecraft-centered $\hat{R}
\hat{H}$ frame. In this frame, the $\hat{R}$ direction is the radial direction from the \hat{\theta} \hat{H}$ frame. In this frame, the $\hat{R}$ direction is the radial
center of the primary to the center of the spacecraft. The $\hat{H}$ hat is direction from the center of the primary to the center of the spacecraft. The $\hat{H}$
perpendicular to this, in the direction of orbital momentum (out-of-plane) and the hat is perpendicular to this, in the direction of orbital momentum (out-of-plane) and
$\hat{\theta}$ direction completes the right-handed orthonormal frame. the $\hat{\theta}$ direction completes the right-handed orthonormal triad.
This frame is useful because, for a given orbit, especially a nearly circular one, the This frame is useful because, for a given orbit, especially a nearly circular one, the
$\hat{\theta}$ direction is nearly aligned with the velocity direction for that orbit at $\hat{\theta}$ direction is nearly aligned with the velocity direction for that orbit at
that moment. This allows us to define a set of two angles, which we'll call $\alpha$ and that moment. This allows us to define a set of two angles, which we'll call $\alpha$ and
$\beta$, to represent the in and out of plane pointing direction of the thrusters. This $\beta$, to represent the in and out of plane pointing direction of the thrusters. This
convention is useful because a $(0,0)$ set represents a thrust force more or less convention is useful because, in a near-circular path, a $(0,0)$ set represents a thrust
directly in line with the direction of the velocity, a commonly useful thrusting force more or less directly in line with the direction of the velocity, a commonly
direction for most effectively increasing (or decreasing if negative) the angular useful thrusting direction for most effectively increasing (or decreasing if negative)
momentum and orbital energy of the trajectory. the angular momentum and orbital energy of the trajectory.
Using these conventions, we can then redefine our thrust vector in terms of the $\alpha$ Using these conventions, we can then redefine our thrust vector in terms of the $\alpha$
and $\beta$ angles in the chosen frame: and $\beta$ angles in the chosen frame:
\begin{align} \begin{align}
F_r &= F \cos(\beta) \sin (\alpha) \\ F_r &= F \cos(\beta) \sin (\alpha) \\
F_\theta &= F \cos(\beta) \cos (\alpha) \\ F_\theta &= F \cos(\beta) \cos (\alpha) \\
@@ -792,12 +713,12 @@
\subsubsection{Thrust Magnitude} \subsubsection{Thrust Magnitude}
However, there is actually another variable that can be varied by the majority of There is another variable that can be varied by the majority of electric thrusters.
electric thrusters. Either by controlling the input power of the thruster or the duty Either by controlling the input power of the thruster or the duty cycle, the thrust
cycle, the thrust magnitude can also be varied, limited by the maximum thrust available magnitude can also be varied, limited by the maximum thrust available to the thruster.
to the thruster. Not all control laws allow for this fine-tuned control of the thruster. Not all control laws allow for this fine-tuned control of the thruster.
The algorithm used in this thesis does vary the magnitude of the thrust control. In The approach used in this thesis does vary the magnitude of the thrust control. In
certain cases it actually can be useful to have some fine-tuned control over the certain cases it actually can be useful to have some fine-tuned control over the
magnitude of the thrust. Since the optimization in this algorithm is automatic, it is magnitude of the thrust. Since the optimization in this algorithm is automatic, it is
relatively straightforward to consider the control thrust as a 3-dimensional vector in relatively straightforward to consider the control thrust as a 3-dimensional vector in
@@ -818,21 +739,14 @@
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{fig/low_efficiency} \includegraphics[width=\textwidth]{LaTeX/fig/low_efficiency}
\caption{Graphic of an orbit-raising with a low efficiency cutoff} \caption{Graphic of an orbit-raising with a low efficiency cutoff}
\label{low_efficiency_fig} \label{low_efficiency_fig}
\end{figure} \end{figure}
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{fig/high_efficiency} \includegraphics[width=\textwidth]{LaTeX/fig/high_efficiency}
\caption{Graphic of an orbit-raising with a high efficiency cutoff} \caption{Graphic of an orbit-raising with a high efficiency cutoff}
\label{high_efficiency_fig} \label{high_efficiency_fig}
\end{figure} \end{figure}
All of 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.

56
LaTeX/trajectory_optimization.tex
View File

@@ -6,10 +6,11 @@
highly non-linear, unpredictable systems such as this. The field that developed to highly non-linear, unpredictable systems such as this. The field that developed to
approach this problem is known as Non-Linear Programming (NLP) Optimization. approach this problem is known as Non-Linear Programming (NLP) Optimization.
A Non-Linear Programming Problem is defined by an attempt to optimize a function A Non-Linear Programming Problem involves finding a solution that optimizes a function
$f(\vec{x})$, subject to constraints $\vec{g}(\vec{x}) \le 0$ and $\vec{h}(\vec{x}) = 0$ $f(\vec{x})$, subject to constraints $\vec{g}(\vec{x}) \le 0$ and $\vec{h}(\vec{x}) = 0$
where $n$ is a positive integer, $x$ is any subset of $R^n$, $g$ and $h$ can be vector where $n$ is a positive integer, $x$ is any subset of $R^n$, $g$ and $h$ can be vector
valued functions of any size, and at least one of $f$, $g$, and $h$ must be non-linear. valued functions of any size, and at least one of $f$, $\vec{g}$, and $\vec{h}$ must be
non-linear.
There are, however, two categories of approaches to solving an NLP. The first category, 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 indirect methods, involve declaring a set of necessary and/or sufficient conditions for
@@ -20,10 +21,10 @@
The other category is the direct methods. In a direct optimization problem, the cost The other category is the direct methods. In a direct optimization problem, the cost
function itself provides a value that an iterative numerical optimizer can measure function itself provides a value that an iterative numerical optimizer can measure
itself against. The optimal solution is then found by varying the inputs $x$ until the itself against. The optimal solution is then found by varying the inputs $\vec{x}$ until
cost function is reduced to a minimum value, often determined by its derivative the cost function is reduced to a minimum value, often determined by its derivative
jacobian. A number of tools have been developed to optimize NLPs via this direct method jacobian. A number of tools have been developed to formulate NLPs for optimization via
in the general case. this direct method in the general case.
Both of these methods have been applied to the problem of low-thrust interplanetary Both of these methods have been applied to the problem of low-thrust interplanetary
trajectory optimization \cite{Casalino2007IndirectOM} to find local optima over trajectory optimization \cite{Casalino2007IndirectOM} to find local optima over
@@ -40,7 +41,7 @@
Therefore, a direct optimization method was leveraged by transcribing the problem into Therefore, a direct optimization method was leveraged by transcribing the problem into
an NLP and using IPOPT to find the local minima. an NLP and using IPOPT to find the local minima.
\subsubsection{Non-Linear Solvers} \subsection{Non-Linear Solvers}
One of the most common packages for the optimization of NLP problems is One of the most common packages for the optimization of NLP problems is
SNOPT\cite{gill2005snopt}, which is a proprietary package written primarily using a SNOPT\cite{gill2005snopt}, which is a proprietary package written primarily using a
@@ -48,7 +49,7 @@
University. It uses a sparse sequential quadratic programming algorithm as its University. It uses a sparse sequential quadratic programming algorithm as its
back-end optimization scheme. back-end optimization scheme.
Another common NLP optimization packages (and the one used in this implementation) Another common NLP optimization package (and the one used in this implementation)
is the Interior Point Optimizer or IPOPT\cite{wachter2006implementation}. It uses is the Interior Point Optimizer or IPOPT\cite{wachter2006implementation}. It uses
an Interior Point Linesearch Filter Method and was developed as an open-source an Interior Point Linesearch Filter Method and was developed as an open-source
project by the organization COIN-OR under the Eclipse Public License. project by the organization COIN-OR under the Eclipse Public License.
@@ -63,7 +64,7 @@
libraries that port these are quite modular in the sense that multiple algorithms can be libraries that port these are quite modular in the sense that multiple algorithms can be
tested without changing much source code. tested without changing much source code.
\subsubsection{Interior Point Linesearch Method} \subsection{Interior Point Linesearch Method}
As mentioned above, this project utilized IPOPT which leveraged an Interior Point 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 Linesearch method. A linesearch algorithm is one which attempts to find the optimum
@@ -74,12 +75,7 @@
step the initial guess, now labeled $x_{k+1}$ after the addition of the ``step'' 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. 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 \subsection{Shooting Schemes for Solving a Two-Point Boundary Value Problem}
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.
\subsubsection{Shooting Schemes for Solving a Two-Point Boundary Value Problem}
One straightforward approach to trajectory corrections is a single shooting One straightforward approach to trajectory corrections is a single shooting
algorithm, which propagates a state, given some control variables forward in time to algorithm, which propagates a state, given some control variables forward in time to
@@ -87,39 +83,30 @@
iterative process, using the correction scheme, until the target state and the iterative process, using the correction scheme, until the target state and the
propagated state matches. propagated state matches.
As an example, we can consider the Two-Point Boundary Value Problem (TPBVP) defined As an example, we can consider the one-dimensional Two-Point Boundary Value Problem
by: (TPBVP) defined by:
\begin{equation} \begin{equation}
y''(t) = f(t, y(t), y'(t)), y(t_0) = y_0, y(t_f) = y_f y''(t) = f(t, y(t), y'(t)), y(t_0) = y_0, y(t_f) = y_f
\end{equation} \end{equation}
\noindent
We can then redefine the problem as an initial-value problem: We can then redefine the problem as an initial-value problem:
\begin{equation} \begin{equation}
y''(t) = f(t, y(t), y'(t)), y(t_0) = y_0, y'(t_0) = x y''(t) = f(t, y(t), y'(t)), y(t_0) = y_0, y'(t_0) = \dot{y}_0
\end{equation} \end{equation}
\noindent
With $y(t,x)$ as a solution to that problem. Furthermore, if $y(t_f, x) = y_f$, then With $y(t,x)$ as a solution to that problem. Furthermore, if $y(t_f, x) = y_f$, then
the solution to the initial-value problem is also the solution to the TPBVP as well. the solution to the initial-value problem is also the solution to the TPBVP as well.
Therefore, we can use a root-finding algorithm, such as the bisection method, Therefore, we can use a root-finding algorithm, such as the bisection method,
Newton's Method, or even Laguerre's method, to find the roots of: Newton's Method, or even Laguerre's method, to find the roots of:
\begin{equation} \begin{equation}
F(x) = y(t_f, x) - y_f F(x) = y(t_f, x) - y_f
\end{equation} \end{equation}
\noindent
To find the solution to the IVP at $x_0$, $y(t_f, x_0)$ which also provides a To find the solution to the IVP at $x_0$, $y(t_f, x_0)$ which also provides a
solution to the TPBVP. This technique for solving a Two-Point Boundary Value solution to the TPBVP. This technique for solving a Two-Point Boundary Value
Problem can be visualized in Figure~\ref{single_shoot_fig}. Problem can be visualized in Figure~\ref{single_shoot_fig}.
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{fig/single_shoot} \includegraphics[width=\textwidth]{LaTeX/fig/single_shoot}
\caption{Visualization of a single shooting technique over a trajectory arc} \caption{Single shooting over a trajectory arc}
\label{single_shoot_fig} \label{single_shoot_fig}
\end{figure} \end{figure}
@@ -138,24 +125,23 @@
each of these points we can then define a separate control, which may include the each of these points we can then define a separate control, which may include the
states themselves. The end state of each arc and the beginning state of the next states themselves. The end state of each arc and the beginning state of the next
must then be equal for a valid arc (with the exception of velocity discontinuities must then be equal for a valid arc (with the exception of velocity discontinuities
if allowed for maneuvers at that point), as well as the final state matching the if allowed for maneuvers or gravity assists at that point), as well as the final
target final state. state matching the target final state.
\begin{figure}[H] \begin{figure}[H]
\centering \centering
\includegraphics[width=\textwidth]{fig/multiple_shoot} \includegraphics[width=\textwidth]{LaTeX/fig/multiple_shoot}
\caption{Visualization of a multiple shooting technique over a trajectory arc} \caption{Visualization of a multiple shooting technique over a trajectory arc}
\label{multiple_shoot_fig} \label{multiple_shoot_fig}
\end{figure} \end{figure}
In this example, it can be seen that there are now more constraints (places where In this example, it can be seen that there are now more constraints (places where
the states need to match up, creating an $x_{error}$ term) as well as control the states need to match up, creating an $\vec{x}_{error}$ term) as well as control
variables (the $\Delta V$ terms in the figure). This technique actually lends itself variables (the $\Delta V$ terms in the figure). This technique actually lends itself
very well to low-thrust arcs and, in fact, Sims-Flanagan Transcribed low-thrust arcs very well to low-thrust arcs and, in fact, Sims-Flanagan Transcribed low-thrust arcs
in particular, because there actually are control thrusts to be optimized at a in particular, because there actually are control thrusts to be optimized at a
variety of different points along the orbit. This is, however, not an exhaustive variety of different points along the orbit. This is, however, not an exhaustive
description of ways that multiple shooting can be used to optimize a trajectory, description of ways that multiple shooting can be used to optimize a trajectory.
simply the most convenient for low-thrust arcs.
\section{Monotonic Basin Hopping Algorithms} \section{Monotonic Basin Hopping Algorithms}

85
Makefile
View File

@@ -1,68 +1,39 @@
OPTIONS = markdown+yaml_metadata_block+smart SRC_DIR = LaTeX/
THESIS_SRC_NAMES = thesis.tex thesis.bib approach.tex conclusion.tex introduction.tex \
NOTES = $(wildcard prelim_notes/*.md) results.tex trajectory_design.tex trajectory_optimization.tex
NOTES_PDFS = $(patsubst %.md,%.pdf,$(NOTES)) THESIS_SRC = $(addprefix $(SRC_DIR)/,$(THESIS_SRC_NAMES))
THESIS_PRES_SRC_NAMES = presentation.tex
THESIS = LaTeX/thesis.tex THESIS_PRES_SRC = $(addprefix $(SRC_DIR)/,$(THESIS_PRES_SRC_NAMES))
THESIS_PRES = presentation.pdf
THESIS_PDF = thesis.pdf THESIS_PDF = thesis.pdf
BUILD_DIR = temp/
all: $(THESIS_PDF) $(NOTES_PDFS) all: $(THESIS_PDF) $(THESIS_PRES)
$(NOTES_PDFS): $(NOTES) $(BUILD_DIR):
pandoc \ mkdir -p $(BUILD_DIR)
--variable mainfont="Roboto" \
--variable monofont="Fira Code" \
--variable fontsize=11pt \
--variable geometry:"top=1in, bottom=1in, left=1in, right=1in" \
--variable geometry:letterpaper \
-f markdown $< \
-o $@
thesis_pdf/: $(BUILD_DIR)/$(THESIS_PDF): $(BUILD_DIR) $(THESIS_SRC)
mkdir $@ xelatex --output-directory $(BUILD_DIR) $(SRC_DIR)/thesis
bibtex $(BUILD_DIR)/thesis
xelatex --output-directory $(BUILD_DIR) $(SRC_DIR)/thesis
xelatex --output-directory $(BUILD_DIR) $(SRC_DIR)/thesis
$(THESIS_PDF): $(THESIS) LaTeX/thesis.bib $(BUILD_DIR)/$(THESIS_PRES): $(BUILD_DIR) $(THESIS_PRES_SRC)
mkdir -p temp xelatex --output-directory $(BUILD_DIR) $(SRC_DIR)/presentation
cp -r LaTeX/fig .
cp -r LaTeX/flowcharts . $(THESIS_PDF): $(BUILD_DIR)/$(THESIS_PDF)
cp -r LaTeX/thesis.tex . cp $(BUILD_DIR)/thesis.pdf .
cp -r LaTeX/macros.tex .
cp -r LaTeX/thesis.bib . $(THESIS_PRES): $(BUILD_DIR)/$(THESIS_PRES)
cp -r LaTeX/thesis.cls . cp $(BUILD_DIR)/presentation.pdf .
xelatex thesis
bibtex thesis
xelatex thesis
xelatex thesis
rm -rf fig
rm -rf flowcharts
cp thesis.pdf temp/.
rm -rf macros.tex
rm -rf thesis.*
cp temp/thesis.pdf .
rm -rf temp
.PHONY: clean revise .PHONY: clean revise
clean: clean:
rm -rf $(THESIS_PDF) rm -rf $(THESIS_PDF)
rm -rf $(NOTES_PDFS) rm -rf $(THESIS_PRES)
rm -rf thesis_pdf rm -rf $(BUILD_DIR)
revise: $(THESIS) LaTeX/thesis.bib final: $(THESIS_PDF) $(THESIS_PRES)
mkdir -p temp rm -rf $(BUILD_DIR)
cp -r LaTeX/fig .
cp -r LaTeX/flowcharts .
cp -r LaTeX/thesis.tex .
cp -r LaTeX/macros.tex .
cp -r LaTeX/thesis.bib .
cp -r LaTeX/thesis.cls .
xelatex thesis
bibtex thesis
xelatex thesis
xelatex thesis
rm -rf fig
rm -rf flowcharts
cp thesis.pdf temp/.
rm -rf macros.tex
rm -rf thesis.*
cp temp/thesis.pdf .

3
best_missions
View File

@@ -1,3 +0,0 @@
./archive/EMMJS_2021-12-03T17:38:37.546/mission
./archive/EMS_2021-12-12T16:01:04.514/mission
./archive/EVMS_2021-12-14T16:31:21.927/mission

15
email_to_commitee.md
View File

@@ -1,15 +0,0 @@
Dr. Scheeres,
I'm currently working on finishing up my Master's Thesis. The topic will be a routine for optimizing
Interplanetary Low-Thrust Trajectories, combining approaches for flyby optimization and low-thrust
trajectory arc optimization. I think that your research interests would make you an ideal fit for a
place on my committee, so I would like to ask if you have any availability to be on my committee
this semester? Currently, I am targeting finish up sometime around March 12th, so it would probably
be a defense in late March or early April, but I have a little room for flexibility if that doesn't
work.
If you think you may be able to sit on my committee then just let me know and I can begin the
process of determining everyone's availabilities and finding the most appropriate defense date.
Thanks!
Connor Johnstone

0
ipopt.out
View File

102
prelim_notes/paper_notes.md
View File

@@ -1,102 +0,0 @@
# Notes on Research Papers
## Neural Networks in Time-Optimal Low-Thrust Interplanetary Transfers
- For the most part this paper isn't *that* relevant.
- It seems to be a low-fidelity method, or at least, it's fidelity is kind of hard to pin down,
since it uses a neural net.
- However, the neural net is an interesting concept.
- And in fact, this paper advocates for a concept of using the neural net as a predictor. So, say,
given a leg of a particular journey, as in the Englander paper, could I use this technique rather
than Sims-Flanagan transcription in the outer loop?
- Something to consider as an alternative to the method proposed in Englander
- Could also be used to generate initial guesses for the single-shooting methods that I'll have to
use if I use an indirect optimization method.
## Design and optimization of interplanetary low-thrust trajectory with planetary aerogravity-assist maneuver
- I didn't realize that this paper is specifically talking about *aero*gravity assists. I don't
think that level of complication is necessary for this paper. I'm going to stick with gravity
assists that don't get into atmospheric effects.
## Orbital and Angular Motion Construction for Low Thrust Interplanetary Flight
- This one actually isn't even about optimization. Not relevant
- To be honest, I don't actually understand it very well anyway
## A Rapid Estimation for Interplanetary Low-Thrust Trajectories Using Support Vector Regression
- This is another machine-learning approach
- It uses a different approach that I'm not that familiar with (Support Vector Regression)
- It looks like this is a form of regression similar to linear regression, I suppose being used
by the machine learning algorithm for predicting optimal trajectories
- However, it seems like everything that applies in the neural net paper probably apply here as well
- This could be an alternative for predicting optimal trajectories over certain legs of the journey
## Automated Solution of the Low-Thrust Interplanetary Trajectory Problem
- This is the Englander paper I mentioned earlier. It seems highly relevant as they're essentially
doing what I'd like to do: producing an automated method for high-level interplanetary low-thrust
mission design including the flybys.
- I need to look up MALTO and GALLOP (established tools that do this)
- This paper also open sources it's code: available [here](https://opensource.gsfc.nasa.gov/projects/emtg/index.php)
- This paper formulates the problem as a *Hybrid Optimal Control Problem*. This requires some
further research by me, but from the paper it seems to be a way of optimizing two seperables
subproblems where one of the subproblems exists as a sub-loop of the other problem (for instance
optimizing flybys in a high-level loop and particular planet-to-planet trajectories as a sub-loop of
that problem). This apparently works because the "outer-loop" uses discrete variables while the
"inner-loop" uses continuous variables.
- The outer loop is based on the "null-gene" transciption from another Englander paper and uses a
genetic algorithm.
- I'm not going to go too deep here into the details of the GA. But there's an entire paper on
it
- The paper does mention that it lends itself well to parallelization, which is true. Kubernetes
cluster?
- The inner loop uses Sim-Flanagan transcription combined with Monotonic Basin Hopping, a method for
single-shooting without initial guesses
- Sims-Flanagan transcription is where you discretize the flight arc into many smaller time
steps and the thrust applied is approximated as an impulsive thrust in the middle of each time
step.
- SFT is considered to be a "medium-fidelity" approach
- For this solver, the trajectory betweens these points is produced as a solution to Kepler's
problem, which is basically just the analytical solution to the 2BP, so that no derivatives or
orbit-propagation is needed, for speed.
- One thing I noted about this approach is that it doesn't seem to include a possibility for
"coasting arcs" or throttling anything less than 100% (though the modeling of what 100% means is
quite thorough), so perhaps we're missing some fidelity there?
- I think SNOPT is used to optimize these "inner-loops". This should be pretty fast since it
just uses Kepler's eq
- The MBH method eliminates the need to solve Lambert's problem for initial guesses. This allows
for a more robust analysis of the search base (if there are global optima further from the local
optima near lambert's solution) but might be slower? I'm not sure.
- The technique is kind of weird, but I suppose it works.
- This paper uses a hierarchy of events starting with the overall *mission*, which separates into
*journies*, which, in the example I'm pursuing will be Earth -> Neptune and Neptune -> Earth (if
applicable, but probably just the first one). Then these journies are further divided into
*phases*, which include each planet -> planet leg. The number of phases and the identities of the
planets are chosen by the algorithm.
- The paper goes into some length to determine what the launch C3, propulsion, power, and ephemeris
modeling are. This is all very useful, but as far as I can tell it's pretty typical, so I won't note
too much about it.
- However, it does mention that SPICE presented some challenges for using a preferable method of
parallelization. As an alternative, the paper mentions that FIRE could be used instead for
ephemeris. Which might be worth looking into.
- The paper also include pseudocode, which is nice
## Automated Mission Planning via Evolutionary Algorithms
- This is another Englander paper that gives more details on the outer-loop GA. Useful for details.
## Multi-Objective Low-Thrust Interplanetary Trajectory Optimization Based on Generalized Logarithmic Spirals
- This is the first paper I looked at. It's actually quite similar to the Englander paper, but I
think not quite as good
- Again, it formulates the problem as an HOCP.
- However, the "inner-loop" for this problem is an optimization of generalized logarithmic spirals.
I don't think this is a very high fidelity method.
- The outer step uses collocation and an NLP optimizer (looks like it actually might just feed the
guesses into GALLOP, which I assume uses SNOPT, though, to be honest, I can't find much on it from
a quick search)
- I'm leaning toward Englander's approach over this one, perhaps with an alternative being to use
one of the machine-learning approaches from above for the inner loop instead

37
prelim_notes/plan.md
View File

@@ -1,37 +0,0 @@
# Plan After Literature Review
After reading the papers in the Google Drive (see `paper_notes.md`) I've come up with the following
plan:
I think that I'd like to follow an approach similar to what I saw in Englander and Morante. However,
I think it makes more sense to follow the Englander approach for the specific optimizers being used
in the inner and outer loops. Specifically this means:
- Set up the problem as a Hybrid Optimal Control Problem (HOCP) with an inner and an outer loop
- The outer loop will determine the number and identities of the flybys and optimize using a Binary
Genetic Algorithm described by Englander
- The inner loop in Englander uses Sims-Flanagan Transcription optimized using monotonic basin
hopping. This seems like a good approach. I'd like to use this approach, but also consider using,
as an alternative method for comparison, one of the machine-learning algorithms from the other
papers.
- There are a number of other details including modeling launch C3, power, thrust, and ephemeris.
For all of these I'll use either the exact approach from Englander or a similar approach. There
exists an option to use alternatives to SPICE for ephemeris, but I think the parallelization
problems that SPICE poses can be solved in other ways.
- Specifically, I'd like to build this program using a micro-service architecture. This could
allow for deployment using Kubernetes clusters. This will handle the parallelization (by
running multiple inner loop micro-services at once) and allow for simpler use in production
environments if that's ever needed, as kubernetes has a robust integration with most
web-hosting services. This also allows for flexible scalability if improved speed is ever
needed.
## Open Questions
- Are there other inner and outer loop optimization approaches I should consider?
- I didn't see a ton of information regarding the specific parameters of the orbital flybys outside
of the number and the identities of the planets. I need to figure out whether Englander approaches
that optimization in the inner or the outer loop
- I would like to use a fast language for this, and ideally one that I could benefit from learning
about (by this I mean, preferably not C++, as I've used that before -- though it does integrate
well with tools like SPICE, so it is a good option). Preferably Rust, as it has the fastest server
libraries, which I don't expect to have a huge overhead, but is something to consider.

76
prelim_notes/technical_plan.md
View File

@@ -1,76 +0,0 @@
# Technical Plan
## Notes
- Each of these steps will be completed first in Julia, as I can use that to get
to a working point really quickly. Then, I'll redo the efforts in Rust inside
of docker containers, providing increased robustness and speed. I'll provide
expected time for completion after each item in the form (x/y) where x is the
number of hours I expect to take to write the Julia code and y is the number
of hours I expect to take to write the Rust code.
## Stages of Development
1. First I need to set up the inner loop. The inner loop is a Sims-Flanagan transcription
single-shooting algorithm, optimized using monotonic basin-hopping.
a. I've constructed single-shooting algorithms before, and the Sims-Flanagan
transcription is just a really simple notation for approximating low-thrust
arcs. So first I'll need to set up a very simple Sims-Flanagan single
shooter. I don't expect this to take long. I will write a simple
single-shooting algorithm that takes in a start and end state (for now... in
stage 2b below I'll update this to only need the J200 time (for ephemeris
lookup) and the velocity) and solves Kepler's equation to satisfy the
continuity condition. For the rust code, this will require some considerable
initial setup of the docker containers and such, though. (6/20)
b. Next I'll need to set up the monotonic basin-hopping algorithm. I've
written similar optimization algorithms in my spacecraft trajectory
optimization course, but none that use basin-hopping. However, it's more
or less just a genetic algorithm, so it should be simple to set up in a
very general way. The function I'm optimizing can be seen as a system of
inputs (initial and final states) and outputs (mass used, or some other
more complicated measure of optimality), so I can follow many very
descriptive accounts online of this type of genetic algorithm. (12/18)
2. Then I'll need to set up the outer loop. This will require optimizing over
inputs (selection of flybys) by calling the inner loop to optimize each leg
and summing these by way of some cost function. These flyby selections will
be optimized using the genetic algorithm described in the first Englander
paper.
a. First I'll set up a very simple version of the outer loop, that
doesn't optimize anything. It will take in a selection of flybys and
simply call the inner loop for each of the flybys, outputting some cost
function. This shouldn't take long, but at this point I'm adding time
because the inner loop may take some time to call (8/12).
b. There are also some other inputs to this algorithm, that won't change from
problem to problem, but will take some effort to calculate. These include
initial launch C3, ephemeris, problem constraints, etc. SPICE and basic
modeling can be used for most of this, but I should include some time for
incorporating SPICE via C-bindings, which should be relatively easy in
both Julia and Rust (in fact, a wrapper exists at least in Julia
already). As a less high fidelity backup, I've also modeled ephemeris by
polynomials before, so I can do that again, if I find incorporating SPICE
is taking too long (16/16)
c. Next, I'll write the genetic algorithm in a very general way, as before.
This is a different GA (a Binary Genetic Algorithm, suited for discrete
problems), but the algorithm is described in-depth in the first Englander
paper and I also can refer to the source code from the second Englander
paper. In the first paper, he mentions that he used Matlab's
implementation. It seems Julia also has an implementation, so this should
be simple in Julia, but I will include significant extra time in Rust, in
case that is more difficult (6/32)
d. From there, calling my outer-loop implementation from the GA should be
simple. However, I will include extra time for the rust implementation,
because I will also have to consider spinning up extra kubernetes
instances for each inner and outer loop. (6/32)
3. From here, the work is more or less done. There will be some finalizing to
do, including setting up the Kubernetes cluster for the Rust implementation.
I've not set up a Kubernetes cluster before, but I have set up docker-compose
setups for Rust code, so I would anticipate an extra 10 hours being needed
for transitioning that. I will also have to provide a sort of wrapper
function that parses the inputs and spins up/calls each lower loop of this
optimizer. That should be really simple, but I'm including some slop time in
this category to account for administrative issues, particularly in setting
up the deployment of the cluster (12/32)
## Summary
So, in total, that comes out to 66 hours of coding for the initial Julia implementation and then another 162 hours for the re-implementation in Rust using docker + kubernetes and the deployment of the cluster. At a rate of 20 hours/week, which seems reasonable, as I will have my job to account for, but I am finished with classes, that comes out to roughly 3.5 weeks to finish the first implementation and then 8 weeks to finalize. I think, if anything, that's very conservative, particularly on the Rust side, as I've got a decent amount of experience in web-based deployment, so I shouldn't have as much trouble with Kubernetes and Docker as I've outlined. But it's also worth noting that the Julia implementation is probably perfectly proficient for a Master's Thesis. If the Rust work is too much, then I can write the thesis on the Julia implementation.

17
readme.md
View File

@@ -2,10 +2,7 @@
This will be a repository for my code while I work on my thesis. The general idea is to generate a This will be a repository for my code while I work on my thesis. The general idea is to generate a
method for automatically producing optimal low-thrust trajectories (including optimizations of the method for automatically producing optimal low-thrust trajectories (including optimizations of the
planetary flybys) for a sample mission to Neptune. planetary flybys) for a sample mission to Saturn.
For now, it will just be a location for storing my notes as I come up with a plan. I will update
this readme as I flush out the plan a little better and then as I begin producing code.
## Dependencies ## Dependencies
@@ -34,8 +31,14 @@ Pkg.build()
To produce the pdfs, simply run: To produce the pdfs, simply run:
```bash ```bash
make make
``` ```
### Julia Code
They will go in the build directory
To produce final pdfs:
```bash
make final
```