Our project focused on finding the optimal path a space probe could take while traveling in the inner solar system. The model we created accounts for the gravity that large celestial objects produce, as well as their orbits. The inspiration for our problem comes from the Nasa Missions Voyager and Mariner. The Voyager missions launched in 1977 (“Voyager - Planetary,” 2015) and the Mariner Missions occurred throughout the 1960s and early 1970s ("Mariner Missions," n.d.). These missions involved space probes traveling throughout the inner and outer Solar System, visiting multiple planets and studying each with various onboard instruments. Here's a link to some facts about the Mariner missions, and another link about the Voyager missions.
For our project we will be finding an optimized route for a probe to follow starting at Earth and subsequently visiting Mars, Venus, and finally Mercury. Our model differs from the Mariner and Voyager missions in that we are not basing our assumptions on the mechanics and limitations of a traditional chemical rocket. Rather, the probe will use an ion thruster for propulsion, thus limiting its max thrust output to less than 0.25 newtons but allowing for a much greater amount of burn time. In addition, the payload of the probe will be nothing but a small camera. We decided early on that trying to visit each of the planets for any length of time would result in too complicated of a calculation, so instead we want to conduct a "flyby" and simply take a few nice photographs, thus allowing for a lower mass. In total, the mass of the probe will be about 100 kg, which is a reasonable estimate as modern ion thrusters are lightweight.
Now that ion thrusters are becoming ever smaller, longer-lasting, and powerful, more missions could be executed in the inner solar system and beyond. Our problem represents an example of a mission that could be executed with current day technology, but it is interesting to consider missions of much greater scale. Perhaps in the not-so-distant future a probe could be launched that visits all the planets of our solar system for a fraction of what it would cost with chemical propellants.
The data we use in our model comes from Jet Propulsion Library’s Horizons database (JPL Solar System Dynamics, n.d.) . This database contains ephemerides (tabulated movements) for various celestial bodies. We were able to generate cartesian coordinates for the position of the solar system starting on May 9th, 2016. The data uses a reference frame with the Sun at the origin.
We will be showcasing the best route for a space probe leaving Earth and visiting other inner solar system planets. First will be a mathematical explanation of our model. Following that will be a toy model to demonstrate the principles being applied. After that we will move into our more complex models and find an optimized route through the inner Solar System. Finally, we will discuss and interpret our results.
Our model is based in Newtonian physics. Technically, the only "real" variable we have control over is the 2-vector representing thrust of the probe. However, to simplify the constraints of our model, we chose to have variables for velocity and position in addition to thrust (each of these being 2-vectors as well).
As far as constraints are concerned, we have several that are applied at every time step to insure our probe follows the basic laws of kinematics.
In addition, we have several constraints that are applied only once.
Our model does not require it, but for our purposes we set initial position to be close to the starting planet and the final position to be the destination planet. By "close" to the starting planet we mean an offset was applied to prevent the gravity of the planet from being too strong for the probe to overcome. Additionally, the initial velocity was specified to match the starting planet's velocity (with respect to the sun). This does not take into account orbital velocity, but we chose to ignore this as it does not greatly impact the final solution given the distance at which we start our probe from the planet.
In our static models as well as our intermediate model, we are either trying to minimize fuel usage or journey duration. However, for our final model, we formulated the objective function to minimize fuel usage and applied $L_2$-regularization to smooth the x- and y-components of the thrust. This regularization makes our solution more feasible because in real life a rocket cannot fire in an arbitrary direction instantly—it takes time to pivot.
Additionally, with our objective function we assumed that fuel use as a function of thrust is quadratic in nature. Theoretically, by conservation of momentum it should be linear—when exhaust velocity is assumed to be constant, doubling the thrust means doubling the momentum of the exhaust and therefore doubling the mass/fuel loss. However, because our thrust can take on both positive and negative values, we would have to take the norm of the 2-vector thrust at every time step and add it up to get the total fuel usage. Instead we chose to square each component of thrust and add them together to obtain the magnitude squared, and it is the sum of these squared magnitudes that we use as our objective value. Making this assumption allows us to code the solution as well as simplifying our model (no square roots are involved).
Formulated in this manner, the only classification applicable to our particular model is NLP, a non-linear non-convex program. However, two details to note are that a). the problem itself is not necessarily non-convex, just our model of it, and b). all of our constraints could be part of a convex model if our gravity system were simpler. With five different gravity sinks, four of which are moving, the solution space becomes more convoluted and we can no longer assume it is convex. This removes any possiblity of efficient solving using specialized algorithms, and also forces us to keep our number of variables fairly small. However, by using a couple of tricks detailed later, we managed to overcome this and solve the problem on a large scale (~30,000 variables).
This section is meant to faciliate understanding of the final model in the same way we found the it—starting with the basics and working our way up the ladder of complexity.
To better illustrate the concepts and model we are applying, following are a couple of simplified models optimizing a trajectory through a static gravity system. When trying to optimize such a trajectory, one has to choose what they want to optimize. In our case, we narrowed the list down to whether to find the quickest route or the most fuel-efficient route, so we have provided examples illustrating the pros and cons of both. Below is a contour plot displaying the strength of gravity in the 2-sink gravity system we are using as an example. Before proceeding any further, make a guess as to what you think the optimal trajectory is to get from the start to the destination and then later compare your prediction to what is calculated.
using PyPlot
G = .05 # toy gravitational constant
sink1Pos = [8 1] # location of gravity sink 1 (x,y)
sink2Pos = [2 -1] # location of gravity sink 2 (x,y)
initPos = [0 0] # initial position (x,y)
destPos = [10 0] # destination position (x,y)
figure(figsize=(14,8))
buffer = .0000001 # removes hole in contour plot
x = -2:.1:12; # x and y axis limits
y = -4:.1:4;
#***** GENERATE GRAVITY DATA *****#
X = zeros(length(x), length(y))
Y = zeros(length(x), length(y))
for x_ in 1:length(x)
for y_ in 1:length(y)
X[x_,y_] = x[x_]
Y[x_,y_] = y[y_]
end
end
Z = (((G*(-X + sink1Pos[1])./((-X + sink1Pos[1]).^2 + (-Y + sink1Pos[2]).^2).^1.5 +
G*(-X + sink2Pos[1])./((-X + sink2Pos[1]).^2 + (-Y + sink2Pos[2]).^2).^1.5)).^2 +
((G*(-Y + sink1Pos[2])./((-X + sink1Pos[1]).^2 + (-Y + sink1Pos[2]).^2).^1.5 +
G*(-Y + sink2Pos[2])./((-X + sink2Pos[1]).^2 + (-Y + sink2Pos[2]).^2).^1.5)).^2) + buffer
Z = log(Z)
#***** PLOT DATA *****#
contourf(X, Y, Z, linspace(-17,2, 100))
plot(sink1Pos[1], sink1Pos[2], "mo", sink2Pos[1], sink2Pos[2], "mo")
plot(initPos[1], initPos[2], "kx")
plot(destPos[1], destPos[2], "kx")
legend(["Gravity Sinks"], "upper left")
text(initPos[1], initPos[2], " Initial Position")
text(destPos[1], destPos[2], " Destination")
text(8.5,-3.1,"Weaker Gravity")
text(2.5,-1.5,"Stronger Gravity")
text(4.1,0,"Equilibrium Point")
axis([-2, 12, -4, 4])
title("Static Gravity Environment")
xlabel("x []")
ylabel("y []")
;
One thing to keep in mind with our models, both the time- and fuel-optimal ones, is that we are constraining the probe to arrive at a specific time, i.e. removing the variability of arrival time. The journey is constrained to take the entries of durationSecs as its duration, so the solver is free to take as roundabout a way as it pleases to the goal so long as it minimizes fuel consumption and arrives on time. This can sometimes result in routes that take some interesting detours, as evidenced below.
durationSecs = [60,120,180] # array storing durations of each journey
G = .05 # toy gravitational constant
sink1Pos = [8 1] # location of gravity sink 1 (x,y)
sink2Pos = [2 -1] # location of gravity sink 2 (x,y)
thrustLimit = .01 # limit on magnitude of thrust vector
maxDist = .05 # max distance between destination and probe at end of journey
initPos = [0 0] # initial position (x,y)
initVelocity = [0 0] # final position (x,y)
destPos = [10 0] # destination position (x,y)
;
using JuMP, Ipopt, PyPlot
lims = round(durationSecs) # number of seconds accounted for in problem
buffer = .0000001 # makes contour plot look nicer
for n in lims
n = Int64(n)
m = Model(solver = IpoptSolver(print_level=0))
#***** DEFINE VARIABLES FOR PROBE POSITIONS, VELOCITIES, AND THRUSTS *****#
@variable(m, AX[1:n,1:2])
@variable(m, AV[1:n,1:2])
@variable(m, AU[1:n,1:2])
#***** DYNAMICS CONSTRAINTS *****#
for i in 1:n-1
# current position is previous + velocity
@constraint(m, AX[i+1,:] .== AX[i, :] + AV[i, :])
# current velocity is previous + thrust + gravity
@NLconstraint(m, AV[i+1, 1] == AV[i, 1] + AU[i, 1] +
G*(-AX[i, 1] + sink1Pos[1])/((-AX[i, 1] + sink1Pos[1])^2 + (-AX[i, 2] + sink1Pos[2])^2)^1.5 +
G*(-AX[i, 1] + sink2Pos[1])/((-AX[i, 1] + sink2Pos[1])^2 + (-AX[i, 2] + sink2Pos[2])^2)^1.5)
@NLconstraint(m, AV[i+1, 2] == AV[i, 2] + AU[i, 2] +
G*(-AX[i, 2] + sink1Pos[2])/((-AX[i, 1] + sink1Pos[1])^2 + (-AX[i, 2] + sink1Pos[2])^2)^1.5 +
G*(-AX[i, 2] + sink2Pos[2])/((-AX[i, 1] + sink2Pos[1])^2 + (-AX[i, 2] + sink2Pos[2])^2)^1.5)
end
#***** KEEP THRUST BELOW CERTAIN VALUE *****#
@constraint(m, (AU[:, 1].^2 + AU[:, 2].^2) .<= thrustLimit^2)
#***** START CONSTRAINTS *****#
@constraint(m, AX[1, :] .== initPos)
@constraint(m, AV[1, :] .== initVelocity)
#***** RENDEZVOUS CONSTRAINT *****#
@constraint(m, sum((AX[n, :] - destPos).^2) <= maxDist^2)
#***** MINIMIZE FUEL USAGE *****#
@objective(m, Min, sum(AU[:,1].^2 + AU[:,2].^2))
solve(m)
AX = getvalue(AX); AU = getvalue(AU)
x = -2:.1:12;
y = -4:.1:4;
X = zeros(length(x), length(y))
Y = zeros(length(x), length(y))
for x_ in 1:length(x)
for y_ in 1:length(y)
X[x_,y_] = x[x_]
Y[x_,y_] = y[y_]
end
end
Z = (((G*(-X + sink1Pos[1])./((-X + sink1Pos[1]).^2 + (-Y + sink1Pos[2]).^2).^1.5 +
G*(-X + sink2Pos[1])./((-X + sink2Pos[1]).^2 + (-Y + sink2Pos[2]).^2).^1.5)).^2 +
((G*(-Y + sink1Pos[2])./((-X + sink1Pos[1]).^2 + (-Y + sink1Pos[2]).^2).^1.5 +
G*(-Y + sink2Pos[2])./((-X + sink2Pos[1]).^2 + (-Y + sink2Pos[2]).^2).^1.5)).^2) + buffer
Z = log(Z)
figure(figsize=(14,8))
contourf(X, Y, Z, linspace(-17,2, 100))
plot(AX[:, 1], AX[:, 2], "k-", sink1Pos[1], sink1Pos[2], "mo", sink2Pos[1], sink2Pos[2], "mo")
plot(initPos[1], initPos[2], "kx")
plot(destPos[1], destPos[2], "kx")
grid("on")
axis([-2, 12, -4, 4])
title(string("Minimize Fuel Use: ", n, " seconds"))
text(-1.5,-3.5,string("Objective value is: ", getobjectivevalue(m)))
xlabel("x []")
ylabel("y []")
end
Notice that just by changing a single thing, the amount of time taken for the journey, the most fuel-efficient route can change drastically. The optimal route for $t = 60 \text{ seconds}$ is what you probably predicted. Such a route strikes a good balance between the amount of fuel used and the amount of time "wasted" on detours. An interesting thing to notice about these plots is the objective values of each. The route for $t = 60$ has a value of $\approx0.001952$ (no meaningful units) while the route for $t = 180$ has a value of $\approx0.001697$, indicating that the longer path is actually more fuel-efficient.
This can be explained by the way the $t = 180$ route uses gravity. It rides the gravity well of the first sink in the same way a penny rides a spiral wishing well, circling rather than resisting the pull "downward." After it reaches a certain point, it uses its thrusters to give it a little push out of the well and then flies along the gravitational equilibrium before taking the plunge and riding the other well to its destination.
The $t = 60$ route, on the other hand, takes a more direct approach, following an $S$-curve to utilize gravity to a certain extent but relying more heavily on thrusters to get the probe where it needs to be quicker. In a sense it fights through the waves rather than ride them, and this is evident in the plots above.
In the more complicated, dynamic gravity systems we use for optimizing a probe's trajectory in the inner solar system, striking a balance between time and fuel factors heavily into the way we formulate our model. Another tool we used to perform this balancing act was the minimization of journey duration.
For this model, on the other hand, we are trying to get the probe to its destination as fast as possible. To be clear, we would like to minimize t where time t is the first time at which the distance between the probe and the destination is less than the maxDist parameter. However, there is not a simple, standard way to constrain t thusly, and we were forced to improvise a little to optimize for trip duration.
In the end, to minimize the time taken we ended up constructing an objective function that strongly penalizes being far away from the destination, encouraging solutions that move toward the destination rapidly. We implemented this by squaring the distance between the destination and the probe and adding it up for all the time steps.
$$ f_{\text{obj}} = \sum_{i=1}^{n}(\mathbf{x}_i-\mathbf{x_{final}})\cdot(\mathbf{x}_i-\mathbf{x_{final}}) $$durationSecs = [75,60,45] # array storing durations of each journey
G = .05 # toy gravitational constant
sink1Pos = [8 1] # location of gravity sink 1 (x,y)
sink2Pos = [2 -1] # location of gravity sink 2 (x,y)
thrustLimit = .01 # limit on magnitude of thrust vector
maxDist = .05 # max distance between destination and probe after journey
initPos = [0 0] # initial position (x,y)
initVelocity = [0 0] # final position (x,y)
destPos = [10 0] # destination position (x,y)
;
using JuMP, Ipopt, PyPlot
lims = round(durationSecs) # number of seconds accounted for in problem
buffer = .0000001 # makes contour plot look nicer
for n in lims
n = Int64(n)
m = Model(solver = IpoptSolver(print_level=0))
#***** DEFINE VARIABLES FOR PROBE POSITIONS, VELOCITIES, AND THRUSTS *****#
@variable(m, AX[1:n,1:2])
@variable(m, AV[1:n,1:2])
@variable(m, AU[1:n,1:2])
#***** DYNAMICS CONSTRAINTS *****#
for i in 1:n-1
# current position is previous + velocity
@constraint(m, AX[i+1,:] .== AX[i, :] + AV[i, :])
# current velocity is previous + thrust + gravity
@NLconstraint(m, AV[i+1, 1] == AV[i, 1] + AU[i, 1] +
G*(-AX[i, 1] + sink1Pos[1])/((-AX[i, 1] + sink1Pos[1])^2 + (-AX[i, 2] + sink1Pos[2])^2)^1.5 +
G*(-AX[i, 1] + sink2Pos[1])/((-AX[i, 1] + sink2Pos[1])^2 + (-AX[i, 2] + sink2Pos[2])^2)^1.5)
@NLconstraint(m, AV[i+1, 2] == AV[i, 2] + AU[i, 2] +
G*(-AX[i, 2] + sink1Pos[2])/((-AX[i, 1] + sink1Pos[1])^2 + (-AX[i, 2] + sink1Pos[2])^2)^1.5 +
G*(-AX[i, 2] + sink2Pos[2])/((-AX[i, 1] + sink2Pos[1])^2 + (-AX[i, 2] + sink2Pos[2])^2)^1.5)
end
#***** KEEP THRUST BELOW CERTAIN VALUE *****#
@constraint(m, (AU[:, 1].^2 + AU[:, 2].^2) .<= thrustLimit^2)
#***** START CONSTRAINTS *****#
@constraint(m, AX[1, :] .== initPos)
@constraint(m, AV[1, :] .== initVelocity)
#***** RENDEZVOUS CONSTRAINT *****#
@constraint(m, sum((AX[n, :] - destPos).^2) <= maxDist^2)
#***** MINIMIZE TIME (APPROX.) *****#
@objective(m, Min, sum((AX[:,1]-destPos[1]).^2 + (AX[:,2]-destPos[2]).^2))
solve(m)
AX = getvalue(AX)
x = -2:.1:12;
y = -4:.1:4;
X = zeros(length(x), length(y))
Y = zeros(length(x), length(y))
for x_ in 1:length(x)
for y_ in 1:length(y)
X[x_,y_] = x[x_]
Y[x_,y_] = y[y_]
end
end
Z = (((G*(-X + sink1Pos[1])./((-X + sink1Pos[1]).^2 + (-Y + sink1Pos[2]).^2).^1.5 +
G*(-X + sink2Pos[1])./((-X + sink2Pos[1]).^2 + (-Y + sink2Pos[2]).^2).^1.5)).^2 +
((G*(-Y + sink1Pos[2])./((-X + sink1Pos[1]).^2 + (-Y + sink1Pos[2]).^2).^1.5 +
G*(-Y + sink2Pos[2])./((-X + sink2Pos[1]).^2 + (-Y + sink2Pos[2]).^2).^1.5)).^2) + buffer
Z = log(Z)
figure(figsize=(14,8))
contourf(X, Y, Z, linspace(-17,2, 100))
plot(AX[:, 1], AX[:, 2], "k-", sink1Pos[1], sink1Pos[2], "mo", sink2Pos[1], sink2Pos[2], "mo")
plot(initPos[1], initPos[2], "kx")
plot(destPos[1], destPos[2], "kx")
grid("on")
axis([-2, 12, -4, 4])
title(string("Minimize Duration: ", n, " seconds"))
xlabel("x []")
ylabel("y []")
end
When this objective function is used, changing the duration of the journey does not have as large of an impact on the final trajectory and the routes tend to look more similar, as evidenced above.
Additionally, the times we used were chosen specifically to demonstrate how we went about finding the fastest route from planet to planet in later models. Looking at the trajectory for $t=75 \text{ seconds}$, it is clear that the "best" route it could find is not the fastest possible. There is no need to make the detour that it does at the end, it is clearly just an unintended consequence of the objective function not being perfectly fine-tuned to our needs.
However, by systematically adjusting our $t$ value we can eventually find a reasonable estimate for the fastest possible time subject to our constraints, and in this particular instance this is approximately $t=45 \text{ seconds}$.
Our "Intermediate" model is named thusly because it is just that—an intermediate between the toy models we made to verify our model and the final product that routes a space probe between 4 planets. In it we use actual planetary data we obtained from NASA and the actual empirical values of constants (e.g. Gravitational constant, acceleration conversions, etc.).
This model uses the same objective function as the second proof-of-concept and its goal is to minimize the time to get from one planet to another. We ran the model sequentially through each pair of planets we planned to travel between and tabulated our results for use in the final model. Below is an example of calculating the fastest route for the first leg of the journey.
This section is just that—retrieving data for planetary positions and velocities, and setting the various parameters used in the model. Several things to note about units:
Additionally, it is worth mentioning that the variable startOffset is used to say how long after the start date of our data (May 9th) we wish to begin our model.
using NamedArrays
using JuMP
using PyPlot, Ipopt
# NOTE: Ephemerides must be in same directory as .ipynb file
# get Mercury data
rawMercury = readdlm("ScottPatenaude_Mercury.csv",',')
# get Venus data
rawVenus = readdlm("ScottPatenaude_Venus.csv",',')
# get Earth data
rawEarth = readdlm("ScottPatenaude_Earth.csv",',')
# get Mars data
rawMars = readdlm("ScottPatenaude_Mars.csv",',')
# masses in Kg
sunMass = 1.99e30
mercuryMass = 3.30e23
venusMass = 4.87e24
earthMass = 5.97e24
marsMass = 6.42e23
probeMass = 100
# miscellaneous variables
xCoor = 0 # x-coordinate of Sun
yCoor = 0 # y-coordinate of Sun
G = 2.98e-33 # Gravitational constant in units of newtons*AU^2/Kg^2
timeStep = 1/12 # in units of days
thrustLimit = .25 # in units of newtons
accelConv = .0499 # converts M/s^2 to AU/day^2
minDistance = .005 # how close is considered "visiting" a planet [Au] (2*earth-moon distance)
startOffset = (Int64)(0/timeStep)
journeyLength = 117 # how long journey takes in days (has to be manually adjusted)
dataPoints = Int64(round(journeyLength/timeStep))
# get position data
mercuryData = rawMercury[1:dataPoints,3:8]
venusData = rawVenus[1:dataPoints,3:8]
earthData = rawEarth[1:dataPoints,3:8]
marsData = rawMars[1:dataPoints,3:8]
startData = earthData
destData = marsData
;
This is a non-convex model, so just because a trajectory is a local minimum does not imply it is the global minimum. In order for the solver to succeed in finding a reasonable minimum, it has to be prodded toward the answer. We tried a few approaches to making an initial guess and settled on linearly moving between the initial planet's and destination planet's positions: $$ \mathbf{x}_{i}^{\text{guess}} = (1-\frac{i}{n})*\mathbf{x}_{i}^{\text{start}} + \frac{i}{n}*\mathbf{x}_{i}^{\text{dest}} \qquad \forall i \in {1,\dots,n} $$
We then set the velocities to whatever they need to be to result in the guessed positions. Using the simple kinematic law that $\mathbf{v} = \frac{d\mathbf{x}}{dt}$, we came up with this discrete version: $$ \mathbf{v}_{i}^{\text{guess}} = \tau\left(\mathbf{x}_{i+1}^{\text{guess}} - \mathbf{x}_{i}^{\text{guess}}\right) \qquad \forall i \in {1,\dots,n-1}\\ (\tau \text{ is time step in days}) $$
We also included an offset in the actual code to avoid gravity singularities.
m = Model(solver = IpoptSolver(print_level=0))
#***** DEFINE VARIABLES FOR PROBE POSITIONS, VELOCITIES, AND THRUSTS *****#
@variable(m, AX[1:2,1:dataPoints])
@variable(m, AV[1:2,1:dataPoints])
@variable(m, AU[1:2,1:dataPoints])
#***** SET INITIAL GUESS *****#
offSet = ones(2,1)*.005
initialGuessPos = zeros(2,dataPoints)
initialGuessVel = zeros(2,dataPoints)
# leg to Mars
for i in 1:dataPoints
factor = i/dataPoints
initialGuessPos[:,i] = (earthData[i,1:2]'*(1-factor)+marsData[i,1:2]'*factor + offSet)
end
for i in 1:dataPoints-1
initialGuessVel[:,i] = (initialGuessPos[:,i+1]-initialGuessPos[:,i])*timeStep
end
setvalue(AX, initialGuessPos)
setvalue(AV, initialGuessVel)
In this section we added the constraints to our model as described in the Mathematical Model section. Then we solve the model and display the time at which the probe first gets within maxDist of the destination.
#***** ADD IN CONSTRAINTS FOR PHYSICAL DYNAMICS OF ROCKET TRAVEL AS SPECIFIED IN PROBLEM *****#
for i in 1:dataPoints-1
# current position is previous + velocity*timestep
@constraint(m, AX[:,i+1] .== AX[:,i] + timeStep*AV[:,i])
# current velocity is previous + thrust + gravity
@NLconstraint(m, AV[1, i+1] == AV[1, i] + timeStep*accelConv*(AU[1, i]/probeMass +
sunMass*G*(-AX[1,i] + xCoor)/((-AX[1,i] + xCoor)^2 + (-AX[2,i] + yCoor)^2)^1.5 +
mercuryMass*G*(-AX[1,i] + mercuryData[i, 1])/((-AX[1,i] + mercuryData[i, 1])^2 +
(-AX[2,i] + mercuryData[i, 2])^2)^1.5 +
venusMass*G*(-AX[1,i] + venusData[i, 1])/((-AX[1,i] + venusData[i, 1])^2 + (-AX[2,i] + venusData[i, 2])^2)^1.5 +
earthMass*G*(-AX[1,i] + earthData[i, 1])/((-AX[1,i] + earthData[i, 1])^2 + (-AX[2,i] + earthData[i, 2])^2)^1.5 +
marsMass*G*(-AX[1,i] + marsData[i, 1])/((-AX[1,i] + marsData[i, 1])^2 + (-AX[2,i] + marsData[i, 2])^2)^1.5))# +
@NLconstraint(m, AV[2, i+1] == AV[2, i] + timeStep*accelConv*(AU[2, i]/probeMass +
sunMass*G*(-AX[2,i] + yCoor)/((-AX[1,i] + xCoor)^2 + (-AX[2,i] + yCoor)^2)^1.5 +
mercuryMass*G*(-AX[2,i] + mercuryData[i, 2])/((-AX[1,i] + mercuryData[i, 1])^2 +
(-AX[2,i] + mercuryData[i, 2])^2)^1.5 +
venusMass*G*(-AX[2,i] + venusData[i, 2])/((-AX[1,i] + venusData[i, 1])^2 + (-AX[2,i] + venusData[i, 2])^2)^1.5 +
earthMass*G*(-AX[2,i] + earthData[i, 2])/((-AX[1,i] + earthData[i, 1])^2 + (-AX[2,i] + earthData[i, 2])^2)^1.5 +
marsMass*G*(-AX[2,i] + marsData[i, 2])/((-AX[1,i] + marsData[i, 1])^2 + (-AX[2,i] + marsData[i, 2])^2)^1.5))
end
#***** KEEP THRUST BELOW CERTAIN VALUE *****#
@constraint(m, (AU[1,:][:].^2+AU[2,:][:].^2) .<= thrustLimit^2)
#***** INITIAL CONSTRAINTS FOR PROBE *****#
@constraint(m, AX[:,1] .== [startData[1,1] - .005, startData[1,2]])
@constraint(m, AV[:,1] .== [startData[1,4], startData[1,5]])
#***** RENDEZVOUS CONSTRAINTS FOR PROBE *****#
@constraint(m, (AX[1, dataPoints]-destData[dataPoints, 1])^2 + (AX[2, dataPoints]-destData[dataPoints, 2])^2 <= .02^2)
@objective(m, Min, sum((AX[1,:]-destData[:,1]').^2 + (AX[2,:]-destData[:,2]').^2))
#***** SOLVE *****#
@time(solve(m))
AX = getValue(AX)
mercuryDistance = AX[:,:] - mercuryData[:,1:2]'
venusDistance = AX[:,:] - venusData[:,1:2]'
marsDistance = AX[:,:] - marsData[:,1:2]'
firstOccurence = -1
for i in 1:dataPoints
dist = (marsDistance[1,i].^2+marsDistance[2,i].^2).^.5
if dist <= minDistance && firstOccurence == -1
firstOccurence = i
end
end
if firstOccurence != -1
println("Takes approx. ", round(firstOccurence*timeStep*1000)/1000, " days")
println(rawEarth[startOffset+1, 2], " to", rawEarth[startOffset+firstOccurence, 2])
else
println("Did not reach destination.")
end
1.781441 seconds (1.83 M allocations: 82.395 MB, 1.23% gc time) Takes approx. 113.917 days A.D. 2016-May-09 00:00:00.0000 to A.D. 2016-Aug-30 20:00:00.0000
#***** PLOT TRAJECTORIES *****#
figure(figsize=(10,7))
plot(AX[1,:][:],AX[2,:][:], "g-")
plot(initialGuessPos[1,:][:],initialGuessPos[2,:][:], "c-.")
plot(xCoor, yCoor, "oy")
plot(mercuryData[:,1], mercuryData[:,2], "b--",
venusData[:,1], venusData[:,2], "r--",
earthData[:,1], earthData[:,2], "m--",
marsData[:,1], marsData[:,2], "k--")
axis("equal")
grid("on")
title("Trajectories")
xlabel("x [AU]")
ylabel("y [AU]")
legend(["Probe's trajectory", "Initial Guess", "Sun", "Mercury", "Venus", "Earth", "Mars"], loc="upper left",
fontsize="12");
This plot represents a significant step toward our goal of traversing the inner solar system. Using our model we have found a plausible trajectory that minimizes the duration of the probe's journey through a dynamical gravity system. By iterating this solution over each pair of planets between which we wish to travel, we can find the shortest journey times. Then, we can apply these data as constraints in the final model in order to optimize a trip that visits several planets. Below is the data we found and applied to our final model.
Start | Destination | Time [days] | Departure | Arrival |
---|---|---|---|---|
Earth | Mars | 113.917 | May 9, 2016 00:00:00 | Aug 30, 2016 20:00:00 |
Mars | Venus | 170.667 | Aug 30, 2016 22:00:00 | Feb 17, 2017 12:00:00 |
Venus | Mercury | 119.25 | Feb 17, 2017 14:00:00 | Jun 16, 2017 18:00:00 |
Our final model is the culmination of everything we learned from the previous models. Using the tabulated data found with the intermediate model as additional constraints, we are now able to optimize a journey that visits the four inner planets of our solar system. The additional constraints are included in the same way the destination constraint is in the Mathematical Model section, namely:
$$ (\mathbf{x}_{i_w} - \mathbf{x}_{\text{waypoint}}) \cdot (\mathbf{x}_{i_w} - \mathbf{x}_{\text{waypoint}}) - d_{max}^2 \le 0\\ $$for each waypoint (flyby of planet).
This section is nearly identical to the corresponding one in the Intermediate Model; additionally, there is now a Lagrange multiplier for use in the regularization of the objective function.
using NamedArrays
using JuMP
using PyPlot, Ipopt
# NOTE: Ephemerides must be in same directory as .ipynb file
# get Mercury data
rawMercury = readdlm("ScottPatenaude_Mercury.csv",',')
# get Venus data
rawVenus = readdlm("ScottPatenaude_Venus.csv",',')
# get Earth data
rawEarth = readdlm("ScottPatenaude_Earth.csv",',')
# get Mars data
rawMars = readdlm("ScottPatenaude_Mars.csv",',')
# masses in Kg
sunMass = 1.99e30
mercuryMass = 3.30e23
venusMass = 4.87e24
earthMass = 5.97e24
marsMass = 6.42e23
probeMass = 100
# miscellaneous variables
xCoor = 0 # x-coordinate of Sun
yCoor = 0 # y-coordinate of Sun
G = 2.98e-33 # Gravitational constant in units of newtons*AU^2/Kg^2
timeStep = 1/12 # in units of days
thrustLimit = .25 # in units of newtons
accelConv = .0499 # converts M/s^2 to AU/day^2
minDistance = .005 # how close is considered "visiting" a planet [Au] (2*earth-moon distance)
λ = 100 #trade-off between minimizing thrust curve and smoothing it
startOffset = (Int64)(0/timeStep)
journeyLength = 403.834 # amount of time journey takes in days
#***** WAYPOINTS AS FOUND IN INTERMEDIATE MODEL *****#
marsWay = Int64(round(113.917/timeStep))
venusWay = Int64(round(284.584/timeStep))
mercuryWay = Int64(round(403.834/timeStep))
dataPoints = Int64(round(journeyLength/timeStep))
# get position/velocity data
mercuryData = rawMercury[1:dataPoints,3:8]
venusData = rawVenus[1:dataPoints,3:8]
earthData = rawEarth[1:dataPoints,3:8]
marsData = rawMars[1:dataPoints,3:8]
startData = earthData
destData = mercuryData
;
Again, this is similar to the corresponding section in the Intermediate Model. However, the initial guess now has to take into account traveling to three planets as opposed to only one. We account for this by breaking the journey down into its "legs" and making a guess for each.
m = Model(solver = IpoptSolver(print_level=0))
#***** DEFINE VARIABLES FOR PROBE POSITIONS, VELOCITIES, AND THRUSTS *****#
@variable(m, AX[1:2,1:dataPoints])
@variable(m, AV[1:2,1:dataPoints])
@variable(m, AU[1:2,1:dataPoints])
#***** SET INITIAL GUESS *****#
offSet = ones(2,1)*.005
initialGuessPos = zeros(2,dataPoints)
initialGuessVel = zeros(2,dataPoints)
# leg to Mars
for i in 1:marsWay
factor = i/marsWay
initialGuessPos[:,i] = (earthData[i,1:2]'*(1-factor)+marsData[i,1:2]'*factor + offSet)
end
for i in 1:marsWay-1
initialGuessVel[:,i] = (initialGuessPos[:,i+1]-initialGuessPos[:,i])*timeStep
end
# leg to Venus
off = marsWay
for i in marsWay:venusWay
factor = (i-off)/(venusWay-off)
initialGuessPos[:,i] = (marsData[i,1:2]'*(1-factor)+venusData[i,1:2]'*factor + offSet)
end
for i in marsWay:venusWay-1
initialGuessVel[:,i] = (initialGuessPos[:,i+1]-initialGuessPos[:,i])*timeStep
end
# leg to Mercury
off = venusWay
for i in venusWay:dataPoints
factor = (i-off)/(mercuryWay-off)
initialGuessPos[:,i] = (venusData[i,1:2]'*(1-factor)+mercuryData[i,1:2]'*factor + offSet)
end
for i in venusWay:dataPoints-1
initialGuessVel[:,i] = (initialGuessPos[:,i+1]-initialGuessPos[:,i])*timeStep
end
setvalue(AX, initialGuessPos)
setvalue(AV, initialGuessVel)
For this section, instead of having a single destination constraint as in the Intermediate Model, there are three separate waypoint constraints slated to occur at the specific times found a priori. Additionally we have added the new objective function that minimizes fuel use as well as smoothing the thrust using $L_2$-regularization.
#***** ADD IN CONSTRAINTS FOR PHYSICAL DYNAMICS OF ROCKET TRAVEL AS SPECIFIED IN PROBLEM *****#
for i in 1:dataPoints-1
# current position is previous + velocity*timestep
@constraint(m, AX[:,i+1] .== AX[:,i] + timeStep*AV[:,i])
# current velocity is previous + thrust + gravity
@NLconstraint(m, AV[1, i+1] == AV[1, i] + timeStep*accelConv*(AU[1, i]/probeMass +
sunMass*G*(-AX[1,i] + xCoor)/((-AX[1,i] + xCoor)^2 + (-AX[2,i] + yCoor)^2)^1.5 +
mercuryMass*G*(-AX[1,i] + mercuryData[i, 1])/((-AX[1,i] + mercuryData[i, 1])^2 +
(-AX[2,i] + mercuryData[i, 2])^2)^1.5 +
venusMass*G*(-AX[1,i] + venusData[i, 1])/((-AX[1,i] + venusData[i, 1])^2 + (-AX[2,i] + venusData[i, 2])^2)^1.5 +
earthMass*G*(-AX[1,i] + earthData[i, 1])/((-AX[1,i] + earthData[i, 1])^2 + (-AX[2,i] + earthData[i, 2])^2)^1.5 +
marsMass*G*(-AX[1,i] + marsData[i, 1])/((-AX[1,i] + marsData[i, 1])^2 + (-AX[2,i] + marsData[i, 2])^2)^1.5))# +
@NLconstraint(m, AV[2, i+1] == AV[2, i] + timeStep*accelConv*(AU[2, i]/probeMass +
sunMass*G*(-AX[2,i] + yCoor)/((-AX[1,i] + xCoor)^2 + (-AX[2,i] + yCoor)^2)^1.5 +
mercuryMass*G*(-AX[2,i] + mercuryData[i, 2])/((-AX[1,i] + mercuryData[i, 1])^2 +
(-AX[2,i] + mercuryData[i, 2])^2)^1.5 +
venusMass*G*(-AX[2,i] + venusData[i, 2])/((-AX[1,i] + venusData[i, 1])^2 + (-AX[2,i] + venusData[i, 2])^2)^1.5 +
earthMass*G*(-AX[2,i] + earthData[i, 2])/((-AX[1,i] + earthData[i, 1])^2 + (-AX[2,i] + earthData[i, 2])^2)^1.5 +
marsMass*G*(-AX[2,i] + marsData[i, 2])/((-AX[1,i] + marsData[i, 1])^2 + (-AX[2,i] + marsData[i, 2])^2)^1.5))
end
#***** KEEP THRUST BELOW CERTAIN VALUE *****#
@constraint(m, (AU[1,:][:].^2+AU[2,:][:].^2) .<= thrustLimit^2)
#***** INITIAL CONSTRAINTS FOR PROBE *****#
@constraint(m, AX[:,1] .== [startData[1,1] - .005, startData[1,2]])
@constraint(m, AV[:,1] .== [startData[1,4], startData[1,5]])
#***** WAYPOINT CONSTRAINTS FOR PROBE *****#
# Mars
@constraint(m, (AX[1, marsWay]-marsData[marsWay, 1])^2 +
(AX[2, marsWay]-marsData[marsWay, 2])^2 <= minDistance^2)
# Venus
@constraint(m, (AX[1, venusWay]-venusData[venusWay, 1])^2 +
(AX[2, venusWay]-venusData[venusWay, 2])^2 <= minDistance^2)
# Mercury
@constraint(m, (AX[1, mercuryWay]-mercuryData[mercuryWay, 1])^2 +
(AX[2, mercuryWay]-mercuryData[mercuryWay, 2])^2 <= minDistance^2)
@objective(m, Min, λ*sum((AU[1,1:dataPoints-1]-AU[1,2:dataPoints]).^2 + (AU[2,1:dataPoints-1]-AU[2,2:dataPoints]).^2) +
sum((AU[1,:].^2 + AU[2,:].^2)))
#***** SOLVE *****#
@time(solve(m))
AX = getvalue(AX)
AV = getvalue(AV)
AU = getvalue(AU)
mercuryDistance = AX[:,:] - mercuryData[:,1:2]'
venusDistance = AX[:,:] - venusData[:,1:2]'
earthDistance = AX[:,:] - earthData[:,1:2]'
marsDistance = AX[:,:] - marsData[:,1:2]'
firstOccurence = -1
for i in 1:dataPoints
dist = (mercuryDistance[1,i].^2+mercuryDistance[2,i].^2).^.5
if dist <= minDistance && firstOccurence == -1
firstOccurence = i
end
end
endDist = (mercuryDistance[1,dataPoints].^2+mercuryDistance[2,dataPoints].^2).^.5
if firstOccurence != -1
println("Takes approx. ", round(firstOccurence*timeStep*1000)/1000, " days")
println(rawEarth[startOffset+1, 2], " to", rawEarth[startOffset+firstOccurence, 2])
end
9.874584 seconds (8.07 M allocations: 313.417 MB, 2.62% gc time) Takes approx. 403.333 days A.D. 2016-May-09 00:00:00.0000 to A.D. 2017-Jun-16 06:00:00.0000
The route our model found takes $\approx403$ days to traverse the inner planets. Leaving on May 9th, 2016 our probe will reach Mercury, its final destination, on June 16, 2017.
using NamedArrays
using JuMP
using PyPlot, Ipopt
#***** PLOT TRAJECTORIES *****#
figure(figsize=(10,7))
plot(AX[1,:][:],AX[2,:][:], "g-")
plot(initialGuessPos[1,:][:],initialGuessPos[2,:][:], "c-.")
plot(xCoor, yCoor, "oy")
plot(mercuryData[:,1], mercuryData[:,2], "b--",
venusData[:,1], venusData[:,2], "r--",
earthData[:,1], earthData[:,2], "m--",
marsData[:,1], marsData[:,2], "k--")
plot(AX[1,marsWay],AX[2,marsWay], "m>")
plot(AX[1,venusWay],AX[2,venusWay], "m>")
plot(AX[1,mercuryWay],AX[2,mercuryWay], "m>")
axis("equal")
grid("on")
title("Trajectories")
xlabel("x [AU]")
ylabel("y [AU]")
legend(["Probe's trajectory", "Initial Guess", "Sun", "Mercury", "Venus", "Earth", "Mars", "Waypoints"], loc="upper left",
fontsize="12")
;
When plotted, the proposed trajectory for the probe looks reasonable. It bears a fair resemblance to the initial guess, which is to be expected when dealing with a non-convex model like ours. There are likely local minima located all over the solution space, so our result is tailored to look as it does.
figure(figsize=(12,5))
plot(collect(1:dataPoints)*timeStep, sqrt(AU[1,:][:].^2+AU[2,:][:].^2), "g-")
plot(collect(1:dataPoints)*timeStep, -AU[1,:][:], "r--")
plot(collect(1:dataPoints)*timeStep, -AU[2,:][:], "b--")
plot(collect(1:dataPoints)*timeStep, thrustLimit*ones(dataPoints), "c-.")
plot(marsWay*timeStep, sqrt(AU[1,marsWay]^2+AU[2,marsWay]^2), "m>")
plot(venusWay*timeStep, sqrt(AU[1,venusWay]^2+AU[2,venusWay]^2), "m>")
plot(mercuryWay*timeStep, sqrt(AU[1,mercuryWay]^2+AU[2,mercuryWay]^2), "m>")
axis([0, dataPoints*timeStep, -1.1*thrustLimit, 1.1*thrustLimit])
grid("on")
title("Thrust vs. Time")
xlabel("Time [days]")
ylabel("Thrust [Newtons]")
legend(["Total Thrust", "x-Thrust", "y-Thrust", "Thrust Limit", "Waypoints"], loc="lower left", fontsize="12")
;
The first thing to notice about this plot is that the only time our thrust is at the limit is near the beginning. This is when the probe is making its way away from the Sun. After reaching Mars, the thrust never has to rise as high again because it can act more as a guide and let gravity do the work, simply "falling with style" to the other planets.
A more subtle thing to see is that the thrust curve has local maxima located close to waypoints. This could have several explanations. One is that when near the waypoints, there is more gravity that necessitates burning more fuel. However, as shown in the Gravity section, the gravity of the planets, no matter the stage of the journey, is orders of magnitude less than the sun and does not impact the solution much.
Another explanation is that upon its destination, the probe fires its thrusters extra high to change its course to the next one. However, this explanation is more applicable to a human pilot. He or she would think of the journey in three separate parts and have to fire their thrusters to move in the right direction at the start of each leg, whereas when planning the journey as we are this is not necessary. Instead the solver could account for the directions needed for each waypoint and perhaps construct a more contiguous thrust curve. We suspect the reason for these local maxima lies in the way the solver itself works—we give it an initial guess consisting of 3 parts that are independently derived, so it gives us a solution that looks like it has three separate parts.
figure(figsize=(12,5))
plot(collect(1:dataPoints)*timeStep, sqrt(AV[1,:][:].^2+AV[2,:][:].^2)*1.731e6, "g-")
plot(marsWay*timeStep, sqrt(AV[1,marsWay]^2+AV[2,marsWay]^2)*1.731e6, "m>")
plot(venusWay*timeStep, sqrt(AV[1,venusWay]^2+AV[2,venusWay]^2)*1.731e6, "m>")
plot(mercuryWay*timeStep, sqrt(AV[1,mercuryWay]^2+AV[2,mercuryWay]^2)*1.731e6, "m>")
grid("on")
title("Velocity vs. Time")
xlabel("Time [days]")
ylabel("Velocity [m/s]")
legend(["Probe velocity", "Waypoints"], loc="upper left", fontsize="12")
;
The reason we include this plot is to perform a sanity-check of sorts on the physics of our problem. We wanted to insure two things; namely, insure classical kinematics are followed (e.g. Newtonian mechanics, Kepler's laws), and make sure we do not venture into the realm of special relativity ($v << c$).
As far as classical physics is concerned, this plot makes total sense. A good way to think about this is with Kepler's second law:
$$ \frac{d}{dt}\left( \mathbf{r}\times\mathbf{v} \right) = 0 $$Also known as the law of equal areas, this states that the area of the sector traced out by the probe with respect to the Sun per unit time does not change. When further away from the Sun, the probe should move more slowly than when it is closer. Indeed, this is the case for our probe as evident when comparing this plot to the Trajectories plot.
Additionally, it is clear from this plot that the velocity never exceeds 70,000 m/s. In special relativity, there is a number referred to as the $\gamma$-factor that at low speeds is close to 1. When this is the case, relativistic physics reduces to classical physics. Calculating $\gamma$ with our maximum speed yields:
$$ \gamma = \frac{1}{\sqrt{1-\frac{v^2}{c^2}}}\Bigg|_{v=70,000\text{m/s}} \approx 1.000000027 $$Therefore we do not need to worry about special relativity.
figure(figsize=(12,5))
plot(collect(1:dataPoints)*timeStep, ((mercuryDistance[1,:].^2+mercuryDistance[2,:].^2).^.5)', "b-")
plot(collect(1:dataPoints)*timeStep, ((venusDistance[1,:].^2+venusDistance[2,:].^2).^.5)', "r-")
plot(collect(1:dataPoints)*timeStep, ((earthDistance[1,:].^2+earthDistance[2,:].^2).^.5)', "m-")
plot(collect(1:dataPoints)*timeStep, ((marsDistance[1,:].^2+marsDistance[2,:].^2).^.5)', "k-")
plot(collect(1:dataPoints)*timeStep, ones(dataPoints)*6e-5, "c--")
plot(collect(1:dataPoints)*timeStep, ones(dataPoints)*minDistance, "y--")
plot(mercuryWay*timeStep, ((mercuryDistance[1,mercuryWay].^2+mercuryDistance[2,mercuryWay].^2).^.5)', "m>")
plot(venusWay*timeStep, ((venusDistance[1,venusWay].^2+venusDistance[2,venusWay].^2).^.5)', "m>")
plot(marsWay*timeStep, ((marsDistance[1,marsWay].^2+marsDistance[2,marsWay].^2).^.5)', "m>")
grid("on")
title("Distance from source vs. Time")
xlabel("Time [days]")
ylabel("Distance from source [AU]")
yscale("log")
legend(["Mercury", "Venus", "Earth", "Mars", "minDist", "maxDist"])
;
This plot is another sanity check used to make sure our constraints are being met. As shown, each of our waypoints gets close enough to be considered "visiting" each planet.
Also, nowhere in our journey do we get so close to the planet that the physics of the situation breaks down. We are modeling the planets as point-masses, so when the probe gets too close our time steps of 2 hours would not be sufficient to handle the changing accelerations. When we tried constraining the probe to not get within a certain distance of each planet, we greatly complicated the model and even with good initial guesses it often could not find a solution. This is to be expected because as it is formulated currently the set of allowable positions is nearly all of $\mathbb{R}^2$ (see the Mathematical Model section), while with the extra constraints there are holes that make the set more non-convex.
One more thing to note is how sharp the waypoints are. This is indicative of the fact that they are true flybys, too short for much to be done besides snap a few pictures. For Earth the distance drops off much more gradually because one of our initial constraints was to match the velocity of Earth as well as its location, so it takes more time to diverge.
figure(figsize=(12,5))
plot(collect(1:dataPoints)*timeStep,
(sunMass*probeMass*G*1*((-AX[1,:] + xCoor).^2 + (-AX[2,:] + yCoor).^2).^-.5)',"y-")
plot(collect(1:dataPoints)*timeStep,
(mercuryMass*probeMass*G*1*((-AX[1,:] + mercuryData[:,1]').^2 + (-AX[2,:] + mercuryData[:,2]').^2).^-.5)',"b-")
plot(collect(1:dataPoints)*timeStep,
(venusMass*probeMass*G*1*((-AX[1,:] + venusData[:,1]').^2 + (-AX[2,:] + venusData[:,2]').^2).^-.5)',"r-")
plot(collect(1:dataPoints)*timeStep,
(earthMass*probeMass*G*1*((-AX[1,:] + earthData[:,1]').^2 + (-AX[2,:] + earthData[:,2]').^2).^-.5)',"m-")
plot(collect(1:dataPoints)*timeStep,
(marsMass*probeMass*G*1*((-AX[1,:] + marsData[:,1]').^2 + (-AX[2,:] + marsData[:,2]').^2).^-.5)',"k-")
plot(marsWay*timeStep,
(marsMass*probeMass*G*1*((-AX[1,marsWay] + marsData[marsWay,1]').^2 +
(-AX[2,marsWay] + marsData[marsWay,2]').^2).^-.5)',"m>")
plot(venusWay*timeStep,
(venusMass*probeMass*G*1*((-AX[1,venusWay] + venusData[venusWay,1]').^2 +
(-AX[2,venusWay] + venusData[venusWay,2]').^2).^-.5)',"m>")
plot(mercuryWay*timeStep,
(mercuryMass*probeMass*G*1*((-AX[1,mercuryWay] + mercuryData[mercuryWay,1]').^2 +
(-AX[2,mercuryWay] + mercuryData[mercuryWay,2]').^2).^-.5)',"m>")
yscale("log")
grid("on")
title("Force of gravity vs. Time")
xlabel("Time [days]")
ylabel("Force [Newtons]")
legend(["Sun", "Mercury", "Venus", "Earth", "Mars", "Waypoints"], loc="upper left", fontsize="12")
;
This plot is similar to the Distances one in that it shows nice spikes right around the waypoints. In fact, one can transform one plot into the others using Newton's law of gravitation. However, we included it because we feel it does a nice job showing one simple fact—the gravity of the Sun is orders of magnitude above the gravity of the planets.
When we first formulated our model, we assumed we had to include the gravity of all the planets in order to make our solution realistic. Further, we had hopes of utilizing gravity-assists (slingshots) in the same way NASA space probes do sometimes in order to make our journey more fuel-efficient.
However, what we did not count on was lack of pull from other planets. We tried running our model without the gravity of other planets, and when we did the speed with which the solver finds a solution tripled while the solution itself was almost exactly the same. If we wanted to use gravity assists in the future, it would most likely require a combination of using a larger planet like Jupiter or Saturn, getting closer to the planets than we are, and shortening our time steps.
To summarize, the main limitations we have run into are:
If this problem were to be pursued further, there are several directions one could take. An interesting extension to our project could be limiting when we can use our thrusters. For instance, we could limit the number of discrete burns that a thruster can use. Our model currently assumes that a thruster can produce thrust for the duration of its trip. Limiting its thrust to only certain times would better simulate how a rocket would travel through space in using traditional chemical rockets. We expect that this extension would produce significantly different results than our current model.
An additional extension could be incorporating mass loss into the model. This would be more applicable to a chemical rocket than a ion thruster, as an ion thruster is very efficient with its fuel usage. A chemical rocket burning through its fuel will result in a dramatic change in mass during flight, and thus require less thrust to produce the same acceleration.
In the end, by formulating a model based upon Newtonian physics, we successfully found a route for a space probe equipped with an ion drive that visits all four inner planets. Due to limitations of our model and the nature of the problem, we cannot state that this is the absolute fastest or most fuel-efficient route; however, using heuristic techniques a good upper bound on the absolute best was calculated.
To determine the "truly best" route, a complete reformulation of our model would be required. One model we considered early on involved binary variables indicating whether the probe was visiting a planet or not, which would allow for variable times of arrival at each waypoint. We scrapped this idea because it is too complicated to be solved by any solver we have in our limited toolbox. It is not unreasonable, however, to say that such a model could be constructed and solved with the right solver. We acted according to the wise words of one Theodore Roosevelt: "Do what you can, with what you have, where you are."
Angrum, Andrea, Enrique Medina, and Daniel Sedlacko. "Voyager - The Interstellar Mission." Voyager - The Interstellar Mission. Jet Propulsion Laboratory, n.d. Web. 08 May 2016.
"Mariner Missions." NASA Science. NASA, n.d. Web. 05 May 2016.
Standish, Myles, Robert Jacobson, Steve Ostro, Alan Chamberlin, Don Yeomans, Jon Giorgini, and Paul Chodas. "HORIZONS Web-Interface." HORIZONS Web-Interface. JPL Solar System Dynamics, n.d. Web. 06 May 2016.
"Voyager - Planetary Voyage." Jet Propulsion Library. NASA, 17 Feb. 2015. Web. 05 May 2016.