July 2014 

NASA/TM–2014-218509 

Design of the VISITOR Tool: A Versatile 

ImpulSive Interplanetary Trajectory OptimizeR 

Luca Corpaccioli and Harry Linskens 

Science and Technology Corporation, Hampton, Virginia 

David R. Komar 

Langley Research Center, Hampton, Virginia 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page 

at http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443-757-5803 

• Phone the NASA STI Information Desk at 

443-757-5802 

• Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

July 2014 

NASA/TM–2014-218509 

Design of the VISITOR Tool: A Versatile 

ImpulSive Interplanetary Trajectory OptimizeR 

Luca Corpaccioli and Harry Linskens 

Science and Technology Corporation, Hampton, Virginia 

David R. Komar 

Langley Research Center, Hampton, Virginia 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 

official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics 

and Space Administration.
1

ABSTRACT

The design of trajectories for interplanetary missions represents one of the most 

complex and important problems to solve during conceptual space mission design. To 

facilitate conceptual mission sizing activities, it is essential to obtain sufficiently 

accurate trajectories in a fast and repeatable manner. To this end, the VISITOR tool was 

developed. This tool modularly augments a patched conic MGA-1DSM model with a mass 

model, launch window analysis, and the ability to simulate more realistic arrival and 

departure operations. This was implemented in MATLAB, exploiting the built-in 

optimization tools and vector analysis routines. The chosen optimization strategy uses a 

grid search and pattern search, an iterative variable grid method. A genetic algorithm can 

be selectively used to improve search space pruning, at the cost of losing the 

repeatability of the results and increased computation time. The tool was validated 

against seven flown missions: the average total mission ΔV offset from the nominal 

trajectory was 9.1%, which was reduced to 7.3% when using the genetic algorithm at the 

cost of an increase in computation time by a factor 5.7. It was found that VISITOR was 

well-suited for the conceptual design of interplanetary trajectories, while also facilitating 

future improvements due to its modular structure.

List of Symbols

Symbol Definition Unit Subscript Definition

V∞ Excess Velocity km/s ib Inbound

α Hyperbolic Turn Angle rad ob Outbound

r Radius km tot Total

μ Gravitational Parameter km3/s2p Periapsis

ΔV Impulsive Velocity Change km/s a Apoapsis

V Velocity km/s GA Gravity Assist

i Inclination rad esc Escape

a Semi-Major Axis km hyp Hyperbolic Orbit

φ Rotation Angle rad park Parking Orbit

Ω Right Ascension Of The Ascending Node rad RGA Resonant Gravity Assist

ω Argument Of Periapsis rad helio Heliocentric

n Mean Motion rad/s sun Of The Sun

R Planet Radius km pl Of The Planet

e Eccentricity [-] inert Inert Stage

J2 J2 Spherical Harmonic Coefficient [-] prop Propellant

P Orbital Period s pld Effective Payload

j Integer Number Of Spacecraft Revolutions Around Sun [-] pay Mission Payload

k Integer Number Of Planet Revolutions Around Sun [-] 0 Initial

M Mass kg i Generic Counter

Isp Specific Impulse s j Generic Counter

ζ Tankage Mass Fraction [-] D Departure

ε Inert Mass Fraction [-] A Arrival

MJD Modified Julian Date days I Interplanetary

TOF Time-Of-Flight days grid Grid Search

Δdate Difference In Date days
2

1. Introduction

The problem of trajectory optimization is one of the central topics in space mission design. 

This is especially true when considering interplanetary flight, in which it becomes one of the 

central drivers for the mission. Additionally, this is a problem with high degrees of complexity, in 

which there is no “easy” analytical way to approximate trajectories, even at a conceptual phase. 

Said problem almost always involves the optimization of a multivariate functions with several 

constraints.

This paper focuses on the need for sizing interplanetary trajectories rather than developing 

the best optimization possible. Here, the user is not necessarily after the absolute best possible 

trajectory but seeks to compute multiple cases rapidly, valuing computational efficiency over 

fidelity. Furthermore, for a sizing tool, it is essential to be able to obtain repeatable results such 

that the user can, for example, investigate the sensitivity of certain variables. This is a rather 

stringent requirement, which is almost never fulfilled when considering all modern efforts done 

in this area [1].

Furthermore, since a sizing tool should serve many purposes, it is important to obtain a 

modular program which is easily usable and modifiable by many without too much programming 

knowledge. Such capability would allow it to interface easily with other tools, and therefore 

serve a broader purpose.

This paper describes how the VISITOR tool works, and how it answers this specific need. It 

is divided in three main parts. The first shows how the trajectory aspects are modelled; then 

moves onto the optimization strategy used, and finally ends by benchmarking the results against 

several flown missions to ensure the program produces realistic trajectories.
3

2. Trajectory model

This section describes how each routine of the program is modeled, including how solution 

methods were implemented. First, the basic multiple gravity assist (MGA) patched conic 

formulation is described, which serves as a basic structure to the other subroutines. Then, the 

various additions to such basic model, which serve to more accurately simulate reality, are 

detailed. These are the inclusion of departure and arrival orbits in planetary space, the presence 

of potential deep space maneuvers (DSM’s), the ability to perform more complex transfer orbits

(multiple revolutions, resonant orbit gravity assist), the ability to rendezvous with asteroids, and 

a mass model which allows computing and updating spacecraft mass throughout the mission as 

elements are spent or added.

2.1 MGA Model

An interplanetary trajectory is defined here by a series of Keplerian arcs connecting together 

relevant nodes around the Sun. For now, such nodes are the planets themselves; however this 

is not always the case. Such arcs are initially determined by selecting a planetary sequence and 

the times of flight between each planet. The times of flight specify the date of each planetary 

encounter; in turn, a planetary ephemeris [2] [3] is used to determine the exact coordinates of 

the planets in three-dimensional space. The location of these nodes and the time of flight 

between them fully determines Lambert’s problem [4], which is then solved using Gooding’s 

method [5]. Thus, one obtains a series of connecting Keplerian arcs which the spacecraft must 

follow. Note that multiple-revolution transfers are also possible.

2.1.1 Patched Conic and Gravity Assist Model

After having determined the necessary transfer orbits to fly a determined trajectory, the first 

step is to patch them together at each of the nodes. This is done using the patched conic 

approximation [4]. The required change in velocity to go from one orbit to another is applied 

instantaneously as a powered gravity assist (GA) at the periapsis of each hyperbolic fly-by. 

The known variables are therefore the incoming and outgoing heliocentric velocities at each 

node necessary to fly the Keplerian arcs. These in turn translate to an inbound and outbound 

excess velocity vectors V∞,ib and V∞,ob within the planet’s sphere of influence. The unknowns to 

be calculated for such maneuvers are the hyperbolic periapsis radius and the ΔV magnitude of 

the powered assist. Together, these two variables must ensure that when the spacecraft arrives 

with V∞,ib, it departs with V∞,ob. 

First, the angle between V∞,ib and V∞,ob is determined, which is the required turn angle αtot to 

be provided by the hyperbola. Thus Equation (1) can be written as shown below. 

 (1)

In Equation (1), αib is the turn angle provided by the first half of the hyperbola, going from 

V∞,ib to periapsis, and αob is the turn angle of the second half of the hyperbola, going from 

periapsis to V∞,ob. The simple relation for hyperbolic turn angle is then substituted in Equation 

(1), yielding Equation (2) for total turning angle of the powered gravity assist. 

 

(

 

 

 

 | |

 

 )

 

 

 

(

 

 

 

 | |

 

 )

 

 

(2)
4

The transcendental Equation (2) is then numerically solved for hyperbolic periapsis. Note 

that since the outgoing and incoming V∞ are potentially different in magnitude, two distinct 

hyperbola legs will be flown, with a common periapsis. Also note that, based on the planetary 

body in question, there is an absolute minimum periapsis radius that can be flown due to 

practical constraints (atmospheres and solid surfaces), hence anything below this limit makes 

the trajectory impossible to fly. 

Having determined the periapsis radius, the last step is to find the required ΔV of the 

powered assist. This is readily computed by subtracting the magnitudes of periapsis velocities of 

the two hyperbola legs as described in Equation (3).

 √ 

 

 √ 

 

 

(3)

Thus, after having selected a planetary combination, and having computed all the periapsis 

radii and powered assist magnitudes, the Keplerian arcs described in Section 2.1 can be flown 

(provided the radii are above the pre-set minima).

2.1.2 Deep Space Manuever (DSM) Model

DSMs can significantly improve a trajectory’s feasibility; in fact, several past missions have 

flown an MGA-1DSM trajectory, wherein one DSM was given along certain arcs [6] [7]. It is 

therefore important to be able to accurately simulate this. A schematic of the model used is 

shown in Figure 1.

Figure 1: Schematic of DSM model

At first, an initial guess for a time of flight is used to compute the unperturbed trajectory with 

no DSM, corresponding to the gray trajectory in Figure 1. Next, a point along this trajectory is 

chosen based on a percentage of the time of flight, known as the search point. For example, if 

the time percentage is set to 40% and time of flight is 300 days, the search point is set to the 

location on the unperturbed trajectory after 120 days, indicated in Figure 1 by the gray square.

At this point, bounds are placed around the search point, which form a “cube” around it 

(although in the planar figure they are shown as a square).
5

Now the DSM computations can begin: a point is chosen inside this search space, referred 

to here as the DSM point, and Lambert’s problem is solved, connecting it to the departing planet 

in the time allocated by the time percentage (the 120 days). Next Lambert’s problem is solved 

again, connecting the DSM point with the arrival planet in the remaining time of flight (in this 

case 180 days). These two Keplerian arcs are patched by a DSM burn to form a new trajectory,

shown in black in Figure 1. The results would then have to be further optimized by choosing 

better times of flight, DSM search point along the trajectory, and point inside the search space. 

A point of note is that DSMs are generally modeled by varying the departure velocity from a 

planet rather than varying the mid-point [8]. The method described in this paper was chosen to

allow better estimation of the exact position of the DSM through linear constraints (see Section 

4.3), which in turn can better simulate mission restrictions. For example, the spacecraft may

want to avoid certain regions (asteroid presence, proximity to the sun etc...), which is easily 

done with such a formulation. Furthermore, this implementation is seldom used, and therefore it 

can allow for future verification with more popular simulation methods.

 

2.2 Planetary Operations

One of the major additions to the MGA model described in Section 2.1 is the analysis of the 

trajectory within planetocentric space upon departure and arrival. Doing this allows the full ΔV 

requirement of the mission to be determined, which can then be used for more accurate vehicle 

sizing. In order to do this, the periapsis distance rp, apoapsis distance ra, and inclination i of the 

departing or arriving orbit must be specified by the user (also known as parking and target orbit 

respectively). Choosing this incomplete set of orbital elements allows the user to define the 

most important characteristics of the orbit for the mission at hand, while allowing simple 

analytical methods to be used to determine burn magnitudes.

A number of different types of departure and arrival procedures can be analyzed, which will 

have a different impact on the final cost function. These procedures are summarized in the 

following subsections. For all cases, it is assumed that burns to change velocity are purely 

tangential burns. Burns to change inclination will contain a component out of the orbital plane. It 

should be noted that for sections 2.2.1 and 2.2.2 only the departure from a parking orbit to a 

hyperbolic trajectory is described, as the arrival is simply the reverse of this process.

2.2.1 One-Burn Transfers

The simplest form of planetary operations is to perform a one-burn transfer. During such a 

maneuver, the spacecraft transitions directly from the parking orbit to the hyperbolic escape 

trajectory, using only a single tangential ΔV at periapsis. In applying this form of analysis, the 

orientation of the escape trajectory with respect to the target orbit is disregarded, and it is thus 

assumed that the orbits are correctly aligned. As such, one-burn transfers are only realistic 

forms of analysis if there are no requirements on the orientation of the parking or target. The ΔV 

required for the one-burn transfer is given by the vis-viva equation, shown in Equation (4) [4]. In 

general, this equation will hold for all planar velocity change maneuvers.

 √ (

 

 

) (

 

 ) √ (

 

 

) (

 

 ) (4)

2.2.2 Bi-elliptic Transfers

In order to simulate the departure from an orbit with a specified inclination around a 

specified planet, bi-elliptic transfers can be used. A bi-elliptic transfer will transition from the 
6

parking orbit to the escape trajectory or vice versa with four burns, through an intermediate orbit 

with a high apoapsis [4]. Performing the inclination change at a high altitude will dramatically 

reduce the amount of ΔV required. The burn sequence is as follows: burn at periapsis of parking 

orbit to raise apoapsis, burn at apoapsis to set departure hyperbolic periapsis, burn at equatorial 

crossing to set polar inclination, and burn at new periapsis to enter hyperbolic escape. Each of 

these will be explained in detail.

The amount of ΔV required to rotate the orbit over an arbitrary angle φ can be determined 

using Equation (5) [4].

 ( ⁄ ) (5)

The optimal location to change the inclination of the orbit is at apoapsis, as the velocity is

lowest at this point. However, since inclination changes must be performed at equatorial 

crossing, this optimal solution can only be achieved if the line of apsides lies in the equatorial 

plane [4]. In this case, the second and third burns occur simultaneously, thus becoming a three 

burn transfer.

If this is not achievable (argument of periapsis is not 0), it is desirable to approach this 

condition as much as possible by choosing the orbit orientation such the velocity is as low as 

possible at the equatorial crossing. This, in turn, implies that the crossing must occur at the 

highest possible altitude. The latter is achieved by setting the inclination of the departure 

hyperbolic orbit to 90°. Then, the inclination change burn is performed upon crossing the 

equatorial plane, which is visualized in Figure 2.

Figure 2: Schematic of suboptimal bi-elliptic inclination change

In Figure 2, a departure condition is considered in which the declination δvinf of the V∞,ob

vector and the periapsis altitude, necessitate a bi-elliptic transfer. Initially, the spacecraft is in a 

parking orbit with specified rp, ra, and i. At periapsis of the parking orbit, a burn is performed to 

raise the apoapsis, which corresponds to point A in Figure 1. This puts the spacecraft in a highly 

elliptical clockwise orbit in the same plane as the parking orbit, shown by the solid line (the 

parking orbit is not shown). Note that the line of apsides does not lie in the equatorial plane.
7

At point B, a burn is performed to change the periapsis altitude, if the parking and hyperbolic 

periapsis altitudes are different. Then, at the equatorial crossing near apoapsis, at point C, the 

third burn is performed to rotate the orbit around the line of nodes, changing the inclination to 

90°, resulting in the orbit shown by the dashed line. Finally, at point D, the departure burn is 

performed. From this inclination, any δvinf can be targeted. There is a dependence on the 

argument of periapsis, so it is assumed that the argument of periapsis of the initial orbit is set 

such that this maneuver is possible.

While this solution is suboptimal, it allows the burn magnitudes to be found using consistent

analytical relations, maintaining the fast operating speed of the tool. Upon arrival the same 

method is used, but the steps are performed in reverse.

2.2.3 Apo-Twist Transfers

When a round trip mission is simulated (an arrival followed by a departure), the orientation 

of the parking orbit around the target planet becomes important for determining accurate ΔV 

values. It is desirable to line up the inbound and outbound asymptotes during the stay, such that

V∞,ib and V∞,ob can be achieved in minimal ΔV. A method to do this is based on a single rotation 

about the line of apsides to align an inbound and an outbound escape trajectory, known as apotwist [9]. This rotation is then accomplished by a single burn at apoapsis at the correct time 

during the stay. 

This method uses B-plane targeting to align the periapsis vectors of the inbound and 

outbound hyperbola legs, for a set rp and ra of a parking orbit, such that only one tangential burn 

at arrival and one tangential burn at departure is necessary in addition to the apo-twist 

maneuver itself. The angle between the angular momentum vectors of the resulting hyperbolic 

trajectories then forms the angle to rotate at apoapsis. The ΔV required for this rotation can then 

be determined using Equation (5). Thus, by choosing the correct B-plane arrival and departure 

point, given a target parking orbit, it is possible to align the inbound and outbound excess 

velocities with three burns.

Additionally, precession due to the J2 effect is taken into account. This secularly affects the 

right ascension of the ascending node Ω and argument of periapsis ω, the behavior of which are 

governed by Equation (6) [9]. This solution methodology requires the use of numerical solvers, 

and is also not guaranteed to find a solution. However, to disregard the J2 effect in this case is 

not acceptable, and thus the effect thereof is included in the simulation. The full implementation 

of apo-twist is too complex to be fully detailed here; however it follows the method described in

Landau [9], the source of which is in the bibliography.

 ̇ 

 

 

 

 

 

 

 (  ) 

 ( )

 ̇ 

 

 

 

 

 

 

 (  ) 

( 

 )

(6)

2.2.4 Aerocapture

It is also possible to consider an aerocapture maneuver upon arrival. This allows the 

spacecraft to reach the target apoapsis using aerodynamic drag instead of a rocket burn. At 

present, no tool is integrated within VISITOR to allow for the analysis of such maneuvers. 

However, a switch exists in the code which can be used to incorporate an external aerocapture 

tool. This switch simply sets the cost function for the arrival phase to zero.
8

2.2.5 Low-Thrust Departure Asymptote Alignment

An option also exists to use low-thrust propulsion to improve the energy state of the parking 

orbit around the departure planet, and to align this with the outbound hyperbola. As with 

aerocapture and EDL, the analysis for this is independent of the variables used in the evaluation 

of the interplanetary, allowing this analysis to be performed post-optimization. Also, it has not 

been integrated into VISITOR itself, although the resulting orbit can be input as a parking orbit 

at departure. Furthermore, as for apo-twist, the full description is too long to be included here. It 

follows the methods by Gefert and Hack [10], the source of which is in the appendix.

The method used is based on analytical in-plane steering angle solutions to four subsequent 

control laws [10]. First, the semi-major axis of the initial orbit is raised through tangential 

thrusting. Second, the argument of periapsis is rotated to align the orbit properly. Then, the 

apoapsis is raised while keeping the periapsis constant. Finally, the eccentricity is maximized, 

while respecting a minimum periapsis distance. For each phase, the corresponding control law, 

as a function of true anomaly, is applied to the steering angle, and the equations of motion are 

numerically integrated. Both eclipse and precession due to J2 can be taken into account.

Since the trajectory must be numerically integrated, this analysis is rather slow. However, it 

does represent a mass-efficient method for improving the energy state upon departure, and 

therefore represents a valuable sizing tool.

2.4 Resonant Gravity Assist Model

A special class of GA trajectories occurs when the initial and target planets are the same 

and the phase angle between them is exactly 2π. At this point, Lambert’s problem does not 

have a unique solution, as there are infinitely many solutions connecting the same point in a 

given period. However, such trajectories have used on a number of previous missions (e.g. the 

Galileo [11], Rosetta [6], and MESSENGER [12] missions), so it was chosen to include the 

possibility for this special case in the overall analysis.

Two possible trajectories are considered: Gravity Assist–Resonant Gravity Assist (GARGA) 

trajectories and Departure–Resonant Gravity Assist (DRGA) trajectories. Both suffer from the 

same limitation in the solution to Lambert’s problem, but require slightly different solution 

methodologies. These will be discussed in the subsequent subsections.

2.4.1 Gravity Assist–Resonant Gravity Assist (GARGA) Trajectories

The basic problem of GARGA trajectories is illustrated in Figure 3. In essence, it is desired 

to change the heliocentric velocity from Vib to Vob. However, this might not be possible in a 

single pass, if the turn angle required to do so is too large. Thus, an intermediate resonance 

trajectory is flown, splitting the required velocity change over two GA maneuvers.

The initial heliocentric velocity results in an initial inbound excess velocity at the planet 

(V∞,ib), as seen in Figure 3. During the first GA, the excess velocity vector is changed from V∞,ib

to V∞,RGA. Since the planet is encountered at the same position in the resonant orbit at both GA 

maneuvers, the heliocentric velocity vector of the spacecraft at the beginning and at the end of 

the resonance trajectory will be the same. Thus, after the resonance orbit (an integer amount of 

planet periods later), the inbound excess velocity at the planet will also be V∞,RGA. The second 

GA maneuver then changes this velocity to V∞,ob. Note that the two periapsis altitudes rp,ib and 

rp,ob are both unknown, and can be different.
9

Figure 3: Schematic of RGA model. Left is planetocentric space, and right is heliocentric space.

The magnitude of the heliocentric velocity which must result from V∞,RGA is also known. This 

results from the fact that the resonant orbit has a specified semi-major axis helio dependent on 

the user-input resonant ratio j:k, with j being the number of revolutions performed by the 

spacecraft and k being the number of revolutions performed by the planet. This semi-major axis 

is given in Equation (7), with Ppl being the orbital period of the planet

 √ (

 

 

)

 

(7)

By extension, the excess velocity vector V∞,RGA will also be identical at both planetary 

encounters, due to the patched conic approximation. Thus, for a certain choice of V∞,RGA, it is 

possible to solve for the periapsis distance rp,ib and rp,ob of each GA maneuver. This is done in 

the same fashion as described in Section 2.1.1. A numerical method is then used to solve for 

the V∞,RGA vector, such that the total ΔV for both maneuvers is minimized.

2.4.2 Departure–Resonant Gravity Assist (DRGA) Trajectories

The basic problem for DRGA trajectories is essentially identical to that of GARGA 

trajectories: the unknown excess velocity vector V∞,RGA must be found such that the ΔV of the 

maneuvers is minimized. However, since DRGA trajectories start with a departure phase, there 

is no inbound heliocentric velocity Vib. As such, DGRA trajectories will minimize the magnitude 

of the GA burn, as a departure burn is generally provided by a separate upper stage and would 

thus provide unrealistic numbers for sizing the spacecraft itself. Naturally, although this burn is 

not kept into account in finding optimal resonant orbit, it is instead taken into account in the total 

cost function.

2.5 Rendezvous Model

In the context of this paper, rendezvous refers to the heliocentric position and velocity of the 

spacecraft matching that of a target body. While this can in principle be performed for any 

celestial body, including the major planets, it is most useful when applied to celestial bodies with 

a negligible gravity field, i.e. asteroids or comets. Examples of missions flown to comets or 

asteroids include Rosetta or NEAR Shoemaker [13].
10

The actual analysis of a rendezvous trajectory is straightforward: given the departure and 

arrival epochs, and the state of the initial planet and the target, a Lambert solution is found 

between the two. To match the velocity, a burn is performed at arrival with the object, the vector 

of which is simply given by the difference between the final velocity in the transfer orbit and the 

velocity of the target. This burn is then included in the optimization total cost function which will 

be discussed in Section’s 4 and 5 below.

Thus, the most important element in the rendezvous model is the ability to retrieve the state 

of the target body. To do this, the JPL HORIZONS ephemeris computation is used [3], together 

with the SPICE library [2], in order to read out the data. The binary ephemeris files for minor 

solar system bodies which are read by SPICE can be generated for free on-line [14].

2.6 Mass Model

An important improvement over most trajectory optimization tools is the integration of a 

mass model in VISITOR. This allows it not only to minimize for total mission ΔV, but it can also 

minimize initial mass for a certain amount of payload. The basic mass model considers a single 

stack of i high-thrust propulsion stages and j payload elements. For each stage, the specific 

impulse Isp and the inert mass ratio ε, defined as Minert/(Minert+Mprop), is specified by the user, with 

ε being either constant, or a function of the propellant mass. Per payload element, the user must 

also specify when it is added to or jettisoned from the stack, allowing elements to be on- or offloaded mid-mission. Finally, it is specified which stage performs which burns. Note that one 

stage can perform multiple burns, but the stages must be used sequentially.

The mass iteration for a single stack element is detailed in Figure 4. In this figure, a sample 

case is taken where a single stage performs three burns. The third burn is evaluated first, and 

the corresponding propellant mass Mprop,3 and inert mass Minert,3 are found with the prescribed 

payload mass Mpay,3. Then, the second burn is evaluated, using Mprop,3 summed with the new 

prescribed payload mass Mpay,2 as the total effective payload mass Mpld,2. A new inert mass 

Minert,2 is found, and if this is different from Minert,3, the found value is fed back to the iteration loop 

of the third burn. Once the two inert masses match, the analysis is stepped forward to the first 

burn. This process repeats across all stages until final convergence. If a multi-stage stack is 

considered, the resulting total mass forms the payload for the next stage.

Note that the tankage mass fraction ζ, defined as Minert/Mprop, is used instead of ε for more 

convenient notation; the conversion between the two is defined in Equation (8).

 ⁄( ) (8)

While the mass model described is quite basic, it has been designed to be easily replaced 

later on with other, more detailed mass model, once these become available. The importance is 

that it is included within the optimization routine, and thus the final solution will be mass-optimal.
11

Figure 4: Diagram of mass model iteration for a three-burn, one-stage spacecraft.
12

3. Choice of programming language

This section discusses the chosen programming language, and evaluates its usefulness 

when compared to other options. The choice was made to write VISITOR in MATLAB, which 

has implications for both computation speed and development time. First, previous choices for 

programming languages used in solving the interplanetary trajectory problem will be discussed, 

followed by arguments for the choice of using MATLAB. Finally, some examples will be given on 

its practical use.

3.1 Previous Efforts

Many interplanetary trajectory optimization tools found in the literature use the C++ or 

FORTRAN languages to model the problem. Traditionally, these languages allow the 

computations to be performed more quickly than with MATLAB [15]. However, some aspects of 

these languages may make them less attractive for rapid software development. The main 

consideration is the lack of integrated routines within the standard packages and libraries, 

necessitating that many of the tools used must be coded by the user. This implies that a larger 

skill set may be required to be able to write the tool efficiently. Also, the structuring of code 

requires a lot of user specification in order to run properly, e.g. every variable must be 

predefined as a certain type.

3.2 Advantages of MATLAB

MATLAB has a number of advantages over the commonly used C++ and FORTRAN 

languages. It is a so-called fourth-generation language, which adds an extra layer of abstraction 

above the level used in third-generation languages, of which C++ and FORTRAN are examples 

[15]. Thus, structuring of code requires less user specification, and standard packages and 

libraries contain many integrated routines. The resulting ease of programming significantly 

reduces development time and guarantees robustness, since many basic routines are 

preprogrammed and thus thoroughly verified. This leaves a lot more time for post-evaluation 

and further improvements, rather than writing the actual tool. Furthermore, these routines are 

available on almost all copies of MATLAB, whereas other C++ and FORTRAN libraries are often 

not readily available.

As previously mentioned, MATLAB contains many integrated routines. Examples of routines 

used in VISITOR and other general advantages are summarized below.

 It is an engineering tool, meaning it comes equipped with a large number of useful 

basic functions. Examples include numerical solvers (e.g. “fsolve.m”), and efficient 

integrators (e.g. “ode45.m”). 

 Its vector-based operations aid in many large computations, as well as helping to 

avoid slow loops. Furthermore, it possesses a large number of useful vector and 

matrix functions, such as the ability to reorder and reshape vectors (e.g. “min.m”, 

“reshape.m”). This also helps when dealing with frame transformations and vector 

rotations, for example, the command to convert a rotation matrix to axis-angle

formulation (e.g. “vrrotvec2mat.m”).

 There are many packages that come with it, which further enhance engineering 

capability. The most useful one considered here is the Optimization package. This 

comes with several optimizers, such as genetic algorithms and gradient-based 

methods. These are extremely complex routines which could take significant 

amounts of time to replicate (in fact, there is an entire field in engineering dedicated 
13

to this subject). The user can thus harness this entire body of knowledge with simple 

function calls (e.g. “ga.m”, “fmincon.m”). 

 Ample shareware availability. Even if MATLAB does not have a pre-made function, it 

is often the case that freeware is available to solve a particular problem (even in the 

case of astrodynamics). Examples include Lambert problem solvers, Kepler equation 

solvers and Cartesian to Keplerian converters. These do need to be verified, 

however they are often efficient, and again save programming efforts.

 MATLAB is extremely popular, and is available to anybody who purchases a license. 

This means that VISITOR could be run by a large number of engineers which 

increases its user base, and possibilities of improvement. 

 The SPICE library is easily available in MATLAB through MICE routines [2], which 

makes implementation of planetary ephemerides and asteroids extremely easy, as 

well as simplifying the process of heliocentric to planetocentric frame transformation.

 MATLAB interfaces well with other popular programs; examples include Microsoft 

Excel, Visual Basic for Applications (VBA) and Satellite ToolKit (STK). Furthermore, 

it is possible to easily convert MATLAB script to C script via the MATLAB Coder [16].

It must be noted that most capabilities summarized above could be implemented in other 

languages than MATLAB. However, the routines that are built into MATLAB have been 

developed and improved by experts in the field, and therefore can be said to be extremely 

robust and efficient. This speeds up VISITOR remarkably, as well as significantly cutting on 

development time.
14

4. Problem Formulation

In this section, the method in which the interplanetary trajectory problem was formulated will 

be detailed. The formulation described is then used in the numerical analysis of the problem, 

which will be discussed further in Section 5. First, the simplified and full problem formulations 

are explained followed by a summary of the constraints used to bound the search space.

4.1 Simplified Problem Formulation

The simplified problem formulation only considers the basic interplanetary trajectory, and is 

set up as a multidimensional static optimization problem. This formulation allows for the analysis 

of the MGA model, along with the possibility for RGA trajectories. As such, it excludes all forms 

of planetary operations and DSMs, and it is thus requires to use the excess velocities in the cost 

function, along with GA burns. It is possible to specify which planets’ excess velocity should be 

taken into account in the cost function.

A result of this formulation is a problem of which the dimension is relatively low: the control 

vector xs consists only of the initial date (in MJD) and the time-of-flight (in days) for each leg. 

Thus, a relatively complex trajectory such as the EVVEJS trajectory flown by Cassini can be 

analyzed using only six variables. The structure of xs is shown for a trajectory with i legs in 

Equation (9).

 ̅ [ ]

 

(9)

4.2 Full Problem Formulation

With the full problem, all elements detailed in Section 2 are incorporated in the problem 

formulation. Doing so results in a need to expand the control vector, increasing the 

dimensionality of the problem.

First, the time-of-flight for each leg is split up into three separate time-of-flights: one 

allocated to the planetary departure phase (TD), one allocated to the interplanetary phase (TI), 

and one allocated to the planetary arrival phase (TA). This allows more favorable geometries in 

heliocentric space to be traded off against more favorable geometries in planetocentric space. 

For certain legs, the departure or the arrival time-of-flight is set to zero, for instance during GA 

maneuvers or one-burn departures or arrivals.

Second, to account for DSMs, four variables are added to the control vector per DSM 

performed: the time at which the DSM occurs in its leg (TDSM) and the three Cartesian 

coordinates of the DSM location.

The final structure of the control vector of the full problem xfis shown in Equation (10), for 

an example with i legs and j DSMs.

 ̅ 

[ 

 ]

 

(10)

As a result of the use of the full problem formulation, a number of different cost functions 

can be selected. As with the simplified problem, it is possible to formulate the cost function in 

terms of excess velocities along with GA and DSM burns. Alternatively, planetary departure and 

arrival delta-V’s can be used instead of excess velocities. Finally, it is possible to use the burn 

magnitudes and the mass model to express the cost function in terms of initial stack mass.
15

4.3 Constraints on Problem

Any static optimization problem ( ̅) can contain particular constraints. These are 

functions of the control variables that must be satisfied. In turn, these are divided into equality 

constraints ( ̅) and inequality constraints ( ̅) . A further distinction can be made, 

wherein they can be said to be linear or non-linear, which occurs if the functions ( ̅) or ( ̅)

are either linear or not. This is an important distinction to make, as numerical optimizers have 

particular difficulty in enforcing nonlinear constraints, especially when the function is complex 

and not easily differentiable. It is crucial to the speed of a program to be able to find fast 

solutions to such constraints. 

4.3.1 Penalty functions

A popular method of dealing with nonlinear constraints is through penalty functions [17]. 

This is when the optimizer does not try to enforce the nonlinear constraint function, but rather it 

evaluates it, and determines whether the current combination of control variables is violating the 

constraint or not. If it is, the cost function is “artificially” increased by applying a penalty term. 

This will force the optimizer to shy away from that combination of state variables and seek an 

area where the constraints are no longer violated. 

Throughout the program, linear penalty functions are used to replace nonlinear constraints. 

Given the problem ( ̅), with nonlinear constraint ( ̅) , a candidate solution is used to 

evaluate ( ̅). If this is found to be outside the constraint, the optimization problem changes to 

 ( ( ̅) 

), where the constants k1 and k2 are to be selected. Usually they have to be 

large enough (in the order of 103-104) such that no solution outside the constraint could possibly 

be better than a solution inside. Thus if the cost function appears to be extremely large (<103) it 

means a nonlinear constraint is violated, and the trajectory is infeasible.

There are some disadvantages in using this method. First, the function cannot distinguish 

when a nonlinear constraint is violated, so it will evaluate more points than it needs to, and the 

user must decide whether an optimizer is within constraints or not. In addition, by changing the 

cost function upon violation of the constraints, large discontinuities are added to the problem, 

which makes the function more “jagged” and hard to differentiate. 

4.3.2 Constraints on the simplified problem

All the constraints used on the simplified problem are found in Table 1. Three different types 

of constraints are considered, each of which can fall in three separate categories. The types of 

constraint considered are bounds, linear constraints, and nonlinear constraints. Bound 

constraints are a special type of linear constraint where . 

The three categories which can be defined are ceiling, floor, or two-side. A ceiling indicates 

a situation in which the function is not allowed to be greater than the constraint. Conversely, a 

floor indicates that the function is not allowed to be smaller than the constraint. Finally, a twoside indicates that both a ceiling and a floor are used. Note that bound constraints are 

necessary to some optimizers, as these require a limited search range (see Section 5).
16

Table 1: Active constraints on simplified problem

Constraint name Explanation Type

Departure date Set constraints on desired departure date. Bound, Two-side

Times of flight Set constraints on desired time of flight for 

each leg

Bound, Two-side

Minimum 

heliocentric distance

Spacecraft must not come closer than a 

particular limit to the Sun*

Nonlinear, Floor

Minimum flyby 

periapsis

Constraint on how close spacecraft can 

perform a planet flyby 

Nonlinear, Floor

Excess velocity Bound excess velocity around any planet* Nonlinear, Two-side

Max total flight time Ensure the total flight time of all the legs 

combined is lower than a threshold*

Linear, Ceiling

* constraint can be turned off by user

4.3.3 Constraints on the full problem

The full problem imposes a few additional constraints to those used for the simplified 

problem. These additional constrains are shown in Table 2. Note that the constraints from Table 

1 are also still valid. Naturally the first two constraints in Table 2 are only used if the trajectory 

uses either DSMs or apo-twists. 

Table 2: Added active constraints on full problem

Constraint name Explanation Type

Heliocentric DSM distance Impose a minimum heliocentric 

distance that the spacecraft must keep 

during a leg with a DSM*

Nonlinear, floor

Apo-twist Alignment of periapsis vectors in apotwist. The constraint is violated if the 

difference between the periapsis 

vectors is more 50km*.

Nonlinear, ceiling

Total leg time-of-flight The sum of the departure, 

interplanetary, and arrival phases 

must match the assigned total leg 

time-of-flight.

Linear, ceiling

* constraint can be turned off by user 
17

5. Optimization Strategy

This section describes the optimization strategy used for the program. As already noted in 

Section 4.1 and 4.2, the complete trajectory analysis can be reduced to a static optimization 

problem (a simplified and a full model). This section will explain which optimizers have been 

found to yield the best results, and how they are used. It is assumed that the reader is familiar 

with typical optimization techniques; therefore this paper will not go in much detail explaining

how these work.

5.1 Grid Search

A grid search is perhaps the simplest type of optimization. It involves discretizing each 

variable of the control vector into nodes, usually regularly spaced from the lower bound to the 

upper bound. This gives the resemblance to a grid (or hyper grid, if the problem is 

multidimensional) [18]. Then, each possible combination of variable is computed, and the lowest 

cost function value is taken as best.

5.1.1 Implementation within the program

MATLAB does not possess built-in grid search routines. However, its capability to perform 

fast vector operations was used to accelerate the process. Within the program, the user 

specifies the grid resolution desired for each of the variables (departure date and times of flight). 

The program then creates N vectors (where N is the number of variables of the problem), each 

containing the regularly spaced grid points. 

An example is shown in Equation (11), where the control vector is made up of two variables: 

departure date (MJD0), and one time of flight (TOF1). Furthermore, bounds are given to each of 

the variables (note that for departure date, the MJD formulation is used).

 ̅ [ ]

 < 

 < 

(11)

Assuming the grid resolution is 100 days for the departure date, and 50 days for the time of 

flight, the grid vectors shown in Equation (12) are obtained

 [

 

 

 

] [

 

 

 

]

(12)

At this point, the program makes use of MATLAB’s in built command “combvec.m”. This tool 

combines any number of vectors and returns a matrix with every possible combination of the 

elements (referred to here as the combinations matrix). In this example, the output would be 

Equation (13).

 [

 

 ]

(13)

As can be seen, the columns provide all the combinations of control vectors to be computed 

in the grid search. The advantage of using this method is that MATLAB is able to very quickly 
18

generate this combinations matrix, which can contain hundreds of thousands of columns. This 

circumvents the problem of having to use many nested loops, which are notoriously slow to run.

The final step is simply to take each column and run it through the function, and determine 

the best point(s). This will only require one loop.

5.1.2 Notes on the grid search

Despite being a global optimization strategy, it is well known that a grid search is rather slow 

to run, since it has to go through every possible combination of variables. This also means that 

it suffers from the classic curse of dimensionality [19], where the number of combinations rises 

exponentially with the number of variables. Thus, the more Keplerian arcs a problem has, the 

longer the grid search will take without sacrificing grid resolution. It is therefore advisable to 

reduce the number of variables as much as possible when using this technique.

A useful application of a grid search is that rather than using it for finding the global optima, 

a coarse grid search can instead be used to prune the search space for favorable regions, 

which can be used as initial guesses by other optimizers. Furthermore, since there are no 

random elements in the algorithm, it will always produce equal results with the same inputs, 

thus satisfying the criterion of repeatability. 

A final advantage is that the user knows beforehand the inputs that the grid search will run 

through the cost function (these are the columns of the combinations matrix, e.g. in equation 

(13)), and therefore the user can choose to skip some, based on external criteria. For example, 

the user can impose a maximum total flight time, as in Table 1. Based on this, the user can 

eliminate some combinations before evaluating the function, which can save computation time. 

5.2 Pattern Search

Pattern search is a very powerful built-in MATLAB optimizer, called “patternsearch.m” [20]. It 

is essentially an iterative adaptive grid search which is also capable of finding global minima. It 

requires an initial guess to start, and employs the following algorithm: 

For a problem ( ̅), an initial guess ̅ and a grid size for every variable ̅ is defined.

Then, the following iteration is performed. Note that the grid size for each variable will be used 

both as a positive step and a negative step.

1. Evaluate ( ̅ 

)

2. Evaluate ( ̅ ̅)

3. If ( ̅ ̅ )< ( ̅ 

), set ̅ ̅ ̅, increase ̅

4. If ( ̅ ̅ ) ( ̅ 

) , decrease ̅

5. If ̅ , stop, minimum is found. Else, repeat from step 1.

Note that this is a simplified version of the pseudocode; more detailed information is publicly 

available [20]. Also note that it is capable of handling all types of constraints.

5.2.1 Implementation within the program

Since pattern search relies on an initial guess, is to be used at least after another preoptimizer (such as the grid search), which is to give it its best solution. 

The tolerance for which to consider the problem converged is set to relatively high (tol=10-4). 

This implies the program stops once the solutions are within an accuracy of 10-4days (~9 

seconds), or 100 km for DSM point. This is considered sufficient for a solution of the trajectory.
19

5.2.2 Notes on pattern search

Pattern search has the two major advantages needed at this stage: it is a semi-global 

optimizer, and it produces repeatable results, since there are no random elements. Being a 

standard MATLAB function it is also easily implemented, and, like discussed in Section 3: it has 

already been very efficiently programmed by experts in the field.

It has been found to be very suitable for the problems at hand. It is capable of moving 

through the search space reasonably well with low number of function evaluations, and since it 

does not contain derivative based decisions, it can easily deal with strong gradients associated 

with penalty functions. Both are important, since there are a lot of nonlinear constraints and 

many minima in the search space. 

5.3 Genetic Algorithm

This paper will only outline the general principle of genetic algorithms, since this technique is 

well known in the optimization world, and many sources exist to learn more about this [21]. It is 

a heuristic global optimization method which evolves an initial population (a number of control 

vector combinations) through several generations (iterations), each time changing the 

population according to particular rules. In the case of genetic algorithms, the rule mimics 

natural selection, wherein members of the population are able to semi-randomly “reproduce”

and yield new members of population which are combinations of the parents. Members who 

produce lower cost function values are given advantages, and others are eliminated. Other 

elements are introduced, like migration, mutation and crossover; this adds another element of 

randomness to the algorithm.

5.3.1 Implementation within the program

The MATLAB built-in genetic algorithm “ga.m” is used [22]. It is set up to take a predefined 

population from the best results of the grid search optimization. If the grid search cannot provide 

the initial population desired by the GA, the population size is adjusted to the number available.

It has been found that a high population, low generation algorithm tends to yield the best 

results (usually in the order of 200 members for 10-15 generations is sufficient). Since this 

algorithm also does not need to cycle through very many combinations, it can deal with 

relatively larger problems. 

5.3.2 Notes on genetic algorithm

It has been found that the genetic algorithm is capable of exploring more regions of the 

search space that hadn’t been previously discovered by grid search (due to poor resolution) or 

pattern search (due to poor starting condition). Thus it is capable of finding more solutions. 

However it does not yield repeatable results, which are sometimes critical. 

Since it does not require an initial guess, it can be the only algorithm used in an optimization 

strategy. However, it has been found that even a very rough grid search to get a few points 

significantly helps the runs. Similar to the pattern search, it is capable of handling penalty 

function regions very well, since the strong gradients involved are not used by the algorithm.

The GA tends to be rather slow, which is typical of heuristic algorithms. However, since it is 

a built-in routine, it is robust and contains many customizable options which the user can 

choose to change depending on the problem [22]. It has also been observed that rather than 

using it to find the absolute global optima, a secondary use is to navigate the search space 
20

better, looking for more favorable regions. Other optimizers can then take the best members of 

the population and use faster methods to converge.

5.4 Other Optimizers

It is worth noting that there are many other optimizers available. Due to the reasons 

mentioned in Section 3, only the ones pre-programmed by MATLAB have been investigated, 

and all are inferior to the optimizer’s discussed in this paper. This might not always be the case

for any optimization problem, as the performance of optimizers is often dependent on the 

problem at hand.

An unusual choice was that no gradient based methods have been used. This was a 

deliberate decision for a number of reasons. It was found that the cost function is extremely 

“jagged”, and thus contains many local minima and maxima, which in turn poses a limit for any 

gradient based method. Furthermore, penalty functions produce very large gradients (in fact, 

since the cost function “changes” when a nonlinear constraint is violated, the gradient is 

‘infinite’). This causes many problems for gradient based methods, which rely on smooth 

functions. In the end, it was found that such methods could at best provide extremely modest 

improvements, at the expense of computation time. Hence they were eliminated entirely in favor 

of non-gradient based solvers.

Due to the modular structure of the code, further developers of the program will be able to 

easily add new optimization techniques, even outside the scope of MATLAB’s built-in

algorithms.

5.5 Final Optimization Strategy

Having described the optimizers to be used, the final step is to determine the order in which 

to use them. The following order has proven to be the best strategy.

Grid search – This is the starting optimizer. It has been found that even a coarse grid tends 

to help any further optimization attempt made. It is set as the starter, as it will sample the entire 

space evenly without any preference. However, due to the “curse of dimensionality”, the grid 

search only optimizes the simplified problem, in order to achieve an acceptable grid resolution 

with shorter run times.

Genetic algorithm – This is the second optimizer in line, which the user can choose to 

deactivate. The algorithm will take the best N points from the grid search (to be determined by 

the user) as initial population (if there are fewer grid points than desired population, the 

population size is adjusted). This algorithm runs with the full problem, and takes the initial guess 

of DSM position in the middle of the unperturbed trajectory found by the grid search.

This algorithm has the advantage of being able to explore the search space better, and 

therefore it can find more local optima, however it does not yield repeatable results with equal 

inputs. Therefore, the user might want to deactivate it depending on what the program is used 

for. Generally, if one wants to perform sizing studies, repeatable results are necessary. The 

typical case is if the user wants to investigate the effect of changing a particular variable. In that 

case, the entire program needs to be run multiple times, each time with slightly different inputs. 

If the genetic algorithm is active, the user cannot tell whether the difference in results is due to 

the heuristic nature of the algorithm, or due to the different inputs. If on the other hand, the user 

wishes to find the absolute best point, it is useful to activate it.

Pattern search – This is the final optimization algorithm to be applied. It will use the best 

result from any previous optimizers as a starting point. It is rather fast in converging, and can 
21

reach ‘basin bottoms’ rather quickly. The full problem is used with this algorithm. Due to its 

proven effectiveness with this type of functions, this optimizer is always left active. Even in the

situation where it is not strictly necessary (because previous optimizers have already reached 

the global optimum) the pattern search does not slow down the program, as convergence will be 

much faster.
22

6. Results and Discussion

To validate the VISITOR tool, previously flown missions were used as test cases. The 

mission planetary sequence for each of these missions was input in VISITOR, along with ranges 

around the flown initial date and time-of-flights. The resulting optimal trajectory was then 

compared to the actual trajectory. The results of the analysis for seven flown interplanetary 

missions are shown in Section 6.1. Then, the results will be discussed in Section 6.2.

6.1 Case Analysis Results

To perform a preliminary evaluation of the performance of VISITOR, the mission sequence 

of a large number of flown missions were used. The selected missions are shown in Table 3, 

including which of the capabilities of VISITOR will be tested. Note that only missions with at 

least one GA maneuver were chosen, as direct trajectories do not allow for much variation. Also 

listed are in which leg(s) or encounter(s) these elements will occur. The bracketed numbers 

refer to asteroid and comet identification numbers.

Table 3: Summary of test case missions

Mission Planet sequence Planet 

ops

DSM RGA MultirevAsteroid

Voyager 2 [23] EJSUN - - - - -

Galileo [11] EVE(951)E(243)J J – arrival - - Flyby, flyby

Cassini [24] EVVEJS S – arrival VV - - -

NEAR 

Shoemaker 

[13]

E(253)E(433)(433) - (253)E - - Flyby, 

rendezvous

Rosetta [6] EEME(2867)E(21)(67P) - (21)(67P) EE EM Flyby, flyby, 

rendezvous

New 

Horizons [25]

EJP - - - - -

Juno [7] [26] EEJ J – arrival EE - EE -

It must be noted that strictly speaking Galileo flew an RGA trajectory between the two Earth 

GA maneuvers. However, as it performed a flyby of 951 Gaspra in this leg, the RGA capability 

of VISITOR will not be tested. The same holds for Rosetta, which visited 2867 Šteins during an 

RGA trajectory with Earth. In addition, the original trajectory for NEAR Shoemaker did not plan 

for two encounters with 433 Eros, but due to an engine failure, an extra leg was necessary to 

complete the mission [13].

The case analysis was performed in two ways: first, the exact encounter dates were input in 

VISITOR without allowing these to vary, in order to confirm that the actual trajectory could be 

reconstructed to a sufficient degree, this was done without optimizers. Second, bounds were 

increased and the optimization routines were used to find a solution with the same mission 

sequence. These two analyses will be discussed in the following subsections. Note that both 

analyses allow VISITOR to optimize for DSM position and timing, as it is not possible to fix 

these values in the user input.
23

6.1.1 Reproduction of Flown Trajectories

To reproduce the flown trajectories, the exact encounter dates were entered in VISITOR, 

and were not allowed to vary. If the mission used DSMs, the position and timing thereof were 

still allowed to be changed by the pattern search routine. No genetic algorithm was used for this 

analysis. As the computed trajectories have identical encounter dates to the flown trajectories, 

the computed excess velocity at Earth, DSM burn magnitudes (if applicable), and insertion burn 

magnitudes (if applicable) are compared to the actual values. The results of this analysis are 

shown in Table 4. Note that the computed DSM burn magnitudes include the GA burns that are 

affected by the DSM. This analysis provides benchmark values for the subsequent optimization 

analyses.

Table 4: Benchmark results for reproduction of flown trajectories

Mission V∞ at Earth [km/s] ΔV at DSM [km/s] ΔV at insertion [km/s]

Actual Computed Actual Computed Actual Computed

Voyager 2 [27] 10.119 10.128 - - - -

Galileo [28] [11] 3.61 3.933 - - 0.805[1]0.742

Cassini [24] [29] 4.254 4.043 0.466 0.451 0.613 0.587

NEAR 

Shoemaker [13]

5.089 4.896 0.269[2]0.489 0.010 0.794

0.932[2]0.540

Rosetta [6] 3.547 3.978 0.789[2]0.900 0.794 0.797

New 

Horizons [25]

12.806 12.587 - - - -

Juno [7] [26] 5.550 5.536 0.776[3]0.759 ?

[4]

0.475

[1]Value includes the effect of a GA maneuver around Io

[2]Value does not include numerous small trajectory correction maneuvers (TCM’s)

[3]Combined value of two DSM’s in quick succession

[4]

Juno has not yet arrived at Jupiter, value unknown

6.1.2 Optimization of Flown Mission Sequences

To test the optimization capabilities of VISITOR, the launch and encounter dates were given 

windows, within which these variables could be freely varied by the optimization routines. The 

bounds used consisted of the initial date ±30 days, and the window size around the time-of-flight 

of each leg of the nominal trajectory dependent on the target: ±30 days was used for the inner 

planets and most asteroids, and ±100 days was used for the outer planets and 67P. The 

progressively larger windows allow larger angular variations of bodies which have longer orbits.

For each case, the same optimization strategy will be used. A relatively quick grid search will 

generate the 200 best combinations, which will be used as the initial population for a genetic 

algorithm run. Three grid points per variable will be used in the case of missions with more than 

four legs, and five grid points per variable will be used in the other cases. The genetic algorithm 

will then run with a population of 200 individuals for 50 generations, to allow decent 

convergence while avoiding excessive run times. The best result from the genetic algorithm is 

then fed into pattern search to obtain a final solution, with a limit set of 1000 function 

evaluations. 

The results of this analysis for the seven test cases can be seen in Table 5. As with Table 4, 

the excess velocity at Earth, DSM (plus corresponding GA’s) burn magnitudes, and insertion 

burns will be compared, along with total mission ΔV. In this comparison, the computed results 

from Table 4 are treated as the “nominal” trajectory. In addition, the mean of the absolute value 
24

of the difference between encounter dates is included per mission, to provide an indication of 

how close the optimized result is to the nominal result. Also shown is the runtime of the code, 

performed on an Intel® CoreTM i5-2520M 2.50 GHz CPU with 4.00 GB of RAM. Note that 

MATLAB will only use 25% of the CPU per core at any given time. In this case, only one core is 

used, as no parallel processing capability has yet been added.

For comparison, Table 6 contains the results for identical runs, but without the genetic 

algorithm such that repeatable results are obtained. Thus, the single best result from the grid 

search will be used as a starting point for pattern search, a technique designated by GSPS. 

These results are compared to those which include the genetic algorithm, treating the optimized 

results as the new nominal trajectory. Thus, the difference between the encounter dates is also 

measured with respect to the optimized solution.

Note that all evaluations assume a 300 km altitude circular parking orbit around the Earth as 

an initial orbit, for computing total mission ΔV. This was chosen as it is a good approximation of 

actual launch practices, and since more accurate data was not available.

Table 5: Results of full optimization compared to the nominal trajectory

Mission Runtime

[s]

Mean 

|Δdate| 

[days]

Total mission ΔV 

[km/s]

V∞ at Earth [km/s] ΔV at DSM [km/s] ΔV at insertion 

[km/s]

Nominal Optimized Nominal Optimized Nominal Optimized Nominal Optimized

Voyager 2 813.4 59.2 7.255 6.691 10.128 9.406 - - - -

Galileo 682 9.3 5.279 4.682 3.933 3.866 - - 0.742 0.721

Cassini 409.2 3.2 5.124 5.346 4.043 3.890 0.451 0.883 0.587 0.591

NEAR 

Shoemaker 

375.3 24.6 6.071 7.122 4.896 4.938 0.489 1.591 0.794 0.789

0.540 0.478

Rosetta 13,442.1 22.9 6.086 5.560 3.978 3.793 0.900 0.087 0.797 0.744

New 

Horizons 

142.9 38 8.970 9.023 12.587 12.703 - - - -

Juno 123.5 2.7 5.756 5.688 5.536 5.676 0.759 0.663 0.475 0.438

Table 6: Results of the GSPS method compared to the results of the full optimization

Mission Runtime

[s]

Mean 

|Δdate| 

[days]

Total mission ΔV 

[km/s]

V∞ at Earth [km/s] ΔV at DSM [km/s] ΔV at insertion [km/s]

Optimized GSPS Optimized GSPS Optimized GSPS Optimized GSPS

Voyager 2 256.8 0.2 6.691 6.691 9.406 9.406 - - - -

Galileo 137.3 3.4 4.682 5.365 3.866 3.933 - - 0.721 0.724

Cassini 104.4 26.6 5.346 5.184 3.890 2.947 0.883 1.025 0.591 0.553

NEAR 

Shoemaker 

100.7 5 7.122 7.285 4.938 4.853 1.591 1.924 0.789 0.729

0.478 0.403

Rosetta 3,873.5 3.4 5.560 6.541 3.793 3.902 0.087 0.614 0.744 0.665

New 

Horizons 

10.4 0.7 9.023 9.035 12.703 12.710 - - - -

Juno 18.1 33 5.688 7.169 5.676 5.715 0.663 2.126 0.438 0.439

The average ΔV offset from the nominal result to the optimized and GSPS results across all 

missions was 8.2%. For the optimized results, this value is slightly lower, 7.3%, and the GSPS 

solutions have a slightly higher average offset of 9.1%. The average increase in run time when 
25

using the full optimization rather than the GSPS technique is a factor 5.7. Note that in a few 

cases, particularly Juno, the GSPS solution is markedly worse than the optimized solution. 

6.2 Discussion

From the results shown in Section 6.1, a number of differences can be noted. The first part 

of the discussion will be devoted to the reproduction of the flown trajectories, and the second 

part will relate to the optimization process.

6.2.1 Discussion on Trajectory Reproduction

Based on the results from Table 4, VISITOR is able to relatively accurately reproduce most 

trajectories, with the notable exception of NEAR Shoemaker. Particularly with respect to the 

insertion burn, the difference is very significant, and while the DSM burns have the correct order 

of magnitude, the computed values show major discrepancies with the actual burns.

These large errors can have a number of different causes. First, the trajectory flown by 

NEAR Shoemaker has numerous planned TCM’s [13], which cannot be simulated by VISITOR. 

While each of these maneuvers is small, their added effect can have a significant influence on 

the orbit. 

Second, and perhaps most significantly, NEAR Shoemaker was never intended to first fly by 

433 Eros before waiting more than a year to actually rendezvous with the asteroid [13]. 

Originally, the first encounter with 433 Eros would have been the orbit insertion. However, due 

to an engine misfire, NEAR Shoemaker was not able to perform this maneuver, and the extra 

revolution around the Sun and eventual rendezvous was implemented to salvage the mission. 

This could explain the large discrepancy between the magnitudes of the actual and computed 

second DSM burn, as this burn would have been suboptimal. However, when combined with the 

ΔV required for orbit insertion, the total expenditure of this last phase of the trajectory of the real 

mission is still less than computed by VISITOR.

To explain the large difference in insertion ΔV, it is possible be that after the second DSM, 

NEAR Shoemaker remained close enough to 433 Eros to be influenced by the latter’s gravity, 

as the DSM reduced the relative velocity to near-zero [13]. Again, this is something that 

VISITOR is unable to simulate, as it only uses the patched conic approximation.

Apart from the large discrepancies when discussing NEAR Shoemaker, the other missions 

examined show good result reproduction, matching all real values to within a few hundred m/s. 

This shows that VISITOR is able to reasonably accurately reproduce the results of flown 

missions, which validates the tool for high-level sizing purposes.

6.2.2 Case Discussion on Trajectory Optimization

A number of remarks can be made with regards to the trajectory optimization performance of 

VISITOR. First, the full optimization sequence will be considered, which includes the genetic 

algorithm pre-optimization. The most important remark to make is that, as with all heuristic 

methods, genetic algorithms incorporate an element of randomness [21], and can thus have a 

significant impact on the solution to the problem, especially one as irregular as interplanetary 

trajectory optimization. Ideally, if ran for a very long time (thousands of individuals for hundreds 

of generations), the solution should converge to the global optimum, but there is no guarantee 

of this.

With this in mind, the results from Table 5 and Table 6 can be examined. Each mission will 

be discussed individually below. A general discussion of the performance of VISITOR’s 

trajectory optimization will be done in Section 6.2.3.
26

Voyager 2 – The optimized solution for Voyager 2 shows some improvement over the 

nominal case, with respect to reduction in excess velocity and total cost function. The reason for 

this is simply due to the added windows around the time of flights, and VISITOR makes liberal 

use of this: The launch date is delayed by 12 days, with the subsequent Jupiter, Saturn, Uranus, 

and Neptune encounters occurring 24, 50, 87, and 123 days later than nominal, respectively. In 

reality, Voyager 2 was most likely subject to launch constraints, as Voyager 1 was to be 

launched 12 days after Voyager 2 in order to reach Titan [30].

As a result of the relatively simple, continuously outwards trajectory of Voyager 2, turning off 

the genetic algorithm does not change the result significantly, only delaying the launch by an 

extra day. However, it does significantly reduce run time of the code, as expected, a trend which 

will be seen throughout these cases.

Galileo – The optimized solution for Galileo also shows some improvement over the 

nominal trajectory with respect to total mission ΔV. The main difference occurs in the final 

mission leg from 243 Ida to Jupiter, which takes 53 days longer than the nominal result. For the 

other encounters, the maximum deviation from the nominal result is 3 days. The result of the 

longer final mission leg is a better flyby geometry of 243 Ida, resulting in a reduction in flyby ΔV 

by more than 200 m/s. The slight variations in the other dates also reduce flyby ΔV’s at the 

other encounters. It is possible that the final leg of Galileo’s mission was subject to additional

constraints, necessitating it to arrive at Jupiter earlier. Also, the slight variations in the other 

dates can be explained by the errors introduced by the patched conic approximation, resulting in 

a slightly different trajectory being optimal. 

However, when switching off the genetic algorithm, the solution becomes markedly worse, 

requiring more total mission ΔV than the nominal case. While the excess velocity and insertion 

ΔV are good matches for the optimized solution, the GA geometry is not as good when only 

using GSPS, resulting in more ΔV being expended at flybys. It is possible that the optimum 

solution lies in a very narrow range, requiring the better search capability of a genetic algorithm.

Cassini – The optimized solution for Cassini is slightly worse in terms of total mission ΔV 

than the nominal one, mostly resulting from a larger DSM ΔV. For the rest, the encounter dates 

are a very close match for the nominal trajectory. The fact that the DSM is suboptimal can be 

contributed to the fact that the genetic algorithm was not able to properly reproduce its 

positioning and timing.

However, the GSPS solution actually performs better than the optimized solution. While this 

seems counter-intuitive, it is most likely a result of the inherent random component of the 

genetic algorithm: it could be possible that the genetic algorithm is “stuck” in a wide local 

optimum, from which pattern search is then unable to “escape". The GSPS solution also shows 

quite a large difference in encounter dates with respect to the optimized solution: both the 

departure and arrival dates are significantly delayed, by 33 and 84 days respectively, while the 

other dates are fairly close to the optimized solution. The DSM performance is still a bit 

inaccurate.
27

NEAR Shoemaker – The optimized solution for NEAR Shoemaker is also slightly worse 

than nominal in terms of total mission ΔV. As with Cassini, most of this difference is accounted 

for by the first DSM, which expends more than 1 km/s of ΔV more than the nominal case. A 

probable cause is the large variations in encounter dates with respect to the nominal solution: 

while the optimized departure date is 15 days earlier than nominal, the final arrival date is 55 

days later than nominal.

The GSPS solution has similar performance to the optimized solution, with an average 

deviation of 5 days at each encounter. Total mission ΔV values are also very similar, implying 

that the nominal solution lies within a very narrow range of the search space.

Rosetta – The optimized solution for Rosetta is a bit better than the nominal solution, 

despite it being the most complex mission analyzed in terms of capabilities tested. The main 

source of the improvement in total mission ΔV is the reduction of the DSM burn by more than 

800 m/s, which in turn is caused by the large differences in encounter dates: The optimized 

solution departs 30 days before and arrives 87 days after the nominal solution. As with Galileo 

and Voyager, Rosetta was likely subject to similar launch or mission-related constraints which 

are not taken into account here. Also note the long run time, which is caused by the analysis of 

the DGRA trajectory flown between launch and the first flyby.

When analyzing the GSPS results, a worse result than both the nominal and optimized 

cases is observed. The encounter dates are close to the optimized trajectory, except for the 

arrival date which has been pushed back by 12 days. The result of this is that the DSM ΔV is 

once again significantly larger.

New Horizons – The optimized solution for New Horizons is very similar to the nominal 

case, with the largest difference being the arrival date, which is 109 days earlier than the 

nominal case. The GSPS solution reflects this, being almost exactly the same as the optimized 

case. While 109 days difference seems significant, the target of this leg is Pluto, which moves a 

very small amount in this length of time. Therefore, the total cost function will not change much 

over this range. 

Juno – The optimized result for Juno shows a slight improvement over the nominal case, 

caused by smaller DSM and insertion ΔV’s. The encounter dates are also very close to the 

nominal solution. However, the GSPS solution is significantly worse, primarily due to a much 

higher DSM ΔV. The likely cause is that the best point from the grid solution does not allow 

pattern search to “escape” from that local minimum, while the genetic algorithm can more easily 

cover a larger region of the search space.

6.2.3 General Discussion on Trajectory Optimization

A number of trends can be discerned based on the results obtained. First, it is clear that the 

function describing interplanetary trajectories is highly irregular, as it is very rare for two different 

methods to converge to the same minimum. Even if the trajectory is only a few days from the 

nominal result, as with Cassini, the total mission ΔV will still be quite different. While not 

unexpected, this does highlight the challenges of finding good interplanetary trajectories.

Second, in general, the use of a genetic algorithm in the optimization scheme can improve 

the total mission ΔV, although it is possible that the convergence is to a suboptimal solution, as 

with Cassini and NEAR Shoemaker. As such, there is a certain element of luck involved, 

dependent on the fineness of the initial grid search and the amount of generations and initial 

population in the genetic algorithm. However, given enough time, heuristic methods will in 

general converge in a stable manner, and for this reason such methods are often used in 
28

solving interplanetary trajectory problems [1]. A disadvantage of this is that it is difficult to obtain 

repeatability without excessively long run times.

Third, it is apparent that if the absolute global optimum is desired, the GSPS method is not 

well suited for the task. The main disadvantage lies in the fact that pattern search is only given a 

single initial point, and while it can fairly effectively examine the local search space, it is not 

capable of covering the entire search space. Thus, if the grid search returns a result which is in 

the wrong area of the search space, it is very difficult for pattern search to converge to the 

global optimum with a very irregular function. 

However, the GSPS method is repeatable: the same inputs will always yield the same 

outputs. Also, the grids used in the sample cases were relatively coarse (3 or 5 points per 

variable), so increasing the grid fineness might improve the solution. Finally, the simpler the 

trajectory, such as Voyager 2 or New Horizons, the more regular the cost function, and 

therefore pattern search is more easily able to find the “global” minimum.

Overall it can be said that the optimizers are generally able to find close trajectories to the 

nominal. Note also that computation times were purposely kept as low as possible to allow 

VISTOR to serve as a proper computation tool. For example, reproducing Cassini with a 

FORTRAN genetic algorithm to high fidelity can take several hours [31], whereas it took 

VISITOR a little under 7 minutes. Thus the goal of obtaining fast results for fast sizing purposes 

has also been achieved.

7. Conclusions

The VISITOR tool was developed to solve the impulsive interplanetary trajectory problem. 

For a given mission sequence, it uses a multiple gravity assist model based on the patched 

conic approximation, with the ability to analyze deep space maneuvers, asteroid encounters, 

resonant gravity assists, and a number of different types of planetary operations based on real 

planet orientations. The cost function can then be formulated in terms of excess velocity, burn 

ΔV, or initial stack mass for a given final payload. This cost function is then minimized using a 

hybrid optimization approach, consisting of an initial coarse grid search over a simplified 

problem formulation, an optional genetic algorithm pre-optimization over the full problem, and a 

final optimization using pattern search over the full problem.

To validate the results produced by VISITOR, seven flown interplanetary missions were 

used as test cases, being Voyager 2, Galileo, Cassini, NEAR Shoemaker, Rosetta, New 

Horizons, and Juno. First, benchmark values were produced by fixing the encounter dates of 

these missions in VISITOR and comparing the results to the actual missions. Then, two 

optimization strategies were compared: a full optimization using the genetic algorithm preoptimization, and a fast optimization which did not employ the genetic algorithm. The average 

offset between the benchmark values and the optimized results in total mission ΔV was 8.2%.

The full optimization achieved lower values than the fast optimization, 7.3% as opposed to 

9.1%, at the cost of longer run times by an average factor of 5.7. These results prove that 

VISITOR is able to find solutions with realistic values for total mission ΔV.

A number of recommendations for future VISITOR development can be made, covering the 

areas of modeling and optimization. With respect to modeling, the most important 

recommendation is the inclusion of the ability to analyze low-thrust MGA trajectories. Also, the 

inclusion of more detailed planetary operation analysis, such as the simulation of large planetary 

moons, aerocapture, EDL, ascent analysis, and spacecraft rendezvous, would result in a more 

realistic simulation. Finally, more intricate mass models than the simple iterative method 

presented here could also be implemented.

With respect to optimization, new optimization strategies should be investigated, along with 

fine-tuning the parameters of each implemented method. Through the inclusion of more 
29

advanced heuristic methods, such as evolutionary branching or monotonic basin hopping, the 

overall global optimization performance could be improved. 

As a final point, it should be mentioned that it would be more appropriate to validate the 

performance of VISITOR with cases also using an externally-developed simplified model, rather 

than flown missions which can be subject to other constraints than cannot be considered by 

VISITOR.
30

Bibliography

[

1] 

D. Izzo and T. Vinko, "ACT / ESA - Informatics and Applied Mathematics - GTOP 

Database," ESA, [Online]. Available: http://www.esa.int/gsp/ACT/inf/projects/gtop/gtop.html. 

[Accessed 16 May 2014].

[

2] 

"SPICE," NASA, 31 3 2014. [Online]. Available: http://naif.jpl.nasa.gov/naif/. [Accessed 

5 5 2014].

[

3] 

"Solar System Dynamics," NASA - Jet Propulsion Laboratory, [Online]. Available: 

http://ssd.jpl.nasa.gov/?ephemerides. [Accessed 5 5 2014].

[

4] 

Vallado, Fundamentals of Astrodynamics and Applications, El Segundo, CA: Microcosm 

Press & Kluwer Academic Publishers, 1999. 

[

5] 

R. H. Gooding, "A procedure for the solution of Lambert's orbital boundary-value 

problem," Celestial Mechanics and Dynamical Astronomy, p. 48:145, 1990. 

[

6] 

J. Schoenmaekers and R. Bauske, "Re-Design of the Rosetta Mission For Launch in 

2004," in Proceedings of the 18th International Symposium on Space Flight Dynamics, 

Munich, Germany, 2004. 

[

7] 

T. D. Kowalkowski, J. R. Johannesen and T. Lam, "Launch Period Development for the 

Juno Mission to Jupiter," in AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 

Honolulu, HI, 2008. 

[

8] 

J. Englander, "Automated Trajectory Planning for Multiple-flyby Interplanetary 

Missions," University of Illinois, Urbana, 2013.

[

9] 

D. F. Landau, J. M. Longuski and P. A. Penzo, "Method for Parking-Orbit Reorientation 

for Human Missions to Mars," Journal of Spacecraft and Rockets, vol. 42, no. 3, pp. 517-

522, 2005. 

[

10] 

L. P. Gefert and K. J. Hack, "Low-Thrust Control Law Development for Transfer from 

Low Earth Orbits to High Energy Elliptical Parking Orbits," in AAS/AIAA Astrodynamics 

Specialist Conference, Girdwood, AK, 1999. 

[

11] 

L. A. D'Amario, L. E. Bright and A. A. Wolf, "Galileo Trajectory Design," Space Science 

Reviews, vol. 60, no. 1-4, pp. 23-78, 1992. 

[

12] 

J. V. McAdams, D. W. Dunham, R. W. Farquhar, A. H. Taylor and B. G. Williams, 

"Trajectory Design and Maneuver Strategy for the MESSENGER Mission to Mercury," 

Journal of Spacecraft and Rockets, vol. 43, no. 5, pp. 1054-1064, 2006. 

[

13] 

D. M. Dunham, J. V. McAdams and R. W. Farquhar, "NEAR Mission Design," John 

Hopkins APL Technical Digest, vol. 23, no. 1, pp. 18-33, 2002. 

[

14] 

JPL SSD Group, "Asteroid & Comet SPK File Generation Request," Jet Propulsion 

Laboratory, [Online]. Available: http://ssd.jpl.nasa.gov/x/spk.html. [Accessed 18 March 

2014].

[

15] 

D. F. Rico, H. H. Sayani and R. F. Field, "History of Computers, Electronic Commerce 

and Agile Methods," in Advances in Computers, Elsevier Science Publishing Co Inc, 2014, 

pp. Volume 73, pages 1-55.

[

16] 

MathWorks, "Quickly convert MATLAB code to C with MATLAB Coder.," MathWorks, 

[Online]. Available: http://www.mathworks.com/programs/nrd/matlab-convert-to-c-trialrequest.html?ref=ggl&s_eid=ppc_2790. [Accessed 15 05 2014].

[

17] 

A. E. Smith and D. W. Coit, "Penalty Functions," Oxford University Press and Institute of 

Physics Publishing, Oxford, 1996.
31

[

18] 

J. Pintér, "Global Optimization," MathWorld--A Wolfram Web Resource, [Online]. 

Available: http://mathworld.wolfram.com/GlobalOptimization.html. [Accessed 13 05 2014].

[

19] 

R. E. Bellman, Dynamic programming, Princeton: Princeton University Press, 1957. 

[

20] 

MathWorks, "Direct Search," MathWorks, [Online]. Available: 

http://www.mathworks.com/help/gads/direct-search.html. [Accessed 14 05 2014].

[

21] 

Y. Davidor, H. Schwefel and R. Männer, Genetic algorithms with multi-parent 

recombination, Springer, 1994. 

[

22] 

MathWorks, "Genetic Algorithm," [Online]. Available: 

http://www.mathworks.com/help/gads/genetic-algorithm.html. [Accessed 14 05 2014].

[

23] 

R. P. Rudd, J. C. Hall and G. L. Spradlin, "The Voyager Interstellar Mission," Acta 

Astronautica, vol. 40, no. 2-8, pp. 383-396, 1997. 

[

24] 

F. Peralta and S. Flanagan, "Cassini Interplanetary Trajectory Design," Control 

Engineering Practice, vol. 3, no. 11, pp. 1603-1610, 1995. 

[

25] 

Y. Guo and R. W. Farquhar, "New Horizons Mission Design," Space Science Reviews, 

vol. 140, no. 1-4, pp. 49-74, 2008. 

[

26] 

T. Greicius and B. Dunbar, "NASA - Juno Spacecraft Cruise Trajectory," NASA, 28 July 

2011. [Online]. Available: 

http://www.nasa.gov/mission_pages/juno/multimedia/Juno20110727-3.html. [Accessed 14 

May 2014].

[

27]

R. L. Heacock, "The Voyager Spacecraft," in James Watt International Gold Medal 

Lecture, London, UK, 1980. 

[

28] 

H. F. Meissinger, S. Dawson and J. R. Wertz, "A Low-Cost Modified Launch Mode for 

High-C~ 3 Interplanetary Missions," Advances in the Astronautical Sciences, vol. 97, pp. 

1693-1708, 1998. 

[

29] 

A. Wessen, K. Munsell, E. Piazza and A. Yu, "Cassini Solstice Mission: Saturn Arrival --

A Guide to Saturn Orbit Insertion," NASA, [Online]. Available: 

http://saturn.jpl.nasa.gov/photos/imagedetails/index.cfm?imageId=794. [Accessed 15 May 

2014].

[

30] 

C. E. Kohlhase and P. A. Penzo, "Voyager Mission Description," Space Science 

Reviews, vol. 21, no. 2, pp. 77-101, 1977. 

[

31] 

B. K. B. W. Sam Wagner, "GPU Accelerated Genetic Algorithm for Multiple GravityAssist and Impulsive dV Maneuvers," in AIAA/AAS Astrodynamics Specialist Conference, 

Minneapolis, Minnesota, 2012. 

[

32] 

D. Eagle, "Modeling Aero-assist Flight Mechanics with MATLAB - File Exchange -

MATLAB Central," The MathWorks, Inc., 19 August 2013. [Online]. Available: 

http://www.mathworks.com/matlabcentral/fileexchange/43121-modeling-aero-assist-flightmechanics-with-matlab. [Accessed 19 May 2014].
32

Appendix

Here the full results from the validation cases are shown. This is the raw data obtained from 

VISITOR. For each of the seven missions analyzed, three tables will be shown: the first shows 

the nominal case where the dates were hard coded into the program, then the case where the 

genetic algorithm was not used, and lastly the case where the genetic algorithm is used.

NEAR-Showmaker:
33

Juno:

Rosetta

:
34

New Horizons:

Voyager 2:
35

Galileo:

Cassini:
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

Design of the VISITOR Tool: A Versatile ImpulSive Interplanetary 

Trajectory OptimizeR 

5a. CONTRACT NUMBER

 6. AUTHOR(S)

 Corpaccioli, Luca; Linskens, Harry; Komar, David R.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20442

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 13

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

The design of trajectories for interplanetary missions represents one of the most complex and important problems to solve during conceptual 

space mission design. To facilitate conceptual mission sizing activities, it is essential to obtain sufficiently accurate trajectories in a fast and 

repeatable manner. To this end, the VISITOR tool was developed. This tool modularly augments a patched conic MGA-1DSM model with a 

mass model, launch window analysis, and the ability to simulate more realistic arrival and departure operations. This was implemented in 

MATLAB, exploiting the built-in optimization tools and vector analysis routines. The chosen optimization strategy uses a grid search and 

pattern search, an iterative variable grid method. A genetic algorithm can be selectively used to improve search space pruning, at the cost of 

losing the repeatability of the results and increased computation time. The tool was validated against seven flown missions: the average total 

mission ¿V offset from the nominal trajectory was 9.1%, which was reduced to 7.3% when using the genetic algorithm at the cost of an 

increase in computation time by a factor 5.7. 

15. SUBJECT TERMS

Astrodynamics; Deep space manuever; Gravity assist; Interplanetary; Lambert; Mission design; Optimization

18. NUMBER

 OF 

 PAGES

19b. TELEPHONE NUMBER (Include area code)

40 (443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

662122.04.07.16 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2014-218509

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 07 - 2014