thesis/LaTeX/trajectory_design.tex

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

		\subsection{Kepler's Laws}

			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{Kepler's Equation}

			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\cite{vallado2001fundamentals}, $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}

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

			\noindent
			Where $E$, notably, is the eccentric anomaly, or the angular distance from the
			periapsis to the vertical projection of the spacecraft on the auxiliary circle. 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}) for period of time since the spacecraft left periapsis:

			\begin{equation}
				\frac{\Delta t}{T} = \frac{t_2 - t_{peri}}{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 in order to use this equation to propagate a
			spacecraft.

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

		For this thesis, the algorithm used to solve Kepler's equation was the general numeric
		root-finding scheme first developed by LaGuerre in the 1800s and first applied to
		Kepler's equation by Bruce Conway in 1985\cite{laguerre_conway}. In his paper, Conway
		makes a compelling argument for utilizing the less common LaGuerre method over higher
		order Newton or Newton-Raphson methods.

		The Newton-Raphson methods, while found to generally have quite impressive convergence
		rates (generally successfully solving Kepler's equation correctly within 5 iterations),
		were prone to failures in convergence given certain specific initial conditions.
		Therefore LaGuerre's algorithm is proposed as an alternative.

		The algorithm can be relatively easily derived by examining the polynomial equation with
		$m$ roots:

		\begin{equation}
			g(x) = (x - x_1) (x - x_2) ... ( x - x_m)
		\end{equation}

		\noindent
		We can then generate some useful convenience functions as:

		\begin{align}
			\ln|g(x)| &= \ln|(x - x_1)| + \ln|(x - x_2)| + ... + \ln|( x - x_m)| \\
			\frac{d\ln|g(x)|}{dx} &= \frac{1}{x - x_1} + \frac{1}{x - x_2} + ... + \frac{1}{x -
			x_m} = G_1(x)
		\end{align}

		and

		\begin{align}
			\frac{-d^2\ln|g(x)|}{dx^2} &= \frac{1}{(x - x_1)^2} + \frac{1}{(x - x_2)^2} + ... +
			\frac{1}{(x - x_m)^2} = G_2(x)
		\end{align}

		Now we define the targeted root as $x_1$ and make the approximation that all of the
		other roots are equidistant from the targeted root, which means:

		\begin{equation}
			x - x_i = b, i=2,3,...,m
		\end{equation}

		\noindent
		We can then rewrite $G_1$ and $G_2$ as:

		\begin{align}
			G_1 &= \frac{1}{a} + \frac{n-1}{b} \\
			G_2 &= \frac{1}{a^2} + \frac{n-1}{b^2}
		\end{align}

		\noindent
		Which may be solved for $a$ in terms of $G_1$, $G_2$:

		\begin{equation}
			a = \frac{n}{G_1 \pm \sqrt{(n-1)(nG_2 - G_1^2)}}
		\end{equation}

		\noindent
		With corresponding iteration function:

		\begin{equation}
			x_{i+1} = x_i - \frac{n g(x_i)}{g'(x_i) \pm \sqrt{(n-1)^2 f'(x_i)^2 - n (n-1) f(x_i)
			f''(x_i)}}
		\end{equation}

		This iteration scheme can be shown to be globally convergent, regardless of the initial
		guess. More relevantly, Conway also showed that the application of this method to
		Kepler's equation was shown to converge with similar speed to many of the best common
		higher order Newton-Raphson solvers. However, LaGuerre's method was also found to be
		incredibly robust, converging to the correct value for every one of Conway's 500,000
		tests. Because of this robustness, it is very useful for propagating spacecraft states.

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