From 5c0232dbd7e0be38e2e3bf9a8ce2528ce2af3c4d Mon Sep 17 00:00:00 2001 From: Connor Date: Sat, 8 Jan 2022 19:21:58 -0700 Subject: [PATCH] Finished section 2! --- LaTeX/flowcharts/kepler.svg | 535 ++++++++++++++++++++++++++++++++++++ LaTeX/thesis.tex | 174 ++++++++++-- julia/find_neptune.jl | 2 +- 3 files changed, 692 insertions(+), 19 deletions(-) create mode 100644 LaTeX/flowcharts/kepler.svg diff --git a/LaTeX/flowcharts/kepler.svg b/LaTeX/flowcharts/kepler.svg new file mode 100644 index 0000000..25ab028 --- /dev/null +++ b/LaTeX/flowcharts/kepler.svg @@ -0,0 +1,535 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + a + + ae + acos(E) + E + υ + + rcos(υ) + r + rsin(υ) + k + k2 + B + C + Primary + P + O + B' + a + b + + asin(E) + + diff --git a/LaTeX/thesis.tex b/LaTeX/thesis.tex index c8194df..8df1e3a 100644 --- a/LaTeX/thesis.tex +++ b/LaTeX/thesis.tex @@ -130,6 +130,8 @@ Monotonic Basin Hopping} \newpage \tableofcontents + + \listoffigures \newpage @@ -230,37 +232,173 @@ Monotonic Basin Hopping} Where $\mu = G m_1$ is the specific gravitational parameter for our primary body of interest. - \subsubsection{Kepler's Equations} + \subsubsection{Kepler's Laws and Equations} - Now that we've fully qualified the forces acting within the Two Body Problem, we can note - that the Problem is actually analytically solvable in the case when the position of the - spacecraft and the $\mu$ value of the primary body are known. This can be easily observed by - noting that there are three one-dimensional equations (one for each component of the - three-dimensional position) and three unknowns (the three components of the second - derivative of the position). + % TODO: Can I segue better from 2BP to Keplerian geometry? - Therefore, we can use this analytically solvable force model to model the spacecraft's - motion in time, a more useful re-interpretation of the equations of motion. + Now that we've fully qualified the forces acting within the Two Body Problem, we can concern + ourselves with more practical applications of this 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 position) 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: + + \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 body at one of the foci''. + \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. + \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: $T = 2 + \pi \sqrt{\frac{a^3}{\mu}}$ where $T$ is the period and $a$ is the semi-major axis. + \end{enumerate} \subsection{Analytical Solutions to Kepler's Equations} - Discuss how, since the 2BP is analytically solvable, there exists algorithms for solving - these 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. 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: + + \begin{equation}\label{swept} + \frac{\Delta t}{T} = \frac{k}{\pi a b} + \end{equation} + + 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} + \centering + \includegraphics[width=0.8\textwidth]{LaTeX/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 $OB'C$, 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( \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{align} + M &= E - e \sin E \\ + &= \sqrt{\frac{\mu}{a^3}} \Delta t + \end{align} + + 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} - Step through the LaGuerre-Conway algorithm, as it's used in this implementation. + For this application, I used an algorithm known as the LaGuerre-Conway algorithm, which was + presented in 1986 as a faster 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 the code will be present + in the Appendix. \subsection{Non-Linear Problem Optimization} - Discuss the concept of a non-linear problem as it related to the study of 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. In these cases, generally + speaking, direct methods are impossible, so a number of tools have been developed to optimize + NLPs in the general case, via indirect, iterative methods. \subsubsection{Non-Linear Solvers} - Mention common non-linear solvers (Ipopt, SNOPT, etc.). + 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. - \subsubsection{Newton Trust-Region Method} - Discuss how the Newton Trust-Region method works. + One of the most common of these NLP optimizers is SNOPT, 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. 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. \subsubsection{Linesearch Method} - I may take this section out, because I'm not currently using a linesearch. But I would cover - the additions of linesearch methods. + 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 automatic differentiation 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. \section{Low-Thrust Considerations} \label{low_thrust} Highlight the differences between high and low-thrust mission profiles. diff --git a/julia/find_neptune.jl b/julia/find_neptune.jl index 58bebc9..491809d 100644 --- a/julia/find_neptune.jl +++ b/julia/find_neptune.jl @@ -55,7 +55,7 @@ end log(_::Nothing, _::Vector{Mission}, _::Vector{Body}) = println("No Mission Found...") -planets = [Earth, Venus, Venus Jupiter, Saturn] +planets = [Earth, Venus, Venus, Jupiter, Saturn] best, archive = mbh(planets) log(best, archive, planets)