\documentclass[defaultstyle,11pt]{thesis}

\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{float}
\usepackage{xfrac}

\title{Designing Optimal Low-Thrust Interplanetary Trajectories Utilizing Monotonic Basin Hopping}
\author{Richard C.}{Johnstone}
\otherdegrees{B.S., Unviersity of Kentucky, Mechanical Engineering, 2016 \\
			  B.S., University of Kentucky, Physics, 2016}
\degree{Master of Science}{M.S., Aerospace Engineering}
\dept{Department of}{Aerospace Engineering}
\advisor{Prof.}{Natasha Bosanac}
\reader{Kathryn Davis}
\readerThree{Daniel Scheeres}

\abstract{  \OnePageChapter

	There are a variety of approaches to finding and optimizing low-thrust trajectories in
	interplanetary space. This thesis analyzes one such approach, namely the application of a
	Monotonic Basin Hopping (MBH) algorithm to a series of Sims-Flanagan transcribed trajectory arcs
	and its applications in a multiple-shooting non-linear solver for the purpose of finding valid
	low-thrust trajectories between planets given poor initial conditions. These valid arcs are then
	fed into the MBH algorithm, which combines them in order to find and optimize interplanetary
	trajectories, given a set of flyby planets. This allows for a fairly rapid searching of a very
	large solution space of low-thrust profiles via a medium fidelity inner-loop solver and a
	well-suited optimization routine. The trajectories found by this method can then be optimized
	further by feeding the solutions back, once again, into the non-linear solver, this time
	allowing the solver to perform optimization.

}

\dedication[Dedication]{
	Dedicated to some people.
}

\acknowledgements{	\OnePageChapter
	This will be an acknowledgement.
}

\emptyLoT

\begin{document}

	\input macros.tex

	\chapter{Introduction}

		Continuous low-thrust arcs utilizing technologies such as Ion propulsion, Hall thrusters, and
		others can be a powerful tool in the design of interplanetary space missions. They tend to be
		particularly suited to missions which require very high total change in velocity ($\Delta
		V$) values and take place over a particularly long duration. Traditional impulsive thrusting
		techniques can achieve these changes in velocity, but typically have a far lower specific
		impulse and, as such, are much less fuel efficient, costing the mission valuable
		financial resources that could instead be used for science. Because of their inherently high
		specific impulse (and thus efficiency), low-thrust techniques are well-suited to interplanetary
		missions.

		The first such mission by NASA to use an electric ion-thruster for an interplanetary mission
		was the Deep Space 1 mission\cite{brophy2002}, which tested the ``new'' technology, first
		appearing as a concept in science fiction stories of the early 1900's and first tested
		successfully during NASA's Space Electric Rocket Test (SERT) mission of
		1964\cite{cybulski1965results}, on an interplanetary mission for the first time. The Ion
		thruster used on Deep Space 1 allowed the mission to rendezvous with both an asteroid (9969
		Braille) and later with a comet (Borrelly), when the technologies being tested, such as the
		ion thruster, proved robust enough and efficient enough to allow for two mission extensions.

		After this initial successful test, ion thrusters and other forms of low-thrust electric
		propulsion have been used in a variety of missions. The NASA Dawn \cite{rayman2006dawn}
		spacecraft in 2015 became the first spacecraft to successfully orbit two planetary bodies,
		thanks in large part to the efficiency of its ion propulsion system. Also notable is the
		joint ESA and JAXA spacecraft Bepi-Colombo\cite{benkhoff2010bepicolombo}, which was launched
		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.

		In general, whenever an interplanetary trajectory is posed, it is advisable to
		first explore the possibility of low-thrust technologies. In an interplanetary mission, the
		primary downside to low-thrust orbits (that they require significant time to achieve large
		$\Delta V$ changes) is made irrelevant by the fact that interplanetary trajectories take such a
		long time as a matter of course.

		Another technique often leveraged by interplanetary trajectory designers is the gravity assist.
		Gravity assists cleverly utilize the inertia of a large planetary body to ``slingshot'' a
		spacecraft, modifying the direction of its velocity with respect to the central body, the Sun.
		This technique lends itself very well to impulsive trajectories. The gravity assist maneuver
		itself can be modeled very effectively by an impulsive maneuver with certain constraints, placed
		right at the moment of closest approach to the (flyby) target body. Because of this,
		optimization with impulsive trajectories and gravity assists are common.

		However, there is no physical reason why low-thrust trajectories can't also incorporate gravity
		assists. The optimization problem simply becomes much more complicated. The separate problems of
		optimizing flyby parameters (planet, flyby date, etc.) and optimizing the low-thrust control
		arcs don't combine very easily. This concept has been explored heavily by Dr. Jacob
		Englander \cite{englander2014tuning}, \cite{englander2017automated},
		\cite{englander2012automated} recently in an effort to develop a generalized and automated
		routine for producing unconstrained, globally optimal trajectories for realistic
		interplanetary mission development that utilizes both planetary flybys and efficient
		low-thrust electric propulsion techniques.

		This thesis will attempt to recreate the techniques developed by Dr Englander in his 2012
		paper, which explored a Hybrid Optimal Control Algorithm for the optimization of
		interplanetary missions\cite{englander2012automated}. Several changes have been made to the
		approach presented in that paper. Most notably, the Evolutionary Algorithm used to automate
		flyby selection has been removed in favor of manual designer flyby selection and the
		addition of a Lambert's Solver was used to improved robustness during the Monotonic Basin
		Hopping step. However, in large, the intention of this thesis is to explore the capabilities
		of this algorithm as an initial design tool for an example mission to Saturn in the near
		future.

		This thesis will explore these concepts in a number of different sections. Section
		\ref{traj_opt} will explore the basic principles of trajectory optimization in a manner agnostic
		to the differences between continuous low-thrust and impulsive high-thrust techniques. Section
		\ref{low_thrust} will then delve into the different aspects to consider when optimizing a low
		thrust mission profile over an impulsive one. Section \ref{interplanetary} provides more detail
		on the interplanetary considerations, including force models and gravity assists. Section
		\ref{algorithm} will cover the implementation details of the optimization algorithm
		developed for this paper. Finally, section \ref{results} will explore the results of some
		hypothetical missions to Saturn.

	\chapter{Trajectory Optimization} \label{traj_opt}

		Trajectory optimization is concerned with finding a narrow solution (namely, optimizing a
		spaceflight trajectory to an end state) with a wide range of possible techniques,
		approaches, and controls making up a very large solution space. In this section, the
		foundations for direct optimization of these sorts of problems will be explored by first
		introducing the Two-Body Problem, then an algorithm for directly solving for states in that
		system, then exploring approaches to Non-Linear Problem (NLP) solving in general and how
		they apply to spaceflight trajectories.

		\section{The Two-Body Problem}
			The motion of a spacecraft in space is governed by a large number of forces. When planning and
			designing a spacecraft trajectory, we often want to use the most complete (and often complex)
			model of these forces that is available. However, in the process of designing these
			trajectories, we often have to compute the path of the spacecraft many hundreds, thousands, or
			even millions of times. Utilizing very high-fidelity force models that account for aerodynamic
			pressures, solar radiation pressures, multi-body effects, and many others may be infeasible
			for the method being used if the computations take too long.

			Therefore, a common approach (and the one utilized in this implementation) is to first look
			simply at the single largest force governing the spacecraft in motion, the gravitational force
			due to the primary body around which it is orbiting. This can provide an excellent
			low-to-medium fidelity model that can be extremely useful in categorizing the optimization
			space as quickly as possible. In many cases, including the algorithm used in this paper, it is
			unlikely that local cost-function minima would be missed due to the lack of fidelity of the
			Two Body Problem.

			In order to explore the Two Body Problem, we must first examine the full set of
			assumptions associated with the force model\cite{vallado2001fundamentals}. Firstly, we
			are only concerned with the nominative two bodies: the spacecraft and the planetary body
			around which it is orbiting. Secondly, both of these bodies are modeled as simple point
			masses. This removes the need to account for non-uniform densities and asymmetry. The
			third assumption is that the mass of the spacecraft ($m_2$) is much much smaller than
			the mass of the planetary body ($m_1$) and enough so as to be considered negligible. The
			only force acting on this system is then the force of gravity that the primary body
			enacts upon the secondary. Lastly, we'll assume a fixed inertial frame. This isn't
			necessary for the formulation of a solution, but will simplify the derivation.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.65\textwidth]{fig/2bp}
				\caption{Figure representing the positions of the bodies relative to each other and
				the center of mass in the two body problem}
				\label{2bp_fig}
			\end{figure}

			Reducing the system to two point masses with a single gravitational force acting between them
			(and only in one direction) we can model the opposite forces on the two bodies as:

			\begin{align}
				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|}
			\end{align}

			Where $\vec{r}$ is the position of the spacecraft relative to the primary body,
			$\vec{r}_1$ is the position of the primary body relative to the center of the inertial
			frame, and $\vec{r}_2$ is the position of the spacecraft relative to the center 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
			then determine the acceleration of the spacecraft relative to the primary planet of the
			system:

			\begin{equation}
				\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|}
			\end{equation}

			Due to our assumption that the mass of the spacecraft is significantly smaller than the
			mass of the primary body ($m_1 >> m_2$) we can simplify the problem by removing the
			negligible $m_2$ term. We can also introduce, for convenience, a number $\mu$ which
			represents, for a given planet, the universal gravitational parameter multiplied by the
			mass of the planet. Doing so and simplifying produces:

			\begin{equation}
				\ddot{\vec{r}} = - \frac{\mu}{r^2} \hat{r}
			\end{equation}

			Where $\mu = G m_1$ is the specific gravitational parameter for our primary body of interest.

			\subsection{Kepler's Laws and Equations}

				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,
				firstly, that the spacecraft's position and velocity (given an initial position and velocity
				and of course the $\mu$ value of the primary body) is actually analytically solvable for all
				future points in time. This can be easily observed by noting that there are three
				one-dimensional equations (one for each component of the three-dimensional space) and
				three unknowns (the three components of the second derivative of the position).

				In the early 1600s, Johannes Kepler produced just such a solution, by taking advantages of
				what is also known as ``Kepler's Laws'' which are\cite{murray1999solar}:

				\begin{enumerate}
						\item Each planet's orbit is an ellipse with the Sun at one of the foci. This can be
							expanded to any orbit by re-wording as ``all orbital paths follow a conic section
							(circle, ellipse, parabola, or hyperbola) with a primary mass at one of the foci''.

							Specifically the path of the orbit follows the trajectory equation:

							\begin{equation}
								r = \frac{\sfrac{h^2}{\mu}}{1 + e \cos(\theta)}
							\end{equation}

							Where $h$ is the angular momentum of the satellite, $e$ is the
							eccentricity of the orbit, and $\theta$ is the true anomaly, or simply
							the angular distance the satellite has traversed along the orbit path.

						\item The area swept out by the imaginary line connecting the primary and secondary
							bodies increases linearly with respect to time. This implies that the magnitude of the
							orbital speed is not constant. For the moment, we'll just take this
							value to be a constant:

							\begin{equation}\label{swept}
								\frac{\Delta t}{T} = \frac{k}{\pi a b}
							\end{equation}

							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
							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
							axis of the orbit, regardless of eccentricity. Specifically, the relationship is:

							\begin{equation}
								T = 2 \pi \sqrt{\frac{a^3}{\mu}}
							\end{equation}

							Where $T$ is the period and $a$ is the semi-major axis.
				\end{enumerate}

		\section{Analytical Solutions to Kepler's Equations}

				Kepler was able to produce an equation to represent the angular displacement of an
				orbiting body around a primary body as a function of time, which we'll derive now
				for the elliptical case\cite{vallado2001fundamentals}. Since the total area of an
				ellipse is the product of $\pi$, the semi-major axis, and the semi-minor axis ($\pi
				a b$), we can relate (by Kepler's second law) the area swept out by an orbit as a
				function of time, as we did in Equation~\ref{swept}.

				This leaves just one unknown variable $k$, which we can determine through use of the
				geometric auxiliary circle, which is a circle with radius equal to the ellipse's semi-major
				axis and center directly between the two foci, as in Figure~\ref{aux_circ}.

				\begin{figure}[H]
					\centering
					\includegraphics[width=0.8\textwidth]{fig/kepler}
					\caption{Geometric Representation of Auxiliary Circle}\label{aux_circ}
				\end{figure}

				In order to find the area swept by the spacecraft, $k$, we can take advantage of the fact
				that that area is the triangle $k_1$ subtracted from the elliptical segment $PCB$:

				\begin{equation}\label{areas_eq}
					k = area(seg_{PCB}) - area(k_1)
				\end{equation}

				Where the area of the triangle $k_1$ can be found easily using geometric formulae:

				\begin{align}
					area(k_1) &= \frac{1}{2} \left( ae - a \cos E \right) \left( \frac{b}{a} a \sin E \right) \\
										&= \frac{ab}{2} \left(e \sin E - \cos E \sin E \right)
				\end{align}

				Now we can find the area for the elliptical segment $PCB$ by first finding the circular
				segment $POB'$, subtracting the triangle $COB'$, then applying the fact that an ellipse is
				merely a vertical scaling of a circle by the amount $\frac{b}{a}$.

				\begin{align}
					area(PCB) &= \frac{b}{a} \left( area(POB') - area(COB') \right) \\
					&= \frac{b}{a} \left( \frac{a^2 E}{2} - \frac{1}{2} \left( a \cos E \right)
					\left( a \sin E \right) \right) \\
					&= \frac{abE}{2} - \frac{ab}{2} \left( \cos E \sin E \right) \\
					&= \frac{ab}{2} \left( E - \cos E \sin E \right)
				\end{align}

				By substituting the two areas back into Equation~\ref{areas_eq} we can get the $k$ area
				swept out by the spacecraft:

				\begin{equation}
					k = \frac{ab}{2} \left( E - e \sin E \right)
				\end{equation}

				Which we can then substitute back into the equation for the swept area as a function of
				time (Equation~\ref{swept}):

				\begin{equation}
					\frac{\Delta t}{T} = \frac{E - e \sin E}{2 \pi}
				\end{equation}

				Which is, effectively, Kepler's equation. It is commonly known by a different form:

				\begin{equation}
					M = \sqrt{\frac{\mu}{a^3}} \Delta t = E - e \sin E
				\end{equation}

				Where we've defined the mean anomaly as $M$ and used the fact that $T =
				\sqrt{\frac{a^3}{\mu}}$. This provides us a useful relationship between Eccentric Anomaly
				($E$) which can be related to spacecraft position, and time, but we still need a useful
				algorithm for solving this equation.

			\subsection{LaGuerre-Conway Algorithm}\label{laguerre}

				For this application, I used an algorithm known as the LaGuerre-Conway
				algorithm\cite{laguerre_conway}, which was presented in 1986 as a faster and more
				robust algorithm for directly solving Kepler's equation and has been in use in many
				applications since. This algorithm is known for its convergence robustness and also
				its speed of convergence when compared to higher order Newton methods.

				This thesis will omit a step-through of the algorithm itself, but psuedo-code for
				the algorithm will be discussed briefly in Section~\ref{conway_pseudocode}.

		\section{Non-Linear Problem Optimization}

			Now we can consider the formulation of the problem in a more useful way. For instance, given a
			desired final state in position and velocity we can relatively easily determine the initial
			state necessary to end up at that desired state over a pre-defined period of time by solving
			Kepler's equation. In fact, this is often how impulsive trajectories are calculated since,
			other than the impulsive thrusting event itself, the trajectory is entirely natural.

			However, often in trajectory design we want to consider a number of other inputs. For
			instance, a low thrust profile, a planetary flyby, the effects of rotating a solar panel on
			solar radiation pressure, etc. Once these inputs have been accepted as part of the model, the
			system is generally no longer analytically solvable, or, if it is, is too complex to calculate
			directly.

			Therefore an approach is needed, in trajectory optimization and many other fields, to optimize
			highly non-linear, unpredictable systems such as this. The field that developed to approach
			this problem is known as Non-Linear Problem (NLP) Optimization.

			There are, however, two categories of approaches to solving an NLP. The first category,
			indirect methods, involve declaring a set of necessary and/or sufficient conditions for declaring
			the solution optimal. These conditions then allow the non-linear problem (generally) to be
			reformulated as a two point boundary value problem. Solving this boundary value problem can
			provide a control law for the optimal path. Indirect approaches for spacecraft trajectory
			optimization have given us the Primer Vector Theory\cite{jezewski1975primer}.

			The other category is the direct methods. In a direct optimization problem, the cost
			function itself is calculated to provide the optimal solution. The problem is usually
			thought of as a collection of dynamics and controls. Then these controls can be modified
			to minimize the cost function. A number of tools have been developed to optimize NLPs
			via this direct method in the general case. For this particular problem, direct
			approaches were used as the low-thrust interplanetary system dynamics adds too much
			complexity to quickly optimize indirectly and the individual optimization routines
			needed to proceed as quickly as possible.

			\subsection{Non-Linear Solvers}
				For these types of non-linear, constrained problems, a number of tools have been developed
				that act as frameworks for applying a large number of different algorithms. This allows for
				simple testing of many different algorithms to find what works best for the nuances of the
				problem in question.

				One of the most common of these NLP optimizers is SNOPT\cite{gill2005snopt}, which
				is a proprietary package written primarily using a number of Fortran libraries by
				the Systems Optimization Laboratory at Stanford University. It uses a sparse
				sequential quadratic programming approach.

				Another common NLP optimization packages (and the one used in this implementation)
				is the Interior Point Optimizer or IPOPT\cite{wachter2006implementation}. It can be
				used in much the same way as SNOPT and uses an Interior Point Linesearch Filter
				Method and was developed as an open-source project by the organization COIN-OR under
				the Eclipse Public License.

				Both of these methods utilize similar approaches to solve general constrained non-linear
				problems iteratively. Both of them can make heavy use of derivative Jacobians and Hessians
				to improve the convergence speed and both have been ported for use in a number of
				programming languages, including in Julia, which was used for this project.

				This is by no means an exhaustive list, as there are a number of other optimization
				libraries that utilize a massive number of different algorithms. For the most part, the
				libraries that port these are quite modular in the sense that multiple algorithms can be
				tested without changing much source code.

			\subsection{Linesearch Method}

				As mentioned above, this project utilized IPOPT which leveraged an Interior Point
				Linesearch method. A linesearch algorithm is one which attempts to find the optimum
				of a non-linear problem by first taking an initial guess $x_k$. The algorithm then
				determines a step direction (in this case through the use of either automatic
				differentiation or finite differencing to calculate the derivatives of the
				non-linear problem) and a step length. The linesearch algorithm then continues to
				step the initial guess, now labeled $x_{k+1}$ after the addition of the ``step''
				vector and iterates this process until predefined termination conditions are met.

				In this case, the IPOPT algorithm was used, not as an optimizer, but as a solver. For
				reasons that will be explained in the algorithm description in Section~\ref{algorithm} it
				was sufficient merely that the non-linear constraints were met, therefore optimization (in
				the particular step in which IPOPT was used) was unnecessary.

			\subsection{Multiple-Shooting Algorithms}

				Now that we have software defined to optimize non-linear problems, what remains is
				determining the most effective way to define the problem itself. The most simple
				form of a trajectory optimization might employ a single shooting algorithm, which
				propagates a state, given some control variables forward in time to the epoch of
				interest. The controls over this time period are then modified in an iterative
				process, using the NLP optimizer, until the target state and the propagated state
				matches. This technique can be visualized in Figure~\ref{single_shoot_fig}.

				\begin{figure}[H]
					\centering
					\includegraphics[width=\textwidth]{fig/single_shoot}
					\caption{Visualization of a single shooting technique over a trajectory arc}
					\label{single_shoot_fig}
				\end{figure}

				In this example, the initial trajectory is the green arc, which contains a certain
				control thrust $\Delta V_{init}$ and is propagated for a certain amount of time and
				results in the end state $x_{init}$. The target state $x_{final}$ can be achieved by
				varying the control and propagating forward in time until this final state is
				achieved. This type of shooting algorithm can be quite useful for simple cases such
				as this one.

				However, some problems require the use of a more flexible algorithm. In these cases,
				sometimes a multiple-shooting algorithm can provide that flexibility and allow the
				NLP solver to find the optimal control faster. In a multiple shooting algorithm,
				rather than having a single target point at which the propagated state is compared,
				the target orbit is broken down into multiple arcs, then end of each of which can be
				seen as a separate target. At each of these points we can then define a separate
				control. The end state of each arc and the beginning state of the next must then be
				equal for a valid arc, as well as the final state matching the target final state.
				This changes the problem to have far more constraints, but also increased freedom
				due to having more control variables.

				\begin{figure}[H]
					\centering
					\includegraphics[width=\textwidth]{fig/multiple_shoot}
					\caption{Visualization of a multiple shooting technique over a trajectory arc}
					\label{multiple_shoot_fig}
				\end{figure}

				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
				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
				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
				description of ways that multiple shooting can be used to optimize a trajectory,
				simply the most convenient for low-thrust arcs.

	\chapter{Low-Thrust Considerations} \label{low_thrust}

		Thus far, the techniques that have been discussed can be equally useful for both impulsive and
		continuous thrust mission profiles. In this section, we'll discuss the intricacies of continuous
		low-thrust trajectories in particular. There are many methods for optimizing such profiles and
		we'll briefly discuss the difference between a direct and indirect optimization of a low-thrust
		trajectory as well as introduce the concept of a control law and the notation used in this
		thesis for modelling low-thrust trajectories more simply.

		\section{Low-Thrust Control Laws}

			In determining a low-thrust arc, a number of variables must be accounted for and, ideally,
			optimized. Generally speaking, this means that a control law must be determined for the
			thruster. This control law functions in exactly the same way that an impulsive thrust
			control law might function. However, instead of determining the proper moments at which
			to thrust, a low-thrust control law must determine the appropriate direction, magnitude,
			and presence of a thrust at each point along its continuous orbit.

			\subsection{Angle of Thrust}

			Firstly, we can examine the most important quality of the low-thrust control law, the
			direction at which to point the thrusters while they are on. 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 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} \hat{H}$ frame. In this
			frame, the $\hat{R}$ direction is the radial direction from the center of the primary to
			the center of the spacecraft. The $\hat{H}$ hat is perpendicular to this, in the
			direction of orbital momentum (out-of-plane) and the $\hat{\theta}$ direction completes
			the right-handed orthonormal frame.

			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
			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
			convention is useful because a $(0,0)$ set represents a thrust force more or less
			directly in line with the direction of the velocity, a commonly useful thrusting
			direction for most effectively increasing (or decreasing if negative) the angular
			momentum and orbital energy of the trajectory.

			Therefore, at each point, the first controls of a control-law, whichever frame or
			convention is used to define them, need to represent a direction in 3-dimensional space
			that the force of the thrusters will be applied.

		\subsection{Thrust Magnitude}

			However, there is actually another variable that can be varied by the majority of
			electric thrusters. Either by controlling the input power of the thruster or the duty
			cycle, the thrust magnitude can also be varied in the direction of thrust, limited by
			the maximum thrust available to the thruster. Not all control laws allow for this
			fine-tuned control of the thruster. Generally speaking, it's most efficient either to
			thrust or not to thrust. Therefore, controlling the thrust magnitude may provide too
			much complexity at too little benefit.

			The algorithm used in this thesis, however, does allow the magnitude of the thrust
			control to be varied. In 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 relatively straightforward to consider the control thrust as a
			3-dimensional vector in space limited in magnitude by the maximum thrust, which allows
			for that increased flexibility.

		\subsection{Thrust Presence}

			The alternative to this approach of modifying the thrust magnitude, is simply to modify
			the presence or absence of thrust. At certain points along an arc, the efficiency of
			thrusting, even in the most advantageous direction, may be such that a thrust is
			undesirable (in that it will lower the overall efficiency of the mission too much) or,
			in fact, be actively harmful.

			For instance, we can consider the case of a simple orbit raising. Given an initial orbit
			with some eccentricity and some semi-major axis, we can define a new orbit that we'd
			like to achieve that simply has a higher semi-major axis value, regardless of the
			eccentricity of the new orbit. It is well known by analysis of the famous Hohmann
			Transfer\cite{hohmann1960attainability}, that thrusting for orbit raising is most
			effective near the periapsis of an orbit, where changes in velocity will have a higher
			impact on total orbital energy. Therefore, for a given low-thrust control law that
			allows for the presence or absence of thrusting at different efficiency cutoffs, we can
			easily come up with two different orbits, each of which achieve the same semi-major
			axis, but in two different ways at two different rates, both in time and fuel use, as
			can be seen in Figures~\ref{low_efficiency_fig} and \ref{high_efficiency_fig}.

			\begin{figure}[H]
				\centering
				\includegraphics[width=\textwidth]{fig/low_efficiency}
				\caption{Graphic of an orbit-raising with a low efficiency cutoff}
				\label{low_efficiency_fig}
			\end{figure}

			\begin{figure}[H]
				\centering
				\includegraphics[width=\textwidth]{fig/high_efficiency}
				\caption{Graphic of an orbit-raising with a high efficiency cutoff}
				\label{high_efficiency_fig}
			\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.

		\section{Direct vs Indirect Optimization}

			As previously mentioned, there are two different approaches to optimizing non-linear
			problems such as trajectory optimizations in interplanetary space. These methods are the
			direct method, in which a cost function is developed and used by numerical root-finding
			schemes to drive cost to the nearest local minimum, and the indirect method, in which a
			set of sufficient and necessary conditions are developed that constrain the optimal
			solution and used to solve a boundary-value problem to find the optimal solution.

			Both of these methods have been applied to the problem of low-thrust interplanetary
			trajectory optimization \cite{Casalino2007IndirectOM}. The common opinion of the
			difference between these two methods is that the indirect methods are more difficult to
			converge and require a better initial guess than the direct methods. However, they also
			require less parameters to describe the trajectory, since the solution of a boundary
			value problem doesn't require discretization of the control states.

			In this implementation, robustness is incredibly valuable, as the Monotonic Basin
			Hopping algorithm is leveraged to attempt to find all minima basins in the solution
			space by ``hopping'' around with different initial guesses. Since these initial guesses
			are not guaranteed to be close to any particular valid trajectory, it is important that
			the optimization routine be robust to poor initial guesses. Therefore, a direct
			optimization method was leveraged by transcribing the problem into an NLP and using
			IPOPT to find the local minima.

		\section{Sims-Flanagan Transcription}

			The major problem with optimizing low thrust paths is that the control law must necessarily be
			continuous. Also, since indirect optimization approaches are, in the context of
			interplanetary trajectories including flybys, quite difficult the problem must
			necessarily be reformulated as a discrete one in order to apply a direct approach. Therefore,
			this thesis chose to use a model well suited for discretizing low-thrust paths: the
			Sims-Flanagan transcription (SFT)\cite{sims1999preliminary}.

			The SFT is actually quite a simple method for discretizing low-thrust arcs. First the
			continuous arc is subdivided into a number ($N$) of individual consistent timesteps of length
			$\frac{tof}{N}$. The control thrust is then applied at the center of each of these time
			steps. This approach can be seen visualized in Figure~\ref{sft_fig}.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.6\textwidth]{fig/sft}
				\caption{Example of an orbit raising using the Sims-Flanagan Transcription with 7
				Sub-Trajectories}
				\label{sft_fig}
			\end{figure}

			Using the SFT, it is relatively straightforward to propagate a state (in the context of the
			Two-Body Problem) that utilizes a continuous low-thrust control, without the need for
			computationally expensive numeric integration algorithms, by simply solving Kepler's equation
			(using the LaGuerre-Conway algorithm introduced in Section~\ref{laguerre}) $N$ times. This
			greatly reduces the computation complexity, which is particularly useful for cases in which
			low-thrust trajectories need to be calculated many millions of times, as is the case in this
			thesis. The fidelity of the model can also be easily fine-tuned. By simply increasing the
			number of sub-arcs, one can rapidly approach a fidelity equal to a continuous low-thrust
			trajectory within the Two-Body Problem, with only linearly-increasing computation time.

	\chapter{Interplanetary Trajectory Considerations} \label{interplanetary}

		The question of interplanetary travel opens up a host of additional new complexities. While
		optimizations for simple single-body trajectories are far from simple, it can at least be
		said that the assumptions of the Two Body Problem remain fairly valid. In interplanetary
		travel, the primary body most responsible for gravitational forces might be a number of
		different bodies, dependent on the phase of the mission. In fact, at some points along the
		trajectory, there may not be a ``primary'' body, but instead a number of different forces of
		roughly equal magnitude vying for ``primary'' status.

		In the ideal case, every relevant body would be considered as an ``n-body'' perturbation
		during the entire trajectory. For some approaches, this method is sufficient and preferred.
		However, for most uses, a more efficient model is necessary. The method of patched conics
		can be applied in this case to simplify the model.

		Interplanetary travel does not merely complicate trajectory optimization. The increased
		complexity of the search space also opens up new opportunities for orbit strategies. The
		primary strategy investigated by this thesis will be the gravity assist, a technique for
		utilizing the gravitational energy of a planet to modify the direction of solar velocity.

		Finally, the concept of multiple gravity assists and the techniques used to visualize the
		space in which we might accomplish stringing together multiple flybys will be analyzed. A
		number of tools have been developed to assist mission designers in manually visualizing the
		search space, but some of these tools can also be leveraged by the automated optimization
		algorithm.

		\section{Launch Considerations}

			Before considering the dynamics and techniques that interplanetary travel imposes upon
			the trajectory optimization problem we must first concern ourself with getting to
			interplanetary space. Generally speaking, interplanetary trajectories require a lot of
			orbital energy and the simplest and quickest way to impart orbital energy to a satellite
			is by using the entirety of the launch energy that a launch vehicle can provide.

			In practice, this value, for a particular mission, is actually determined as a parameter
			of the mission trajectory to be optimized. The excess velocity at infinity of the
			hyperbolic orbit of the spacecraft that leaves the Earth can be used to derive the
			launch energy. This is usually qualified as the quantity $C_3$, which is actually double
			the kinetic orbital energy with respect to the Sun, or simply the square of the excess
			hyperbolic velocity at infinity\cite{wie1998space}.

			This algorithm and many others will take, essentially for granted, that the initial
			orbit at the beginning of the mission 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 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 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$ to.

			A similar approach is taken at the end of the mission. This algorithm, and many others,
			doesn't attempt to exactly match the velocity of the planet at the end of the mission.
			Instead, the excess hyperbolic velocity is also treated as a parameter that can be
			minimized by the cost function. If a mission is to then end in insertion, a portion of
			the mass budget can then be used for an impulsive thrust engine, which can provide a
			final insertion burn at the end of the mission. This approach also allows flexibility
			for missions that might end in a flyby rather than insertion.

		\section{Patched Conics}

			The first hurdle to deal with in interplanetary space is the problem of reconciling
			Two-Body dynamics with the presence of multiple and varying planetary bodies. The most
			common method for approaching this is the method of patched
			conics\cite{bate2020fundamentals}. In this model, we break the interplanetary trajectory
			up into a series of smaller sub-trajectories. During each of these sub-trajectories, a
			single primary is considered to be responsible for the trajectory of the orbit, via the
			Two-Body problem.

			The transition point can be calculated a variety of ways. The most typical method is to
			calculate the gravitational force due to the two bodies separately, via the Two-Body
			models. Whichever primary is a larger influence on the motion of the spacecraft is
			considered to be the primary at that moment. In other words, the spacecraft, at that
			epoch, is within the Sphere of Influence of that primary. Generally for missions in this
			Solar System, the spacecraft is either within the Sphere of Influence of a planetary
			body or the Sun. However, there are points in the Solar System where the gravitational
			influence of two planetary bodies are roughly equivalent to each other and to the
			influence of the Sun. These are considered LaGrange points\cite{euler1767motu}, but are
			beyond the scope of this initial analysis of interplanetary mission feasibility.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.8\textwidth]{fig/patched_conics}
				\caption{Patched Conics Example Figure}
				\label{patched_conics_fig}
			\end{figure}

			This effectively breaks the trajectory into a series of orbits defined by the Two-Body
			problem (conics), patched together by distinct transition points. These transition
			points occur along the spheres of influence of the planets nearest to the spacecraft.
			Generally speaking, for the orbits handled by this algorithm, the speeds involved are
			enough that the orbits are always elliptical around the Sun and hyperbolic around the
			planets.

		\section{Gravity Assist Maneuvers}

			As previously mentioned, there are methods for utilizing the orbital energy of the other
			planets in the Solar System. This is achieved via a technique known as a Gravity Assist,
			or a Gravity Flyby. During a gravity assist, the spacecraft enters into the
			gravitational sphere of influence of the planet and, because of its excess velocity,
			proceeds to exit the sphere of influence. Relative to the planet, the speed of the
			spacecraft increases as it approaches, then decreases as it departs. From the
			perspective of the planet, the velocity of the spacecraft is unchanged. However, the
			planet is also orbiting the Sun.

			From the perspective of a Sun-centered frame, though, this is effectively an elastic
			collision. The overall momentum remains the same, with the spacecraft either gaining or
			losing some in the process (dependent on the directions of travel). The planet also
			loses or gains momentum enough to maintain the overall system momentum, but this amount
			is negligible compared to the total momentum of the planet. The overall effect is that
			the spacecraft arrives at the planet from one direction and, because of the influence of
			the planet, leaves in a different direction\cite{negri2020historical}.

			This effect can be used strategically. The ``bend'' due to the flyby is actually
			tunable via the exact placement of the fly-by in the b-frame, or the frame centered at
			the planet, from the perspective of the spacecraft at $v_\infty$. By modifying the
			turning angle of this bend. In doing so, one can effectively achieve a (restricted) free
			impulsive thrust event.

		\section{Flyby Periapsis}

			Now that we understand gravity assists, the natural question is then how to leverage
			them for achieving certain velocity changes. This can be achieved via a technique called
			``B-Plane Targeting''\cite{cho2017b}. But first, we must consider mathematically the
			effect that a gravity flyby can have on the velocity of a spacecraft as it orbits the
			Sun. Specifically, we can determine the turning angle of the bend mentioned in the
			previous section, given an excess hyperbolic velocity entering the planet's sphere of
			influence ($v_{\infty, in}$) and a target excess hyperbolic velocity as the spacecraft
			leaves the sphere of influence ($v_{\infty, out}$):

			\begin{equation}
				\delta = \arccos \left( \frac{v_{\infty,in} \cdot v_{\infty,out}}{|v_{\infty,in}|
				|v_{\infty,out}|} \right)
			\end{equation}

			From this turning angle, we can also determine, importantly, the periapsis of the flyby
			that we must target in order to achieve the required turning angle. The actual location
			of the flyby point can also be determined by B-Plane Targeting, but this technique was
			not necessary in this implementation as a preliminary feasibility tool, and so is beyond
			the scope of this thesis. The periapsis of the flyby, however, can provide a useful
			check on what turning angles are possible for a given flyby, since the periapsis:

			\begin{equation}
				r_p = \frac{\mu}{v_\infty^2} \left[ \frac{1}{\sin\left(\frac{\delta}{2}\right)} - 1 \right]
			\end{equation}

			Cannot be lower than some safe value that accounts for the radius of the planet and
			perhaps its atmosphere if applicable.

		\section{Multiple Gravity Assist Techniques}

			Now that we can leverage gravity flybys for their change in velocity direction, the
			final remaining question is that of stringing together flybys. How, for instance, do we
			know which planets can provide feasible flyby opportunities given a known hyperbolic
			energy leaving the previous planet?

			\subsection{Lambert's Problem}

				The answer comes from the application of the solution to the problem, posed by
				Johann Lambert in the 18th century\cite{blanchard1969unified}, of how to determine,
				given an initial position, a final position (the ephemeris of the two planets), and
				a time of flight between the two positions, what velocity was necessary to connect
				the two states.

				The actual numerical solution to this boundary value problem is not important to
				include here, but there have been a large number of algorithms written to solve
				Lambert's problem quickly and robustly for given inputs\cite{jordan1964application}.

			\subsection{Planetary Ephemeris}

				Applying Lambert's problem to interplanetary travel requires knowing the positions
				of the planets in the inertial reference frame at a specific epoch. Fortunately,
				many packages have been developed for this purpose. The most commonly used for this
				is the SPICE package, developed by NASA in the 1980's. This software package, which
				has ports that are widely available in a number of languages, including Julia,
				contains many useful functions for astrodynamics. 

				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
				ecliptic plane J2000 frame could be easily determined. A method for quickly
				determining the ephemeris using a polynomial fit was also employed as an option for
				faster ephemeris-finding, but ultimately not used.

			\subsection{Porkchop Plots}

				Armed with a routine for quickly determining the outgoing velocity necessary to
				reach a planet at a given time, as well as the ephemeris of the planets in question
				at any given time, one can produce a grid of departure and arrival times between two
				planetary encounters. Within this grid, one can then plot a variety of useful
				values.

				The solution to Lambert's equation provides both the velocity vectors at departure
				and the velocity vectors at arrival. Often, these will be overlayed on the gridded
				time plots, as normalized values, or sometimes converted to characteristic energy
				$C_3$. This ``porkchop plot'' allows for a quick and concise view of what orbital
				energies are required to reach a planet at a given time from a given location, as
				well as an idea of what outgoing velocities one can expect.

				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
				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
				might be possible for a given timeline, spacecraft state, and planet selection.

				%TODO: Create my own porkchop plot
				\begin{figure}[H]
					\centering
					\includegraphics[width=\textwidth]{fig/porkchop}
					\caption{A sample porkchop plot of an Earth-Mars transfer}
					\label{porkchop}
				\end{figure}

				However, this is an impulsive thrust-centered approach. The solution to Lambert's
				problem assumes a natural trajectory. However, to the low-thrust designer, this is
				needlessly limiting. A natural trajectory is unnecessary when the trajectory can be
				modified by a continuous thrust profile along the arc. Therefore, for the hybrid problem
				of optimizing both flyby selection and thrust profiles, porkchop plots are less helpful,
				and an algorithmic approach is preferred.

	\chapter{Algorithm Overview} \label{algorithm}

		In this section, we will review the actual execution of the algorithm developed. As an
		overview, the routine was developed to enable the determination of an optimized spacecraft
		trajectory from the selection of some very basic mission parameters. Those parameters
		include:

		\begin{itemize}
			\item Spacecraft dry mass
			\item Thruster Specific Impulse
			\item Thruster Maximum Thrusting Force
			\item Thruster Duty Cycle Percentage
			\item Number of Thruster on Spacecraft
			\item Total Starting Weight of the Spacecraft
			\item A Maximum Acceptable $V_\infty$ at arrival and $C_3$ at launch
			\item The Launch Window Timing and the Latest Arrival
			\item A cost function relating the mass usage, $v_\infty$ at arrival, and $C_3$ at
				launch to a cost
			\item A list of flyby planets starting with Earth and ending with the destination
		\end{itemize}

		Which allows for extremely automated optimization of the trajectory, while still providing
		the mission designer with the flexibility to choose the particular flyby planets to
		investigate.

		This is achieved via an optimal control problem in which the ``inner loop'' is a
		non-linear programming problem to determine the optimal low-thrust control law and flyby
		parameters given a suitable initial guess. Then an ``outer loop'' monotonic basin hopping
		algorithm is used to traverse the search space and more carefully optimize the solutions
		found by the inner loop.

		\section{Trajectory Composition}

			In this thesis, a specific nomenclature will be adopted to define the stages of an
			interplanetary mission in order to standardize the discussion about which aspects of the
			software affect which phases of the mission.

			Overall, a mission is considered to be the entire overall trajectory. In the context of
			this software procedure, a mission is taken to always begin at the Earth, with some
			initial launch C3 intended to be provided by an external launch vehicle. This C3 is not
			fully specified by the mission designer, but instead is optimized as a part of the
			overall cost function (and normalized by a designer-specified maximum allowable value).

			This overall mission can then be broken down into a variable number of ``phases''
			defined as beginning at one planetary body with some excess hyperbolic velocity and
			ending at another. The first phase of the mission is from the Earth to the first flyby
			planet. The final phase is from the last flyby planet to the planet of interest.

			Each of these phases are then connected by a flyby event at the boundary. Each flyby
			event must satisfy the following conditions:

			\begin{enumerate}
				\item The planet at the end of one phase must match the planet at the beginning of
					the next phase.
				\item The magnitude of the excess hyperbolic velocity coming into the planet (at the
					end of the previous phase) must equal the magnitude of the excess hyperbolic
					velocity leaving the planet (at the beginning of the next phase).
				\item The flyby ``turning angle'' must be such that the craft maintains a safe
					minimum altitude above the surface or atmosphere of the flyby planet.
			\end{enumerate}

			These conditions then effectively stitch the separate mission phases into a single
			coherent mission, allowing for the optimization of both individual phases and the entire
			mission as a whole. This nomenclature is similar to the nomenclature adopted by Jacob
			Englander in his Hybrid Optimal Control Problem paper, but does not allow for missions
			with multiple targets, simplifying the optimization.

		\section{Inner Loop Implementation}\label{inner_loop_section}

			The optimization routine can be reasonable separated into two separate ``loops'' wherein
			the first loop is used, given an initial guess, to find valid trajectories within the
			region of the initial guess and submit the best. The outer loop is then used to traverse
			the search space and supply the initial loop with a number of well chosen initial
			guesses.

			Figure~\ref{nlp} provides an overview of the process of breaking a mission guess down
			into an NLP, but there are essentially three primary routines involved in the inner
			loop. A given state is propagated forward using the LaGuerre-Conway Kepler solution
			algorithm, which itself is used to generate powered trajectory arcs via the
			Sims-Flanagan transcribed propagator. Finally, these powered arcs are connected via a
			multiple-shooting non-linear optimization problem. The trajectories describing each
			phase complete one ``Mission Guess'' which is fed to the non-linear solver to generate
			one valid trajectory within the vicinity of the original Mission Guess.

			\begin{figure}[H]
				\centering
				\includegraphics[width=\textwidth]{flowcharts/nlp}
				\caption{A flowchart of the Non-Linear Problem Solving Formulation}
				\label{nlp}
			\end{figure}

			\subsection{LaGuerre-Conway Kepler Solver}\label{conway_pseudocode}

				The most basic building block of any trajectory is a physical model for simulating
				natural trajectories from one point forward in time. The approach taken by this
				paper uses the solution to Kepler's equation put forward by
				Conway\cite{laguerre_conway} in 1986 in order to provide simple and very
				processor-efficient propagation without the use of integration. The code logic
				itself is actually quite simple, providing an approach similar to the Newton-Raphson
				approach for finding the roots of the Battin form of Kepler's equation.

				The following pseudo-code outlines the approach taken for the elliptical case. The
				approach is quite similar when $a<0$:

				% TODO: Some symbols here aren't recognized by the font
				\begin{singlespacing}
				\begin{verbatim}
					i = 0
					# First declare some useful variables from the state
					sig0 = (position ⋅ velocity) / √(mu)
					a = 1 / ( 2/norm(position) - norm(velocity)^2/mu )
					coeff = 1 - norm(position)/a

					# This loop is essentially a second-order Newton solver for ΔE
					ΔM = ΔE_new = √(mu/a^3) * time
					ΔE = 1000
					while abs(ΔE - ΔE_new) > 1e-10
						ΔE = ΔE_new
						F = ΔE - ΔM + sig0 / √(a) * (1-cos(ΔE)) - coeff * sin(ΔE)
						dF = 1 + sig0 / √(a) * sin(ΔE) - coeff * cos(ΔE)
						d2F = sig0 / √(a) * cos(ΔE) + coeff * sin(ΔE)
						ΔE_new = ΔE - n*F / ( dF + sign(dF) * √(abs((n-1)^2*dF^2 - n*(n-1)*F*d2F )))
						i += 1
					end

					# ΔE can then be used to determine the F/Ft and G/Gt coefficients
					F = 1 - a/norm(position) * (1-cos(ΔE))
					G = a * sig0/ √(mu) * (1-cos(ΔE)) + norm(position) * √(a) / √(μ) * sin(ΔE)
					r = a + (norm(position) - a) * cos(ΔE) + sig0 * √(a) * sin(ΔE)
					Ft = -√(a)*√(mu) / (r*norm(position)) * sin(ΔE)
					Gt = 1 - a/r * (1-cos(ΔE))

					# Which provide transformations from the original position and velocity to the
					# final
					final_position = F*position + G*velocity
					final_velocity = Ft*position + Gt*velocity
				\end{verbatim}
				\end{singlespacing}

				This approach was validated by generating known good orbits in the 2 Body Problem.
				For example, from the orbital parameters of a certain state, the orbital period can
				be determined. If the system is then propagated for an integer multiple of the orbit
				period, the state should remain exactly the same as it began. In
				Figure~\ref{laguerre_plot} an example of such an orbit is provided.

				\begin{figure}[H]
					\centering
					\includegraphics[width=\textwidth]{fig/laguerre_plot}
					\caption{Example of a natural trajectory propagated via the Laguerre-Conway
					approach to solving Kepler's Problem}
					\label{laguerre_plot}
				\end{figure}
				% TODO: Consider adding a paragraph about the improvements in processor time

			\subsection{Sims-Flanagan Propagator}

				Until this point, we've not yet discussed how best to model the low-thrust
				trajectory arcs themselves. The Laguerre-Conway algorithm efficiently determines
				natural trajectories given an initial state, but it still remains, given a control
				law, that we'd like to determine the trajectory of a system with continuous input
				thrust.

				For this, we leverage the Sims-Flanagan transcription mentioned earlier. This allows
				us to break a single phase into a number of ($n$) different arcs. At the center of
				each of these arcs we can place a small impulsive burn, scaled appropriately for the
				thruster configured on the spacecraft of interest. Therefore, for any given phase,
				we actually split the trajectory into $2n$ sub-trajectories, with $n$ scaled
				impulsive thrust events. As $n$ is increased, the trajectory becomes increasingly
				accurate as a model of low-thrust propulsion in the 2BP. This allows the mission
				designer to trade-off speed of propagation and the fidelity of the results quite
				effectively.

				\begin{figure}[H]
					\centering
					\includegraphics[width=\textwidth]{fig/spiral_plot}
					\caption{An example trajectory showing that classic continuous-thrust orbit
					shapes, such as this orbit spiral, are easily achievable using a Sims-Flanagan
					model}
					\label{sft_plot}
				\end{figure}

				Figure~\ref{sft_plot} shows that the Sims-Flanagan transcription model can be used
				to effectively model these types of orbit trajectories. In fact, the Sims-Flanagan
				model is capable of modeling nearly any low-thrust trajectory with a sufficiently
				high number of $n$ samples.

				Finally, it should be noted that, in any proper propagation scheme, mass should be
				decremented proportionally to the thrust used. The Sims-Flanagan Transcription
				assumes that the thrust for a given sub-trajectory is constant across the entirety
				of that sub-trajectory. Therefore, the mass used by that particular thrusting event
				can be determined by knowledge of the percentage of maximum thrust being provided
				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
				efficiency)\cite{sutton2016rocket}:

				\begin{equation}
					\Delta m = \Delta t \frac{f d}{I_{sp} g_0}
				\end{equation}

				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
				Earth.

			\subsection{Non-Linear Problem Solver}

				Now that we have the basic building blocks of a continuous-thrust trajectory, we can
				leverage one of the many non-linear optimization packages to find solutions near to
				a (proposed) trajectory. This trajectory need not be valid.

				For the purposes of discussion in this Section, we will assume that the inner-loop
				algorithm starts with just such a ''Mission Guess``, which represents the proposed
				trajectory. However, we'll briefly mention what quantities are needed for this
				input:

				A Mission Guess object contains:
				\begin{singlespacing}
				\begin{itemize}
					\item The spacecraft and thruster parameters for the mission
					\item A launch date
					\item A launch $v_\infty$ vector representing excess Earth velocity
					\item For each phase of the mission:
						\begin{itemize}
							\item The planet that the spacecraft will encounter (either flyby or
								complete the mission) at the end of the phase
							\item The $v_{\infty,out}$ vector representing excess velocity at the
								planetary flyby (or launch if phase 1) at the beginning of the phase
							\item The $v_{\infty,in}$ vector representing excess velocity at the
								planetary flyby (or completion of mission) at the end of the phase
							\item The time of flight for the phase
							\item The unit-thrust profile in a sun-fixed frame represented by a
								series of vectors with each element ranging from 0 to 1.
						\end{itemize}
				\end{itemize}
				\end{singlespacing}

				From this information, as can be seen in Figure~\ref{nlp}, we can formulate the
				mission in terms of a non-linear problem. Specifically, the Mission Guess object can
				be represented as a vector, $x$, the propagation function as a function $F$, and the
				constraints as another function $G$ such that $G(x) = \vec{0}$.

				This is a format that we can apply directly to the IPOPT solver, which Julia (the
				programming language used) can utilize via bindings supplied by the SNOW.jl
				package\cite{snow}.

				IPOPT also requires the derivatives of both the $F$ and $G$ functions in the
				formulation above. Generally speaking, a project designer has two options for
				determining derivatives. The first option is to analytically determine the
				derivatives, guaranteeing accuracy, but requiring processor time if determined
				algorithmically and sometimes simply impossible or mathematically very rigorous to
				determine manually. The second option is to numerically derive the derivatives,
				using a technique such as finite differencing. This limits the accuracy, but can be
				faster than algorithmic symbolic manipulation and doesn't require rigorous manual
				derivations.

				However, the Julia language has an excellent interface to a new technique, known as
				automatic differentiation\cite{RevelsLubinPapamarkou2016}. Automatic differentiation
				takes a slightly different approach to numerical derivation. It takes advantage of
				the fact that any algorithmic function, no matter how complicated, can be broken
				down into a series of smaller arithmetic functions, down to the level of simple
				arithmetic. Since all of these simple arithmetic functions have a known derivative,
				we can define a new datatype that carries through the function both the float and a
				second number representing the derivative. Then, by applying (to the derivative) the
				chain rule for every minute arithmetic function derivative as that arithmetic
				function is applied to the main float value, the derivative can be determined,
				accurate to the machine precision of the float type being used, with a processing
				equivalent of two function calls (this of course depends on the simplicity of the
				chained derivatives compared to the function pieces themselves). Generally speaking
				this is much faster than the three or more function calls necessary for accurate
				finite differencing and removes the need for the designer to tweak the epsilon value
				in order to achieve maximum precision.

		\section{Outer Loop Implementation}

			Now we have the tools in place for, given a potential ''mission guess`` in the
			vicinity of a valid guess, attempting to find a valid and optimal solution in that
			vicinity. Now what remains is to develop a routine for efficiently generating these
			random mission guesses in such a way that thoroughly searches the entirety of the
			solution space with enough granularity that all spaces are considered by the inner loop
			solver.

			Once that has been accomplished, all that remains is an ''outer loop`` that can generate
			new guesses or perturb existing valid missions as needed in order to drill down into a
			specific local minimum. In this thesis, that is accomplished through the use of a
			Monotonic Basin Hopping algorithm. This will be described more thoroughly in
			Section~\ref{mbh_subsection}, but Figure~\ref{mbh_flow} outlines the process steps of
			the algorithm.

				\begin{figure}[H]
					\centering
					\includegraphics[width=\textwidth]{flowcharts/mbh}
					\caption{A flowchart visualizing the steps in the monotonic basin hopping
					algorithm}
					\label{mbh_flow}
				\end{figure}

			\subsection{Random Mission Generation}\label{random_gen_section}

				At a basic level, the algorithm needs to produce a mission guess (represented by all
				of the values described in Section~\ref{inner_loop_section}) that contains random
				values within reasonable bounds in the space. This leaves a number of variables open
				to for implementation. For instance, it remains to be determined which distribution
				function to use for the random values over each of those variables, which bounds to
				use, as well as the possibilities for any improvements to a purely random search.

				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
				based propagator. For this implementation, that was chosen to be 20, based upon a
				number of tests in which the calculation time for the propagation was compared
				against the accuracy of a much higher $n$ value for some known thrust controls, such
				as a simple spiral orbit trajectory. This value of 20 tends to perform well and
				provide reasonable accuracy, without producing too many variables for the NLP
				optimizer to control for (since the impulsive thrust at the center of each of the
				sub-trajectories is a control variable). This leaves some room for future
				improvements, as will be discussed in Section~\ref{improvement_section}.

				The bounds for the launch date are provided by the user in the form of a launch
				window, so the initial launch date is just chosen as a standard random value from a
				uniform distribution within those bounds.

				A unit launch direction is then also chosen as a 3-length vector of uniform random
				numbers and normalized. This vector is then multiplied by a uniform random number
				between 0 and the root of the maximum launch $C_3$ specified by the user to generate
				an initial $\vec{v}_\infty$ vector at launch.

				Next, the times of flight of each phase of the mission is then decided. Since launch
				date has already been selected, the maximum time of flight can be calculated by
				subtracting the launch date from the latest arrival date provided by the mission
				designer. Then, each leg is chosen from a uniform distribution with bounds currently
				set to a minimum flight time of 30 days and a maximum of 70\% of the maximum time of
				flight. These leg flight times are then iteratively re-generated until the total
				time of flight (represented by the sum of the leg flight times) is less than the
				maximum time of flight. This allows for a lot of flexibility in the leg flight
				times, but does tend toward more often producing longer missions, particularly for
				missions with more flybys.

				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
				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
				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
				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
				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
				available in Section~\ref{mbh_subsection}.

				However, for the truly random mission guess, there are still the $v_\infty$ values
				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
				velocity was calculated. However, instead of multiplying the randomly generate unit
				direction vector by a random number between 0 and the square root of the maximum
				launch $C_3$, bounds of 0 and 10 kilometers per second are used instead, to provide
				realistic flyby values.

				The outgoing excess hyperbolic velocity at infinity is then calculated by again
				choosing a uniform random unit direction vector, then by multiplying this value by
				the magnitude of the incoming $v_{\infty}$ since this is a constraint of a
				non-powered flyby.

				From these two velocity vectors the turning angle, and thus the periapsis of the
				flyby, can then be calculated by the following equations:

				\begin{align}
					\delta &= \arccos \left( \frac{\vec{v}_{\infty,in} \cdot
					\vec{v}_{\infty,out}}{|v_{\infty,in}| \cdot {|v_{\infty,out}}|} \right) \\
					r_p &= \frac{\mu}{\vec{v}_{\infty,in} \cdot \vec{v}_{\infty,out}} \cdot \left(
					\frac{1}{\sin(\delta/2)} - 1 \right)
				\end{align}

				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.
				Since the thrust is defined as a 3-vector of values between -1 and 1 representing
				some percentage of the full thrust producible by the spacecraft thrusters in that
				direction, the initial thrust controls can then be generated as a $20 \times 3$
				matrix of uniform random numbers within that bound.

			\subsection{Monotonic Basin Hopping}\label{mbh_subsection}

				Now that a generator has been developed for mission guesses, we can build the
				monotonic basin hopping algorithm. Since the implementation of the MBH algorithm
				used in this paper differs from the standard implementation, the standard version
				won't be described here. Instead, the variation used in this paper, with some
				performance improvements, will be considered.

				The aim of a monotonic basin hopping algorithm is to provide an efficient method for
				completely traversing a large search space and providing many seed values within the
				space for an ''inner loop`` solver or optimizer. These solutions are then perturbed
				slightly, in order to provide higher fidelity searching in the space near valid
				solutions in order to fully explore the vicinity of discovered local minima. This
				makes it an excellent algorithm for problems with a large search space, including
				several clusters of local minima, such as this application.

				The algorithm contains two loops, the size of each of which can be independently
				modified (generally by specifying a ''patience value``, or number of loops to
				perform, for each) to account for trade-offs between accuracy and performance depending on
				mission needs and the unique qualities of a certain search space.

				The first loop, the ''search loop``, first calls the random mission generator. This
				generator produces two random missions as described in
				Section~\ref{random_gen_section} that differ only in that one contains random flyby
				velocities and control thrusts and the other contains Lambert's-solved flyby
				velocities and zero control thrusts. For each of these guesses, the NLP solver is
				called. If either of these mission guesses have converged onto a valid solution, the
				lower loop, the ''drill loop`` is entered for the valid solution. After the
				convergence checks and potentially drill loops are performed, if a valid solution
				has been found, this solution is stored in an archive. If the solution found is
				better than the current best solution in the archive (as determined by a
				user-provided cost function of fuel usage, $C_3$ at launch, and $v-\infty$ at
				arrival) then the new solution replaces the current best solution and the loop is
				repeated. Taken by itself, the search loop should quickly generate enough random
				mission guesses to find all ''basins`` or areas in the solution space with valid
				trajectories, but never attempts to more thoroughly explore the space around valid
				solutions within these basins.

				The drill loop, then, is used for this purpose. For the first step of the drill
				loop, the current solution is saved as the ''basin solution``. If it's better than
				the current best, it also replaces the current best solution. Then, until the
				stopping condition has been met (generally when the ''drill counter`` has reached
				the ''drill patience`` value) the current solution is perturbed slightly by adding
				or subtracting a small random value to the components of the mission.

				The performance of this perturbation in terms of more quickly converging upon the
				true minimum of that particular basin, as described in detail by
				Englander\cite{englander2014tuning}, is highly dependent on the distribution
				function used for producing these random perturbations. While the intuitive choice
				of a simple Gaussian distribution would make sense to use, it has been found that a
				long-tailed distribution, such as a Cauchy distribution or a Pareto distribution is
				more robust in terms of well chose boundary conditions and initial seed solutions as
				well as more performant in time required to converge upon the minimum for that basin.

				Because of this, the perturbation used in this implementation follows a
				bi-directional, long-tailed Pareto distribution generated by the following
				probability density function:

				\begin{equation}
					1 +
					\left[ \frac{s}{\epsilon} \right] \cdot
					\left[ \frac{\alpha - 1}{\frac{\epsilon}{\epsilon + r}^{-\alpha}} \right]
				\end{equation}

				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
				distributed random array with dimension equal to the perturbed variable and bounds
				of 0 and 1, $\epsilon$ is a small value (nominally set to $1e-10$), and $\alpha$ is
				a tuning parameter to determine the size of the tails and width of the distribution
				set to $1.01$, but easily tunable.

				The perturbation function then steps through each parameter of the mission,
				generating a new guess with the parameters modified by the Pareto distribution.
				After this perturbation, the NLP solver is then called again to find a valid
				solution in the vicinity of this new guess. If the solution is better than the
				current basin solution, it replaces that value and the drill counter is reset to
				zero. If it is better than the current total best, it replaces that value as well.
				Otherwise, the drill counter increments and the process is repeated. Therefore, the
				drill patience allows the mission designer to determine a maximum number of
				iterations to perform without improvement in a row before ending the drill loop.
				This process can be repeated essentially ''search patience`` number of times in
				order to fully traverse all basins.

	\chapter{Results Analysis} \label{results}

		The algorithm described in this thesis is quite flexible in its design and could be used as
		a tool for a mission designer on a variety of different mission types. However, to consider
		a relatively simple but representative mission design objective, a sample mission to Saturn
		was investigated.

		Ultimately, two optimized trajectories were selected. The results of those trajectories can
		be found in Table~\ref{results_table} below:

		\begin{table}[h!]
		\begin{small}
		\centering
		\begin{tabular}{ | c c c c c c | }
			\hline
			\bfseries Flyby Selection &
			\bfseries Launch Date &
			\bfseries Mission Length &
			\bfseries Launch $C_3$ &
			\bfseries Arrival $V_\infty$ &
			\bfseries Fuel Usage \\
			&  & (years) & $\left( \frac{km}{s} \right)^2$ & ($\frac{km}{s}$) & (kg) \\
			\hline
			EMS & 2024-06-27 & 7.9844 & 60.41025 & 5.816058 & 446.9227 \\
			EMJS & 2023-11-08 & 14.1072 & 40.43862 & 3.477395 & 530.6683 \\
			\hline
		\end{tabular}
		\end{small}
		\caption{Comparison of the two most optimal trajectories}
		\label{results_table}
		\end{table}

		\section{Mission Constraints}

			The sample mission was 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
			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
			initial mission designer (though not necessarily fixed in the design either) and is that
			of the spacecraft parameters. The application accepts as input a spacecraft object
			containing: the dry mass of the craft, the fuel mass at launch, the number of onboard
			thrusters, 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
			a fuel mass of 3300 kilograms. This was chosen in order to have an overall mass roughly
			in the same zone as that of the Cassini spacecraft, which launched with 5712 kilograms
			of total mass, with the fuel accounting for 2978 of those kilograms\cite{cassini}. The
			dry mass of the craft 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
			delivered dry mass to Saturn could be thought of as a metric of success, without
			discounting mission that may have delivered just under whatever more realistic dry mass
			one might set, 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
			250 millinewtons, and a 100\% duty cycle. This puts the thruster roughly in line with
			having an array of three NSTAR ion thrusters, which were used on the Dawn and Deep Space
			1 missions\cite{polk2001performance}.

			Also of relevance to the mission were the maximum $C_3$ at launch and $v_\infty$ at
			arrival values. In order to not exclude the possibility of a non-capture flyby mission,
			it was decided to not include the arrival $v_\infty$ term in the cost function and,
			because of this, the maximum value was set to be extremely high at 500 kilometers per
			second, in order to fully explore the space. In practice, though, the algorithm only
			looks at flybys below 10 kilometers per second in magnitude. The maximum launch $C_3$
			energy was set conservatively to 200 kilometers per second squared. This is upper limit
			is only possible, for the given start mass, using a heavy launch system such as the
			SLS\cite{stough2021nasa} or the comparable SpaceX Starship, though at values below about
			half of this maximum, it begins to become possible to use existing launch solutions.

			Finally, the mission is meant to represent a near future mission. Therefore the launch
			window was set to allow for a launch in any day in 2023 or 2024 and a maximum total time
			of flight of 20 years. This is longer than most typical Saturn missions, but allows for
			some creative trajectories for higher efficiency.

			It should be noted that each of these trajectories was found using an $n$ value of 20 as
			mentioned previously, but in post-processing, the trajectory was refined to utilize a
			slightly higher fidelity model that uses 60 sub-trajectories per orbit. This serves to
			provide better plots for display, higher fidelity analyses, as well as to highlight the
			efficacy of the lower fidelity method. Orbits can be found quickly in the lower fidelity
			model and easily refined later by re-running the NLP solver at a higher $n$ value.

			\subsection{Cost Function}

				Each mission optimization also allows for the definition of a cost function. This
				cost function accepts as inputs all parameters of the mission, the maximum $C_3$ at
				launch and the maximum excess hyperbolic velocity at arrival.

				The cost function used for this mission first generated normalized values for fuel
				usage and launch energy. The fuel usage number is determined by dividing the fuel
				used by the mass at launch and the $C_3$ number is determined by dividing the $C_3$
				at launch by the maximum allowed. These two numbers are then weighted, with the fuel
				usage value getting a weight of three and the launch energy value getting a weight
				of one. The values are summed and returned as the cost value.

			\subsection{Flybys Analyzed}

				Since the algorithm itself makes no decisions on the actual choice of flybys, that
				leaves the mission designer to determine which flyby planets would make good
				potential candidates. A mission designer can then re-run the algorithm for each of
				these flyby plans and determine which optimized trajectories best fit the needs of
				the mission.

				For this particular mission scenario, the following flyby profiles were
				investigated:

				\begin{itemize}
					\item EJS
					\item EMJS
					\item EMMJS
					\item EMS
					\item ES
					\item EVMJS
					\item EVMS
					\item EVVJS
				\end{itemize}

		\section{Faster, Less Efficient Trajectory}

			In order to showcase the flexibility of the optimization algorithm (and the chosen cost
			function), two different missions were chosen to highlight. One of these missions is a
			slower, more efficient trajectory more typical of common low-thrust trajectories. The
			other is a faster trajectory, quite close to a natural trajectory, but utilizing more
			launch energy to arrive at the planet.

			It is the faster trajectory that we'll analyze first. Most interesting about this
			particular trajectory is that it's actually quite efficient despite its speed, in
			contrast to the usual dichotomy of low-thrust travel. The cost function used for this
			analysis did not include the time of flight as a component of the overall cost, and yet
			this trajectory still managed to be the lowest cost trajectory of all trajectories found
			by the algorithm.

			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 launch energy required to effectively used the gravity assist with
			Mars at this time is quite high. The $C_3$ value was found to be $60.4102$ kilometers
			per second squared. However, 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 the in-plane thrust direction was neither
			zero nor $\pi$, implying that these thrusts were steering thrusts rather than
			momentum-increasing thrusts.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.9\textwidth]{fig/EMS_plot}
				\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
				in which the spacecraft approached them}
				\label{ems}
			\end{figure}

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.9\textwidth]{fig/EMS_plot_noplanets}
				\caption{Another depiction of the EMS trajectory, without the planetary ephemeris
				arcs}
				\label{ems_nop}
			\end{figure}

			The second and final leg of this trip exits the Mars flyby and, initially burns quite
			heavily along the velocity vector in order to increase it's semi-major axis. After an
			initial period of thrusting, though, the spacecraft effectively coasts with minor
			adjustments until its rendezvous with Saturn just four and a half years later in June of
			2032. The arrival $v_\infty$ is not particularly small, at $5.816058$ kilometers per
			second, but this is to be expected as the arrival excess velocity was not considered as
			a part of the cost function. If capture was not the final intention of the mission, this
			may be of little concern. Otherwise, the low fuel usage of $446.92$ kilograms for a
			$3500$ kilogram launch mass leaves much margin for a large impulsive thrust to enter
			into a capture orbit at Saturn.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.9\textwidth]{fig/EMS_thrust_mag}
				\caption{The magnitude of the unit thrust vector over time for the EMS trajectory}
				\label{ems_mag}
			\end{figure}

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.9\textwidth]{fig/EMS_thrust_components}
				\caption{The inertial x, y, and z components of the unit thrust vector over time for
				the EMS trajectory}
				\label{ems_components}
			\end{figure}

			In this case the algorithm effectively realized that a higher-powered launch from
			the Earth, then a natural coasting arc to Mars flyby would provide the spacecraft with
			enough velocity that a short but efficient powered-arc to Saturn was possible with
			effective thrusting. It also determined that the most effective way to achieve this
			flyby was to increase orbital energy in the beginning of the arc, when increasing
			the semi-major axis value is most efficient. All of these concepts are known to skilled
			mission designers, but finding a trajectory that combined all of these concepts would
			have required much time-consuming analysis of porkchop plots and combinations of
			mission-design techniques. This approach is far more automatic than the traditional
			approach.

			The final quality to note with this trajectory is that it shows a tangible benefit of
			the addition of the Lambert's solver in the monotonic basin hopping algorithm. Since the
			initial arc is almost entirely natural, with very little thrust, it is extremely likely
			that the trajectory was found in the Lambert's Solution half of the MBH algorithm
			procedure.

		\section{Slower, More Efficient Trajectory}

			Next we'll analyze the nominally second-best trajectory. While the cost function
			provided to the algorithm can be a useful tool for narrowing down the field of search
			results, it can also be very useful to explore options that may or may not be of similar
			"efficiency" in terms of the cost function, but beneficial for other reasons. By
			outputting many different optimal trajectories, the MBH algorithm can allow for this
			type of mission design flexibility.

			To highlight the flexibility, a second trajectory has been selected, which has nearly
			equal value by the cost function, coming in slightly lower. However, this trajectory
			appears to offer some benefits to the mission designer who would like to capture into
			the gravitational field of Saturn or minimize launch energy requirements, perhaps for a
			smaller mission, at the expense of increased speed.

			The first leg of this three-leg trajectory is quite similar to the first leg of the
			previous trajectory. However, this time the launch energy is considerably lower, with a
			$C_3$ value of only $40.4386$ kilometer per second squared. Rather than employ an almost
			entirely natural coasting arc to Mars, however, this trajectory performs some thrusting
			almost entirely in the velocity direction, increasing its orbital energy in order to
			achieve the same Mars rendezvous. In this case, the launch was a bit earlier, occurring
			in November of 2023, with the Mars flyby occurring in mid-April of 2026. This will prove
			to be helpful in comparison with the other result, as this mission profile is much
			longer.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.9\textwidth]{fig/EMJS_plot}
				\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
				the phase in which the spacecraft approached them}
				\label{emjs}
			\end{figure}

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.9\textwidth]{fig/EMJS_plot_noplanets}
				\caption{Another depiction of the EMJS trajectory, without the planetary ephemeris
				arcs}
				\label{emjs_nop}
			\end{figure}

			The second phase of this trajectory also functions quite similarly to the second phase
			of the previous trajectory. In this case, there is a little bit more thrusting necessary
			simply for steering to the Jupiter flyby than was necessary for Saturn rendezvous in the
			previous trajectory. However, most of this thrusting is for orbit raising in the
			beginning of the phase, very similarly to the previous result. In this trajectory, the
			Jupiter flyby occurs late July of 2029.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.9\textwidth]{fig/EMJS_thrust_mag}
				\caption{The magnitude of the unit thrust vector over time for the EMJS trajectory}
				\label{emjs_mag}
			\end{figure}

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.9\textwidth]{fig/EMJS_thrust_components}
				\caption{The inertial x, y, and z components of the unit thrust vector over time for
				the EMJS trajectory}
				\label{emjs_components}
			\end{figure}

			Finally, this mission also has a third phase. The Jupiter flyby provides quite a strong
			$\Delta V$ for the spacecraft, allowing the following phase to largely be a coasting arc
			to Saturn almost one revolution later. During the most efficient part of the journey,
			some thrust in the velocity direction accounts for a little bit of orbit-raising, but
			the phase is largely natural. Because of this long coasting period, the mission length
			increases considerably during this leg, arriving at Saturn in December of 2037, over 8
			years after the Jupiter flyby.

			However, there are many advantages to this approach relative to the other trajectory.
			While the fuel use is also slightly higher at $530.668$ kilograms, plenty of payload
			mass is still capable of delivery into the vicinity of Saturn. Also, it should be noted
			that the incoming excess hyperbolic velocity at arrival to Saturn is significantly
			lower, at only $3.4774$ kilometers per second, meaning that less of the delivered
			payload mass would need to be taken up by impulsive thrusters and fuel for Saturn orbit
			capture, should the mission designer desire this.

			Also, as mentioned before, the launch energy requirements are quite a bit lower. Having
			a second mission trajectory capable of launching on a smaller vehicle could be valuable
			to a mission designer presenting possibilities. According to an analysis of the Delta IV
			and Atlas V launch configurations\cite{c3capabilities} in Figure~\ref{c3}, this
			reduction of $C_3$ from around 60 to around 40 brings the sample mission to just within
			range of both the Delta IV Heavy and the Atlas V in its largest configuration, neither
			of which are possible for the other result, meaning that either different launch
			vehicles must be found or mission specifications must change.

			\begin{figure}[H]
				\centering
				\includegraphics[width=\textwidth]{fig/c3}
				\caption{Plot of Delta IV and Atlas V launch vehicle capabilities as they relate to
				payload mass}
				\label{c3}
			\end{figure}

	\chapter{Conclusion} \label{conclusion}

		\section{Overview of Results}

			A mission designer's job is quite a difficult one and it can be very useful to have
			tools to automate some of the more complex analysis. This paper attempted to explore one
			such tool, meant for automating the initial analysis and discovery of useful
			interplanetary, low-thrust trajectories including the difficult task of optimizing the
			flyby parameters. This makes the mission designer's job significantly simpler in that
			they can simply explore a number of different flyby selection options in order to get a
			good understanding of the mission scope and search space for a given spacecraft, launch
			window, and target.

			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
			resulting trajectories have trade-offs in mission length, launch energy, fuel usage, and
			more. However, both results should be considered very useful low-thrust trajectories in
			comparison to other missions that have launched on similar interplanetary trajectories,
			using both impulsive and low-thrust arcs with planetary flybys. Each of these missions
			should be feasible or nearly feasible (feasible with some modifications) using existing
			launch vehicle and certainly even larger missions should be reasonable with advances in
			launch capabilities currently being explored.

		\section{Recommendations for Future Work}\label{improvement_section}

			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 papers\cite{englander2014tuning}\cite{englander2017automated}
			\cite{englander2012automated} in which he explored the hybrid optimal control problem of
			multi-objective low-thrust interplanetary trajectories.

			In light of this, there are a number of additional approaches that Englander took in
			preparing his algorithm that were not implemented here in favor of reducing complexity
			and time constraints. For instance, many of the Englander papers explore the concept of
			an outer loop that utilizes a genetic algorithm to compare many different flyby planet
			choice against each other. This would create a truly automated approach to low-thrust
			interplanetary mission planning. However, a requirement of this approach is that the
			monotonic basin hopping algorithm algorithm must converge on optimal solutions very
			quickly. Englander typically runs his for 20 minutes each for evolutionary fitness
			evaluation, which is over an order of magnitude faster than the implementation in this
			paper to achieve satisfactory results.

			Further improvements to performance stem from the field of computer science. An
			evolutionary algorithm such as the one proposed by Englander would benefit from high
			levels of parallelization. Therefore, it would be worth considering a GPU-accelerated or
			even cluster-computing capable implementation of the monotonic basin hopping algorithm.
			These cluster computing concepts scale very well with new cloud infrastructures such as
			that provided by AWS or DigitalOcean.

			Finally, the monotonic basin hopping algorithm as currently written provides no
			guarantees of actual global optimization. Generally optimization is achieved by running
			the algorithm until it fails to produce newer, better trajectories for a sufficiently
			long time. But it would be worth investigating the robustness of the NLP solver as well
			as the robustness of the MBH algorithm basin drilling procedures in order to quantify
			the search granularity needed to completely traverse the search space. From this
			information, a new MBH algorithm could be written that is guaranteed to explore the
			entire space.

	\bibliographystyle{plain}
	\nocite{*}
	\bibliography{thesis}

	\appendix

\end{document}
