diff --git a/LaTeX/flowcharts/mbh.png b/LaTeX/flowcharts/mbh.png
new file mode 100644
index 0000000..2e272c3
Binary files /dev/null and b/LaTeX/flowcharts/mbh.png differ
diff --git a/LaTeX/flowcharts/mbh.svg b/LaTeX/flowcharts/mbh.svg
new file mode 100644
index 0000000..7fe629f
--- /dev/null
+++ b/LaTeX/flowcharts/mbh.svg
@@ -0,0 +1,915 @@
+
+
+
+
diff --git a/LaTeX/flowcharts/nlp.png b/LaTeX/flowcharts/nlp.png
new file mode 100644
index 0000000..9afc1ae
Binary files /dev/null and b/LaTeX/flowcharts/nlp.png differ
diff --git a/LaTeX/flowcharts/nlp.svg b/LaTeX/flowcharts/nlp.svg
new file mode 100644
index 0000000..5f556aa
--- /dev/null
+++ b/LaTeX/flowcharts/nlp.svg
@@ -0,0 +1,1194 @@
+
+
+
+
diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl
index 7ca2c0d..7f51b72 100644
--- a/julia/src/Thesis.jl
+++ b/julia/src/Thesis.jl
@@ -25,8 +25,8 @@ module Thesis
include("./plotting.jl")
include("./inner_loop/laguerre-conway.jl")
include("./inner_loop/propagator.jl")
- include("./inner_loop/phase.jl")
- include("./inner_loop/monotonic_basin_hopping.jl")
+ include("./inner_loop/nlp_solver.jl")
+ # include("./inner_loop/monotonic_basin_hopping.jl")
# include("./outer_loop.jl")
end
diff --git a/julia/src/errors.jl b/julia/src/errors.jl
index 4b0a6b6..0b4f9ad 100644
--- a/julia/src/errors.jl
+++ b/julia/src/errors.jl
@@ -26,3 +26,8 @@ struct Mass_Error{T} <: Exception where T <: Real
mass::T
end
Base.showerror(io::IO, e::Mass_Error) = print(io, "mass (", e.mass, ") got too low in propagation")
+
+struct Mission_Creation_Error{T} <: Exception where T <: Real
+ num_entries::T
+end
+Base.showerror(io::IO, e::Mission_Creation_Error) = print(io, "Tried to create a mission with $(e.num_entries) entries")
diff --git a/julia/src/inner_loop/nlp_solver.jl b/julia/src/inner_loop/nlp_solver.jl
new file mode 100644
index 0000000..eace95e
--- /dev/null
+++ b/julia/src/inner_loop/nlp_solver.jl
@@ -0,0 +1,110 @@
+using SNOW
+
+export solve_mission
+
+struct Result
+ converged::Bool
+ info::Symbol
+ sol::Matrix{Float64}
+end
+
+"""
+This function will take an initial guess to a mission profile and leverage Ipopt to solve the mission
+"""
+function solve_mission( guess::Mission_Guess;
+ tol=1e-8,
+ verbose::Bool=false )
+ # First we need to define a function that propagates the mission all the way through
+ # The function takes in a constraints vector (g) and an input vector (x) and spits out
+ # the objective function. In my case, the objective will be constant.
+ #
+
+ # First we define our starting point
+ x0 = Vector(guess)
+ n = size(guess.phases[1].thrust_profile)[1]
+ flybys = [ p.planet for p in guess.phases ]
+
+ # And our NLP
+ function optimizer!(g,x)
+ # Establish initial conditions
+ mass = guess.start_mass
+ v∞_out = x[2:4]
+ current_planet = Earth
+ launch_date = Dates.unix2datetime(x[1])
+ time = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS"))
+ start = state(current_planet, time, v∞_out, mass)
+
+ # Now, for each phase we must require:
+ # - That the ending state matches (unless final leg)
+ # - That the v∞_out == v∞_in (unless final leg)
+ # - That the minimum height is acceptable (unless final leg)
+ # - That the ending position matches (if final leg)
+ # - That the ending mass is acceptable (if final leg)
+ i = 0
+ for flyby in flybys
+ phase_params = x[ 5 + (3n+7) * i : 5 + (3n+7) * (i+1) - 1 ]
+ v∞_in = phase_params[1:3]
+ v∞_out = phase_params[4:6]
+ tof = phase_params[7]
+ thrusts = reshape(phase_params[8:end], (n,3))
+ current_planet = flyby
+ time += tof
+ goal = state(current_planet, time, v∞_in)
+ final = prop(reshape(thrusts, (n,3)), start, guess.sc, tof)[2]
+ if flyby != flybys[end]
+ g[1+8i:6+8i] .= final[1:6] .- goal[1:6]
+ g[7+8i] = norm(v∞_out) - norm(v∞_in)
+ δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
+ g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
+ else
+ g[8length(flybys)+1:8length(flybys)+3] .= final[1:3] .- goal[1:3]
+ g[8length(flybys)+4] = final[7]
+ end
+ start = state(current_planet, time, v∞_out, final[7])
+ i += 1
+ return 1.0
+ end
+ end
+ g = zeros(8length(flybys)+4)
+ optimizer!(g,x0)
+ return g
+
+
+
+ # n = size(x0)[1]
+
+ # function f!(g,x)
+ # try
+ # g[1:6] .= prop(reshape(x,(n,3)), start, craft, tof, primary)[2][1:6]
+ # return 1.
+ # catch e
+ # if isa(e,LaGuerreConway_Error) || isa(e, Mass_Error) || isa(e, PropOne_Error)
+ # return 10_000.
+ # else
+ # rethrow
+ # end
+ # end
+ # end
+
+ # lower_x = -1 * ones(3n)
+ # upper_x = ones(3n)
+ # num_constraints = 6
+ # filename = "ipopt_"*string(rand(UInt))
+ # ipopt_options = Dict("constr_viol_tol" => tol,
+ # "acceptable_constr_viol_tol" => 1e-4,
+ # "max_iter" => 10_000,
+ # "max_cpu_time" => 30.,
+ # "print_level" => 0,
+ # "output_file" => filename)
+ # options = Options(solver=IPOPT(ipopt_options),
+ # derivatives=ForwardAD())
+ # x, _, info = minimize(f!, vec(x0), num_constraints, lower_x, upper_x, final[1:6], final[1:6], options)
+ # rm(filename)
+ # if info in [:Solve_Succeeded, :Solved_To_Acceptable_Level]
+ # return Result(true, info, reshape(x,(n,3)))
+ # else
+ # if verbose println(info) end
+ # return Result(false, info, x0)
+ # end
+
+end
diff --git a/julia/src/inner_loop/phase.jl b/julia/src/inner_loop/phase.jl
deleted file mode 100644
index 09f9aad..0000000
--- a/julia/src/inner_loop/phase.jl
+++ /dev/null
@@ -1,60 +0,0 @@
-using SNOW
-
-export solve_phase
-
-struct Result
- converged::Bool
- info::Symbol
- sol::Matrix{Float64}
-end
-
-"""
-This function will take a single phase (so an initial state, and a final state) and an initial guess
-to the thrust profile and use an NLP solver to find the nearest thrust profile to that initial guess
-that satisfies the final state condition
-"""
-function solve_phase( start::Vector{Float64},
- final::Vector{Float64},
- craft::Sc,
- tof::Float64,
- x0::Matrix{Float64},
- primary::Body=Sun;
- tol=1e-8,
- verbose::Bool=false )
- n = size(x0)[1]
-
- function f!(g,x)
- try
- g[1:6] .= prop(reshape(x,(n,3)), start, craft, tof, primary)[2][1:6]
- return 1.
- catch e
- if isa(e,LaGuerreConway_Error) || isa(e, Mass_Error) || isa(e, PropOne_Error)
- return 10_000.
- else
- rethrow
- end
- end
- end
-
- lower_x = -1 * ones(3n)
- upper_x = ones(3n)
- num_constraints = 6
- filename = "ipopt_"*string(rand(UInt))
- ipopt_options = Dict("constr_viol_tol" => tol,
- "acceptable_constr_viol_tol" => 1e-4,
- "max_iter" => 10_000,
- "max_cpu_time" => 30.,
- "print_level" => 0,
- "output_file" => filename)
- options = Options(solver=IPOPT(ipopt_options),
- derivatives=ForwardAD())
- x, _, info = minimize(f!, vec(x0), num_constraints, lower_x, upper_x, final[1:6], final[1:6], options)
- rm(filename)
- if info in [:Solve_Succeeded, :Solved_To_Acceptable_Level]
- return Result(true, info, reshape(x,(n,3)))
- else
- if verbose println(info) end
- return Result(false, info, x0)
- end
-
-end
diff --git a/julia/src/mission.jl b/julia/src/mission.jl
index 250154b..65ab7c5 100644
--- a/julia/src/mission.jl
+++ b/julia/src/mission.jl
@@ -1,6 +1,7 @@
export Phase, Mission_Guess, Mission, Bad_Mission
export test_phase1, test_phase2
export test_mg, test_mission
+export Vector
struct Phase
planet::Body
@@ -55,6 +56,51 @@ function Mission_Guess(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float
Mission_Guess(sc, mass, date, v∞, phases, false)
end
+"""
+This is the opposite of the function below. It takes a vector and any other necessary information and outputs a mission guess
+"""
+function Mission_Guess(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body})
+ # Variable mission params
+ launch_date = Dates.unix2datetime(x[1])
+ launch_v∞ = x[2:4]
+
+ # Try to intelligently determine n
+ three_n = ((length(x)-4)/length(flybys)) - 7
+ three_n % 3 == 0 || throw(Mission_Creation_Error(three_n))
+ n::Int = three_n/3
+
+ # Build the phases
+ i = 0
+ phases = Vector{Phase}()
+ for flyby in flybys
+ phase_params = x[ 5 + (3n+7) * i : 5 + (3n+7) * (i+1) - 1 ]
+ v∞_in = phase_params[1:3]
+ v∞_out = phase_params[4:6]
+ tof = phase_params[7]
+ thrusts = reshape(phase_params[8:end], (n,3))
+ push!(phases, Phase(flyby, v∞_in, v∞_out, tof, thrusts))
+ i += 1
+ end
+ return Mission_Guess(sc, mass, launch_date, launch_v∞, phases, false)
+end
+
+"""
+This is a function to convert a mission guess into an nlp-ready vector. It doesn't contain all
+information from the guess though.
+"""
+function Base.Vector(g::Mission_Guess)
+ result = Vector{Float64}()
+ push!(result, Dates.datetime2unix(g.launch_date))
+ push!(result, g.launch_v∞...)
+ for phase in g.phases
+ push!(result,phase.v∞_in...)
+ push!(result,phase.v∞_out...)
+ push!(result,phase.tof)
+ push!(result,phase.thrust_profile...)
+ end
+ return result
+end
+
const test_mg = Mission_Guess(bepi, 12_000., DateTime(1992,11,19), [-3.4,1.2,0.1], test_phases)
struct Mission
diff --git a/julia/test/inner_loop/nlp_solver.jl b/julia/test/inner_loop/nlp_solver.jl
new file mode 100644
index 0000000..49787f1
--- /dev/null
+++ b/julia/test/inner_loop/nlp_solver.jl
@@ -0,0 +1,36 @@
+@testset "Phase" begin
+
+ println("Testing NLP solver")
+
+ # We'll start by testing the mission_guess -> vector function
+ vec = Vector(test_mg)
+ @test typeof(vec) == Vector{Float64}
+
+ # Now we go in the other direction
+ flybys = [ p.planet for p in test_mg.phases ]
+ guess = Mission_Guess(vec, test_mg.sc, test_mg.start_mass, flybys)
+ @test typeof(guess) == Mission_Guess
+ @test guess.sc == test_mg.sc
+ @test guess.start_mass == test_mg.start_mass
+ @test guess.launch_date == test_mg.launch_date
+ @test guess.launch_v∞ == test_mg.launch_v∞
+ @test guess.converged == test_mg.converged
+ for i in 1:length(guess.phases)
+ @test guess.phases[i].planet == test_mg.phases[i].planet
+ @test guess.phases[i].v∞_in == test_mg.phases[i].v∞_in
+ @test guess.phases[i].v∞_out == test_mg.phases[i].v∞_out
+ @test guess.phases[i].tof == test_mg.phases[i].tof
+ @test guess.phases[i].thrust_profile == test_mg.phases[i].thrust_profile
+ end
+
+ # Now we test an example run of the basic "inner function"
+ leave, arrive = DateTime(1992,11,19), DateTime(1993,4,1)
+ earth_state = state(Earth, leave)
+ venus_state = state(Venus, arrive)
+ v∞_out, v∞_in, tof = Thesis.lamberts(Earth, Venus, leave, arrive)
+ phase = Phase(Venus, v∞_in, v∞_in, tof, zeros(20,3))
+ guess = Mission_Guess(bepi, 12_000., leave, v∞_out, [phase])
+ g = solve_mission(guess)
+ println(g)
+
+end
diff --git a/julia/test/inner_loop/phase.jl b/julia/test/inner_loop/phase.jl
deleted file mode 100644
index 3f428d3..0000000
--- a/julia/test/inner_loop/phase.jl
+++ /dev/null
@@ -1,74 +0,0 @@
-@testset "Phase" begin
-
- println("Testing NLP solver")
-
- using SNOW, PlotlyJS
-
- # First we'll test an Earth case
- # Initial Setup
- T = rand( 8hour : second : day )
- revs = 8
- n = 10 * revs
- start_mass = 10_000.
-
- # A simple orbit raising
- start = gen_orbit(T, start_mass, Earth)
- thrust = spiral(0.9, n, start, test_sc, revs*T, Earth)
- final = prop(thrust, start, test_sc, revs*T, Earth)[2]
- new_T = 2π * √( xyz_to_oe(final, Earth.μ)[1]^3 / Earth.μ )
-
- # This should be close enough to converge, but Earth orbits are picky
- thrust_guess = spiral(0.89, n, start, test_sc, revs*T, Earth)
- guess_final = prop(thrust_guess, start, test_sc, revs*T, Earth)[2]
- result = solve_phase(start, final, test_sc, revs*T, thrust_guess, Earth, verbose=true)
- calc_path, calc_final = prop(result.sol, start, test_sc, revs*T, Earth, interpolate=true)
-
- # Test
- @test norm(guess_final[1:6] - final[1:6]) > 1e-4
- @test result.converged
- @test norm(calc_final[1:6] - final[1:6]) < 1e-4
-
- # Plot
- paths = Pathlist()
- push!(paths, prop(start, T, Earth),
- calc_path,
- prop(calc_final, T, Earth),
- prop(final, T, Earth) )
- fig = plot_orbits(paths, Earth,
- labels=["init", "transit", "post-transit", "final"],
- colors=["#FFF","#F44","#4F4","#44F"])
- savefig(fig, "../plots/nlp_test_earth.html")
-
- # Now the sun case
- T = rand( 0.5year : hour : 1.5year )
- tof = 0.7T
- n = 20
- start_mass = 12_000.
-
- # A simple orbit raising again, to keep things simple
- start = gen_orbit(T, start_mass)
- thrust = spiral(0.6, n, start, bepi, tof)
- final = prop(thrust, start, bepi, tof)[2]
- new_T = 2π * √( xyz_to_oe(final, Sun.μ)[1]^3 / Sun.μ )
-
- # This should be close enough to converge
- thrust_guess = spiral(0.55, n, start, bepi, tof)
- guess_final = prop(thrust_guess, start, bepi, tof)[2]
- result = solve_phase(start, final, bepi, tof, thrust_guess, verbose=true)
- calc_path, calc_final = prop(result.sol, start, bepi, tof, interpolate=true)
-
- # Test
- @test norm(guess_final[1:6] - final[1:6]) > 1e-4
- @test result.converged
- @test norm(calc_final[1:6] - final[1:6]) < 1e-4
-
- # Plot
- paths = Pathlist()
- push!(paths, prop(start, T), calc_path, prop(calc_final, T), prop(final, T) )
- fig = plot_orbits(paths,
- labels=["init", "transit", "post-transit", "final"],
- colors=["#FFF","#F44","#4F4","#44F"])
- savefig(fig, "../plots/nlp_test_sun.html")
-
-
-end
diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl
index 46ca670..841f503 100644
--- a/julia/test/runtests.jl
+++ b/julia/test/runtests.jl
@@ -14,9 +14,9 @@ catch
end
@testset "All Tests" begin
- include("plotting.jl")
- include("inner_loop/laguerre-conway.jl")
- include("inner_loop/propagator.jl")
- include("inner_loop/phase.jl")
- include("inner_loop/monotonic_basin_hopping.jl")
+ # include("plotting.jl")
+ # include("inner_loop/laguerre-conway.jl")
+ # include("inner_loop/propagator.jl")
+ include("inner_loop/nlp_solver.jl")
+ # include("inner_loop/monotonic_basin_hopping.jl")
end