Fully working mbh for two phase missions!

.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

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")
 

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...")

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)
+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)
-  ϵ = 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)
+  else
     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)).^(-α)
+  return 1 .+ (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
   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

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
 

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

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

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
 
-  # # Then the inner loop builder function
-  println("Testing inner loop solver function")
-  mission = Thesis.inner_loop_solve(test_mg)
-  @test !mission.converged
+  # Mission Parameters that won't change (they're very lenient)
+  sc, fuel = bepi, 3_600.
+  c3, v∞ = 100., 20.
 
-  # 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
+  # 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)
   
-  # 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

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),

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

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