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

Succesful MBH runs, with some improvements

Browse Source
This commit is contained in:
Connor
2021-10-12 21:49:46 -06:00
parent c7913d2944
commit bb78578d9a
19 changed files with 426 additions and 146 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
julia/src/errors.jl
View File

@@ -6,8 +6,10 @@ struct ΔVsize_Error <: Exception
end
Base.showerror(io::IO, e::ΔVsize_Error) = print(io, "ΔV was too big: $(e.ΔV)")
struct HitPlanet_Error <: Exception end
Base.showerror(io::IO, e::HitPlanet_Error) = print(io, "spacecraft hit the planet...")
struct HitPlanet_Error <: Exception
peri::Float64
end
Base.showerror(io::IO, e::HitPlanet_Error) = print(io, "spacecraft hit the planet..., periapsis: $(e.peri)")
struct Convergence_Error <: Exception end
Base.showerror(io::IO, e::Convergence_Error) = print(io, "NLP solver didn't converge...")

31
julia/src/inner_loop/monotonic_basin_hopping.jl
View File

@@ -51,10 +51,9 @@ function mission_guess( flybys::Vector{Body},
sc::Sc,
start_mass::Float64,
launch_window::Tuple{DateTime, DateTime},
max_C3_out::Float64,
max_C3::Float64,
max_v∞_in_mag::Float64,
latest_arrival::DateTime,
primary::Body=Sun )
latest_arrival::DateTime)
mission_guess = Bad_Mission("Keep trying to generate a guess")
while mission_guess == Bad_Mission("Keep trying to generate a guess")
# TODO: Eventually I can calculate n more intelligently
@@ -62,8 +61,9 @@ function mission_guess( flybys::Vector{Body},
# Determine the launch conditions
launch_date = rand(launch_window...)
launch_v∞_normalized = rand(-1:0.0001:1, 3)
launch_v∞ = rand(0:0.0001:√max_C3_out) * launch_v∞_normalized/norm(launch_v∞_normalized)
launch_v∞_unit = rand(-1:0.0001:1, 3)
launch_v∞_normalized = launch_v∞_unit/norm(launch_v∞_unit)
launch_v∞ = √max_C3 * rand() * launch_v∞_normalized
# Determine the leg lengths
num_phases = length(flybys) - 1
@@ -126,14 +126,14 @@ function mbh( flybys::Vector{Body},
max_C3::Float64,
max_v∞::Float64,
latest_arrival::DateTime,
cost_fn::Function,
primary::Body=Sun;
cost_fn::Function;
search_patience::Int=10_000,
drill_patience::Int=50,
verbose::Bool=false)
verbose::Bool=false,
test::Bool=false )
# Convenience Functions
random_guess() = mission_guess(flybys,sc,start_mass,launch_window,max_C3,max_v∞,latest_arrival,primary)
random_guess() = mission_guess(flybys,sc,start_mass,launch_window,max_C3,max_v∞,latest_arrival)
solve(g::Mission_Guess) = solve_mission(g, launch_window, latest_arrival)
cost(m::Mission) = cost_fn(m, max_C3, max_v∞)
cost(_::Nothing) = Inf
@@ -143,13 +143,16 @@ function mbh( flybys::Vector{Body},
drill_count = 0
archive = Vector{Mission}()
function log()
sct = search_count
dct = drill_count
lar = length(archive)
flyby_abbreviation = join([ f.name[1] for f in flybys ])
print("\r ")
print("\rsearch: $(search_count) drill: $(drill_count) archive: $(length(archive))\t")
print("\rMBH\t$(flyby_abbreviation)\tsearch: $(sct) drill: $(dct) archive: $(lar)\t")
end
# The main loop
x_current = nothing
if verbose println("Starting Monotonic Basin Hopper") end
while search_count < search_patience
# Intialize an x_star, if it doesn't converge, hop on to the next basin
@@ -188,6 +191,12 @@ function mbh( flybys::Vector{Body},
x_current in archive || push!(archive, x_current)
# If in test mode, we don't need to actually optimize. Just grab the first valid basin-best
if test
println()
return x_current, [ x_current ]
end
end
if verbose println() end

59
julia/src/inner_loop/nlp_solver.jl
View File

@@ -2,23 +2,27 @@ using SNOW
export solve_mission
struct Result
converged::Bool
info::Symbol
sol::Mission_Guess
end
const pos_scale = ones(3)
# const vel_scale = 100. * ones(3)
const vel_scale = ones(3)
# const v∞_scale = norm(vel_scale)
const v∞_scale = 1.
const state_scale = [pos_scale; vel_scale ]
function constraint_bounds(guess::Mission_Guess)
low_constraint = Vector{Float64}()
high_constraint = Vector{Float64}()
for phase in guess.phases
state_constraint = [100., 100., 100., 0.005, 0.005, 0.005] .* state_scale
push!(low_constraint, (-1 * state_constraint)... )
push!(high_constraint, state_constraint... )
if phase != guess.phases[end]
push!(low_constraint, -100., -100., -100., -0.005, -0.005, -0.005, -30., 100.)
push!(high_constraint, 100., 100., 100., 0.005, 0.005, 0.005, 30., 100_000.)
push!(low_constraint, -0.005*v∞_scale, 100.)
push!(high_constraint, 0.005*v∞_scale, 1e7)
else
push!(low_constraint, -100., -100., -100., 0.)
push!(high_constraint, 100., 100., 100., guess.start_mass - guess.sc.dry_mass)
push!(low_constraint, 0.)
push!(high_constraint, guess.start_mass - guess.sc.dry_mass)
end
end
return low_constraint, high_constraint
@@ -33,7 +37,7 @@ guess otherwise
function solve_mission( guess::Mission_Guess,
launch_window::Tuple{DateTime,DateTime},
latest_arrival::DateTime;
tol=1e-10,
tol=1e-12,
verbose::Bool=false,
print_level=0 )
@@ -45,64 +49,71 @@ function solve_mission( guess::Mission_Guess,
# 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)
start = state(current_planet, time, v∞_out, guess.start_mass)
# Now, for each phase we must require:
# - That the ending state matches (unless final leg)
# - That the ending state matches (all legs)
# - 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
try
for flyby in flybys
# Get the values from the vector
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))
# Propagate
final = prop(thrusts, start, guess.sc, tof)[2]
current_planet = flyby
time += tof
goal = state(current_planet, time, v∞_in)
final = prop(reshape(thrusts, (n,3)), start, guess.sc, tof)[2]
goal = state(current_planet, time, v∞_in)[1:6]
start = state(current_planet, time, v∞_out, final[7])
# Do Checks
g[1+8i:6+8i] .= (final[1:6] .- goal[1:6]) .* state_scale
if flyby != flybys[end]
g[1+8i:6+8i] .= (final[1:6] .- goal[1:6]) .* [1., 1., 1., 10_000., 10_000., 10_000.]
g[7+8i] = (norm(v∞_out) - norm(v∞_in)) * 10_000.
g[7+8i] = (norm(v∞_out) - norm(v∞_in)) * v∞_scale
δ = 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)+1:8*(length(flybys)-1)+3] .= final[1:3] .- goal[1:3]
g[8*(length(flybys)-1)+4] = final[7] - guess.sc.dry_mass
g[8*(length(flybys)-1)+7] = final[7] - guess.sc.dry_mass
end
start = state(current_planet, time, v∞_out, final[7])
mass = final[7]
i += 1
end
return 1.0
catch e
if isa(e,LaGuerreConway_Error)
g[1:8*(length(flybys)-1)+4] .= 1e10
g[1:8*(length(flybys)-1)+7] .= 1e10
return 1e10
else
rethrow()
end
end
end
max_time = Dates.datetime2unix(latest_arrival) - Dates.datetime2unix(launch_window[1])
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) + 4
num_constraints = 8*(length(guess.phases)-1) + 7
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" => 60.,
"max_cpu_time" => 5. * length(guess.phases),
"print_level" => print_level)
options = Options(solver=IPOPT(ipopt_options), derivatives=ForwardFD())

2
julia/src/inner_loop/propagator.jl
View File

@@ -24,7 +24,7 @@ function prop_one(ΔV_unit::Vector{<:Real},
ΔV = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time, 0., state[7]) * ΔV_unit
halfway = laguerre_conway(state, time/2, primary) + [zeros(3); ΔV]
final = laguerre_conway(halfway, time/2, primary)
return [final; state[7] - craft.mass_flow_rate*norm(ΔV_unit)*time]
return [final; state[7] - mfr(craft)*norm(ΔV_unit)*time]
end

96
julia/src/mission.jl
View File

@@ -25,7 +25,7 @@ function mass_consumption(sc::Sc, phase::Phase)
for i in 1:size(phase.thrust_profile,1)
weighted_thrusting_time += norm(phase.thrust_profile[i,:]) * phase.tof/n
end
return weighted_thrusting_time*sc.mass_flow_rate
return weighted_thrusting_time*mfr(sc)
end
struct Mission_Guess
@@ -133,45 +133,82 @@ function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, p
force=false)
# First do some checks to make sure that it's valid
if !force
time = date
time = utc2et(Dates.format(date,"yyyy-mm-ddTHH:MM:SS"))
current_planet = Earth
start = state(current_planet, time, v∞, mass)
for phase in phases
#Propagate
final = prop(phase.thrust_profile, start, sc, phase.tof)[2]
mass = final[7]
mass > sc.dry_mass || throw(Mass_Error(mass - mass_used))
current_planet = phase.planet
time += Dates.Second(floor(phase.tof))
start = state(current_planet, time, phase.v∞_out, mass)
p_pos = start[1:3]
abs(norm(final[1:3] - p_pos)) < √30_000. || throw(Planet_Match_Error(final[1:3], p_pos))
v∞_in, v∞_out = phase.v∞_in, phase.v∞_out
time += phase.tof
goal = state(current_planet, time, phase.v∞_in)
start = state(current_planet, time, phase.v∞_out, final[7])
# 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]))
if phase != phases[end]
abs(norm(v∞_in) - norm(v∞_out)) < 0.005 || throw(V∞_Error(v∞_in, v∞_out))
δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
periapsis = (phase.planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
periapsis > 1.1phase.planet.r || throw(HitPlanet_Error())
# 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))
# Check that the sc doesn't hit the planet
δ = acos( ( phase.v∞_in ⋅ phase.v∞_out ) / ( norm(phase.v∞_in) * norm(phase.v∞_out) ) )
periapsis = (current_planet.μ/(phase.v∞_in ⋅ phase.v∞_in)) * ( 1/sin(δ/2) - 1 ) - current_planet.r
periapsis > 90. || throw(HitPlanet_Error(periapsis))
# Check that the spacecraft never thrusted more than 100%
for ΔV in phase.thrust_profile
abs(ΔV) <= 1.0 || throw(ΔVsize_Error(ΔV))
end
end
end
end
Mission(sc, mass, date, v∞, phases, true)
m = Mission(sc, mass, date, v∞, phases, true)
# Check that the fuel never runs out
mass = prop(m)[7]
mass > sc.dry_mass || throw(Mass_Error(mass))
return m
end
function Mission(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(sc, mass, launch_date, launch_v∞, phases)
end
"""
BE CAREFUL!! This just makes a guess converged, whether true or not
"""
function Mission(g::Mission_Guess)
println("Skipping checks...")
return Mission(g.sc, g.start_mass, g.launch_date, g.launch_v∞, g.phases, force=true)
end
function Mission(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body})
guess = Mission_Guess(x, sc, mass, flybys::Vector{Body})
return Mission(guess)
end
struct Bad_Mission
message::String
converged::Bool
@@ -181,23 +218,26 @@ Bad_Mission(s::Symbol) = Bad_Mission(String(s),false)
function Base.write(io::IO, m::Mission)
write(io, m.sc)
write(io, "Launch Mass: $(m.start_mass)\n")
write(io, "Launch Mass: $(m.start_mass) kg\n")
write(io, "Launch Date: $(m.launch_date)\n")
write(io, "Launch V∞: $(m.launch_v∞)\n")
write(io, "Launch V∞: $(m.launch_v∞) km/s\n")
time = m.launch_date
i = 1
for phase in m.phases
time += Dates.Second(floor(phase.tof))
write(io, "Phase $(i):\n")
write(io, "\tPlanet: $(phase.planet.name)\n")
write(io, "\tV∞_in: $(phase.v∞_in)\n")
write(io, "\tV∞_out: $(phase.v∞_out)\n")
write(io, "\ttime of flight: $(phase.tof)\n")
write(io, "\tthrust profile: $(phase.thrust_profile)\n")
write(io, "\tV∞_in: $(phase.v∞_in) km/s\n")
write(io, "\tV∞_out: $(phase.v∞_out) km/s\n")
write(io, "\ttime of flight: $(phase.tof) seconds\n")
write(io, "\tarrival date: $(time)\n")
write(io, "\tthrust profile: $(phase.thrust_profile) %\n")
i += 1
end
write(io, "\n")
write(io, "Mass Used: $(m.start_mass - prop(m)[7])\n")
write(io, "Launch C3: $(m.launch_v∞ ⋅ m.launch_v∞)\n")
write(io, "||V∞_in||: $(norm(m.phases[end].v∞_in))\n")
write(io, "Mass Used: $(m.start_mass - prop(m)[7]) kg\n")
write(io, "Launch C3: $(m.launch_v∞ ⋅ m.launch_v∞) km²/s²\n")
write(io, "||V∞_in||: $(norm(m.phases[end].v∞_in)) km/s\n")
end
function store(m::Union{Mission, Mission_Guess}, filename::AbstractString)

18
julia/src/spacecraft.jl
View File

@@ -1,24 +1,26 @@
export Sc, test_sc, bepi, no_thrust
export Sc, test_sc, bepi, no_thrust, mfr
mutable struct Sc
name::AbstractString
dry_mass::Float64
mass_flow_rate::Float64
isp::Float64
max_thrust::Float64
num_thrusters::Int
duty_cycle::Float64
end
mfr(sc::Sc) = sc.num_thrusters * sc.max_thrust / (sc.isp*0.00981)
function Base.write(io::IO, sc::Sc)
write(io, "Spacecraft: $(sc.name)\n")
write(io, "\tdry_mass: $(sc.dry_mass)\n")
write(io, "\tmass_flow_rate: $(sc.mass_flow_rate)\n")
write(io, "\tmax_thrust: $(sc.max_thrust)\n")
write(io, "\tdry_mass: $(sc.dry_mass) kg\n")
write(io, "\tspecific impulse: $(sc.isp) kg/s\n")
write(io, "\tmax_thrust: $(sc.max_thrust) kN\n")
write(io, "\tnum_thrusters: $(sc.num_thrusters)\n")
write(io, "\tduty_cycle: $(sc.duty_cycle)\n")
end
const test_sc = Sc("test", 8000., 0.00025/(2000*0.00981), 0.00025, 50, 0.9)
const bepi = Sc("bepi", 2000., 2*0.00025/(2800*0.00981), 0.00025, 2, 0.9)
const no_thrust = Sc("no thrust", 0., 0.01, 0., 0, 0.)
const test_sc = Sc("test", 8000., 2000., 0.0005, 50, 0.9)
const bepi = Sc("bepi", 2000., 2800., 0.0005, 2, 0.9)
const no_thrust = Sc("no thrust", 10., 2000., 0., 0, 0.)