connor/thesis
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Maybe nearing the final commit

Browse Source
This commit is contained in:
Connor
2022-04-18 10:59:25 -06:00
parent 0fb875c777
commit 5b6656a9a6
3598 changed files with 79349 additions and 572 deletions
Download Patch File Download Diff File Expand all files Collapse all files

42
julia/cost_table.jl Normal file
View File

@@ -0,0 +1,42 @@
using Thesis
using SPICE
using LinearAlgebra
using Dates
using PlotlyJS
try
furnsh("../../spice_files/naif0012.tls")
furnsh("../../spice_files/de430.bsp")
catch
furnsh("spice_files/naif0012.tls")
furnsh("spice_files/de430.bsp")
end
function cost(m::Union{Mission, Mission_Guess})
norm_mass = (m.start_mass - prop(m)[7]) / m.start_mass
norm_C3 = ( m.launch_v∞ ⋅ m.launch_v∞ ) / 200.
return 3norm_mass + norm_C3
end
function get_id(m::Mission)
letters = "E"
for phase in m.phases letters *= phase.planet.name[1] end
letters *= " "
letters *= string(cost(m))
return letters
end
missions = Vector{Mission}()
folders = readdir("/home/connor/projects/thesis/archive/saturn", join=true)
for folder in folders
try
push!(missions, ingest(folder * "/mission"))
catch e
println("Ingestion failure")
end
end
sorted_missions = sort(missions, by=cost)
@show get_id.(sorted_missions)

36
julia/find_neptune.jl
View File

@@ -1,9 +1,8 @@
include("src/Thesis.jl")
using .Thesis
using Thesis
using Dates
using SPICE
using LinearAlgebra
using PlotlyJS: savefig
using PlotlyJS
try
furnsh("../../spice_files/naif0012.tls")
@@ -25,7 +24,7 @@ function easy_cost(m::Union{Mission, Mission_Guess}, C3::Float64, v∞::Float64)
tof += phase.tof
end
tof = tof/(20Thesis.year)
return 3norm_mass + norm_C3 + tof
return 3norm_mass + norm_C3
end
sc = Sc("mySat", 200., 3200., 0.00025, 1, 1.0)
@@ -47,22 +46,25 @@ Thesis.mbh(fbs) = mbh(fbs,
max_v∞,
latest_arrival,
easy_cost,
search_patience=10_000,
drill_patience=200,
verbose=true)
search_patience=200,
drill_patience=50,
verbose=true,
)
function log(m::Mission, archive::Vector{Mission}, planets::Vector{Body})
function log(archive::Vector{Mission}, planets::Vector{Body})
abbrev = join([ planet.name[1] for planet in planets ])
p = plot(m, title="Best $(abbrev) Trajectory")
savefig(p,"archive/$(abbrev).html")
store(m, "archive/$(abbrev)")
store(archive, "archive/$(abbrev)_archive")
println(length(archive))
for m in archive
println(abbrev)
store(m, "archive/$(abbrev)_$(easy_cost(m, 200., 200.))")
end
end
log(_::Nothing, _::Vector{Mission}, _::Vector{Body}) = println("No Mission Found...")
planets = [Earth, Mars, Mars, Jupiter, Saturn]
best, archive = mbh(planets)
log(best, archive, planets)
println("Complete!")
planets = [Earth, Saturn]
latest_arrival = ideal_launch + Year(5)
best, archive, missions = mbh(planets)
log(missions, planets)
p = scatter(;x=1:length(missions), y=easy_cost.(missions, 200., 200.), mode="markers")
PlotlyJS.plot(p)

52
julia/ingest_mbh.jl Normal file
View File

@@ -0,0 +1,52 @@
using Thesis
using SPICE
using LinearAlgebra
using Dates
using PlotlyJS
try
furnsh("../../spice_files/naif0012.tls")
furnsh("../../spice_files/de430.bsp")
catch
furnsh("spice_files/naif0012.tls")
furnsh("spice_files/de430.bsp")
end
function cost(m::Union{Mission, Mission_Guess})
norm_mass = (m.start_mass - prop(m)[7]) / m.start_mass
norm_C3 = ( m.launch_v∞ ⋅ m.launch_v∞ ) / 200.
return 3norm_mass + norm_C3
end
files = Vector{String}()
contents = readdir("/home/connor/projects/thesis/archive", join=true)
for item in contents
if item[1:39] == "/home/connor/projects/thesis/archive/ES"
push!(files, item)
end
end
sort!(files, by=x->stat(x).ctime)
missions = Vector{Mission}()
for file in files push!(missions, ingest(file)) end
best_cost = Vector{Float64}()
global current_min = cost(missions[1])
for mission in missions
if cost(mission) < current_min global current_min = cost(mission) end
push!(best_cost, current_min)
end
p_dots = scatter(;y = cost.(missions), mode="markers")
p_best = scatter(;y = best_cost)
layout = Layout(;showlegend=false,
title=attr(
font_color="rgb(0,0,0)",
yanchor="top",
y=0.95,
text="Cost Functions Values of Single MBH run for Earth-Saturn Trajectory"
),
yaxis_title = "Cost Function Value",
)
display(PlotlyJS.plot([p_dots, p_best], layout))

24
julia/ingest_mission.jl
View File

@@ -16,14 +16,14 @@ function cost_1(m::Union{Mission, Mission_Guess})
norm_mass = (m.start_mass - prop(m)[7]) / m.start_mass
norm_C3 = ( m.launch_v∞ ⋅ m.launch_v∞ ) / 200.
tof = sum([p.tof for p in m.phases])/(20Thesis.year)
return 2norm_mass + norm_C3 + tof
return 3norm_mass + norm_C3 + tof
end
function cost_2(m::Union{Mission, Mission_Guess})
norm_mass = (m.start_mass - prop(m)[7]) / m.start_mass
norm_C3 = ( m.launch_v∞ ⋅ m.launch_v∞ ) / 200.
tof = sum([p.tof for p in m.phases])/(20Thesis.year)
return 2norm_mass + norm_C3# + tof
return 3norm_mass + norm_C3# + tof
end
function get_id(m::Mission)
@@ -74,16 +74,16 @@ println(get_id.(sorted_missions2[1:5]))
# camera=((-0.2, 0.1, -0.3),(-2.,-2.,1.))),
# "/home/connor/projects/thesis/LaTeX/fig/EMS_plot_noplanets.png",
# width=850, height=400)
savefig(Thesis.plot(long, title="EMJS Mission Profile", mode="light", planet_colors=false,
phase_colors=["#F00","#F0F","#6AF"], markers=false, legend=false,
camera=((-0.2, -0.3, -0.2),(1.75,-1.75,1.5))),
"/home/connor/projects/thesis/LaTeX/fig/EMJS_plot.png",
width=850, height=400)
savefig(Thesis.plot(long, title="EMJS Mission Profile", mode="light", planet_colors=false,
phase_colors=["#F00","#F0F","#6AF"], markers=false, planets=false, legend=false,
camera=((-0.2, -0.3, -0.2),(2.,-2.,1.5))),
"/home/connor/projects/thesis/LaTeX/fig/EMJS_plot_noplanets.png",
width=850, height=400)
# savefig(Thesis.plot(long, title="EMJS Mission Profile", mode="light", planet_colors=false,
# phase_colors=["#F00","#F0F","#6AF"], markers=false, legend=false,
# camera=((-0.2, -0.3, -0.2),(1.75,-1.75,1.5))),
# "/home/connor/projects/thesis/LaTeX/fig/EMJS_plot.png",
# width=850, height=400)
# savefig(Thesis.plot(long, title="EMJS Mission Profile", mode="light", planet_colors=false,
# phase_colors=["#F00","#F0F","#6AF"], markers=false, planets=false, legend=false,
# camera=((-0.2, -0.3, -0.2),(2.,-2.,1.5))),
# "/home/connor/projects/thesis/LaTeX/fig/EMJS_plot_noplanets.png",
# width=850, height=400)
# Display Thrust Plots
# display(plot_thrust(display_mission, title="Thrust Magnitude vs Time", mode="light"))
# display(plot_thrust_components(display_mission, title="Thrust Components vs Time", mode="light"))

61
julia/mars_mbh_analysis.jl Normal file
View File

@@ -0,0 +1,61 @@
using Thesis
using SPICE
using LinearAlgebra
using Dates
using PlotlyJS
try
furnsh("../../spice_files/naif0012.tls")
furnsh("../../spice_files/de430.bsp")
catch
furnsh("spice_files/naif0012.tls")
furnsh("spice_files/de430.bsp")
end
function cost(m::Union{Mission, Mission_Guess})
norm_mass = (m.start_mass - prop(m)[7]) / m.start_mass
norm_C3 = ( m.launch_v∞ ⋅ m.launch_v∞ ) / 200.
return 3norm_mass + norm_C3
end
function get_id(m::Mission)
letters = "E"
for phase in m.phases letters *= phase.planet.name[1] end
letters *= " "
letters *= string(cost(m))
return letters
end
missions = Vector{Mission}()
folders = readdir("/home/connor/projects/thesis/archive/EM_3-29", join=true)
for folder in folders
try
push!(missions, ingest(folder * "/mission"))
catch e
println("Ingestion failure")
end
end
sorted_missions = sort(missions, by=cost)
@show get_id.(missions)
# long, short = improve.(sorted_missions2[1:2], n=7)#, use_cost=true)
best_cost = Vector{Float64}()
for i in 1:length(missions)
push!(best_cost, minimum(cost.(missions[1:i])))
end
p_dots = scatter(;y = cost.(missions), mode="markers")
p_best = scatter(;y = best_cost)
layout = Layout(;showlegend=false,
title=attr(
font_color="rgb(0,0,0)",
yanchor="top",
y=0.95,
text="Cost Functions Values of Single MBH run for Earth-Mars Trajectory"
),
yaxis_title = "Cost Function Value",
)
display(PlotlyJS.plot([p_dots, p_best], layout))

71
julia/src/mbh.jl
View File

@@ -155,13 +155,15 @@ function mbh( flybys::Vector{Body},
search_patience::Int=10_000,
drill_patience::Int=50,
verbose::Bool=false,
test::Bool=false )
test::Bool=false,
)
# Convenience Functions
random_guess() = rand_mission_guess(flybys,sc,start_mass,launch_window,max_C3,max_v∞,latest_arrival)
cost(m::Union{Mission,Mission_Guess}) = cost_fn(m, max_C3, max_v∞)
cost(_::Nothing) = Inf
solve(g::Mission_Guess) = solve_mission(g, launch_window, latest_arrival, max_C3, max_v∞)
solve(g::Mission_Guess) = solve_mission(g, launch_window, latest_arrival, max_C3, max_v∞,
use_cost=true, CPU_multiplier=10., print_level=1)
# Initialize stuff
x_current = nothing
@@ -177,6 +179,8 @@ function mbh( flybys::Vector{Body},
print("\rMBH\t$(flyby_abbreviation)\tsearch: $(sct) drill: $(dct) archive: $(lar)\t")
end
all_missions = Vector{Mission}()
# The main loop
while search_count < search_patience
@@ -187,35 +191,37 @@ function mbh( flybys::Vector{Body},
x_star, x_lambert = solve.(random_guess())
# To improve solution time, we'll try drilling down the lambert's solution first
if x_lambert.converged
x_lambert_best = x_lambert
# If it does converge we check to see if it's better than the current
if cost(x_lambert) < cost(x_current) x_current = x_lambert end
while drill_count < drill_patience
verbose && log()
# If better than the best, then keep it as current
if x_lambert.converged && cost(x_lambert) < cost(x_current)
x_current = x_lambert
x_lambert_best = x_lambert
drill_count = 0
# If better than the best of lambert trajectories so far, keep it as lambert_best
elseif x_lambert.converged && cost(x_lambert) < cost(x_lambert_best)
x_lambert_best = x_lambert
drill_count = 0
# Otherwise, keep drilling
else
drill_count += 1
end
end
drill_count = 0
if x_current ∉ archive
push!(archive, x_current)
record(x_current, flybys)
end
end
# if x_lambert.converged
# x_lambert_best = x_lambert
# # If it does converge we check to see if it's better than the current
# if cost(x_lambert) < cost(x_current) x_current = x_lambert end
# push!(all_missions, x_lambert)
# while drill_count < drill_patience
# verbose && log()
# # If better than the best, then keep it as current
# if x_lambert.converged && cost(x_lambert) < cost(x_current)
# x_current = x_lambert
# x_lambert_best = x_lambert
# drill_count = 0
# # If better than the best of lambert trajectories so far, keep it as lambert_best
# elseif x_lambert.converged && cost(x_lambert) < cost(x_lambert_best)
# x_lambert_best = x_lambert
# drill_count = 0
# # Otherwise, keep drilling
# else
# drill_count += 1
# end
# end
# drill_count = 0
# if x_current ∉ archive
# push!(archive, x_current)
# record(x_current, flybys)
# end
# end
# Afterwards, we'll try a completely random guess
if x_star.converged
push!(all_missions, x_star)
x_basin = x_star
# If it does converge we check to see if it's better than the current
if cost(x_star) < cost(x_current) x_current = x_star end
@@ -226,14 +232,21 @@ function mbh( flybys::Vector{Body},
x_star = solve(perturb(x_basin))
# If better than the best, then keep it as current
if x_star.converged && cost(x_star) < cost(x_current)
println(cost(x_star), '\t', cost(x_basin))
push!(all_missions, x_star)
x_current = x_star
x_basin = x_star
drill_count = 0
# If better than the best in this particular basin, keep it as basin
elseif x_star.converged && cost(x_star) < cost(x_basin)
println(cost(x_star), '\t', cost(x_basin))
push!(all_missions, x_star)
x_basin = x_star
drill_count = 0
# Otherwise, keep drilling
elseif x_star.converged
push!(all_missions, x_star)
drill_count += 1
else
drill_count += 1
end
@@ -254,7 +267,7 @@ function mbh( flybys::Vector{Body},
if verbose println() end
return x_current, archive
return x_current, archive, all_missions
end

21
julia/src/nlp_solver.jl
View File

@@ -8,15 +8,16 @@ function constraint_bounds(guess::Mission_Guess)
for phase in guess.phases
state_constraint = [10_000., 10_000., 10_000., 0.02, 0.02, 0.02]
state_constraint = zeros(6)
if phase != guess.phases[end]
# Constraints on the final state
push!(low_constraint, (-1 * state_constraint)... )
push!(low_constraint, (state_constraint)... )
push!(high_constraint, state_constraint... )
# Constraints on the v∞ and turning angle
push!(low_constraint, -0.02, 100.)
push!(high_constraint, 0.02, 1e7)
push!(low_constraint, 0., 100.)
push!(high_constraint, 0., Inf)
else
push!(low_constraint, (-1 * state_constraint[1:3])... )
push!(low_constraint, (state_constraint[1:3])... )
push!(high_constraint, state_constraint[1:3]... )
end
end
@@ -39,7 +40,7 @@ function solve_mission( guess::Mission_Guess,
latest_arrival::DateTime,
max_C3::Float64,
max_v∞::Float64;
tol=1e-4,
tol=1e-1,
verbose::Bool=false,
print_level=0,
use_cost=false,
@@ -100,7 +101,7 @@ function solve_mission( guess::Mission_Guess,
g[8*(length(flybys)-1)+5] = max_v∞ - norm(v∞_in)
if use_cost
return -final[end]
return 3*(3500-final[end])/3500 + norm(x[2:4])^2 / 200.
else
return 1.0
end
@@ -108,8 +109,8 @@ function solve_mission( guess::Mission_Guess,
catch e
if isa(e,LaGuerreConway_Error) || isa(e, Mass_Error)
g .= 1e10
return 1e10
g .= 1e14
return 1e14
else
rethrow()
end
@@ -124,7 +125,11 @@ function solve_mission( guess::Mission_Guess,
g_low, g_high = constraint_bounds(guess)
ipopt_options = Dict("tol" => tol,
"acceptable_tol" => 100tol,
"constr_viol_tol" => 1e-1,
"acceptable_constr_viol_tol" => 1e3,
"max_cpu_time" => CPU_multiplier * length(guess.phases),
"expect_infeasible_problem" => "yes",
"dual_inf_tol" => 1e3,
"print_level" => print_level)
options = Options(solver=IPOPT(ipopt_options), derivatives=ForwardFD())