\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}