Additional display functions and section 6 draft done

This commit is contained in:
Connor
2022-02-20 22:10:56 -07:00
parent 472662fb60
commit 700f7e28cf
169 changed files with 810 additions and 944 deletions

View File

@@ -20,7 +20,12 @@ function easy_cost(m::Union{Mission, Mission_Guess}, C3::Float64, v∞::Float64)
norm_mass = (m.start_mass - prop(m)[7]) / m.start_mass
norm_C3 = ( m.launch_v∞ m.launch_v∞ ) / C3
norm_v∞ = norm(m.phases[end].v∞_in) / v∞
return 3norm_mass + norm_C3# + norm_v∞
tof = 0.0
for phase in m.phases
tof += phase.tof
end
tof = tof/(20Thesis.year)
return 3norm_mass + norm_C3 + tof
end
sc = Sc("mySat", 200., 3200., 0.00025, 1, 1.0)
@@ -56,7 +61,7 @@ end
log(_::Nothing, _::Vector{Mission}, _::Vector{Body}) = println("No Mission Found...")
planets = [Earth, Venus, Venus, Jupiter, Saturn]
planets = [Earth, Mars, Mars, Jupiter, Saturn]
best, archive = mbh(planets)
log(best, archive, planets)

61
julia/ingest_mission.jl Normal file
View File

@@ -0,0 +1,61 @@
using Thesis
using SPICE
using LinearAlgebra
using Dates
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_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
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
end
function get_id(m::Mission)
letters = "E"
for phase in m.phases
letters *= phase.planet.name[1]
end
letters *= " "
letters *= Dates.format(m.launch_date, "yyyy-mm-dd HH:MM:SS")
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_missions1 = sort(missions, by=cost_1)
sorted_missions2 = sort(missions, by=cost_2)
println(get_id.(sorted_missions1[1:10]))
println(get_id.(sorted_missions2[1:10]))
display_mission = improve(sorted_missions2[1], n=5)
# display(plot(display_mission,
# title="Sample Algorithm Result Mission",
# mode="light",
# planet_colors=false,
# phase_colors=["#F00","#0F0","#00F"],
# markers=false));
store(display_mission, "/home/connor/projects/thesis/archive/best/long_mission")

View File

@@ -7,18 +7,17 @@ function constraint_bounds(guess::Mission_Guess)
high_constraint = Vector{Float64}()
for phase in guess.phases
# Constraints on the final state
state_constraint = [10_000., 10_000., 10_000., 0.05, 0.05, 0.05]
push!(low_constraint, (-1 * state_constraint)... )
push!(high_constraint, state_constraint... )
state_constraint = [10_000., 10_000., 10_000., 0.02, 0.02, 0.02]
if phase != guess.phases[end]
# Constraints on the final state
push!(low_constraint, (-1 * state_constraint)... )
push!(high_constraint, state_constraint... )
# Constraints on the v∞ and turning angle
push!(low_constraint, -0.05, 100.)
push!(high_constraint, 0.05, 1e7)
push!(low_constraint, -0.02, 100.)
push!(high_constraint, 0.02, 1e7)
else
# Constraints on mass
push!(low_constraint, 0.)
push!(high_constraint, guess.start_mass - guess.sc.dry_mass)
push!(low_constraint, (-1 * state_constraint[1:3])... )
push!(high_constraint, state_constraint[1:3]... )
end
end
@@ -40,9 +39,11 @@ function solve_mission( guess::Mission_Guess,
latest_arrival::DateTime,
max_C3::Float64,
max_v∞::Float64;
tol=1e-12,
tol=1e-4,
verbose::Bool=false,
print_level=0)
print_level=0,
use_cost=false,
CPU_multiplier=10.)
# First we define our starting point
x0 = Vector(guess)
@@ -58,9 +59,11 @@ function solve_mission( guess::Mission_Guess,
time = Dates.datetime2julian(launch_date)
start = state(current_planet, time, v∞_out, guess.start_mass)
final = zeros(7)
goal = zeros(6)
v∞_in = zeros(3)
# Now, for each phase we must require:
# - That the ending state matches (all legs)
# - 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)
@@ -84,25 +87,28 @@ function solve_mission( guess::Mission_Guess,
start = state(current_planet, time, v∞_out, final[7])
# Do Checks
g[1+8i:6+8i] .= (final[1:6] .- goal[1:6])
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 ) - current_planet.r
else
g[8*(length(flybys)-1)+7] = final[7] - guess.sc.dry_mass
end
g[8*(length(flybys)-1)+8] = max_C3 - norm(x[2:4])^2
g[8*(length(flybys)-1)+9] = max_v∞ - norm(v∞_in)
i += 1
end
g[8*(length(flybys)-1)+1:8*(length(flybys)-1)+3] .= (final[1:3] .- goal[1:3])
g[8*(length(flybys)-1)+4] = max_C3 - norm(x[2:4])^2
g[8*(length(flybys)-1)+5] = max_v∞ - norm(v∞_in)
return -final[end]
if use_cost
return -final[end]
else
return 1.0
end
catch e
if isa(e,LaGuerreConway_Error) || isa(e, Mass_Error)
g[1:8*(length(flybys)-1)+7] .= 1e10
g .= 1e10
return 1e10
else
rethrow()
@@ -114,13 +120,11 @@ function solve_mission( guess::Mission_Guess,
max_time = (datetime2julian(latest_arrival) - datetime2julian(launch_window[1]))*86400
lower_x = lowest_mission_vector(launch_window, length(guess.phases), n)
upper_x = highest_mission_vector(launch_window, max_time, length(guess.phases), n)
num_constraints = 8*(length(guess.phases)-1) + 9
num_constraints = 8*(length(guess.phases)-1) + 5
g_low, g_high = constraint_bounds(guess)
ipopt_options = Dict("constr_viol_tol" => tol,
"acceptable_constr_viol_tol" => 100tol,
"bound_relax_factor" => 0.,
"max_iter" => 100_000,
"max_cpu_time" => 10. * length(guess.phases),
ipopt_options = Dict("tol" => tol,
"acceptable_tol" => 100tol,
"max_cpu_time" => CPU_multiplier * length(guess.phases),
"print_level" => print_level)
options = Options(solver=IPOPT(ipopt_options), derivatives=ForwardFD())

View File

@@ -13,12 +13,13 @@ end
period(b::Body) = 2π * (b.a^3 / Sun.μ)
function state(p::Body, t::DateTime, add_v∞::AbstractVector=[0., 0., 0.], add_mass::Float64=1e10)
time = Dates.datetime2julian(t)
[ basic_ephemeris(p,time); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
time = utc2et(format(t,"yyyy-mm-ddTHH:MM:SS"))
[ spkssb(p.id, time, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
end
function state(p::Body, t::Real, add_v∞::AbstractVector=[0., 0., 0.], add_mass::Real=1e10)
[ basic_ephemeris(p,t); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
time = utc2et(format(julian2datetime(t),"yyyy-mm-ddTHH:MM:SS"))
[ spkssb(p.id, time, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
end
"""

View File

@@ -41,9 +41,9 @@ struct Mission_Creation_Error{T} <: Exception where T <: Real
end
Base.showerror(io::IO, e::Mission_Creation_Error) = print(io, "Tried to create a mission with $(e.num_entries) entries")
struct Planet_Match_Error{T} <: Exception where T <: Real
state::Vector{T}
planet::Vector{T}
struct Planet_Match_Error <: Exception
state::Vector
planet::Vector
end
Base.showerror(io::IO, e::Planet_Match_Error) = print(io, """"
Spacecraft didn't make it to planet

View File

@@ -1,7 +1,7 @@
export Phase, Mission_Guess, Mission, Bad_Mission
export test_phase1, test_phase2
export test_mg
export store
export store, ingest, improve
export Vector
struct Phase
@@ -146,10 +146,11 @@ function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, p
# Perform checks
# Check that the sc actually makes it to the goal state
abs(norm(final[1:3] - goal[1:3])) < 3e6 || throw(Planet_Match_Error(final[1:3], goal[1:3]))
abs(norm(final[4:6] - goal[4:6])) < 3e-4 || throw(Planet_Match_Error(final[4:6], goal[4:6]))
abs(norm(final[1:3] - goal[1:3])) < 1e5 || throw(Planet_Match_Error(final[1:3], goal[1:3]))
if phase != phases[end]
abs(norm(final[4:6] - goal[4:6])) < 1e-1 || throw(Planet_Match_Error(final[4:6], goal[4:6]))
# Check that the v∞s match
abs(norm(phase.v∞_in) - norm(phase.v∞_out)) < 0.01 || throw(V∞_Error(phase.v∞_in, phase.v∞_out))
@@ -260,3 +261,113 @@ function store(ms::Union{Vector{Mission}, Vector{Mission_Guess}}, filename::Abst
end
end
end
"""
Opposite of store. Ingests from a file. Only works if the mission was valid
"""
function ingest(filename::AbstractString)
phases = Vector{Phase}()
planet = Earth
vinf_in = Vector{Float64}()
vinf_out = Vector{Float64}()
tof = 0.0
thrusts = Vector{Float64}()
dry_mass = 0.0
isp = 0.0
max_thrust = 0.0
num_thrusters = 0
duty_cycle = 0.0
launch_mass = 0.0
launch_date = Dates.now()
launch_vinf = Vector{Float64}()
open(filename) do file
phase = 0
while !eof(file)
line = readline(file)
words = split(line, " ")
if strip(words[1]) == "dry_mass:" && phase == 0
dry_mass = parse(Float64,words[end-1])
elseif length(words) >= 2 && strip.(words[1:2]) == ["specific", "impulse:"] && phase == 0
isp = parse(Float64,words[end-1])
elseif strip(words[1]) == "max_thrust:" && phase == 0
max_thrust = parse(Float64,words[end-1])
elseif strip(words[1]) == "num_thrusters:" && phase == 0
num_thrusters = parse(Int,words[end])
elseif strip(words[1]) == "duty_cycle:" && phase == 0
duty_cycle = parse(Float64,words[end])
elseif length(words) >= 2 && strip.(words[1:2]) == ["Launch", "Mass:"] && phase == 0
launch_mass = parse(Float64,words[end-1])
elseif length(words) >= 2 && strip.(words[1:2]) == ["Launch", "Date:"] && phase == 0
launch_date = DateTime(words[end])
elseif length(words) >= 2 && strip.(words[1:2]) == ["Launch", "V∞:"] && phase == 0
launch_vinf = map(x->parse(Float64, strip(x, [ '[', ']', ',' ])), words[end-3:end-1])
elseif strip(words[1]) == "Phase"
if phase > 0
push!(phases, Phase(planet, vinf_in, vinf_out, tof, thrusts))
end
phase = parse(Int,words[2][1:end-1])
elseif strip(words[1]) == "Planet:"
planet_string = words[end]
if planet_string == "Venus"
planet = Venus
elseif planet_string == "Earth"
planet = Earth
elseif planet_string == "Mars"
planet = Mars
elseif planet_string == "Jupiter"
planet = Jupiter
elseif planet_string == "Saturn"
planet = Saturn
end
elseif strip(words[1]) == "V∞_in:"
vinf_in = map(x->parse(Float64, strip(x, [ '[', ']', ',' ])), words[end-3:end-1])
elseif strip(words[1]) == "V∞_out:"
vinf_out = map(x->parse(Float64, strip(x, [ '[', ']', ',' ])), words[end-3:end-1])
elseif length(words) >= 3 && strip.(words[1:3]) == ["time", "of", "flight:"]
tof = parse(Float64,words[end-1])
elseif length(words) >= 2 && strip.(words[1:2]) == ["thrust", "profile:"]
thrusts = map(x->parse(Float64, strip(x, [ '[', ']', ',', ';' ])), words[3:end-1])
thrusts = reshape(thrusts, (3, Int(length(thrusts)/3)))'
end
end
push!(phases, Phase(planet, vinf_in, vinf_out, tof, thrusts))
end
Mission(Sc("ingested", dry_mass, isp, max_thrust, num_thrusters, duty_cycle),
launch_mass,
launch_date,
launch_vinf,
phases)
end
"""
Increases the fidelity of a mission by using it as a base and resolving with increased number of points
"""
function improve(mission::Mission; n::Int=5)
new_phases = Vector{Phase}()
for phase in mission.phases
new_thrusts = repeat(phase.thrust_profile, inner=(n,1))
push!(new_phases, Phase(phase.planet,
phase.v∞_in,
phase.v∞_out,
phase.tof,
new_thrusts))
end
guess = Mission_Guess(mission.sc,
mission.start_mass,
mission.launch_date,
mission.launch_v∞,
new_phases)
arrival = mission.launch_date + Second(floor(sum([p.tof for p in mission.phases])))
new_mission = solve_mission(
guess,
(mission.launch_date - Dates.Day(7), mission.launch_date + Dates.Day(7)),
arrival + Dates.Day(7),
200.,
20.,
CPU_multiplier = 30.
)
new_mission.converged || error("Mission improvement failed")
return new_mission
end

View File

@@ -42,18 +42,62 @@ Base.getindex(x::Nothing, i::Int) = nothing
"""
Generates a layout that works for most plots. Requires a title and a limit.
"""
function standard_layout(limit::Float64, title::AbstractString)
function standard_layout(limit::Float64, title::AbstractString; mode="dark")
limit *= 1.1
Layout(title=attr(font=attr(color="rgb(250,250,250)"), yanchor="top", y=0.95, text=title),
textfont = attr(color="rgb(250,250,250)"),
margin=attr(b=20, t=0, l=0, r=0),
paper_bgcolor="rgba(5,10,40,1.0)",
plot_bgcolor="rgba(200,200,200,0.01)",
scene = attr(xaxis = attr(autorange = false,range=[-limit,limit],color="rgb(255,255,255)"),
yaxis = attr(autorange = false,range=[-limit,limit],color="rgb(255,255,255)"),
zaxis = attr(autorange = false,range=[-limit,limit],color="rgb(255,255,255)"),
aspectratio=attr(x=1,y=1,z=1),
aspectmode="manual"))
if mode == "dark"
Layout(
title=attr(
font_color="rgb(250,250,250)",
yanchor="top",
y=0.95,
text=title
),
textfont_color = "rgb(250,250,250)",
margin=attr(b=20, t=0, l=0, r=0),
paper_bgcolor="rgba(5,10,40,1.0)",
plot_bgcolor="rgba(200,200,200,0.01)",
scene_xaxis = attr(
autorange = false,
range=[-limit,limit],
color="rgb(255,255,255)"
),
scene_yaxis = attr(
autorange = false,
range=[-limit,limit],
color="rgb(255,255,255)"
),
scene_zaxis = attr(
autorange = false,
range=[-2e8,2e8],
color="rgb(255,255,255)"
),
scene_aspectratio=attr(x=1,y=1,z=1),
scene_aspectmode="manual"
)
elseif mode == "light"
Layout(
title=attr(
font_color="rgb(0,0,0)",
yanchor="top",
y=0.95,
text=title
),
textfont_color = "rgb(0,0,0)",
margin=attr(b=20, t=0, l=0, r=0),
paper_bgcolor="rgba(255,255,255,0.0)",
plot_bgcolor="rgba(255,255,255,0.0)",
scene_aspectmode = "data",
scene_xaxis = attr(
autorange = true,
color="rgb(0,0,0)"
),
scene_yaxis = attr(
autorange = true,
color="rgb(0,0,0)"
),
scene_zaxis_visible = false
)
end
end
"""
@@ -94,13 +138,14 @@ A fundamental function to plot a single orbit path
function plot(path::Vector{Vector{Float64}};
color::AbstractString = random_color(),
title::AbstractString = "Single Orbit Path",
markers=true)
markers=true,
mode="dark")
# Generate the plot
p, limit = gen_plot(path, color=color, markers=markers)
# Now use the standard layout
layout = standard_layout(limit, title)
layout = standard_layout(limit, title, mode=mode)
# Generate the plot
PlotlyJS.plot( PlotlyJS.Plot( p, layout ) )
@@ -118,7 +163,11 @@ An optional prop time is allowed, which will also provide an orbit
function gen_plot(b::Body,
time::Union{DateTime, Nothing}=nothing,
prop_time::Union{Float64, Nothing}=nothing;
scale::Int=100)
scale::Int=100,
color=nothing)
if color==nothing
color = b.line_color
end
# Gotta be a little weird with Sun Scale
if b == Sun && scale > 25 scale = 25 end
@@ -151,7 +200,7 @@ function gen_plot(b::Body,
outputs = [ body_surface ]
if prop_time !== nothing
p, maybe_max = gen_plot(prop(body_state, prop_time), label=b.name, color=b.line_color, markers=false)
p, maybe_max = gen_plot(prop(body_state, prop_time), label=b.name, color=color, markers=false)
push!(outputs, p...)
if maybe_max > max_value max_value = maybe_max end
end
@@ -172,13 +221,14 @@ function plot(b::Body,
time::Union{DateTime, Nothing}=nothing,
prop_time::Union{Float64, Nothing}=nothing;
title::AbstractString="Single Planet Plot",
scale::Int=100)
scale::Int=100,
mode="dark")
# Generate the plot
p, limit = gen_plot(b,time,prop_time,scale=scale)
# Now use the standard layout
layout = standard_layout(limit, title)
layout = standard_layout(limit, title, mode=mode)
# Generate the plot
PlotlyJS.plot( PlotlyJS.Plot( p, layout ) )
@@ -196,7 +246,8 @@ function plot(bodies::Vector{Body},
time::Union{DateTime, Nothing};
orbits::Bool=true,
title::AbstractString="Multiple Planet Plot",
scale::Int=100)
scale::Int=100,
mode="dark")
# Generate the plots
trace, limit = gen_plot(Sun, scale=scale)
@@ -212,7 +263,7 @@ function plot(bodies::Vector{Body},
end
# Now use the standard layout
layout = standard_layout(limit, title)
layout = standard_layout(limit, title, mode=mode)
# Generate the plot
PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
@@ -226,7 +277,8 @@ function plot(paths::Vector{Matrix{Real}};
labels::Union{Vector{String},Nothing}=nothing,
colors::Union{Vector{String},Nothing}=nothing,
markers::Bool=true,
title::String="Spacecraft Position")
title::String="Spacecraft Position",
mode="dark")
# Generate the path traces
color = colors === nothing ? random_color() : colors[1]
@@ -244,20 +296,33 @@ function plot(paths::Vector{Matrix{Real}};
# Determine layout details
# layout = standard_layout(limit, title)
layout = standard_layout(1e9, title)
layout = standard_layout(1e9, title, mode=mode)
# Plot
PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
end
function plot(m::Union{Mission, Mission_Guess}; title::String="Mision Plot")
function plot(m::Union{Mission, Mission_Guess};
title::String="Mission Plot",
mode::String="dark",
planet_colors=true,
phase_colors=nothing,
markers=true)
# First plot the earth
# Then plot phase plus planet phase
time = datetime2julian(m.launch_date)
mass = m.start_mass
current_planet = Earth
earth_trace, limit = gen_plot(current_planet, julian2datetime(time), year)
if planet_colors
earth_trace, limit = gen_plot(current_planet, julian2datetime(time), year)
else
if mode=="dark"
earth_trace, limit = gen_plot(current_planet, julian2datetime(time), year, color="#FFF")
elseif mode=="light"
earth_trace, limit = gen_plot(current_planet, julian2datetime(time), year, color="#000")
end
end
start = state(current_planet, time, m.launch_v∞, mass)
traces = [ earth_trace... ]
@@ -270,12 +335,35 @@ function plot(m::Union{Mission, Mission_Guess}; title::String="Mision Plot")
time += phase.tof/86400
current_planet = phase.planet
start = state(current_planet, time, phase.v∞_out, mass)
path_trace, new_limit = gen_plot(path, label="Phase "*string(i))
if phase_colors == nothing
path_trace, new_limit = gen_plot(path,
label="Phase "*string(i),
markers=markers)
else
path_trace, new_limit = gen_plot(path,
label="Phase "*string(i),
markers=markers,
color=phase_colors[i])
end
push!(traces, path_trace...)
limit = max(limit, new_limit)
# Then the planet
planet_trace, new_limit = gen_plot(current_planet, julian2datetime(time), -phase.tof)
if planet_colors
planet_trace, new_limit = gen_plot(current_planet, julian2datetime(time), -phase.tof)
else
if mode=="dark"
planet_trace, new_limit = gen_plot(current_planet,
julian2datetime(time),
-phase.tof,
color="#FFF")
elseif mode=="light"
planet_trace, new_limit = gen_plot(current_planet,
julian2datetime(time),
-phase.tof,
color="#000")
end
end
push!(traces, planet_trace...)
limit = max(limit, new_limit)
@@ -286,7 +374,7 @@ function plot(m::Union{Mission, Mission_Guess}; title::String="Mision Plot")
push!(traces, gen_plot(Sun)[1]...)
# Determine layout details
layout = standard_layout(limit, title)
layout = standard_layout(limit, title, mode=mode)
# Plot
PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )

View File

@@ -22,11 +22,12 @@ function prop(ΔVs::AbstractArray,
size(ΔVs)[2] == 3 || throw(ΔVsize_Error(size(ΔVs)))
n = size(ΔVs)[1]
output = Array{Real,2}(undef,n,7)
output = Array{Real,2}(undef,n+1,7)
current = copy(state)
full_thrust = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time/n, 0.)
for i in 1:n
output[i,:] = current
halfway = laguerre_conway(current[1:6], time/(2n), primary)
if time > 0
halfway += [zeros(3); full_thrust * ΔVs[i,:]] / current[7]
@@ -36,8 +37,8 @@ function prop(ΔVs::AbstractArray,
halfway += [zeros(3); full_thrust * ΔVs[i,:]] / current[7]
end
current[1:6] = laguerre_conway(halfway, time/(2n), primary)
output[i,:] = current
end
output[end,:] = current
current[7] > craft.dry_mass || throw(Mass_Error(current[7]))
@@ -48,7 +49,7 @@ end
"""
Convenience function for propagating a state with no thrust
"""
prop(x::AbstractArray, t::Real, p::Body=Sun) = prop(zeros(1000,3), x, no_thrust, t, p)
prop(x::AbstractArray, t::Real, p::Body=Sun) = prop(zeros(1000,3), x, no_thrust, t, primary=p)
"""
This is solely for the purposes of getting the final state of a mission or guess