석사학위논문

Modular Design of Trajectories to Sun-Earth L4

2026

1.1Nomenclature and Acronyms/Abbreviations

Vectors with lowercase letters denote normalized versions of their uppercase vectors.

1.2.1Deep-Space Trajectory Design Software

The need for reusable trajectory design software is not unique to Sun-Earth L4. It is a recurring issue in deep-space mission design, where the analyst must explore a large design space before a small number of candidates can be refined with higher-fidelity models. Historically, one influential example was Swingby, a prototype tool developed for gravity-assist trajectory design and launch-window analysis. Carrico et al. [1] described Swingby as combining analytic approximations, multiconic orbit integration, and a Newton-method targeter to design sensitive gravity-assist trajectories quickly. NASA later reported that the Swingby lineage was commercialized as Navigator in 1994 and incorporated into STK/Astrogator in 1996 [2]; a later review of STK/Astrogator similarly traces its roots to Swingby and Navigator [3].

These examples illustrate that trajectory design software occupies more than one position in the design workflow. Tools such as GMAT, FreeFlyer, and STK/Astrogator provide propagation, targeting, optimization, scripting, and higher-fidelity mission analysis capabilities. GMAT is especially important as an open-source, multi-mission space mission analysis tool: NASA describes it as modeling, optimizing, and estimating trajectories from low-Earth orbit to interplanetary and other deep-space missions [4, 5]. FreeFlyer likewise provides high-fidelity astrodynamics functions such as orbit propagation, maneuver planning, optimization, and orbit determination for missions of any orbit regime [6]. STK/Astrogator is a commercial trajectory design and maneuver planning environment that supports development, refinement, and validation of trajectory solutions [3, 7].

However, early deep-space mission design often requires a different layer: a high-level exploration tool that can rapidly generate, filter, and compare families of feasible candidates before the design is committed to detailed numerical propagation. This role is visible in DSTE, which is described by a.i. solutions as a preliminary mission design tool for rapid generation and analysis of complex multibody trajectories, with interactive trajectory maps, Poincare maps, and a workflow for transitioning trajectories to FreeFlyer's ephemeris model [8, 9]. It is also visible in NASA's low-thrust trajectory tool suite, where MALTO is described as the medium-fidelity starting point for nearly all low-thrust preliminary design studies, while other tools such as Mystic and Copernicus serve higher-fidelity optimization roles [10]. These tools show the practical value of high-level design environments, but many such environments are commercial products, released through request-based procedures, or subject to export-control and institutional access constraints [10, 2].

The gap addressed in this thesis is therefore not the absence of propagators or numerical optimizers. Rather, it is the limited availability of source-level, reusable, high-level trajectory design frameworks that can expose the logic of candidate generation itself. The present implementation first targets inclined Sun-Earth L4 because it provides a concrete mission case with multiple competing trajectory families. At the same time, the software structure was designed with broader deep-space use in mind: components such as UniverseComponent and CelestialBody provide a path for extending the model across the spectrum of deep-space settings, including planetary systems, planet-moon systems, libration-point environments, near-Earth objects, asteroids, comets, and other small-body mission targets. Such extensions would still require the corresponding ephemeris, dynamical, and mission-specific search models, but the framework is structured so that those additions can be made by extending the system components and trajectory-arc abstractions rather than rewriting the entire design procedure. The currently empty PropagatedArc class also indicates where a future connection to a high-fidelity propagator such as GMAT could be introduced.

1.2.2Lagrange Points

Lagrange points are equilibrium points under the assumption of a restricted three-body problem. This means that when a body with small mass is placed on those points, it will maintain its position in a frame fixed to other two massive bodies (\(m_1\), \(m_2\)) that are rotating around each other. Five Lagrange points exist for a certain pair of massive bodies, of which L1, L2 and L3 are points of unstable (saddle) equilibrium, and L4 and L5 are of stable equilibrium only if \(\frac{m_1}{m_2} > 24.96\) [11].

1.2.3Sun-Earth L4 and Inclined L4

The Sun-Earth pair satisfies the mass ratio condition required for stability introduced in Section 1.2.2. Thus, the Sun-Earth L4/L5 can provide a stable base for missions of various purposes, such as observing the Sun [12, 13] or providing communication relays to spacecrafts that are obscured by the sun when viewed from Earth [14]. A key motivation for the present inclined-L4 problem is the multi-purpose heliophysics mission concept of Posner et al. [12], which identifies Sun-Earth L4 as a uniquely useful platform for space-weather observations, solar energetic particle warning, inner-heliosphere solar-wind modeling, and dust science. In particular, L4, which is 60\(^\circ\) ahead of Earth on its orbital path around the Sun, provides a continuous view of the solar radiation hemisphere centered near W60, where sources of Earth-affecting solar energetic particles are magnetically well connected to the Earth-Moon system [12, 15].

A spacecraft can not only stay precisely at L4, but also "orbit" around it. That is, it can be in a 1AU-radius, near-circular heliocentric orbit that (i) has some inclination with respect to the ecliptic plane and (ii) passes through L4 every 6 months. Posner et al. [12] proposed that L4 and L5 missions ideally use a 14.5\(^\circ\) inclination with respect to the ecliptic plane, so that the longitude coverage of L4/L5 is complemented by improved visibility of the solar polar regions. Such orbit, hereafter called "inclined L4", enables a view of the high latitude regions of the solar surface without degradation of image quality due to the effect of projection [12, 13].

1.2.4Trajectory Candidates to Inclined L4

Several different ways to reach inclined L4 have been already proposed. The simplest method is a two-chemical-propulsion-based phasing trajectory (hereafter called “Phasing”)[16]. The first burn makes the spacecraft escape from Earth's sphere of influence and inserts it into an inclined phasing orbit with a semi-major axis smaller than that of Earth. After some number of revolutions, the spacecraft will have matched its phase with L4, i.e. 60\(^\circ\) ahead of Earth. And this is when the second burn is applied, to restore the orbit into a near-circular shape and make the spacecraft stay at inclined L4.

A variation to this phasing trajectory can be done by introducing low-thrust propulsion[17]. In this case the launch vehicle delivers only a part of the total inclination required, at the stage of Earth escape and phasing orbit insertion. The rest is done by low-thrust propulsion which takes place once the spacecraft matches its phase with L4. Lesser the inclination of the phasing orbit, smaller the escape energy (C3) that must be delivered by the launch vehicle becomes, thus lowering launch expense. But at the same time the arrival mass at L4 also diminishes, as the spacecraft is now responsible of delivering bigger change in inclination.

Raising orbital inclination of the spacecraft is the most costly part throughout the transfer trajectory in terms of Δv. Given this, an alternative approach to designing the transfer is to leverage the power of gravity assist. In this approach, the spacecraft first goes through a Hohmann-like transfer orbit to Mars, and then performs a sequence of multiple gravity assists around it[18]. The posterior orbit after each gravity assist is formed such that its inclination is gradually raised towards target inclination (10-14.5\(^\circ\)), while at the same time it is in resonance with Mars to enable another rendezvous with it.

After 3-5 gravity assists, the spacecraft will have its inclination (and thus its orbital plane) matched with inclined L4. It now goes through the endgame — a sequence of (i) perihelion matching → (ii) phasing → (iii) aphelion matching to match the orbital shape and phase with inclined L4. Both chemical and low-thrust propulsion can be used in the process of phasing. This kind of strategy was named by the authors as Multiple Mars Gravity-Assist, or MMGA.

As similar to the case of low-thrust phasing trajectory, MMGA unburdens the mission from the need to deliver whole change in inclination solely using the launch vehicle. This again allows use of cheaper launch vehicles with lower C3 capabilities. It was observed that while it is impossible to achieve an inclination of 14.5\(^\circ\) by Phasing using Falcon 9 with recoverable 1st stage, it can be made possible by opting to MMGA.

Table 1.2 summarizes the trajectory-concept space that motivates this thesis. Phasing-only transfers and Mars-based MPGA have already been studied in both impulse-only and low-thrust-augmented forms, and Earth-based candidates have also been reported. However, Venus-based and mixed-planet MPGA concepts remain open in the surveyed literature. This motivates the present work's reconstruction of existing concepts and its extension to a Venus-based MPGA family.

1.2.5Modular Trajectory Design

The software gap discussed in Section 1.2.1 becomes concrete in the inclined-L4 problem. As we have seen in Section 1.2.4, there are different approaches to reach inclined L4, such as Phasing and MMGA. Under each approach there are again different options that we could choose, such as whether and when to use impulsive/low-thrust propulsion, or which resonance ratios to use during a sequence of multiple gravity assists. This leads to a generation of trajectory candidates that counts well over a few thousand.

Comparing those candidates is a painstaking task when each trajectory family is designed and described in a separate software environment, framework, or programming language. As a consequence, the process of assessing and aggregating comparable mission properties such as duration and Δv is highly time-consuming. More importantly for this thesis, the design logic itself becomes difficult to reuse: a new gravity-assist sequence, target body, or final matching condition can require another monolithic script even when most of the underlying design steps are shared.

Inconsistent design of trajectories can not only raise the efforts put into comparing those, but also threaten the credibility of comparison. Simulating a trajectory of a spacecraft is a compute-intensive task, especially when considering all the perturbations that exist in the real world. This includes N-body gravitation from other parts of the Solar System, non-spherical gravity fields and solar radiation pressure. Often, those perturbations are ignored and the motions of spacecrafts are idealized as conic arcs to accelerate the design process. Furthermore, the orbits of planets are sometimes assumed to be circular, to simplify the physical relationship between a spacecraft and a planet during gravity assist[20]. Table 1.3 shows an example of different levels of fidelity that can be employed for orbital simulation.

Trajectories designed under these simplifications later need to be fine-tuned to become an actually workable solution. During this process where the complications are now in action, the timing and vector of burns are adjusted. One might even have to introduce new mid-course burns for correction. The net effect of these in general is a rise in total Δv. Therefore, for fair comparison, it must be guaranteed that all the trajectories were designed under a same level of fidelity in orbital simulation.

This can be achieved by modularizing the design process. Modularization conisists of 3 phases (refer to Figure 1.2 for visual representation):

1. breakdown of trajectories into sequential stages
2. abstraction of stages that are common among different types of trajectories
3. reconstruction of trajectories using abstracted stages

The emphasis is on "abstraction". Common rules and concepts are identified and extracted from multiple stages from different types of trajectories. The stages that used to exist as independent algorithms are now bounded into one, and requires a mere variation of input parameters to be applied to different situations. Consider for example two different subtasks each from Phasing and MMGA, that determines the burn vectors to (i) finally insert the spacecraft to inclined L4 (in Phasing), and (ii) match the aphelion to 1AU as the last part of endgame (in MMGA). Most of the solution process is the same in the two, thus both can be abstracted as a same "deep-space maneuver". The difference in objectives, that is, whether to match the perihelion or aphelion to 1AU, can be left out as an input parameter.

Abstraction greatly reduces the number of points in the whole design process that has to be checked to guarantee consistency and absence of errors, and thus the effort put into it. And the benefits extends beyond this. Suppose one was able to improve a logic to implement some stage. Apply it to the abstracted stage, and the improvement will automatically prevail all over different types of trajectories that adopted the abstraction. Furthermore, if one wants to devise a novel trajectory, say, multiple gravity assist trajectory with respect to a new planet, what they need to do is to just assemble a new sequence of abstracted stages with different inputs.

The framework proposed in this thesis should therefore be understood primarily as a source-level and process-level contribution. It does not by itself reduce the number of dynamically feasible trajectory branches, nor does it claim to replace mission-specific pruning logic. Instead, it reduces the source-level description length needed to express related searches, localizes modifications to configuration or block-level interfaces, and provides a systematic structure in which future high-level search policies can be implemented. Once trajectory generation is represented as operations on reusable blocks and search-tree nodes, additional capabilities such as parallel branch evaluation or more aggressive pruning can be introduced at the framework level rather than being repeatedly embedded inside separate monolithic scripts. Thus, the immediate benefit is improved maintainability and reusability, while reduced wall-clock computation time is an indirect benefit that can be obtained by exploiting the modular search structure.

1.3Research Questions and Objectives of This Thesis

The preceding discussion motivates two research questions.

1. RQ1: Can a framework be built that enables efficient design of trajectories to Sun-Earth L4? The answer requires showing that the seemingly different Phasing, MMGA, and MVGA design procedures can be expressed as compositions of common building blocks. Chapter 2 and Chapter 3 address this question by defining the blocks, implementing them in a MATLAB framework, and constructing multiple trajectory families from the same abstractions.
2. RQ2: What practical utility does such a framework provide for devising and improving trajectory searches? Practical utility in this thesis is interpreted as acceleration in two design activities: devising new trajectory searches, and improving existing trajectory searches through inner-block-level modifications. Chapter 4 addresses this question by comparing code-writing and code-modification workload, by treating recursive search as a code-modification case study for deeper gravity-assist exploration, and by discussing how block-level improvements propagate to multiple trajectory types.

The first objective of this thesis is therefore to show that modularization is possible for existing inclined-L4 transfer concepts, namely Phasing and MMGA. The second objective is to demonstrate why the resulting framework is useful: the same block set is reused to construct a novel Multiple Venus Gravity-Assist (MVGA) transfer, and the same gravity-assist improvements can be applied across multiple trajectory families without rewriting the entire design procedure.

The remainder of this thesis is organized as follows. Chapter 2 defines the five trajectory building blocks, PH, IT, ED, GA, and DSM, and describes their implementation in the MATLAB framework. Chapter 3 uses these blocks to reconstruct Phasing and MMGA, to introduce the MVGA strategy, and to compare the resulting candidates through transfer time, total \(\Delta V\), launch energy, and launch-site considerations. Chapter 4 then evaluates the practical utility of the framework in terms of top-level code workload, modification locality, recursive search depth, and propagation of inner-block improvements. Finally, Chapter 5 summarizes the findings and discusses future extensions such as two-dimensional gravity-assist matching, stronger search pruning, generalized endgame optimization, mid-trajectory maneuvers, mixed-planet sequences, and applications beyond the inclined Sun-Earth L4 case.

2.1Phasing (PH)

The most critical part in designing a transfer trajectory to inclined L4 through phasing is the determination of phasing orbit itself. As L4 is located 60\(^\circ\) ahead of Earth on its orbital path, it is natural for the spacecraft to gradually lead in phase and chase L4 until it gets there. This mandates the phasing orbit to have a smaller period and therefore a smaller semi-major axis than that of Earth.

How many revolutions to stay in the phasing orbit is however still a matter of choice. Without loss of generality, we can compute the required semi-major axis of phasing orbit that makes the spacecraft arrive exactly at L4 after \(N\) revolutions around the Sun.

The spacecraft should make a phase difference of \(\Delta\phi = 60^\circ/N\) per revolution. Then the SMA of phasing orbit is obtained by using Kepler's third law as follows:

Heliocentric speed required to enter a phasing orbit with such \(a_{SC,\ ph}\) \(\left( \vecc{V_{SC,\ ph}} \right)\) can be obtained by using conservation of specific orbital energy:

If the phasing orbit were to be coplanar to that of Earth, the following relation will simply hold:

However for the inclined case of our interest, some trigonometry is involved (refer to Figure 2.3):

where

Therefore

2.2Interplanetary Transfer (IT)

In MPGA, spacecraft should first encounter with a planet to initiate a sequence of multiple gravity assists. The coasting phase between Earth departure and this first encounter is referred to as Interplanetary Transfer (IT) in this thesis. The most efficient way of doing an interplanetary transfer is a Hohmann-like transfer. It is not exactly a Hohmann transfer since the orbits of each planets are neither circular nor coplanar to each other. So instead, a transfer is referred to be "Hohmann-like" when the Earth (at departure)-Sun-planet (at encounter) angle matches 180\(^\circ\) with some margin.

When given a target planet and a time period, the IT block generates a set of "opportunities" for interplanetary transfer by doing the following:

1. Within the period, identify the pairs of Earth departure and planet encounter time that makes the Earth-Sun-planet angle near 180\(^\circ\).
2. Describe each pair in terms of orbital elements by solving Lambert's problem with the initial and final position set to that of Earth and the planet.

The resulting opportunities act as seeds for initiating the chaining of multiple gravity assists in the next step of designing MPGA trajectories, as can be seen in Figure 2.10.

When a launch vehicle is specified, a set of opportunities can further be refined by first filtering out the ones that require more escape energy \(\left( C_3 \right)\) than what the launch vehicle can deliver up to, and then raising the inclination of remaining interplanetary transfer trajectories until they match the maximum deliverable \(C_3\).

2.3Earth Departure (ED)

PH and IT block both produces the required \(\vecc{V_\infty}\) upon Earth escape. Starting from a certain initial condition on/around Earth, the Earth departure (ED) block identifies the appropriate maneuver to make an escape with the \(\vecc{V_\infty}\). In other words, input to the ED block is the outputs from either the PH or IT block. Therefore, ED is computed after the execution of PH and IT even though it is chronologically the initiating block for any types of trajectories. In fact, it is always carried out as the final step, for two reasons: First, as will be reasoned in Section 2.4, not every phasing orbit or transfer opportunity to a planet will turn out to be useful. And second, as a result of the procedure described in this section, ED is guaranteed to produce a feasible solution for any input.

Here, the spacecraft is assumed to start from a circular parking orbit around Earth. However the block can easily be extended to as well support departure from non-circular parking orbits, or even from the launchpad.

The spacecraft will make a single, impulsive burn while in the circular parking orbit to escape Earth. That is, the asymptotic velocity of hyperbolic escape path equals \(\vecc{V_\infty}\). The following explains how to obtain the optimal initial parking orbit, point of burn and the burn vector.

Speed right after burn relative to Earth is related to the speed upon Earth escape, because specific orbital energy is conserved:

The required amount of burn is obtained by subtracting the speed in parking orbit \(\left( V_{SC,\ burn^-} = \sqrt{\frac{\mu_E}{R_{SC,\ burn}}} \right)\) from above:

We now proceed to obtain where and in what direction the burn should happen.

The optimal hyperbolic path of escape after burn should satisfy the following two conditions: (i) its asymptotic velocity must be \(\vecc{V_\infty}\), and (ii) the burn should happen at the perigee of the path, for maximum efficiency. Then the possible location for such burn is not unique, but actually a set of points on a circle (refer to Figure 2.4). If there is a valid hyperbolic path that meets such conditions, produce a set of new paths by revolving the path around an the antipodal line aligned with \(\vecc{V_\infty}\). Then their asymptotic lines are all parallel with the axis, while each are apart from the axis by an identical radius. Such radius is negligible in a heliocentric scope and does not affect phasing significantly. Therefore, this set of new paths obtained is valid as a whole, and their insertion points collectively form a circle (burn circle) whose center is on the antipodal line.

Obtaining a possible point of burn (refer to Figure 2.5) starts from identifying the hyperbolic eccentricity of the escape path, which is related to the parking orbit radius and the speed upon escape[21]:

This is then related to the orientation of the apse line:

Now obtain the point of burn using \(\beta\) on a plane that is spanned by \(\vecc{v_\infty}\) and a perpendicular \(\vecc{v_{\infty\perp}}\):

The plane will become the orbital plane of both parking and escape. Note that the selection of \(\hat{z}\) is purely arbitrary. Choosing some other vector to produce \(\vecc{v_{\infty\perp}}\) will result in a point of burn that is different but still on the burn circle.

Rotate \(\vecc{R_{SC,\ burn}}\) around \(\vecc{V_\infty}\) to obtain the full set of burn circle. To rotate by the angle of \(\theta\), use quaternion as follows:

For a specific point of burn, the velocity direction after burn, which is also the direction of the burn itself, is obtained by exploiting the fact that it should be (i) contained in the orbital plane and (ii) perpendicular to \(\vecc{R_{SC,\ burn}}\) (velocity-wise burn at perigee) at the same time:

Multiply this direction to the speed to obtain the velocity itself:

As we already know the speed in parking orbit:

Furthermore, we can determine from so obtained \(\vecc{R_{SC,\ burn}}\) and \(\vecc{V_{SC,\ burn^-}}\) pairs, a continuous set of parking orbits available, as shown in yellow lines in Figure 2.4. It will be shown in Section 3.1 that while these were computed under the assumption of a two-body problem, they are still robust enough to be applied in a perturbed environment without significant modification. Given this, one can choose among the set a suitable parking orbit that satisfies certain constraints, such as a range of inclination that is deliverable from a specific launch site.

Throughout the rest of this thesis, the altitude of circular parking orbits that are put into the ED block is fixed at 300km.

2.4Gravity Assist (GA)

The primary purpose of performing gravity assists in MPGA is to raise the orbital inclination of the spacecraft. As can be seen from the name itself, MPGA might require a sequence of multiple gravity assists if a single GA is not enough to fully deliver the total desired change in inclination.

In this section, first a single GA block will be implemented, and a discussion on techniques to find a suitable sequence of GA blocks will be followed.

2.4.1Single Gravity Assist

Russel and Ocampo [22] devised a useful geometric plot to represent a planetary flyby: the \(V_\infty\) sphere, as shown in Figure 2.6. A pair of heliocentric velocities of the planet \(\left( \vecc{V_{planet}} \right)\) and the spacecraft \(\left( \vecc{V_{SC,\ in}} \right)\) on encounter uniquely determines a point on the \(V_\infty\) sphere that is centered at the tip of \(\vecc{V_{planet}}\) and has a radius of \(V_\infty\).

A planetary flyby is ideally a hyperbolic motion with respect to the planet. That is, the in and outbound \(V_\infty\) will be the same. This means that the tip of outbound asymptotic velocity vector \(\vecc{V_{\infty,\ out}}\) after flyby will still be on the same sphere.

Strange[20] proposed in his dissertation a spherical coordinate system to describe a point on the sphere. Two angles in this system: pump \(\left( \alpha \right)\) and crank \(\left( \kappa \right)\) are not identical but analogous to the azimuth and elevation in a horizontal coordinate system.

The author further derived a set of formulas to analytically compute the two angles from the orbital elements of the spacecraft, assuming the planet is in circular heliocentric motion. The flow ("forward process") is shown in the form of a directed graph in Figure 2.7. Numbers in parentheses denote the indices of individual formulae as in the dissertation. Vice versa ("backwards process", from pump and crank angles to orbital elements) is possible as well, as shown in Figure 2.8.

The forward and backward processes enable conversion between two different ways that describe the same physical relationship between a spacecraft and a planet upon encounter. Of how much can a flyby change this relationship through a gravity assist is however still in question, and has also been answered by Strange with a process depicted in Figure 2.9 ("dynamics process").

In the dynamics process, \(V_\infty\) and the periapsis distance of flyby \(\left( r_{p,\ fb} \right)\) determines the angle in which the \(\vecc{V_\infty}\) is rotated, or \(\delta\) (bending angle). Then \(\delta\) and the pump angle prior to GA \(\left( \alpha^- \right)\) decides the range of achievable pump angle after GA \(\left( \alpha^+ \right)\). Choosing a value within this range, and the corresponding crank angle after GA \(\left( \kappa^+ \right)\) will be automatically determined.

The GA block is constructed by assembling the forward, dynamics and backward process in order. This will

1. transform the orbital elements before GA in terms of pump and crank angles
2. produce a set of viable pairs of pump and crank angles after GA
3. transform those pairs back to those in terms of orbital elements.

2.4.2Sequential Chaining of Multiple Gravity Assists

Starting from the first encounter with a planet after Earth departure, ensured by the interplanetary transfer (IT) block, the spacecraft will go through a sequence of multiple gravity assists until it achieves an orbital inclination of \(i_{target}\). While chaining GA blocks, the continuity of adjacent GA blocks must be assured. That is, similar to the initial interplanetary transfer, a GA should enable another encounter with a planet for next GA. In missions such as MPGA, where the spacecraft encounters the same planet multiple times, this is equivalent to inserting the spacecraft into a resonant orbit after each GA.

A useful primitive in this process is flyby-matching: find a flyby geometry that makes a scalar post-GA quantity equal to a desired target value. In the one-dimensional version used in the baseline framework, the flyby periapsis radius is fixed to the minimum admissible value and only the rotation angle \(\theta\) of the outgoing \(\vecc{V_\infty}\) is varied. For a target function \(f\) and target value \(z_{target}\), the problem is

where \(\mathcal{T}_{GA}\) denotes the trajectory after applying the GA block with the selected geometry, and \(f_z\) denotes the residual of the scalar target condition. This primitive is used in two separate places: initial matching for resonant flybys, and terminal inclination matching immediately before the MPGA endgame. The current code implements this shared operation as Trajectory.flyby_matching(). The improved two-dimensional formulation and the sampling-based fast rejection used by the implementation are discussed later in Chapter 4; this section only describes the one-dimensional concept inherited from the prior MMGA logic.

A resonance with a given ratio corresponds to a specific semi-major axis. If the spacecraft must return to the same flyby planet after \(N\) periods of the flyby body and \(M\) spacecraft revolutions, the required post-GA semimajor axis is

where \(T_p\) is the orbital period of the flyby planet. The resonant GA block therefore calls flyby-matching with

This SMA-based matching gives an analytic initial resonant branch. A later fine-tuning stage then adjusts the flyby geometry and coast duration using ephemeris-based targeting, so that the propagated spacecraft position actually meets the next encounter body. The framework implements this procedure as Trajectory.flyby_resonant().

A sequence will terminate if it is possible to match the orbital inclination to \(i_{target}\) after last gravity assist. This time the required geometry of flyby is solved in the same manner as before but to achieve \(i_{target}\) instead of resonant semi-major axis.

It is not challenging to find a sequence that eventually matches \(i_{target}\); finding a "good" one is. The difficulty lies in the fact that the flyby geometry that results in a certain semi-major axis or inclination is not unique. The required semi-major axis by itself is not unique as there are numerous resonance ratios to choose from.

The effectiveness of a specific geometry of flyby cannot be fully evaluated until the trajectory is completed, including the solution of endgame. For example, selecting a geometry of flyby that does not raise the inclination the most, can actually enable a shorter resonance ratio on the next flyby and reduce total transfer time, or place the spacecraft into an endgame that requires lower \(\Delta V\).

This forces MPGA to be seen as a global optimization problem. In such approach, a pool of trajectory candidates is populated every time it goes through a flyby by branching to multiple different geometries of flyby. This is depicted in Figure 2.10. Only after solving the endgames for each candidates are the candidates compared in terms of total \(\Delta V\) and transfer time.

2.5Deep Space Maneuver (DSM)

Deep space maneuver in this thesis refers to an impulsive burn performed by a spacecraft when it is far from significant celestial bodies, i.e. Earth, Venus and Mars. DSM is performed at the endgame of both Phasing and MPGA. In Phasing, a single DSM at the end of phasing, i.e. at aphelion, raises the perihelion up to 1AU. This transforms the spacecraft's orbit from elliptical phasing orbit to near-circular inclined L4, thus completing the transfer. In MMGA (MVGA), three DSM blocks are carried out, once at aphelion (perihelion) and then twice at perihelion (aphelion) to match perihelion (aphelion), phase and aphelion (perihelion) with L4 one by one.

One might argue that the use of DSM blocks should not be limited only at the final stages. This is true. Lee et al.[17] designed trajectories to inclined L4 using low-thrust while phasing instead of an impulsive burn upon arrival to stop at L4. It is also possible to perform maneuvers between gravity assists. It is in fact the only way to make an Earth gravity assist right after Earth departure useful. Recall the \(V_\infty\) sphere from Section 2.4. What the spacecraft could get from performing a gravity assist around a planet was mere the rotation of its \(\vecc{V_\infty}\) vector with respect to that planet. Without a DSM, the magnitude of inbound \(\vecc{V_\infty}\) will be the same as what the spacecraft already had when it just departed from Earth. There is no point in rotating it since it could have been made already in the first place on Earth departure.

Adding maneuvers between gravity assists, or having an option to use low-thrust along with chemical propulsion essentially gives more freedom in design and an opportunity to find a better trajectory. However along with freedom comes a greatly expanded search space, and such space cannot be searched thoroughly without a firm understanding of its landscape and rules for pruning. This leads to a whole new field of study on global optimization. Therefore low-thrust propulsion, or DSMs at the middle stages of transfer is out of scope for this thesis, and thus so is the MEGA (transfer to inclined L4 through multiple Earth gravity-assist).

That leaves us only the DSM blocks during endgames as explained in the beginning of this section. In the opening of endgame for either Phasing and MPGA, the orbital plane already has an inclination of \(i_{target}\) and is therefore identical to that of inclined L4. With no need to further change the orbital plane, optimal raising or lowering of an apsis is a relatively simple task; the maneuver must be performed at the other apsis either velocity- or anti-velocity-wise, or in other words the velocity should be scaled. DSM block is implemented with this underlying concept. There are two subtypes of DSM block, which will be explained in the following subsections.

2.5.1Apsis-Matching

First is a simpler subtype among two. The spacecraft is assumed to be already at one apsis, and we know its position \(\vecc{R_{DSM}}\) and velocity \(\vecc{R_{DSM}}\). The required change in speed to match the other apsis to 1AU is to be computed.

Assuming conic motion, it can be analytically driven using the Vis-viva equation:

Or to solve for a perturbation-incorporated environment, the following optimization problem with a nonlinear constraint can be solved:

The \(\Delta V\) factor \(\gamma\) in the range specified above, allows the posterior speed to be more than 0 and less then twice the prior speed. This guarantees that the solver converge to a desired solution no matter which of the two apses we are solving for, as long as for the endgames of Phasing and MPGA.

2.5.2Phase-Matching

In the middle of the endgame for MPGA, between two apsis-matchings, is a step to match phase with L4. At the beginning of this step, the spacecraft just completes a half-revolution coast after first DSM and arrives at the opposite apsis at 1AU. Another DSM is to be applied to insert the spacecraft into phasing orbit. The basic idea is the same as in Phasing: coasting in an orbit with a period different from 1 year will make a shift in phase. However it is a more generalized problem to solve, because this time (i) the initial orbit is in an arbitrary shape with only a guarantee that one of the apses is at 1AU (the other apsis will be at a radius larger (smaller) than 1AU in case of MM(V)GA), and (ii) L4 is not assured to be 60\(^\circ\) ahead, but can be located anywhere from -180\(^\circ\) behind to 180\(^\circ\) ahead in phase.

First the difference in phase between L4 and the spacecraft must be calculated. The spacecraft is at an orbit that is no longer coplanar to the ecliptic plane, but instead has an inclination of \(i_{target}\). In such setting the phase difference cannot be obtained simply by calculating the L4-Sun-spacecraft angle. For example, suppose the spacecraft's orbit is already 1AU-near-circular and its phase also matched with L4. When the spacecraft is at the highest point over the ecliptic plane, the L4-Sun-spacecraft angle will turn out to be \(i_{target}\) even though there is no difference in phase.

The key to calculating the difference correctly is in the use of a common line between two planes, i.e. line of nodes. A cross product of angular momentum vectors of Earth and the spacecraft is in the direction of this line:

Define the absolute phase of Earth as the angle between \(\vecc{L_{LoN}}\) and Earth's position vector \(\vecc{R_E}\):

and don't forget to flip its sign if Earth is behind the ascending node:

Phase of the spacecraft is defined in the same manner. Phase of L4 is just a 60\(^\circ\) addition to that of Earth:

and we finally get the phase difference between L4 and the spacecraft:

We will make the desired phase shift by forming a phasing orbit with a semi-major axis between the value of original orbit and that of inclined L4 (1AU). That is, we will lower or raise the other apsis to match 1AU not at once but twice: by this and final DSM. Phase will be gradually matched in the intermediate period between these two.

The two DSM blocks are performed at the same location and in the same direction. In fact, they partition a fixed amount of total \(\Delta V\) that must be eventually delivered. In other words, the total \(\Delta V\) does not decrease even if we increase the number of revolutions in phasing (\(N\)). Therefore, the smallest \(N\) feasible should be chosen in order to minimize total transfer time without any trade-off.

Maximum possible phase shift per revolution is proportional to the difference in orbital periods of the original spacecraft orbit and L4:

In MMGA, \(T_{SC} > 1yr\) and therefore the phase will shift backwards \(\left( \dot{\phi} < 0 \right)\). Opposite for MVGA.

Dividing the previously obtained phase difference by the negation of above will give the minimum required number of revolutions. Apply the ceiling function to make it integral:

Now the actual shift that has to be made on each revolution:

The corresponding orbital period and semi-major axis while phasing can be obtained using \({\Delta\phi}_{per\ rev}\) in the same manner as in Section 2.1. And again by solving either the Vis-viva equation or the optimization problem introduced in Section 2.5.1 can we compute the required \(\Delta V\) for this DSM. The only modification that has to made is replacing 1AU with the SMA while phasing.

It is worth noting that there exists some corner cases where it is useful to have the semi-major axis outside the range from 1AU to the SMA of original orbit, while phasing. This includes when:

1. \(N\) drops by 1 if the phasing orbit is formed slightly bigger (smaller) than the original orbit.
2. L4 is located right ahead (behind) of the spacecraft in phase, so it can be catched up much faster by chasing (lagging) than by lagging (chasing); requires a phasing orbit with a semi-major axis smaller (bigger) than 1AU.

(for MM(V)GA)

In any of the cases, there will be a little increase in total \(\Delta V\) but a drastic save in transfer time. This transforms the endgame from a problem with a unique solution to a multi-objective optimization problem with its own Pareto front formed. Further, there is more room for optimization if we solve the endgame for MPGA as a whole including all the way from the first DSM, as a global optimization problem. These alternative strategies is however out of scope for this thesis and is left out as future work.

2.6MATLAB Framework

The framework was implemented in MATLAB using a small object-oriented core. The class diagrams in Figure 2.12, Figure 2.13, Figure 2.14 and Figure 2.15 summarize the relationships; the highlighted methods in Figure 2.14 correspond directly to the trajectory building blocks introduced in this chapter.

State is the basic representation of an orbit or instantaneous state. It stores the six Keplerian elements \((a,e,i,\Omega,\omega,\nu)\) with respect to a central body, while Cartesian quantities are exposed as dependent properties. Thus

for elliptic orbits. This keeps the stored state compact while allowing either element-wise or vector-wise calculations inside the blocks.

UniverseComponent identifies a physical or abstract object in the trajectory graph. CelestialBody extends it with mass, radius, and the derived gravitational parameter \(\mu=Gm\). Star and Planet specialize this concept; planets additionally provide ephemeris-based state_at(time) queries and admissible flyby-periapsis bounds.

Arc is the abstract unit connecting two trajectory events. A TransferArc represents finite-time motion between an origin and a target and is specialized into ConicArc and PropagatedArc. The former advances the spacecraft analytically under two-body motion; the latter is reserved for numerically propagated segments. FlybyArc represents an instantaneous planetary encounter and stores the pre- and post-flyby states, from which inbound/outbound \(\vecc{V_\infty}\), turning angle, periapsis radius, and flyby geometry are derived. ImpulseArc stores impulsive maneuvers.

Trajectory is an ordered list of arcs together with common mission constants. It exposes derived mission quantities such as start/end time, current target, total duration, \(C_3\), total \(\Delta V\), and the number of completed flybys. Its methods implement generic block operations: escape_earth() for ED, flyby_matching(), flyby_resonant(), and flyby_targeting() for GA, coast() for transfer propagation, and burn() for impulsive updates analogous to DSM. TrajectoryL4 adds L4-specific phase and circularization utilities. TrajectoryL4Phasing constructs Phasing trajectories, whereas TrajectoryL4Mpga implements MPGA start, recursive search, resonant expansion, inclination matching, and endgame completion.

The remaining helper classes support repeated block operations. FlybyGeometry stores the two flyby design variables, DtRange represents admissible transfer-time windows, Resonance describes resonant ratios, and Porkchop generates interplanetary opportunity maps. This division keeps high-level trajectory scripts short: they specify the target bodies, resonance set, search depth, and stopping criteria, while the block-level algorithms are reused through the class methods.

3.1.1Analytical Computation

A Phasing trajectory is constructed using one phasing (PH), one Earth departure (ED), and one perihelion-matching DSM block. The heliocentric part of the trajectory is first generated by PH and DSM: PH determines the inclined phasing orbit and DSM circularizes it into inclined L4 after the required number of revolutions. The ED block is then applied as a separate post-processing step to find parking-orbit departure candidates that realize the required escape state.

In the MATLAB framework, this separation appears directly at the top level, as shown in Framework Demo 3.1. The phasing construction call produces the heliocentric trajectory, whereas escape_earth() is called afterward to generate Earth-departure candidates.

constants = prepare();

% config: phasing
t_start   = juliandate(2031, 5, 3);
i_target = deg2rad(10);
N_phasing = 1;

% config: Earth departure
h_parking    = 300;
n_candidates = 12;

% phasing and endgame DSM
[trajectory, ~] = TrajectoryL4Phasing.make(t_start, N_phasing, i_target, constants);

% Earth departure
escape_trajectories = trajectory.escape_earth(h_parking, n_candidates);

% final trajectory summary and plot
trajectory.brief();

The internal structure of TrajectoryL4Phasing.make() is shown separately in Framework Demo 3.2 because the method combines two trajectory blocks. The first part solves the PH problem by choosing a velocity scaling that satisfies the phase condition after \(N\) revolutions. The final call to circularize() is the DSM step that restores the orbit to near-circular inclined L4.

function [trajectory, exitflag] = TrajectoryL4Phasing_make( ...
    t_dep, N_phasing, i_target, constants)

state_earth = constants.earth.state_at(t_dep);
    state_minus = make_inclined_circular_state_at_earth( ...
        state_earth, i_target, constants);

    dV_factor_seed = calc_dV_factor_seed(state_earth, N_phasing, constants);
    lb  = max(-0.3, dV_factor_seed - 0.1);
    ub  = min( 0.0, dV_factor_seed + 0.1);
    obj = @(dV_factor) abs(dV_factor);

    [dV_factor_opt, ~, exitflag, ~] = fmincon( ...
        obj, dV_factor_seed, [], [], [], [], lb, ub, ...
        @(dV_factor) constraints_phase_match( ...
            state_minus, t_dep, N_phasing, dV_factor, constants), ...
        constants.options_fmincon);

    state_phasing = TrajectoryL4.scaleVelocity(state_minus, dV_factor_opt);
    dt_phasing = N_phasing * state_phasing.T;

    trajectory = TrajectoryL4Phasing(constants);
    phasingArc = ConicArc(t_dep, constants.earth, state_phasing, ...
        dt_phasing, constants.void, constants.tol_dr);
    phasingArc.description = sprintf("L4 phasing (N=%d)", N_phasing);
    trajectory = trajectory.addArc(phasingArc);

[trajectory, exitflag] = trajectory.circularize();
end

The resulting object exposes common trajectory properties such as trajectory.C3, trajectory.dv, and trajectory.duration. This is the main advantage of reconstructing Phasing through the same arc and trajectory classes later used by MPGA, while keeping the Earth-departure computation as an explicit downstream block.

Table 3.1 shows the analytically driven properties of Phasing trajectories for \(N\) that ranges from \(1\) to \(10\). Total \(\Delta V\) and transfer time shows a clear relationship of trade-off. Transfer time strictly increases as \(N\) increases while total \(\Delta V\) strictly decreases.

3.1.2.1Designing the GMAT Mission

Analytically driven Phasing trajectories were further validated by reconstructing and simulating them numerically on GMAT. Among the trajectories that goes through two revolutions of phasing around the Sun (\(N=2\), \({\Delta t}_{total}=1.83yr\)), three trajectories that each starts from a parking orbit with the lowest, intermediate and the highest inclination were sampled per each \(i_{target}\), for simulation. The GMAT-internal differential corrector and solver were configured to find burn vectors such that by starting from the computed parking orbit and applying the burns, the spacecraft arrives precisely at L4. The analytically driven \(\Delta V\)s were put in as initial guess. The initial version of mission profile on GMAT is shown in Pseudoscript 3.1.

A new frame was defined to check precise arrival of the spacecraft on L4. This frame, denoted as "L4 frame" in the pseudoscript, had the Sun-Earth L4 as its origin and axes parallel to MJ2000Ec. Because of its definition of axes, \(z=0\) in L4 frame meant the intersection point with the ecliptic plane. The spacecraft was deemed to be at L4 if the magnitude of both \(x\) and \(y\) coordinates were under 100km at that moment.

The reason for a short (1 month) propagation after COI is to escape out of the ecliptic plane and nullify \(z=0\) to prevent an instant stopping during the next propagation.

The motion of the spacecraft was propagated in high-fidelity, with in consideration of non-spherical gravitation of the Earth and solar radiation pressure, as well as the gravity from other planets (regarded as point masses).

Initial mission however failed to converge. The solver struggled to find the POI vector that sent the spacecraft precisely to L4 after two revolutions. This led to the addition of midcourse correction burns: MCC1 and MCC2 were added, each after 1 and 1.5 revolutions around the Sun into phasing. The objective of POI was then modified to matching \(i_{target}\) and the SMA required for phasing. MCCs were now responsible of directing the spacecraft precisely to L4.

In addition, SK was introduced to station-keep the spacecraft around L4 after 6 months into mission orbit. The revised mission profile is shown in Pseudoscript 3.2. Underlined are the modified parts.

One month to escape earth was given so that heliocentric SMA and INC could be measured outside the influence of Earth gravity. We can also see more propagations being introduced to nullify \(z=0\) on each stage.

3.1.2.2Simulation Results

The revised mission converged successfully, and the results are shown in Table 3.2. For every sample, required \(\Delta V\)s for POI nearly equaled the theoretical minimum, with maximum overshoots (marked with daggers) of 0.79% and 0.58% for each \(i_{target}\). This was the same for COI, with overshoots all under 0.5%. Not only their magnitude conformed, but examining the individual elements of the burn vectors also revealed that they were mostly aligned to the instantaneous velocity of the spacecraft, as desired.

MCCs showed relatively small magnitudes in range of 0.11-0.16km/s. SKs showed even smaller values of about 3.2-3.3m/s. This agrees with the intended functions of the burns: correction and station-keeping. In fact, the spacecraft was so well-placed on mission orbit such that when experimenting to propagate a decade after COI without any SKs, the spacecraft deviated less than 1\(^\circ\) in phase from L4.

All of the results suggest that the net effect of various perturbations is small, and the analytical procedures implemented in PH, ED and DSM blocks with two-body approximations is valid.

3.2Multiple Planetary Gravity Assist (MPGA)

An MPGA trajectory is constructed using one interplanetary transfer block (IT), one Earth departure block (ED), one or more gravity assist blocks (GA) and three DSM blocks that each matches perihelion, phase and aphelion. MMGA introduced here is mere a dissection of the works done by Lee et al.[18] and reconstruction of them using those blocks. MVGA using Venus instead of Mars is a novel strategy introduced by the author in this thesis, but mostly follows the same sequence of blocks.

Notable changes made are due to the differences between two planets:

1. The altitude of flyby was raised than for Mars due to dense atmosphere.
2. The order of DSM blocks in endgame was reversed; this is because the orbit of Venus is within that of Earth, and the spacecraft orbit after last GA is guaranteed to have a perihelion lower than 1AU. In case of MMGA, it was guaranteed to have an aphelion distant than 1AU.

As already explained in Section 2.4.2, the process of designing an MPGA trajectory must be dealt as a global optimization problem. A single MPGA sequence can be written as a direct composition of blocks, but such code only evaluates one selected branch of the design tree. Framework Demo 3.3 illustrates this fixed-sequence view. The first arc is the IT seed, the repeated flyby calls are GA blocks, and complete() represents the endgame DSM sequence.

constants = prepare();

% config
t_earth_dep = juliandate(2031, 5, 3);
t_mars_arr = juliandate(2032, 3, 7);
dt = (t_mars_arr - t_earth_dep) * constants.day_in_seconds;

% interplanetary transfer
trajectory_0ga = Trajectory.start(t_earth_dep, constants.earth, dt, constants.mars, constants);

% multiple planetary gravity assists
trajectories_1ga = trajectory_0ga.flyby_resonant(3, 4);
trajectories_2ga = trajectories_1ga{1}.flyby_resonant(1, 1);
trajectories_3ga = trajectories_2ga{1}.flyby_resonant(4, 3);
trajectories_4ga = trajectories_3ga{1}.flyby_resonant(1, 1);
trajectories_5ga = trajectories_4ga{1}.flyby_matching(asdf);

% endgame
trajectory_final = trajectories_5ga{1}.complete();

% same Earth departure and briefing as Phasing

In practice, the framework does not commit to one resonance sequence in advance. It starts from an IT root and recursively branches through available GA options until either a terminal inclination match and DSM completion succeeds, or a pruning condition stops the branch. Figure 3.1 visualizes this search as a tree: each node represents a partial trajectory, each GA expansion creates additional children, and each terminal branch is completed by endgame DSM blocks before ED post-processing. Framework Demo 3.4 is the code-level counterpart of this branching process.

% config for interplanetary transfer and recursive search
trajectory_root = Trajectory.start();

function solutions = search_mpga(trajectory, config)
    solutions = {};

    if trajectory.duration > config.duration_max % stop: mission time limit
        return;
    end

    terminal_candidates = match_terminal_inclination(trajectory, config); % terminal flyby
    solutions = [solutions, complete_l4_candidates(terminal_candidates, config)];

    if trajectory.n_flybys >= config.n_flybys_max % stop: maximum recursive depth
        return;
    end

    children = expand_next_flyby(trajectory, config); % branch to next flyby

    for i_child = 1:numel(children)
        child_solutions = search_mpga(children{i_child}, config); % recursion
        solutions = [solutions, child_solutions];
    end
end

solutions = search_mpga(trajectory_root, config);

% same Earth departure-processing as Phasing, per individual solution

Changing this code from MMGA to MVGA primarily changes the target body, the initial transfer opportunity, and the endgame configuration, while the IT, ED, GA, and DSM implementations remain the same. This is the concrete code-level manifestation of the block reconstruction shown in Figure 2.1.

Through such branching, more than 300 and 1,200 feasible trajectories were found each for MMGA and MVGA. They were compared against Phasing trajectories using two Pareto plots: one in the transfer-time-total-\(\Delta V\) plane, and another in the transfer-time-\(C_3\) plane.

The candidates tended to be in groups of similar transfer time but different total \(\Delta V\)s. These are the ones that share the same resonance ratio sequence but differ in the merits of endgames.

Figure 3.2 also shows the Pareto front formed for each set of Phasing, MMGA and MVGA candidates. In the total-\(\Delta V\) comparison, both Mars and Venus gravity-assist trajectories generally require more total \(\Delta V\) than Phasing, although some MVGA candidates reduce transfer time relative to MMGA. This plot alone therefore does not fully expose the launch-side advantage of MPGA.

The \(C_3\) comparison gives a different view. Phasing must acquire the target inclination at Earth departure, so the launch vehicle directly carries much of the plane-change burden. In MPGA, part of this burden is transferred to later planetary flybys. As a result, some MMGA and MVGA candidates can have lower \(C_3\) than Phasing even when their total post-launch \(\Delta V\) is larger. This distinction is important for mission feasibility, because a lower \(C_3\) can permit the use of a less capable and less expensive launch vehicle while leaving the remaining energy change to the spacecraft and gravity-assist sequence.

Launch-site compatibility is another launch-side constraint that should be considered after the \(C_3\) screening. Deep-space missions generally require not only a certain escape-energy level, but also an escape asymptote in a specified inertial direction. The Earth-departure block in Section 2.3 computes a continuous set of parking-orbit and burn-point candidates for a given \(\vecc{V_\infty}\), but a ground launch site can access only the subset that is compatible with its latitude, launch azimuth restrictions, and Earth-rotation phasing. Therefore, for a fixed trajectory candidate, the practical question is often not only whether the escape direction is geometrically possible from a site, but how long and how frequently the site can provide the required launch window.

This effect is especially relevant for launch sites whose useful trajectories are close to polar. Near-polar launch angles can cover a broad range of inertial escape directions as Earth rotates, but the alignment with any particular required \(\vecc{V_\infty}\) direction can be brief. Yoo et al. [23] reported this contrast for Earth-to-Mars ballistic opportunities from Naro Space Center: as summarized in Table 3.3, Naro Space Center produced a much shorter daily window than Kennedy Space Center for the same class of transfer opportunities. The numerical values are not directly reused as constraints in the present L4 search, but they indicate why a launch-site-aware post-processing step is necessary when selecting among Pareto candidates. Among trajectories with comparable total \(\Delta V\) and \(C_3\), a candidate that admits a longer or more frequent launch window can be operationally preferable even if its energy metrics are not strictly minimal.

4.1.1Why Lines of Code Are Used Instead of Man-Hours

Man-hour is the most direct unit of engineering workload, but it is not suitable as the primary metric in this thesis. A fair man-hour comparison would require multiple researchers to independently implement the same trajectory-design tasks using both the legacy script and the proposed framework. The result would depend strongly on each researcher's familiarity with MATLAB, astrodynamics, the legacy code, and the new framework. Such an experiment is outside the scope of this work.

Instead, this thesis uses code volume as a reproducible proxy for the implementation effort associated with devising new trajectory searches. Two quantities are considered: the number of lines that must be written for a search, and the number of lines that must be modified when an existing search is extended to a related case.

4.1.2Source-Level Interpretation of Code Workload

The code-volume comparisons in this chapter are interpreted through source structure, not as a claim that line counts are a universal measure of design difficulty. Three mechanisms are relevant: reuse of common blocks, localization of variation points, and replacement of fixed-depth loop nests by recursive search.

For block reuse, let \(a_{ji}\in\{0,1\}\) indicate whether trajectory family \(j\) uses block \(i\), and let \(\ell_i\) denote the source size of the implementation of that block. In a monolithic implementation, shared logic can be duplicated across families, giving the structural contribution

where \(g_j\) denotes family-specific glue code. In a modular implementation, each used block is implemented once and reused by the families that call it:

The amount of source reduced by reusing blocks is therefore explained by

Here, \(\mathbf{1}(\cdot)\) is the indicator function and \((\cdot)_+=\max(\cdot,0)\). This equation is not used to predict the exact number of physical lines in the present codebase. Rather, it formalizes why adding another trajectory family has a smaller marginal source cost once common operations such as interplanetary transfer construction, gravity-assist matching, and deep-space maneuver targeting are represented as shared blocks.

For modification workload, the relevant concept is variation-point localization. Let \(\theta_q\) be a design choice, such as target body, resonance set, flyby-depth limit, duration limit, or endgame order. If this choice appears in \(o_q\) separate places in a monolithic script, then changing a set of choices \(Q\) requires source edits proportional to

If the same choices are exposed as configuration values or generalized block arguments, the corresponding edit cost is instead proportional to the number of changed choices:

The measured line churn and chunk counts in the following subsections can be read as empirical estimates of these two effects: fewer edited lines indicate a shorter source-level description of the change, while fewer chunks indicate stronger modification locality.

4.1.3Lines of Code That Must Be Written

Table 4.1 compares the top-level code size of the legacy MMGA script and the proposed framework search script. The legacy implementation is a monolithic script that contains opportunity selection, repeated gravity-assist loops, endgame handling, Pareto processing, and plotting in one file. The proposed framework moves the reusable logic into classes and block methods, so the top-level script mostly specifies the initial opportunity and search configuration.

This comparison should not be interpreted as saying that the full framework contains only 59 lines. Rather, it shows the amount of new top-level code needed to express a new MPGA search once the reusable blocks exist. That is the relevant workload for repeated mission-design studies.

4.1.4Lines of Code That Must Be Modified

Absolute code length is sensitive to how much logic is factored into helper functions or classes. A more meaningful test is the effort required to alter an existing MMGA search into a related search. This subsection considers two case studies: diversifying the target flyby planet, and increasing the maximum number of gravity assists. Both are common extensions when devising new MPGA searches from a previous one.

4.1.4.1Case Study 1: Diversifying the Flyby Planet

Table 4.2 compares the number of modified code lines and the number of contiguous change chunks when the flyby planet is changed. Both quantities are computed after excluding blank lines and MATLAB comment lines. For code lines, the reduction proportion is computed from total line churn, meaning additions plus deletions.

The chunk count is included because the same number of modified lines can represent different levels of editing effort. A small number of contiguous chunks indicates that the changes are localized, whereas many small chunks indicate that edits are scattered through the search logic. The measurements in Table 4.2 were obtained with utilities/count_matlab_change_metrics.sh. The reduction proportions are therefore computed as \(1 - C_{proposed}/C_{legacy}\), where \(C\) is total line churn for the code-line columns and chunk count for the chunk columns.

The planet-diversification cases mainly exercise the variation-point model in Equation 4.4 and Equation 4.5. Changing from Mars to another flyby body changes several choices at once: the target body, the initial interplanetary opportunity, admissible resonances, duration limits, and in some cases the endgame order. The endgame change can be written as a composition of primitive maneuver blocks,

where \(\pi\) denotes the order in which the primitive endgame operations are applied. A script that hard-codes each endgame variant must edit or duplicate the expanded sequence when \(\pi\) changes. A generalized implementation instead keeps the primitive blocks \(B_k\) fixed and changes the order or objective parameters passed to them. Thus, the measured reduction should be interpreted as localization of variation points, not as a change in the physical endgame problem.

The MVGA modification is larger in the legacy script because Venus is an inner planet. The endgame order must be reversed relative to the Mars case, while a Jupiter case can preserve the outer-planet apsis-matching order. In the proposed framework, this distinction is localized to the configuration and block parameters. As the set of candidate flyby bodies grows, this gap becomes more important: the legacy approach repeats planet-specific edits inside the search logic, whereas the framework reuses the same block implementations with different inputs. The proposed framework still requires scenario-specific edits such as departure and arrival epochs, target body, resonance set, duration limit, flyby-count limit, and minimum coast time, but these edits remain concentrated in the top-level configuration rather than being interleaved with the search procedure.

4.1.4.2Case Study 2: Increasing the Maximum Gravity-Assist Count

Another modification workload appears when the target trajectory requires more gravity assists than the original search was written to handle. The legacy MMGA script implements repeated gravity assists with manually nested loops. This is manageable when the target inclination can be reached with roughly five flybys, as in the 10-14.5\(^\circ\) cases studied previously. It becomes difficult to maintain when the target inclination is increased and the required number of flybys is not known in advance.

As a stress case, consider a 30\(^\circ\) inclined L4 mission. Such a high-inclination Sun-Earth L4 mission still requires stronger scientific justification, but high-inclination heliocentric missions are not unprecedented; Ulysses, for example, used gravity assist to reach a high solar-latitude orbit. In planet-satellite systems, tour missions can require even longer sequences of gravity assists. For example, the Europa Clipper 21F31 reference tour includes 69 flybys in the Jovian system [24]. This tour is not directly comparable to the present heliocentric L4 design problem, but it illustrates that high-flyby-count trajectory design is a practical maintenance requirement for reusable search software.

The key workload question is not whether a particular high-inclination search succeeds, but how much code must be edited to allow the search to continue to a deeper flyby count. Table 4.3 compares the edits required to increase the maximum flyby depth. The proposed framework changes only top-level configuration values such as the target inclination, mission-duration limit, and n_flybys_max. In contrast, the legacy script must be extended by copying the same manually nested search pattern into deeper hexa and septa branches.

This difference can be expressed more explicitly by writing the MPGA search as a recurrence. Let \(\mathcal{T}\) be the current trajectory node, let \(\mathcal{E}(\mathcal{T})\) be the set of terminal candidates obtained by inclination matching and endgame completion, and let \(\mathcal{B}(\mathcal{T})\) be the set of children produced by one additional admissible gravity-assist expansion. If \(m\) flybys remain before the depth limit is reached, the searched solution set satisfies

A manually nested implementation is a syntactic unrolling of Equation 4.8 up to a fixed depth \(D\). Extending it from depth \(D_0\) to \(D\) therefore requires adding the missing unrolled levels:

where \(\ell_k\) is the source size of the manually written loop pattern for depth \(k\). In the recursive implementation, depth is a stopping parameter rather than a source-code nesting level:

The reduction for a depth extension is therefore

The following table measures this source-level effect. It does not imply that the recursive search visits fewer trajectory nodes than the corresponding unrolled search; both can enumerate the same tree if they use the same branch generator and stopping rules. The difference is that the recursive implementation describes that tree with one transition rule and one depth parameter.

The measured rows use the same utility as Table 4.2. The legacy rows compare five-, six-, and seven-flyby versions of the monolithic MMGA script; the proposed rows use the corresponding examples/mpga/mmga, examples/mpga/mmga-6, and examples/mpga/mmga-7 search scripts.

The 69-flyby row is an extrapolation rather than an actually generated legacy script. It uses Equation 4.9 with the per-depth manual-loop cost calibrated from the observed 5-to-6 and 5-to-7 direct comparisons. For code lines, the direct legacy churn is fitted linearly as \(C_{5\to N}=80+(N-6)(161-80)\), where \(C\) is total line churn after excluding blank and comment lines. This gives \(C_{5\to69}=5183\), reported as approximately \(+5181/-2\) by keeping the two observed replacement lines and assigning the remaining churn to inserted manual-loop code. For chunks, both observed direct comparisons modify the same five regions, so the extrapolated direct 5-to-69 comparison keeps five chunks while allowing the inserted blocks to grow. The proposed framework remains at the same localized configuration edit size because recursive search depth is controlled by a parameter rather than by adding new loop nests.

4.2Improving Existing Trajectories Through Inner-Block Modifications

The second utility of the framework is that improvements to a block can propagate to all trajectory families that use that block. This is different from the workload reductions in the previous section: the point is not to measure the end-to-end time required to improve a trajectory, but to show that the framework provides a natural place to implement an inner-block modification once and reuse it throughout MMGA, MVGA, and later MPGA variants. This section discusses three examples related to GA.

4.2.1Improved One-Dimensional GA Matching

The baseline one-dimensional GA matching problem fixes the flyby periapsis radius to the minimum admissible value and solves only for the approach angle. This formulation is useful because it keeps the branching factor small, but it also makes each attempted resonant or terminal flyby depend on a scalar root solve. The framework implementation improves this block by adding a discrete sampling stage before the local solver. The same samples are used both to reject branches whose target value is outside the observed attainable range and to choose local-solver initial guesses near likely roots.

This is not a new trajectory family by itself. Its value in this chapter is that it is an inner-block improvement: once the GA matching block is changed, the same faster rejection and more accurate root refinement become available to every MPGA search that calls flyby_matching(). The detailed sampling rule, its limitations, and the benchmark values are given in Section B.3 and Section B.5.

4.2.2Two-Dimensional GA Matching

The one-dimensional formulation can return only the roots that appear along the minimum-periapsis cross section. A two-dimensional formulation adds the periapsis radius as a second flyby design variable, so a scalar post-encounter target becomes a contour in the flyby-geometry plane. This changes the role of the GA block: rather than forwarding one or two minimum-periapsis geometries, it can forward a distributed set of feasible resonant or terminal-flyby candidates to the downstream trajectory search.

This improvement is especially consistent with the purpose of the proposed framework. The downstream cost of a gravity-assist branch is not known until subsequent propagation and endgame matching are performed, so preserving multiple local GA geometries can prevent the search from prematurely committing to one branch. The appendix derives the adopted \(r_p\)-\(\theta\) representation, the contour interpretation, and the constrained-Newton contour-sampling method used for the 2DoF prototype; see Section B.1 and Section B.4. The candidate-count comparison is reported in Section B.5.

4.2.3Numerical Fine-Tuning Based on Ephemeris

Both 1DoF and 2DoF GA matching initially solve an analytic post-encounter condition, such as semimajor axis for resonance. This is only an approximation because the actual next encounter depends on ephemeris motion and the finite coast duration. The framework therefore applies a numerical fine-tuning step after the analytic match: the flyby geometry and coast duration are adjusted so that the propagated spacecraft position matches the flyby body within the configured distance tolerance.

Again, the framework-level point is propagation of an inner-block refinement. The analytical matching rule supplies a structured initial branch, and the numerical resonant-flyby routine improves that branch before it is returned to the MPGA search. As a result, the same fidelity correction can be reused by MMGA, MVGA, and later mixed-planet variants without rewriting the top-level recursive search. The residual comparison for the Mars 4:5 resonant-flyby setup is included with the other GA improvement benchmarks in Section B.5.

5.1Conclusion

This thesis addressed the challenge of consistently designing and comparing a large number of transfer trajectories to an inclined Sun-Earth L4 orbit. The key contribution is a modularization of trajectory design based on reusable building blocks. By decomposing transfer strategies into sequential stages and abstracting common stages into PH, IT, ED, GA, and DSM blocks, the proposed framework enables trajectories of different families to be generated, simulated, and assessed under a unified modeling and implementation structure.

Within this framework, the Phasing strategy was derived analytically as a two-impulse transfer composed of Earth escape into an inclined phasing orbit followed by a circularizing endgame maneuver. High-fidelity numerical reconstruction demonstrated that the analytical solutions provide accurate initial designs: the optimized Earth-departure and insertion maneuvers remained close to the theoretical minima, while additional correction and station-keeping requirements were comparatively small. These results support the practical usefulness of the block-level approximations for preliminary design and trade studies.

The same blocks were then used to reconstruct multiple planetary gravity-assist transfers as a global optimization problem. Beyond reproducing the Multiple Mars Gravity-Assist (MMGA) concept in a modular form, this thesis proposed a new Multiple Venus Gravity-Assist (MVGA) strategy by reassembling the same block set with planet-specific constraints and an adapted endgame order. Large candidate sets were generated for both MMGA and MVGA, and Pareto fronts in the transfer-time versus total Δv plane were obtained to quantify the trade-offs among Phasing, MMGA, and MVGA. Overall, the results demonstrate that modular abstraction not only reduces implementation redundancy and improves consistency, but also accelerates exploration of new trajectory concepts by simple recomposition of validated blocks.

The practical utility of the proposed framework was then evaluated from the perspective of design workload and improvement propagation. The top-level framework search required far fewer lines of newly written code than the legacy monolithic MMGA script. Changes to diversify the flyby body and to increase the maximum gravity-assist count were also localized to configuration-level choices rather than repeated manual loop and endgame edits. The gravity-assist block provided a single point where matching improvements, two-dimensional candidate generation, and ephemeris-based fine-tuning can propagate to all MPGA variants. These results answer RQ2 by showing that the proposed framework accelerates both the devising of new trajectory searches and the improvement of existing searches through reusable inner-block modifications.

5.2Future Works

Several extensions can improve both the design space coverage and the realism of the generated solutions.

1. Completion and benchmarking of two-dimensional GA matching. The two-dimensional GA concept has been formulated, but its robustness must be evaluated across flyby bodies, resonance ratios, and target quantities. Future work should define quantitative contour-coverage metrics and determine when the constrained-Newton sampling assumptions are valid.
2. Enhanced global search and pruning. Candidate explosion in multiple gravity-assist design motivates principled pruning rules and search heuristics. Incorporating dominance checks, informed branching policies, and learning-based or surrogate-assisted ranking could preserve Pareto-quality solutions while reducing computational load.
3. Generalized endgame optimization. The endgame was treated as a sequence of DSM subproblems. Solving the endgame as a coupled multi-objective optimization (time versus Δv, including corner-case phasing outside the nominal semi-major-axis bounds) may reveal additional Pareto-superior solutions.
4. Mid-trajectory maneuvers and mixed-planet sequences. Allowing deep-space maneuvers between gravity assists and enabling mixed sequences involving multiple planets can expand attainable geometries (e.g., improving inclination growth or time reduction). This requires tighter coupling between GA and DSM blocks and a broader global optimization formulation.
5. Application to other mission classes. While inclined Sun-Earth L4 is a compelling target, the proposed block-based modularization is mission-agnostic. Applying the framework to other deep-space tour problems can further validate its generality and highlight where additional abstractions are beneficial.

B.1Flyby-Geometry Representation

Once the incoming hyperbolic excess velocity \(\vecc{V_{\infty}^{-}}\) is fixed, an unpowered gravity-assist geometry has two remaining degrees of freedom. The initial GA block in this thesis follows the pump-and-crank-angle description introduced in Section 2.4.1. This representation is convenient for deriving the reachable region on the \(V_\infty\) sphere and for explaining the baseline one-dimensional matching logic inherited from prior MMGA work.

The GA improvement methods in Chapter 4, however, use a different coordinate pair: the periapsis radius \(r_p\) and an approach-angle coordinate \(\theta\). The periapsis radius controls the magnitude of the flyby turn, while \(\theta\) controls the direction in which the outgoing hyperbolic excess velocity is rotated around the incoming asymptote. This change of coordinates was introduced while developing the sampling-based matching improvements because it makes the minimum admissible flyby altitude a direct bound and makes the two-dimensional solution landscape easier to sample. For a flyby body with gravitational parameter \(\mu_B\) and incoming excess-speed magnitude \(v_\infty\), the turning angle associated with a dimensional periapsis radius \(r_p\) is

The implementation often uses a normalized periapsis-radius coordinate for convenience, but the same geometry is obtained after converting back to the corresponding dimensional periapsis radius. This radius-based distance parameter is closely related to a B-plane impact-radius parameter. The reason for using \(r_p\) directly is practical: the lower bound on admissible flyby altitude can be imposed as a simple bound,

As in Section 2.4.2, let \(\mathcal{T}_{GA}(r_p,\theta)\) denote the trajectory after applying the GA block with a selected flyby geometry. Let \(f\) be a scalar target function evaluated on that post-GA trajectory, such as the post-flyby heliocentric semimajor axis for a resonant flyby or the post-flyby inclination for a terminal flyby. The residual of the target condition is

A valid flyby-matching solution satisfies \(f_z\left(\mathcal{T}_{GA}(r_p,\theta)\right)=0\).

The landscape of \(f\left(\mathcal{T}_{GA}(r_p,\theta)\right)\) can be visualized either in B-plane coordinates or in a morphed \(r_p\)-\(\theta\) plane. In the morphed plane, the angular coordinate wraps around at \(\theta=-\pi\) and \(\theta=\pi\), and constant-\(r_p\) sections show how strongly the post-encounter quantity varies with approach direction. The variation is usually largest near \(r_{p,\min}\) because close passages produce the strongest deflection. As \(r_p\) increases, the same general angular structure remains but the flyby perturbation weakens.

This structure makes scalar post-encounter constraints useful for two classes of MPGA subproblems. In a resonant flyby, the next encounter is with the same body, so a two-body approximation reduces the resonance condition to a target post-flyby semimajor axis. In a terminal flyby, the local target may be a final orbital quantity such as inclination. In both cases, the local constraint can be evaluated immediately, while the final trajectory quality still depends on later propagation and endgame matching.

B.2Discrete Sampling for Feasibility Testing

The first use of discrete sampling is branch rejection. Since the largest variation in \(f\left(\mathcal{T}_{GA}(r_p,\theta)\right)\) is expected near the minimum admissible periapsis radius, the algorithm samples the one-dimensional cross section

With uniformly spaced samples \(\{\theta_j\}_{j=1}^{n_{test}}\), the observed attainable range is

If the requested target is outside that sampled range, the branch is rejected before running a local solver:

This rejection rule is deliberately inexpensive. It is also approximate: a coarse angular grid may miss a narrow feasible interval, and the assumption that the largest useful variation is visible near \(r_{p,\min}\) is not a universal proof for all flyby bodies, resonance ratios, and target quantities. In this thesis, the rule is used as a practical screening method for the same Mars 4:5 resonant-flyby setup used in the single-launch single-flyby example: Earth departure on 2031-05-03 and Mars arrival on 2032-03-07.

B.3Improved One-Dimensional Matching

In one-dimensional matching, the flyby periapsis radius is fixed at \(r_{p,\min}\), leaving only the approach angle \(\theta\) to be solved. The target residual is

The discrete samples used for feasibility testing also provide initial guesses for the local solver. If adjacent samples straddle the target,

then the interval \([\theta_j,\theta_{j+1}]\) contains a likely root, assuming no unresolved oscillation between the samples. The implementation uses the corresponding sampled angle or bracket-derived estimate as the initial point for root refinement. Thus, the same low-cost sampling pass supports both fast rejection of infeasible branches and more reliable convergence for feasible branches.

The conceptual difference from an interpolation-only method is that interpolation is not treated as the final answer. The samples identify where a solution is likely to exist, and a local solver then refines the flyby geometry to satisfy the target condition accurately.

B.4Two-Dimensional Contour Sampling

The two-dimensional formulation keeps both flyby variables active. Instead of intersecting a target value with a single \(r_p=r_{p,\min}\) cross section, the solution set is the level contour

Every point on this contour satisfies the same local post-encounter target. However, those points can lead to different downstream trajectories after resonant propagation, later gravity assists, and final endgame matching. The purpose of 2DoF matching is therefore not to select one locally valid flyby immediately, but to retain a distributed set of locally valid candidates for later global comparison.

The prototype uses three steps:

1. Perform the minimum-periapsis sampling test in Section B.2. If the target is outside the sampled range, abort the branch.
2. Place seed points on an ellipse in the morphed \(r_p\)-\(\theta\) plane. The seed set is chosen to span the visible flyby-geometry domain rather than starting all candidates from one constant-\(r_p\) line.
3. Tune each seed with a Newton iteration constrained to move only along the radial line from the morphed-plane origin through that seed.

The ellipse is parameterized by a seed-angle coordinate \(\eta_{seed}\). Given seed-ellipse semi-axes \(a_{seed}\) and \(b_{seed}\), the admissible range is selected so that the lower end of the ellipse remains outside the minimum-periapsis boundary:

The seed parameters are then distributed uniformly,

and mapped to morphed-plane seed coordinates by

For a seed vector \(x_{seed}=[r_{p,seed},\theta_{seed}]^T\) in the morphed plane, the constrained search varies only a scalar radial factor \(\lambda\):

Here, \(\mathcal{T}_{GA}(x)\) is used as shorthand for applying the GA block with the \(r_p\) and \(\theta\) values contained in \(x\). The Newton update is applied to the one-dimensional residual \(f_z(\mathcal{T}_{GA}(x(\lambda)))\):

The derivative can be evaluated by finite difference,

and the update is projected back to the admissible domain if it would violate the minimum-periapsis constraint. The lower bound used for each seed is

The practical motivation for the radial constraint is shown by the Mars 4:5 resonant-flyby case. A first trial placed all seeds at a constant periapsis radius and used an unconstrained Levenberg-Marquardt root solve. Although many seeds converged, the resulting samples collapsed into a narrow angular region because the steepest local gradient was primarily in the \(\theta\) direction.

A second trial used the same unconstrained local solver, but placed seeds along an ellipse so that both \(r_p\) and \(\theta\) varied in the initial population. This improved the spread of converged samples, but the solver still did not cover the full visible extent of the target contour.

The final prototype retained the elliptical seed distribution but constrained each tuning path to the radial line through its seed. This prevents the local search from immediately moving along the same steepest angular direction for many seeds. As a result, the converged samples preserve more of the initial angular diversity and cover a wider portion of the target contour.

The value of \(\Delta\eta\) is not a complete measure of contour quality: it does not measure arc-length uniformity, downstream trajectory diversity, or sensitivity to seed density. It is nevertheless useful for the present comparison because all three trials are applied to the same target landscape. In this controlled case, the constrained-Newton result spans a visibly larger part of the accessible contour than the two unconstrained alternatives.

This 2DoF sampling method should be interpreted as a prototype rather than a fully general theorem. Its behavior depends on the shape of the selected post-encounter landscape. The Mars 4:5 resonant-flyby case used here has a favorable landscape for the adopted coordinate frame, but future work should test other flyby bodies, resonance ratios, and target quantities before claiming general robustness.

B.5GA Improvement Benchmarks

The following benchmarks use the Mars 4:5 resonant-flyby setup from examples/slsf/main.m: Earth departure on 2031-05-03 and Mars arrival on 2032-03-07. They are reported in the appendix because the main claim in Chapter 4 concerns framework-level reuse of block improvements, not the standalone optimization of one GA algorithm.

C.1Constants and Solver Settings

Table C.1 summarizes the constants and default numerical settings used by the MATLAB framework. These values are defined in basics/make_constants.m.

C.2Codebase Reference

Table C.3 lists the main files corresponding to the trajectory blocks and framework components discussed in the thesis. The public MATLAB framework source code is available at https://github.com/GunheeYi/ms-thesis-src-public.