diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 61997c1..2885926 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -14,6 +14,9 @@ unit-test-job: paths: - julia/plots/nlp_test_1_phase.html - julia/plots/nlp_test_2_phase.html + - julia/plots/nlp_test_5_phase.html + - julia/plots/mbh_test_1_phase.html + - julia/plots/mbh_test_2_phase.html - julia/plots/plot_test_mission.html - julia/plots/plot_test_mission_guess.html - julia/plots/plot_test_path.html diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl index 631ef88..2796109 100644 --- a/julia/src/Thesis.jl +++ b/julia/src/Thesis.jl @@ -20,12 +20,12 @@ module Thesis include("./constants.jl") include("./spacecraft.jl") include("./inner_loop/laguerre-conway.jl") - include("./inner_loop/propagator.jl") include("./mission.jl") + include("./inner_loop/propagator.jl") include("./conversions.jl") include("./plotting.jl") include("./inner_loop/nlp_solver.jl") - # include("./inner_loop/monotonic_basin_hopping.jl") + include("./inner_loop/monotonic_basin_hopping.jl") # include("./outer_loop.jl") include("./lamberts.jl") diff --git a/julia/src/errors.jl b/julia/src/errors.jl index ea80e8d..ed60290 100644 --- a/julia/src/errors.jl +++ b/julia/src/errors.jl @@ -1,7 +1,10 @@ struct LaGuerreConway_Error <: Exception end Base.showerror(io::IO, e::LaGuerreConway_Error) = print(io, "LaGuerre-Conway didn't converge") -struct ΔVsize_Error <: Exception end +struct ΔVsize_Error <: Exception + ΔV::Float64 +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...") diff --git a/julia/src/inner_loop/monotonic_basin_hopping.jl b/julia/src/inner_loop/monotonic_basin_hopping.jl index 0eb4478..eb1a694 100644 --- a/julia/src/inner_loop/monotonic_basin_hopping.jl +++ b/julia/src/inner_loop/monotonic_basin_hopping.jl @@ -5,41 +5,22 @@ Generates pareto-distributed random numbers This is usually around one to two percent, but sometimes up to ten or so (fat tails) """ -function pareto(shape::Tuple{Int, Int}, α::Float64=1.01) - s = rand((-1,1), shape) - r = rand(Float64, shape) +function pareto(shape::Union{Tuple{Int, Int}, Int, Nothing}=nothing; α::Float64=1.01) + if shape !== nothing + s = rand((-1,1), shape) + r = rand(Float64, shape) + else + s = rand((-1,1)) + r = rand(Float64) + end ϵ = 1e-10 return 1 .+ (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) end -function pareto(shape::Int, α::Float64=1.01) - s = rand((-1,1), shape) - r = rand(Float64, shape) - ϵ = 1e-10 - return 1 .+ (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) -end - -function pareto(α::Float64=1.01) - s = rand((-1,1)) - r = rand(Float64) - ϵ = 1e-10 - return 1 + (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) -end - -function pareto_add(α::Float64=1.01) - s = rand((-1,1)) - r = rand(Float64) - ϵ = 1e-10 - return (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) -end - """ Returns a random date between two dates """ -function gen_date(date_range::Tuple{DateTime, DateTime}) - l0, lf = date_range - l0 + Millisecond(floor(rand()*(lf-l0).value)) -end +Base.rand(l0::DateTime,lf::DateTime) = l0 + Millisecond(floor(rand()*(lf-l0).value)) """ Perturbs a valid mission with pareto-distributed variables, generating a mission guess @@ -47,7 +28,7 @@ Perturbs a valid mission with pareto-distributed variables, generating a mission function perturb(mission::Mission) mission_guess = Bad_Mission("Starting point") while typeof(mission_guess) == Bad_Mission - new_launch_date = mission.launch_date + Dates.Second(floor(7day * pareto_add())) + new_launch_date = mission.launch_date + Dates.Second(floor(7day * (pareto()-1))) new_launch_v∞ = mission.launch_v∞ .* pareto(3) new_phases = Vector{Phase}() for phase in mission.phases @@ -58,16 +39,7 @@ function perturb(mission::Mission) new_thrust_profile = phase.thrust_profile .* pareto(size(phase.thrust_profile)) push!(new_phases, Phase(phase.planet, new_v∞_in, new_v∞_out, new_tof, new_thrust_profile)) end - # TODO: Mission_Guess.validate() - try - mission_guess = Mission_Guess(mission.sc, mission.start_mass, new_launch_date, new_launch_v∞, new_phases) - catch e - if isa(e, Mass_Error) ||isa(e,V∞_Error) || isa(e,HitPlanet_Error) - continue - else - rethrow() - end - end + mission_guess = Mission_Guess(mission.sc, mission.start_mass, new_launch_date, new_launch_v∞, new_phases) end return mission_guess end @@ -89,7 +61,7 @@ function mission_guess( flybys::Vector{Body}, n = 20 # Determine the launch conditions - launch_date = gen_date(launch_window) + 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) @@ -142,52 +114,6 @@ function mission_guess( flybys::Vector{Body}, return mission_guess end -""" -Sequentially calls the NLP solver to attempt to solve based on the initial guess -""" -function inner_loop_solve(guess::Mission_Guess) - v∞_out = guess.launch_v∞ - mass = guess.start_mass - current_planet = Earth - time = utc2et(Dates.format(guess.launch_date,"yyyy-mm-ddTHH:MM:SS")) - start = state(current_planet, time, v∞_out, mass) - - corrected_phases = Vector{Phase}() - for i in 1:length(guess.phases) - phase = guess.phases[i] - current_planet = phase.planet - time += phase.tof - goal = state(current_planet, time, phase.v∞_in) - result = solve_phase( start, goal, guess.sc, phase.tof, phase.thrust_profile) - result.converged || return Bad_Mission(result.info) # Drop if it's not working - corrected_phase = Phase(phase.planet, phase.v∞_in, phase.v∞_out, phase.tof, result.sol) - push!(corrected_phases, corrected_phase) - mass_used = mass_consumption(guess.sc, corrected_phase) - mass -= mass_used - if i != length(guess.phases) - v∞_out = phase.v∞_out - start = state(current_planet, time, v∞_out, mass) - end - end - - return Mission(guess.sc, guess.start_mass, guess.launch_date, guess.launch_v∞, corrected_phases) -end - -""" -The cost function for the mission -TODO: This will probably move and eventually be passed as an argument -""" -function cost(mission::Mission, max_C3::Float64, max_v∞::Float64) - mass_used = 0.0 - for phase in mission.phases mass_used += mass_used(sc, phase) end - mass_percent = mass_used/mission.sc.dry_mass - C3_percent = ( mission.launch_v∞ ⋅ mission.launch_v∞ ) / max_C3 - v∞_percent = norm(mission.phases[end].v∞_in) / max_v∞ - return 3mass_percent + C3_percent + v∞_percent -end - -cost(_::Bad_Mission) = Inf - """ This is the main monotonic basin hopping function. There's a lot going on here, but the general idea is that hopefully you can provide mission parameters and a list of flybys and get the optimal @@ -200,44 +126,50 @@ function mbh( flybys::Vector{Body}, max_C3::Float64, max_v∞::Float64, latest_arrival::DateTime, + cost_fn::Function, primary::Body=Sun; search_patience::Int=10_000, drill_patience::Int=50, verbose::Bool=false) - # A convenience function + # Convenience Functions random_guess() = mission_guess(flybys,sc,start_mass,launch_window,max_C3,max_v∞,latest_arrival,primary) - cost(m::Mission) = cost(m, max_C3, max_v∞) - status(s,d,l) = print("\r\t", "search: ", s, " drill: ", d, " archive length: ", l, " ") + solve(g::Mission_Guess) = solve_mission(g, launch_window, latest_arrival) + cost(m::Mission) = cost_fn(m, max_C3, max_v∞) + cost(_::Nothing) = Inf # Initialize stuff search_count = 0 - x_current = Bad_Mission("Starting point") + drill_count = 0 archive = Vector{Mission}() + function log() + print("\r ") + print("\rsearch: $(search_count) drill: $(drill_count) archive: $(length(archive))\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 search_count += 1 drill_count = 0 - if verbose status(search_count, drill_count, length(archive)) end - x_star = inner_loop_solve(random_guess()) + verbose && log() + x_star = solve(random_guess()) x_star.converged || continue x_basin = x_star - # If it does, though, we check to see if it's better than the current if cost(x_star) < cost(x_current) x_current = x_star end # Either way, we need to drill into this particular basin, since it's valid while drill_count < drill_patience - if verbose status(search_count, drill_count, length(archive)) end + verbose && log() #Perturb to generate a slightly different x_star - x_star = inner_loop_solve(perturb(x_basin)) + 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) @@ -258,6 +190,8 @@ function mbh( flybys::Vector{Body}, end + if verbose println() end + return x_current, archive end diff --git a/julia/src/inner_loop/nlp_solver.jl b/julia/src/inner_loop/nlp_solver.jl index 82de693..77be0cd 100644 --- a/julia/src/inner_loop/nlp_solver.jl +++ b/julia/src/inner_loop/nlp_solver.jl @@ -31,7 +31,7 @@ NOTE: This function will output a proper mission if it succeeded and just re-ret guess otherwise """ function solve_mission( guess::Mission_Guess, - launch_window::Vector{DateTime}, + launch_window::Tuple{DateTime,DateTime}, latest_arrival::DateTime; tol=1e-10, verbose::Bool=false, @@ -114,7 +114,7 @@ function solve_mission( guess::Mission_Guess, else g = zeros(num_constraints) optimizer!(g,x) - println(g) + if verbose println(g) end return guess end diff --git a/julia/src/inner_loop/propagator.jl b/julia/src/inner_loop/propagator.jl index f179b72..b121cae 100644 --- a/julia/src/inner_loop/propagator.jl +++ b/julia/src/inner_loop/propagator.jl @@ -21,9 +21,6 @@ function prop_one(ΔV_unit::Vector{<:Real}, time::Float64, primary::Body=Sun) - for direction in ΔV_unit - abs(direction) <= 1.0 || throw(PropOne_Error(ΔV_unit)) - end Δ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) @@ -81,3 +78,24 @@ end Convenience function for propagating a state with no thrust """ prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1.], no_thrust, t, p)[1] + +""" +This is solely for the purposes of getting the final state of a mission or guess +""" +function prop(m::Mission) + time = m.launch_date + current_planet = Earth + start = state(current_planet, time, m.launch_v∞, m.start_mass) + + final = zeros(7) + for phase in m.phases + final = prop(phase.thrust_profile, start, m.sc, phase.tof)[2] + mass = final[7] + current_planet = phase.planet + time += Dates.Second(floor(phase.tof)) + start = state(current_planet, time, phase.v∞_out, mass) + end + + return final + +end diff --git a/julia/src/mission.jl b/julia/src/mission.jl index ee607c0..3aeee51 100644 --- a/julia/src/mission.jl +++ b/julia/src/mission.jl @@ -40,13 +40,6 @@ end Constructor for a mission guess. Generally mission guesses are not converged """ function Mission_Guess(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, phases::Vector{Phase}) - # First do some checks to make sure that it's valid - mass_used = 0 - for phase in phases - mass_used += mass_consumption(sc, phase) - mass - mass_used > sc.dry_mass || throw(Mass_Error(mass - mass_used)) - v∞_in, v∞_out = phase.v∞_in, phase.v∞_out - end Mission_Guess(sc, mass, date, v∞, phases, false) end @@ -95,7 +88,7 @@ function Base.Vector(g::Mission_Guess) return result end -function lowest_mission_vector(launch_window::Vector{DateTime}, num_phases::Int, n::Int) +function lowest_mission_vector(launch_window::Tuple{DateTime,DateTime}, num_phases::Int, n::Int) result = Vector{Float64}() push!(result, Dates.datetime2unix(launch_window[1])) push!(result, -10*ones(3)...) @@ -108,7 +101,7 @@ function lowest_mission_vector(launch_window::Vector{DateTime}, num_phases::Int, return result end -function highest_mission_vector(launch_window::Vector{DateTime}, mission_length::Float64, num_phases::Int, n::Int) +function highest_mission_vector(launch_window::Tuple{DateTime,DateTime}, mission_length::Float64, num_phases::Int, n::Int) result = Vector{Float64}() push!(result, Dates.datetime2unix(launch_window[2])) push!(result, 10*ones(3)...) @@ -148,7 +141,8 @@ function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, p mass > sc.dry_mass || throw(Mass_Error(mass - mass_used)) current_planet = phase.planet time += Dates.Second(floor(phase.tof)) - p_pos = state(current_planet, time)[1:3] + 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 if phase != phases[end] @@ -156,6 +150,9 @@ function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, p δ = 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()) + for ΔV in phase.thrust_profile + abs(ΔV) <= 1.0 || throw(ΔVsize_Error(ΔV)) + end end end end diff --git a/julia/test/inner_loop/monotonic_basin_hopping.jl b/julia/test/inner_loop/monotonic_basin_hopping.jl index c1adbd5..55ebd4a 100644 --- a/julia/test/inner_loop/monotonic_basin_hopping.jl +++ b/julia/test/inner_loop/monotonic_basin_hopping.jl @@ -1,58 +1,45 @@ @testset "Monotonic Basin Hopping" begin - using PlotlyJS + using PlotlyJS: savefig println("Testing Monotonic Basin Hopper") - # First we test the random mission guess generator - println("Testing guess generator function") - flybys = [Earth, Venus, Jupiter] - launch_window = ( DateTime(2021,12,25), DateTime(2025,12,25) ) - max_C3 = 10. - max_v∞ = 8. - latest_arrival = DateTime(2030,12,25) - random_guess = Thesis.mission_guess(flybys, bepi, 12_000., launch_window, max_C3, max_v∞, latest_arrival) - @test typeof(random_guess) == Mission_Guess - - # Then the perturb function - println("Testing perturb function") - mission_guess = Thesis.perturb(test_mission) - @test mission_guess.launch_date != test_mission.launch_date - @test mission_guess.launch_v∞ != test_mission.launch_v∞ - for i in 1:2 - @test mission_guess.phases[i].v∞_in != test_mission.phases[i].v∞_in - @test mission_guess.phases[i].v∞_out != test_mission.phases[i].v∞_out - @test mission_guess.phases[i].tof != test_mission.phases[i].tof + """ + The cost function for the mission + """ + function easy_cost(m::Mission, 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∞ end - @test !mission_guess.converged + + # Mission Parameters that won't change (they're very lenient) + sc, fuel = bepi, 3_600. + c3, v∞ = 100., 20. + + # Convenience function for these tests + Thesis.mbh(fbs, lw, la, sp, dp) = mbh(fbs, sc, fuel, lw, c3, v∞, la, easy_cost, + search_patience=sp, drill_patience=dp, verbose=true) - # # Then the inner loop builder function - println("Testing inner loop solver function") - mission = Thesis.inner_loop_solve(test_mg) - @test !mission.converged - - # For the valid case we need to use a lambert's solver - 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, 1.01v∞_in, -1.01v∞_in, tof, 0.01*ones(20,3)) - mission_guess = Mission_Guess(bepi, 12_000., leave, v∞_out, [phase]) - mission = Thesis.inner_loop_solve(mission_guess) - @test mission.converged - if !mission.converged println(mission.message) end - - # Now we're ready to test the whole thing! - println("Testing whole thing") - flybys = [ Earth, Venus ] - start_mass = 12_000. - launch_window = DateTime(1997, 10, 1), DateTime(1997, 10, 31) - max_C3 = 15. - max_v∞ = 20. - latest_arrival = DateTime(2001, 1, 1) - best, archive = mbh(flybys, bepi, start_mass, launch_window, max_C3, max_v∞, latest_arrival, - verbose=true) + # Again, we're going to test a simple case first + # This one seems to converge really easily, so we don't run it much + launch_window = DateTime(1992,11,1), DateTime(1992,12,1) + latest_arrival = DateTime(1993,6,1) + planets = [ Earth, Venus ] + best, archive = mbh(planets, launch_window, latest_arrival, 2, 3) @test typeof(best) == Mission + p = plot(best, title="MBH Test Solution") + savefig(p,"../plots/mbh_test_1_phase.html") + # Now for a more complicated two-phase mission, with a bigger date range + # This is known to have a solution though + planets = [Earth, Venus, Mars] + launch_window = DateTime(2021,6,1), DateTime(2022,6,1) + latest_arrival = DateTime(2024,1,1) + best, archive = mbh(planets, launch_window, latest_arrival, 10, 3) + @test typeof(best) == Mission + p = plot(best, title="MBH Test Solution") + savefig(p,"../plots/mbh_test_2_phase.html") end diff --git a/julia/test/inner_loop/nlp_solver.jl b/julia/test/inner_loop/nlp_solver.jl index 765b9d7..2db034b 100644 --- a/julia/test/inner_loop/nlp_solver.jl +++ b/julia/test/inner_loop/nlp_solver.jl @@ -6,7 +6,7 @@ # Test the optimizer for a one-phase mission # The lambert's solver said this should be pretty valid - launch_window = [DateTime(1992,11,1), DateTime(1992,12,1)] + launch_window = DateTime(1992,11,1), DateTime(1992,12,1) latest_arrival = DateTime(1993,6,1) leave, arrive = DateTime(1992,11,19), DateTime(1993,4,1) test_leave = DateTime(1992,11,12) @@ -25,7 +25,7 @@ # Now we can look at a slightly more complicated trajectory flybys = [Earth, Venus, Mars] - launch_window = [DateTime(2021,10,1), DateTime(2021,12,1)] + launch_window = DateTime(2021,10,1), DateTime(2021,12,1) latest_arrival = DateTime(2023,1,1) dates = [DateTime(2021,11,1), DateTime(2022,3,27), DateTime(2022,8,28)] phases = Vector{Phase}() @@ -45,7 +45,7 @@ # Here is the final, most complicated, trajectory to test flybys = [Earth, Venus, Earth, Mars, Earth, Jupiter] - launch_window = [DateTime(2023,1,1), DateTime(2024,1,1)] + launch_window = DateTime(2023,1,1), DateTime(2024,1,1) latest_arrival = DateTime(2031,1,1) dates = [DateTime(2023,5,23), DateTime(2023,10,21), diff --git a/julia/test/inner_loop/propagator.jl b/julia/test/inner_loop/propagator.jl index 68ef7c4..b0e7b3c 100644 --- a/julia/test/inner_loop/propagator.jl +++ b/julia/test/inner_loop/propagator.jl @@ -18,7 +18,4 @@ state = prop_one([1., 0., 0.], start, bepi, stepsize) @test state[7] == start_mass - bepi.mass_flow_rate*stepsize - # Test that a bad ΔV throws an error - @test_throws Thesis.PropOne_Error prop_one([1.5, 0., 0.], start, bepi, stepsize) - end diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index fba9d59..1f16d54 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -18,5 +18,5 @@ end include("inner_loop/laguerre-conway.jl") include("inner_loop/propagator.jl") include("inner_loop/nlp_solver.jl") - # include("inner_loop/monotonic_basin_hopping.jl") + include("inner_loop/monotonic_basin_hopping.jl") end