\chapter{Interplanetary Trajectory Design} \label{traj_dyn} In order to optimize trajectories in interplanetary space, it is necessary to understand and utilize the dynamical systems that will be acting on the spacecraft throughout the trajectory. This section will explore the models available for understanding the natural motion of a spacecraft in interplanetary space, beginning with the Two Body Problem and a method of patched conics for combining the gravitational effects of multiple primary bodies. \section{System Dynamics} \subsection{The Two-Body Problem} The motion of a spacecraft in space is governed by a complex dynamical environment. However, in the process of designing a trajectory, 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 pressure, solar radiation pressure, multi-body effects, and other forces may be too time intensive for a particular application. Initial surveys of the solution space often don't require such complex models in order to gain valuable insight. Therefore, a common approach (and the one utilized in this implementation) is to first use a lower-fidelity dynamical model that captures only the gravitational force due to the primary body around which the spacecraft is orbiting. This approach can provide an excellent low-to-medium fidelity model that is useful as an underlying model in an algorithm for quickly categorizing a search space for initial mission feasibility explorations. 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 point masses with constant mass. This removes the need to account for non-uniform densities and asymmetry. Finally, for convenience in notation at the end, we'll also assume 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. \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} Under these assumptions, the force acting on the body due to the law of universal gravitation is: \begin{align} F_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\ F_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|} \end{align} And by Newton's second law (force is the product of mass and acceleration), we can derive the following differential equations for $r_1$ and $r_2$: \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 origin 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 planet: \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} Further assuming 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 gravitational parameter $\mu$ which represents the gravity constant for the system about the center of motion ($\mu = G (m_1 + m_2) \approx G m_1$). Doing so and simplifying produces: \begin{equation} \ddot{\vec{r}} = - \frac{\mu}{r^2} \hat{r} \end{equation} \subsubsection{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} \subsection{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. \subsubsection{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{Interplanetary 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. \subsection{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. \subsection{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. \subsection{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. \subsection{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. \subsection{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? \subsubsection{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}. \subsubsection{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. \subsubsection{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. \section{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. \subsection{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. \subsubsection{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. \subsubsection{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. \subsubsection{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. \subsection{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. \subsection{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.