Built the population member gen function

julia/src/Thesis.jl

@@ -20,5 +20,6 @@ module Thesis
   include("./inner_loop/monotonic_basin_hopping.jl")
   include("./inner_loop/phase.jl")
   include("./inner_loop/inner_loop.jl")
+  include("./outer_loop.jl")
 
 end

julia/src/inner_loop/inner_loop.jl

@@ -47,9 +47,9 @@ function inner_loop(launch_date::DateTime,
   thrust_profiles = []
 
   for phase in phases
-    planet1_state = [spkssb(ids[phase.from_planet], time, "J2000"); 0.0]
+    planet1_state = [spkssb(ids[phase.from_planet], time, "ECLIPJ2000"); 0.0]
     time += phase.time_of_flight
-    planet2_state = [spkssb(ids[phase.to_planet], time, "J2000"); 0.0]
+    planet2_state = [spkssb(ids[phase.to_planet], time, "ECLIPJ2000"); 0.0]
     start = planet1_state + [0., 0., 0., phase.v∞_outgoing..., start_mass]
     final = planet2_state + [0., 0., 0., phase.v∞_incoming..., start_mass]
     println(start)

julia/src/outer_loop.jl

@@ -0,0 +1,63 @@
+using Random, Dates
+
+export gen_decision_vector
+
+"""
+Returns a random date between two dates
+"""
+function gen_date(date_range::Vector{DateTime})
+  l0, lf = date_range
+  l0 + Dates.Millisecond(floor(rand()*(lf-l0).value))
+end
+
+"""
+Returns a random amount of time in a range
+"""
+function gen_period(date_range::Vector{DateTime})
+  l0, lf = date_range
+  Dates.Millisecond(floor(rand()*(lf-l0).value))
+end
+
+"""
+So ideally, this should generate a nice random decision vector, given the constraints.
+Everything that you need to produce a vector of phases
+  
+Start with an empty vector of the right size
+You need:
+  - launch_date
+  - 3 components v∞_out for Earth
+  - and then up to four flybys which contain:
+    - a planet
+    - three components v∞_in
+    - turning angle (in ecliptic)
+    - tof to planet
+and finally, the ending planet is held fixed
+"""
+function gen_decision_vector(launch_range::Vector{DateTime},
+                             target::String,
+                             arrival_deadline::DateTime)
+  phases = Vector{Phase}()
+  launch_date = gen_date(launch_range)
+  v∞_out = 20rand(Float64,3) .- 10.
+  
+  # Generate the planets (or null flybys)
+  planets = [ ["Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"];
+              repeat(["None"],8) ] # Just as likely to get a planet as no flyby
+  long_flybys = [ "Earth"; filter(x -> x != "None", rand(planets, 4)); target ]
+  # This will cut the flybys off if the target shows up early
+  flybys = long_flybys[1:findfirst(x->x==target, long_flybys)]
+
+  time = launch_date
+  for i in 1:length(flybys)-1
+    v∞_in = 20rand(Float64,3) .- 10. # Generate the v∞_in components
+    tof = gen_period([time,arrival_deadline])
+    time += tof
+    push!(phases,Phase(flybys[i],flybys[i+1], tof.value/1000, v∞_out, v∞_in))
+    v∞_out_base = rand(Float64,3) .- 0.5
+    v∞_out = norm(v∞_in) * v∞_out_base/norm(v∞_out_base)
+  end
+
+  return phases
+
+end
+

julia/test/inner_loop/inner_loop.jl

@@ -22,9 +22,9 @@
   phase2 = Phase(p1, p2, leg2_tof, v∞s[3], v∞s[4])
 
   # For finding the best trajectories
-  earth_state = [spkssb(Thesis.ids["Earth"], launch_j2000, "J2000"); start_mass]
+  earth_state = [spkssb(Thesis.ids["Earth"], launch_j2000, "ECLIPJ2000"); start_mass]
-  p1_state = [spkssb(Thesis.ids[p1], launch_j2000+leg1_tof, "J2000"); start_mass]
+  p1_state = [spkssb(Thesis.ids[p1], launch_j2000+leg1_tof, "ECLIPJ2000"); start_mass]
-  p2_state = [spkssb(Thesis.ids[p2], launch_j2000+leg1_tof+leg2_tof, "J2000"); start_mass]
+  p2_state = [spkssb(Thesis.ids[p2], launch_j2000+leg1_tof+leg2_tof, "ECLIPJ2000"); start_mass]
   earth = prop(zeros(100,3), earth_state, sc, μs["Sun"], 3600*24*365.)[1]
   p1_path = prop(zeros(100,3), p1_state, sc, μs["Sun"], 3600*24*365*2.)[1]
   p2_path = prop(zeros(100,3), p2_state, sc, μs["Sun"], 3600*24*365*8.)[1]

julia/test/inner_loop/monotonic_basin_hopping.jl

@@ -5,76 +5,78 @@
   println("Testing Monotonic Basin Hopper")
 
   # Initial Setup
-  sc = Sc("test")
+  # sc = Sc("test")
-  a = rand(50_000:1.:100_000)
+  # a = rand(50_000:1.:100_000)
-  e = rand(0.01:0.01:0.5)
+  # e = rand(0.01:0.01:0.5)
-  i = rand(0.01:0.01:π/6)
+  # i = rand(0.01:0.01:π/6)
-  T = 2π*√(a^3/μs["Earth"])
+  # T = 2π*√(a^3/μs["Earth"])
-  prop_time = 0.5T
+  # prop_time = 0.5T
-  n = 20
+  # n = 20
-  start_mass = 10_000.
+  # start_mass = 10_000.
 
-  # A simple orbit raising
+  # # A simple orbit raising
-  start = [oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"]); start_mass]
+  # start = [oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"]); start_mass]
-  Tx, Ty, Tz = conv_T(repeat([0.8], n), repeat([0.], n), repeat([0.], n),
+  # Tx, Ty, Tz = conv_T(repeat([0.8], n), repeat([0.], n), repeat([0.], n),
-                      start,
+                      # start,
-                      sc,
+                      # sc,
-                      prop_time,
+                      # prop_time,
-                      μs["Earth"])
+                      # μs["Earth"])
-  nominal_path, final = prop(hcat(Tx, Ty, Tz), start, sc, μs["Earth"], prop_time)
+  # nominal_path, final = prop(hcat(Tx, Ty, Tz), start, sc, μs["Earth"], prop_time)
-  new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"])
+  # new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"])
 
-  # Find the best solution
+  # # Find the best solution
-  best, archive = mbh(start,
+  # best, archive = mbh(start,
-                      final,
+                      # final,
-                      sc,
+                      # sc,
-                      μs["Earth"],
+                      # μs["Earth"],
-                      0.0,
+                      # 0.0,
-                      prop_time,
+                      # prop_time,
-                      n,
+                      # n,
-                      search_patience_lim=25,
+                      # search_patience_lim=25,
-                      drill_patience_lim=50,
+                      # drill_patience_lim=50,
-                      verbose=true)
+                      # verbose=true)
 
-  # Test and plot
+  # # Test and plot
-  @test best.converged
+  # @test best.converged
-  transit, calc_final = prop(best.zero, start, sc, μs["Earth"], prop_time)
+  # transit, calc_final = prop(best.zero, start, sc, μs["Earth"], prop_time)
-  initial_path = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1]
+  # initial_path = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1]
-  after_transit = prop(zeros((100,3)), calc_final, sc, μs["Earth"], new_T)[1]
+  # after_transit = prop(zeros((100,3)), calc_final, sc, μs["Earth"], new_T)[1]
-  final_path = prop(zeros((100,3)), final, sc, μs["Earth"], new_T)[1]
+  # final_path = prop(zeros((100,3)), final, sc, μs["Earth"], new_T)[1]
-  savefig(plot_orbits([initial_path, nominal_path, final_path],
+  # savefig(plot_orbits([initial_path, nominal_path, final_path],
-                      labels=["initial", "nominal transit", "final"],
+                      # labels=["initial", "nominal transit", "final"],
-                      colors=["#FF4444","#44FF44","#4444FF"]),
+                      # colors=["#FF4444","#44FF44","#4444FF"]),
-                      "../plots/mbh_nominal.html")
+                      # "../plots/mbh_nominal.html")
-  savefig(plot_orbits([initial_path, transit, after_transit, final_path],
+  # savefig(plot_orbits([initial_path, transit, after_transit, final_path],
-                      labels=["initial", "transit", "after transit", "final"],
+                      # labels=["initial", "transit", "after transit", "final"],
-                      colors=["#FFFFFF", "#FF4444","#44FF44","#4444FF"]),
+                      # colors=["#FFFFFF", "#FF4444","#44FF44","#4444FF"]),
-                      "../plots/mbh_best.html")
+                      # "../plots/mbh_best.html")
-  i = 0
+  # i = 0
-  best_mass = calc_final[end]
+  # best_mass = calc_final[end]
-  nominal_mass = final[end]
+  # nominal_mass = final[end]
-  masses = []
+  # masses = []
-  for candidate in archive
+  # for candidate in archive
-    @test candidate.converged
+    # @test candidate.converged
-    path2, calc_final = prop(candidate.zero, start, sc, μs["Earth"], prop_time)
+    # path2, calc_final = prop(candidate.zero, start, sc, μs["Earth"], prop_time)
-    push!(masses, calc_final[end])
+    # push!(masses, calc_final[end])
-    @test norm(calc_final[1:6] - final[1:6]) < 1e-4
+    # @test norm(calc_final[1:6] - final[1:6]) < 1e-4
-  end
+  # end
-  @test best_mass == maximum(masses)
+  # @test best_mass == maximum(masses)
 
-  # This won't always work since the test is reduced in fidelity,
+  # # This won't always work since the test is reduced in fidelity,
-  # but hopefully will usually work:
+  # # but hopefully will usually work:
-  @test (start_mass - best_mass) < 1.1 * (start_mass - nominal_mass)
+  # @test (start_mass - best_mass) < 1.1 * (start_mass - nominal_mass)
 
   # Now let's test a sun case. This should be pretty close to begin with
+  start_mass = 10_000.
   launch_date = DateTime(2016,3,28)
   launch_j2000 = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS"))
-  earth_start = [spkssb(ids["Earth"], launch_j2000, "J2000"); 1e5]
+  earth_start = [spkssb(ids["Earth"], launch_j2000, "ECLIPJ2000"); start_mass]
   earth_speed = earth_start[4:6]
   v∞ = 3.0*earth_speed/norm(earth_speed)
   start = earth_start + [zeros(3); v∞; 0.0]
-  final = [1.62914115303947e8, 1.33709639408102e8, 5.690490452749867e7, -16.298522963602757, 15.193294491415365, 6.154820267250081, 1.0001e8]
   tof = 3600*24*30*10.75
+  mars_state = [spkssb(Thesis.ids["Mars"], launch_j2000+tof, "ECLIPJ2000"); start_mass]
+  final = mars_state + [ zeros(3); [-1.1, -3., -2.6]; 0.0 ]
   a = xyz_to_oe(final, μs["Sun"])[1]
   T = 2π*√(a^3/μs["Sun"])
   n = 20

julia/test/outer_loop.jl

@@ -0,0 +1,17 @@
+@testset "Outer Loop" begin
+
+  using Dates
+
+  println("Testing Genetic Algorithm")
+
+  launch_range = [ DateTime(2016,3,28), DateTime(2019,3,28) ]
+  target = "Saturn"
+  deadline = DateTime(2028,12,31)
+
+  # First let's just test that we can generate a member of the population
+  member = gen_decision_vector(launch_range, target, deadline)
+  println(member)
+
+  @test typeof(member) == Vector{Phase}
+
+end

julia/test/runtests.jl

@@ -21,6 +21,7 @@ end
   include("inner_loop/find_closest.jl")
   include("inner_loop/monotonic_basin_hopping.jl")
   include("inner_loop/inner_loop.jl")
+  include("outer_loop.jl")
 end
 
 print()