This work presents a novel approach for designing preliminary fuel-optimal low-thrust spacecraft trajectories with gravity assists (GAs) by employing Physics-Constrained Neural Networks (PCNNs). The study introduces the Multiple-Leg Trajectory PCNN (MLT-PCNN) which embeds physical dynamics and constraints directly into the neural network architecture and computes planar fuel-optimal low-thrust trajectories with a gravity assist assuming restricted two-body motion. It models a GA as a discontinuous change in the spacecraft's heliocentric velocity vectors, arising from an instantaneous turning of the spacecraft’s hyperbolic excess velocity vectors with respect to the GA body. The model’s performance is demonstrated using Earth–Mars–Ceres and Earth-Venus-Earth-Mars-Jupiter case studies, where MLT-PCNNs are trained with various training schedules using both first-order (Adam) and second-order (L-BFGS) optimizers. Physical accuracy of MLT-PCNN-generated solutions is evaluated by error metrics that measure discrepancies in final position, velocity, and mass, obtained through numerical integration of the MLT-PCNN generated control profile. A control profile correction method based on the state transition matrix validates that only a few grams of additional fuel are required to eliminate residual errors in the final states. The results show that the MLT-PCNN can generate fuel-optimal solutions comparable or superior to state-of-the-art techniques.

The optimization of low-thrust transfer trajectories presents a significant computational challenge due to the vast search space involved. This paper introduces a novel clustering-based successive search space pruning technique to improve low-thrust trajectory optimization. Regions of low ∆V trajectories were extracted from the search space using DBSCAN. Through successive pruning of regions of the search space containing inefficient solutions, effectively zooming in on the most promising areas, the performance of global optimization algorithms was improved. The effectiveness of the search space pruning was evaluated using three interplanetary transfer test cases: Earth-Mars, Earth-Tempel 1, and Earth-Mercury. The performance of three global optimization algorithms, Differential Evolution (DE), Simple Genetic Algorithm (SGA), and Particle Swarm Optimization (PSO) was compared across the original and pruned search spaces. Results show that clustering-based pruning can significantly reduce the search space and subsequently lead to improved optimal transfer trajectories with fewer function evaluations. Differential Evolution outperformed PSO and SGA across all three test cases, demonstrating robustness and efficiency. The results of the optimization highlight the importance of the selection of the population size and number of generations for the optimization algorithms, because insufficient function evaluations resulted in poor convergence for PSO and SGA, including in the reduced search spaces. The results also suggest excessive pruning of the search space can lead to loss of optimal solution. The findings confirmed that clustering-based search space pruning is a promising technique for improving the efficiency of global optimization in low-thrust trajectory design.

With a dedicated mission to the Martian moons, Phobos and Deimos, set to launch soon, there is a growing interest in further exploring these moons using low-thrust propulsion. This paper investigates the trajectory design for a sample-return mission from Deimos using fuel-free solar sail propulsion technology, aiming to maximize operational time near Deimos within a minimum permissible total mission duration. Time-optimal transfers between Earth and Deimos are sought by formulating and solving an optimal control problem using a direct pseudospectral method. Initial guesses for the direct method are generated by considering a patched circular restricted three-body problem (CR3BP) approximation and by searching for heteroclinic-like connections between the Sun-Earth and Mars-Deimos systems using the differential evolution algorithm. The obtained solution, with a maximum mission duration set to eight years based on the insights from initial guess generation, results in an optimal duration of stay at Deimos of 329 days with a mission duration of 7.7 years. Although the patched CR3BP approximation demonstrated valid transfer solutions for this study, it is deemed computationally inefficient for future trajectory designs for similar mission concepts targeting either only Phobos or both moons at once. Nevertheless, the trajectories obtained back and forth from Deimos are sufficiently optimal for a preliminary mission concept and validate the feasibility of achieving such a mission employing a solar sail.

With the development of space research into novel areas, new complex problems arise. The interest in solving space routing problems considering large numbers of targets has recently grown. This paper proposes a novel method to solve the optimal trajectory in such combinatorial space routing problems. This paper focuses on a global optimisation algorithm implemented to solve the problem posed in the 4th Global Trajectory Optimisation Com- petition (GTOC4). The solution is a trajectory of multiple legs, where each leg links two targets and has a specific flight time. To enable the use of the powerful mixed-integer linear problem solver software, Solving Constraint Integer Programs (SCIP), the routing problem concerned with visiting as many target bodies with a predetermined fuel and time budget is split into linear sub-problems. The Fixed Budget sub-problem selects a subset of the given set of targets. The Full Tour sub-problem orders the targets in the subset, and the Fixed Tour sub-problem optimises the flight time for every leg of the given trajectory to find the solution with the lowest total fuel consumption. Each of these sub-problems is formulated in a linear form and is solved using SCIP. The global optimisation algorithm evolves a population where every individual exists of a set of initial guesses for the time of flight values. Analysis shows that initialising this population with a mix of randomly generated individuals and individuals containing a constant value for all entries leads to the fastest convergence towards the optimal solution. In a population of 20, seeding ten individuals is found to be optimal. It is also found that the algorithm performance can be further increased by evolving individuals with infeasible solutions instead of iterating them until a feasible solution is found and eliminating the Full Tour sub-problem. These simplifications allow for an increase in the cost budget multiplier, which leads to finding better objective values without further increasing computational time. The best-performing setup, which uses a cost budget multiplier of 10, can find the optimal solution to the test problem in 100% of the runs, on average in 9 iterations, with a computation time of 5.82 seconds per evaluation. The results show that the global optimisation algorithm produces results that closely match known results for GTOC4 consistently and accurately.

Using low-thrust propulsion for interplanetary space missions has the potential to allow for more payload for the same mass put into orbit compared to what impulsive propulsion would allow for. The disadvantages are found in mission planning, however, as the continuous nature of the thrust yields a more complex problem. One potential solution to help in the early planning and discovery phase of mission design is to employ artificial neural networks (ANNs). This has been done in the past, yet only in a limited capacity. Specifically, the engineering of the feature space used with the neural network has never been investigated. This thesis attempts to provide a first look at the influence of different feature space compositions. This includes the use of nine different state representations but also an analysis of additional values in the feature space. Additionally, the effect of extraneous variables, one of them notably being the target of interest, on the neural network performance is analyzed. The dataset used is generated using indirect optimization, and the case investigated is a set of minimum fuel Earth-Mars transfer trajectories.

Low-thrust spacecraft trajectories are for space exploration what neural networks are for computer science: The vanguard of current trends with a lot of potential. Only recently have the two ideas, trajectory optimization and machine learning, been combined. In all the publications making use of machine learning for low-thrust optimization, a clear gap exists, however: Feature engineering has never been investigated. This thesis attempts to provide a first patch for that gap, limiting itself in scope to interplanetary Earth-Mars trajectories and feedforward neural networks. The data used as the basis to evaluate the performance of a number of factors having a potential influence on the choice of feature is obtained through indirect optimization. A novel method to generate those trajectories is implemented. The data is then used to investigate the effect different targets and network parameters have on the choice of features. On the feature side, a significant number of state representations are analyzed, both in dimensional and nondimensionalized form. Additionally, the feature space is expanded by additional variables, and various transformations are attempted.

Over the course of this work, the importance of properly scaled data has been demonstrated. It is also shown that using Keplerian state and costate as feature and target, respectively, reliably yields good results. When mass is estimated, fuel mass is preferable over total spacecraft mass. Finally, none of the additional parameters or transformations (besides nondimensionalization) attempted resulted in reliable improvements and are thus best avoided.

This thesis proposes an unsupervised Physics-Informed Neural Network (PINN) for solving optimal control problems with the direct method to design and optimize transfer trajectories. The network adheres analytically to boundary conditions and includes the objective fitness as regularization in its loss function. A test scenario of a planar Earth-Mars low-thrust optimal-fuel transfer and rendezvous is chosen. Comprehensive examination of training strategies reveals that convergence is highly dependent on the initialization of the network and that correctly balancing loss terms is essential for navigating the intricate loss landscape. This balance is achieved by carefully selecting loss weights and implementing a refined learning rate schedule. Comparative analysis to hodographic shaping solutions demonstrates that the PINN effectively identifies near-optimal solutions across a wide range of initial and final constraints for the Earth-Mars transfer problem, with a maximum improvement of 4.5 km s⁻¹ and median improvement of 0.55 km s⁻¹. The PCNN shows promise as a preliminary design tool for trajectory optimization in nonlinear dynamics.

Shape-based methods are used in the preliminary optimization of low-thrust trajectories to rapidly search large design spaces and provide initial guesses for higher fidelity methods. The optimization process benefits from the shape-based methods providing initial guesses as quickly and as optimal as possible. This work aims to improve the hodographic shaping method by implementing machine learning (ML) to optimize its free parameters. The addition of free parameters enables a more optimal trajectory, while the ML model reduces the computational effort required to optimize this trajectory. ML incorporated shaping models with 6 and 9 Degrees of Freedom (DoF) and models with different levels of mission generalization are built. The models are tested and compared to the shaping method without ML on different types of single transfers in a grid search and on more complex missions in a genetic algorithm. The ML 6-DoF models can consistently find trajectories that are more optimal compared to the 0-DoF models without ML, by up to 6 km/s for an Earth-Ceres transfer, and they can do so 50-100 times faster than the 6-DoF models without ML. The 9-DoF models produce inconsistent results and can only improve no-ML shaping in areas of low optimality. The ML models significantly accelerate the performance of the genetic algorithm when the relevant transfer bodies are included in the ML training data.

The optimization of interplanetary, low-­thrust trajectories is a computationally expensive aspect of preliminary mission design. To reduce the computational burden associated with it, surrogate models can be used as cheap approximations of the original fitness function. Training the surrogate models in a fully online manner can be done to remove the need of having previously generated datasets, which is another source of computational cost. The Sims­-Flanagan transcription is used to model an Earth-Mars transfer which is optimized through different optimization routines. The development of a C++ library with machine learning tooling was initiated, containing implementations for Generalized Regression Neural Networks (GRNNs) and Radial Basis Function Networks (RBFNs) that are used in global and local surrogates, respectively, having their hyperparameters tuned through cross-­validation. A surrogate model was constructed using Differential Evolution (DE) operators and an uncertainty-based infill criterion for the global search phase, and approximation of the derivative of the original fitness function which is provided to SNOPT (Sparse Nonlinear Optimizer), in the local search phase. An ablation study was performed to assess how each of the components of the surrogate model contribute to the results. It was verified that neither the derivative information nor the local search as a whole led to better results. The surrogate model was also outperformed by the standard optimization strategy found in literature, Monotonic Basin Hopping (MBH). Two new surrogate models incorporating ideas of this strategy were created, with one of them outperforming every other model that was tested. Despite not having performed a full study of the computational effort due to the simulations having been run in a server with a variable load, the new models present better results for similar amounts of fitness function evaluations. A Wilcoxon rank-­sum test was performed to assess whether the results have statistical significance, leading to the conclusion that a surrogate model can be used to improve the optimization of low-­thrust trajectories modeled with the Sims­-Flanagan transcription when inserted in a monotonic basin hopping optimization scheme.

The use of low-thrust propulsion for interplanetary missions requires the implementation of new methods for the preliminary design of their trajectories. This thesis proposes a method using the Monotonic Basin Hopping global optimization algorithm to find feasible trajectories with optimum use of the mass of fuel for the case in which the trajectory is modeled using the Sims-Flanagan transcription method. Due to the large computational time required to find the global optimum, Artificial Neural Networks have been used to predict the objective value and feasibility terms of the local minimum. Therefore, the procedure to set up a working regression Artificial Neural Network is studied as well as its transferability to predict values outside the trained limits and for different missions. In addition to this, the use of pre-training is analyzed to improve the performance of the network without increasing the size of the training database.

Many contemporary interplanetary missions use efficient low-thrust engines to reach the far corners of our Solar System. Their trajectories, however, have proven to be complicated to optimise due to the non-impulsive manoeuvres involved in low-thrust spaceflight. Even though shaping methods have been used extensively to reduce the computational burden, multiple-gravity assists and the presence of constraints create significant computational hurdles. Reducing the number of fitness evaluations during optimisation is one way of speeding up the search and can be done by `pruning away' regions of infeasible trajectories. In this research, we approach this by applying clustering, an unsupervised machine learning approach, to single leg trajectory optimisation problems based on a hodographic shaping trajectory model in combination with restricted two-body dynamics. Through clustering, groups of promising trajectories can be isolated so that unwanted regions can be discarded. Earth --- Mars, Earth --- Venus, and Earth --- 9P/Tempel 1 trajectories are used as test cases and are shown to exhibit periodic behaviour (related to the synodic periods of the departure and target bodies), which enables clustering on grid search-generated datasets. Different clustering algorithms were compared using the Silhouette, Davies-Bouldin, and Calinski-Harabasz internal validation indices. However, traditional clustering algorithms such as (H)DBSCAN, OPTICS, KMeans, and Gaussian Mixture Models, failed to robustly provide clusterings that can be used for pruning, because of the oblong shape of the clusters and the absence of data/noise density differences due to the artificial nature of the problem. Instead, a multimodality-based clustering model called SkinnyDip was found to be much more promising for this task. This algorithm comes with the additional advantage of having very few hyperparameters, eliminating the need for extensive parameter tuning.

Various machine learning algorithms have been applied to find optimal lowthrust
satellite trajectories, however, no fair comparison of their accuracy has been made yet. In this paper, two common and promising supervised machine learning algorithms are compared for their regression capacities:
the artificial neural network and the Gaussian process. A grid search is performed to evaluate the performance of the algorithms on different datasets, varying the input features for training, and the hyperparameters of the algorithms. The best performing Gaussian process and artificial neural network are applied as surrogate in a model that optimizes a trajectory from Earth to Mars using a differential evolution optimization strategy. The performance of the models is evaluated on both the Euclidean distance between the inputs corresponding to the predicted minima of the surrogate method and the nearest local minimum obtained with the shaping method, and the accuracy of the predicted Δ푉 budget required. For both quantities, the Gaussian process outperforms the artificial neural network. The minimum required Δ푉 budget to reach Mars was predicted more accurately, and the duration of the trajectory and the best moment to launch were found more often and accurately as well.

In this thesis, a new method to approximate the cost function of Low-Thrust, Multiple-Gravity-Assist interplanetary trajectories using a Machine Learning surrogate is proposed. This method speeds up the optimization process without fine tuning of the surrogate parameters for every individual case. The computational cost of obtaining training data was identified as the main limitation when using Machine Learning methods for this purpose. Therefore, the surrogate was built with an Online Sequential Extreme Learning Machine Multi-Agent System (OS-ELM-MAS) due to its theoretical good performance when the training data is limited. A high-fidelity global optimization problem was implemented, and a method to include the surrogate during the optimization process was designed. This method does not require specialized optimization algorithms. The parameters that control the interaction between the surrogate and the optimization process were identified and a procedure to obtain the best values was designed and applied. The final results show that the use of the surrogate improves the optimization results when evaluations of the cost function are computationally expensive. However, the values of the parameters that control the interaction between the surrogate and the optimization algorithm had to be carefully selected. The search for a general procedure to obtain these parameters without repeated tests is proposed for future research. Several applications to new optimization problems of the method developed in the thesis are also proposed for future research.

Preliminary design of low-thrust trajectories in the circular restricted three-body problem (CR3BP) frequently relies upon ballistic dynamical structures and optimization algorithms. A fundamental understanding of how these dynamical structures change due to presence of a low-thrust force may lead to trajectories that cannot be obtained otherwise. This paper investigates the effect of a constant low-thrust acceleration on the horizontal Lyapunov (H-L) families in the CR3BP. Families of low-thrust periodic solutions are constructed in vicinity of L1 and L2 using numerical continuation methods. By either varying the Hamiltonian, acceleration magnitude, or acceleration orientation along the solution family, the effect of a low-thrust acceleration on H-L orbits is characterized. Investigating the geometry, bifurcations and hyperbolic unwinding behavior of these families provides insight into the low-thrust periodic solution structure of the Earth-Moon system. The introduction of a constant low-thrust acceleration distorts the geometry of ballistic H-L orbits into ’ear-shaped’ periodic solutions. The bifurcations of the low-thrust periodic solution families imply the existence of low-thrust halo, low-thrust axial, and low-thrust planar double-period families. Finally, low-thrust periodic solutions are identified that possess a higher rate of hyperbolic unwinding behavior than the ballistic L1 and L2 H-L families.

The Kuiper Belt is considered to be formed by remnants of the original Solar System, that is why exploration missions to that region are susceptible of having an enormous scientific impact. Nonetheless, it is far away from Earth, missions to KBO targets a significant challenge. Additionally, current methods to identify trajectories that encounter multiple KBOs are computationally intensive with impractically long run times on the order of months.

The present work deals with applying a rapid low-fidelity technique to identify candidate preliminary KBO sequences, to be used as a starting point for futuremission designers to identify trajectories tomultiple KBOs. This pathfinding approach uses two consecutive grid search types.

First, awell-tuned grid search is implemented using Lambert arcs to reach a first KBO, including the introduction of a new, systematic approach to identify the step sizes in target-body arrival date for the Lambert-based grid search. Then, second and third KBO encounters are assessed from the results of the first grid search and second KBO search respectively. Both second and third KBO searches are performed with a new algorithm consisting of a time grid along with STM linear propagation for maneuver calculation. Finally, an example trajectory resulting from this technique that encounters two KBOs is given as a potential flyable route.

Low-thrust trajectories can benefit the search for propellant-optimal trajectories, but increases in modeling complexity and computational load remain a challenge for efficient mission design and optimization. An approach for developing models utilizing Gaussian Process (GP) regression and classification is proposed to perform computationally efficient optimization while obtaining acceptable accuracies for trajectories based on exponential sinusoid shaping. The goal of this work is to predict a combination of values of input variables which corresponds to a shape-based trajectory with the smallest total velocity increment (∆V) or propellant mass fraction (Jm). A GP classification model is constructed to assess whether a given combination of values of input variables corresponds to a feasible trajectory. GP regression models are developed to predict the total ∆V and Jm corresponding to a combination of shape parameters, which can replace the required integration along the shape. In addition, advanced regression models are developed to predict the target values while requiring only three input parameters, thereby replacing the entire shape computation. In order to develop a GP model that fits the problem at hand, the underlying functions and parameters should be selected rationally. A novel model development approach is proposed to ensure that the mean function, covariance function, likelihood function, inference method, and hyperparameters, which dominate the performance of the models, are chosen rationally in terms of mean absolute percentage error (MAPE) and prediction time. Using this approach, GP models are developed and tested on transfer trajectories from Earth to Mars and Ceres, and from Mars to Earth, and their performance, in terms of MAPE and prediction time, is compared to that of more common optimization techniques in combination with the exponential sinusoid and other shape-based methods. The results demonstrate that the computation time can significantly be reduced while achieving promising MAPE's, especially when the goal is to locate regions of feasible or near-optimal trajectories. The proposed model development procedure is tested for robustness, which provides confidence in the proposed approach. Furthermore, it is found that the models which map three input variables directly to a ∆V or Jm value perform better than the ones trained with shape information, which demonstrates the strength of GP models as applied to low-thrust trajectory optimization.

With all major bodies within the Solar System explored by at least a single fly-by, modern-day missions are becoming increasingly more demanding, up to a point where classical chemical propulsion can no longer supply the required ∆V. Increasingly more is relied upon low-thrust propulsion, characterised by its (very) low thrust force; long continuous thrust arcs, often lasting months at a time; and high specific impulse. To even further increase the possible ∆V budget, thereby allowing more intricate missions, more payload, or a lower transfer time; use is made of low-thrust propulsion combined with gravity assists.

With traditional low-thrust gravity assist optimisation tools heavily relying on astrodynamics and optimal control theory expertise, often requiring an initial guess and heavy modification for each new mission scenario; there is a need for a smart low-thrust gravity assist trajectory optimisation tool. Such a tool should be independent from an initial guess, optimise a trajectory from only a broad description of the mission, and be applicable to a wide variety of mission scenarios.

It was the goal of this thesis to develop such a smart low-thrust gravity assist trajectory optimisation tool, which was found in extending the global low-thrust optimisation software package InTrance. InTrance tackles the problem from the novel perspective of artificial intelligence and machine learning using a method termed evolutionary neurocontrol (ENC), which combines biologically inspired artificial neural networks (ANNs) with evolutionary algorithms (EAs). The internal parameters of the ANN and initial conditions of each phase are optimised by the EA, which then serves as an agent; supplying the spacecraft with a steering strategy at each integration step.

A novel tool capable of optimising the low-thrust gravity assist problem has been developed and has been shown to find more optimal results than available reference trajectories in select cases. Two prominent low-thrust gravity assist missions have been re-optimised: a double-asteroid rendezvous mission with an intermediate Mars gravity assist similar to the Dawn mission, and a low-thrust version of the New Horizons mission to Pluto with an intermediate Jupiter gravity assist. InTrance found a low-thrust New Horizons like trajectory which used 41% less propellant for a similar flight time as the actual mission, even though having a higher drymass and being launched with a lower C3. The results of the Dawn re-optimisation showed excellent agreement with the actual Dawn mission in terms of flight and dwell times, and used significantly less propellant. The efficiency increase due to the inclusion of the gravity assists has furthermore been investigated. The gravity assist at Jupiter in the New Horizons trajectory resulted in a decrease of 3.5% in flight time and over 75% in propellant saving. Dawn’s mission could be flown without a gravity assist w.r.t flight time and propellant usage, however, the inclusion resulted in an increase of 20% in dwell time at Vesta.

Building on recent advances in the fields of low-thrust trajectory optimization based on shaping methods, Artificial Neural Networks, and surrogate models in Evolutionary Algorithms, an investigation into a novel optimization routine is conducted. A flexible Python tool to evaluate linked trajectories in a two-body model based on the hodographic shaping method is implemented and used to develop an evolutionary optimization approach where a Genetic Algorithm is assisted in finding new candidate solutions by a surrogate model. This surrogate is constructed from previous fitness function evaluations using Machine Learning, specifically by training an Artificial Neural Network. After deriving suitable (hyper-)parameters for the Genetic Algorithm and the Artificial Neural Network an experimental investigation into the algorithm's performance is conducted with a focus on the design of the surrogate for low-thrust trajectory problems. Two example problems based on the Dawn trajectory and the GTOC2 problem are studied. The surrogate approach is able to find good new candidate solutions, i.e. solutions that improve the population's overall fitness, especially when the surrogate is designed to approximate the shaping computation. Additionally, the use of a surrogate pretrained on a general data set of low-thrust transfers is tested and found to considerably improve the initial quality of the model, meaning that more good candidate solutions are found early on, accelerating the algorithm's convergence.

The PocketQube is an emerging satellite class, which pushes the miniaturization of space technology beyond the well-established CubeSats, promising rapid design-to-orbit cycles while lowering the cost of accessing space. A showstopper in the success story of nano- and picosatellites are their high mission failure rates. Environmental testing before flight is an effective means to identify design flaws and workmanship errors, and thus improve the chances of mission success. However, inefficiencies in the design process, low-budgets, and stringent schedule requirements often motivate small satellite developers to postpone environmental testing towards the end of the design lifecycle, where recovering for design flaws is inefficient. The project circumstance of nano- and picosatellite missions require, therefore, cost- and time-effective test strategies that facilitate early design evaluation. This research proposes and implements a thermal screening method for PocketQube subsystems to identify temperature hotspots and verify their compliance against operational hardware limits. Key elements of the test method are a thermal IR temperature scan at ambient conditions and an estimation of the worst-case flight temperature using experimentally derived graphs that describe the vacuum heating of thermal hotspots. The study of subsystem layout options complements the screening method by providing solutions to mitigate hotspot overheating. Moreover, the study proposes to lower the required pressure levels for thermal-vacuum testing of PocketQube subsystems. The analysis shows that, due to the small form factor, pressure levels by four orders of magnitude larger than those used in environmental test standards for larger satellites suffice to maintain the resulting temperature errors below 5K on the hot side of the temperature spectrum. Both the screening method and moderation in vacuum requirements contribute to the development of subsystem test methods that match the needs of small satellite developers.

In interplanetary space missions, it is convenient to have a second departure opportunity in case the first is missed. Two distinct approaches to minimizing the maximum of the two Delta-V budgets of such a trajectory pair, are developed. The first (‘a priori’) approach optimizes the variables of both trajectories at once. The second (‘a posteriori’) approach first minimizes Delta-V budgets for a range of discrete departure epochs, and then selects the pair of which the highest Delta-V is minimum. Furthermore, five different pruning and biasing methods are developed, these prove critical for computational efficiency (number of objective function evaluations). Application to three different gravity-assist (and deep space maneuver) trajectories to Saturn, reveals that the a priori approach is more computationally efficient on a trajectory with few variables (3) and that the a posteriori approach is more computationally efficient on a trajectory with many variables (22).

The three-body problem (3BP) formulated by I. Newton has inspired many great mathematicians like L. Euler (1707-1783), J.L. Lagrange (1736-1813), K.G.J. Jacobi (1804-1851), W.R. Hamilton (1805-1865) and J.H. Poincaré (1854-1912) to develop mathematical studies, methods and theories in an attempt to solve this problem. Although these concepts have been crucial to the advancement of the natural sciences, a closed form solution has yet to be found. To aid in solving this problem, a restricted form of the 3BP has been formulated in which the mass of the two primaries greatly exceeds the mass of the third body. Known as the circular restricted three-body problem, this system gives rise to three collinear and two equilateral equilibria referred to as libration points. These locations are surrounded by various families of periodic libration point orbits. The corresponding exotic trajectories are highly non-Keplerian which offer desirable characteristics that have revolutionised space mission design. Crucial to this effort has been the theory on hyperbolic invariant manifolds. These topological structures asymptotically arrive at (stable) or depart from (unstable) a selected target orbit by exploiting the natural dynamics of the system. A connection which joins stable and unstable manifold trajectories with no discrepancy in state space, constitutes to a free (natural) transfer mechanism. In the case that two different equilibria are connected, this path is called a heteroclinic connection. The potential existence of these solutions are of significant scientific interest and are known to exist in the planar case. In this research, the current theory is extended by providing a comprehensive understanding of the phase space of the (spatial) vertical Lyapunov (V-L) orbits and their associated manifolds including the potential existence of heteroclinic connections in the Earth-Moon system.The transition from the planar to the spatial case is associated with a rapid increase in complexity. To reduce the order of the problem, the two hyperbolic manifolds which are to be connected are integrated until they intersect a plane called a Poincaré section. These stopping conditions form a N-1-dimensional subset of the phase space. A completely new effort is aimed at investigating the influence of the orientation of this section on the state vector discrepancy that arises from the connection of hyperbolic manifolds emanating from collinear libration points L1 and L2. The aim of this study is to reveal novel transfer trajectories outside the solution space confined by the traditional locations and orientations of the Poincaré section. The path towards creating these insights consist of three stages, all of which employ a variety of numerical methods. As these steps are sequential, the individual chapters of this thesis cannot be considered as stand-alone work. First of all, the target orbits are generated through the refinement (Differential Correction) and extension (Numerical Continuation) of analytic approximations based on Perturbation Theory. The families of horizontal Lyapunov (H-L), halo, axial, and V-L are analysed to provide an overview of the range of possible solutions. The periodicity of each orbit is verified numerically and the trajectory is validated by studying the so-called monodromy matrix. Secondly, the manifolds are constructed using the eigenspace of this fundamental matrix. Once more a high level of generality is achieved, this time by analysing the sets of manifolds associated with three different families of target orbits at three distinct energy levels. Each member of each manifold is verified through studying deviations in the Jacobi's constant, and validated by using symmetries. Lastly, the stopping conditions for integration of these hyperbolic trajectories are varied to analyse the sensitivity with respect to the orientation of the Poincaré section.This research has revealed that the most optimal connections for V-L orbits are found outside of the traditional stopping conditions. The ramifications of this behaviour are extensive and are observed to arise from the fact that these hyperbolic trajectories curve behind the Moon. Moreover, the orthographic projections display interesting differences in global stability. The manifolds associated with the V-L family retain their shape, whereas those corresponding to the H-L and halo families behave in a chaotic way over time. This enables an accurate approximation over extended integration periods. In addition, the V-L orbits have revealed the largest eigenvalue moduli across all families in both equilibria which indicate a minimum time to unwind from the target orbit. In conclusion, these insights provide a fresh perspective on the range of possible solutions and phenomena that can be further explored. Moreover, new techniques have been employed to ensure the mathematical validity of libration point orbits and their hyperbolic invariant manifolds. As these structures are inherently difficult to verify, this effort forms a valuable contribution to the TU Delft Astrodynamic Toolbox (C++) for future research.