Monotonic Basin Hopping is started

.gitlab-ci.yml

@@ -9,3 +9,8 @@ unit-test-job:
     - apt-get install -y unzip
     - mkdir julia/plots
     - julia --project=julia/Project.toml -E 'using Pkg; Pkg.test()'
+  artifacts:
+    paths:
+      - plots/plot_test.html
+      - plots/find_closest_test.html
+    expire_in: 1 week

julia/src/Thesis.jl

@@ -9,4 +9,5 @@ module Thesis
   include("./laguerre-conway.jl")
   include("./propagator.jl")
   include("./find_closest.jl")
+  include("./monotonic_basin_hopping.jl")
 end

julia/src/find_closest.jl

@@ -1,18 +1,20 @@
 using NLsolve
 
+export nlp_solve
+
 function treat_inputs(x::AbstractVector)
   n::Int = length(x)/3
   reshape(x,(3,n))'
 end
 
-function single_shoot(start::Vector{Float64},
+function nlp_solve(start::Vector{Float64},
-                      final::Vector{Float64},
+                   final::Vector{Float64},
-                      craft::Sc,
+                   craft::Sc,
-                      μ::Float64,
+                   μ::Float64,
-                      t0::Float64,
+                   t0::Float64,
-                      tf::Float64,
+                   tf::Float64,
-                      x0::AbstractVector,
+                   x0::AbstractVector;
-                      tol=1e-6)
+                   tol=1e-6)
 
   n::Int = length(x0)/3
   function f!(F,x)
@@ -20,6 +22,6 @@ function single_shoot(start::Vector{Float64},
     F[7:3n] .= 0.
   end
 
-  return nlsolve(f!, x0, ftol=tol, autodiff=:forward, iterations=10_000)
+  return nlsolve(f!, x0, ftol=tol, autodiff=:forward, iterations=1_000)
 
 end

julia/src/monotonic_basin_hopping.jl

@@ -0,0 +1,56 @@
+function perturb(x::AbstractVector, n::Int)
+  perturb_vector = 0.02 * rand(Float64, (3n)) .- 0.01
+  return x + perturb_vector
+end
+
+function mass_better(x_star::AbstractVector,
+                     x_current::AbstractVector,
+                     start::AbstractVector,
+                     final::AbstractVector,
+                     craft::Sc,
+                     μ::AbstractFloat,
+                     t0::AbstractFloat,
+                     tf::AbstractFloat)
+  mass_star = prop(treat_inputs(x_star), start, craft, μ, tf-t0)[2]
+  mass_current = prop(treat_inputs(x_current), start, craft, μ, tf-t0)[2]
+  return mass_star > mass_current
+end
+
+function mbh(start::AbstractVector,
+             final::AbstractVector,
+             craft::Sc,
+             μ::AbstractFloat,
+             t0::AbstractFloat,
+             tf::AbstractFloat,
+             n::Int,
+             num_iters::Int=10,
+             tol=1e-6)
+  i::Int = 0
+  archive = []
+  x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol)
+  while converged(x_star) == false
+    x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol)
+  end
+
+  x_current = x_star
+  push!(archive, x_current)
+
+  while i < num_iters
+    x_star = nlp_solve(start, final, craft, μ, t0, tf, perturb(x_current.zero,n), tol=tol)
+    if converged(x_star) && mass_better(x_star.zero, x_current.zero, start, final, craft, μ, t0, tf)
+      x_current = x_star
+      push!(archive, x_star)
+    else
+      while converged(x_star) == false
+        x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol)
+      end
+      if mass_better(x_star.zero, x_current.zero, start, final, craft, μ, t0, tf)
+        x_current = x_star
+        push!(archive, x_star)
+      end
+    end
+    i += 1
+  end
+
+  return x_current, archive
+end

julia/test/find_closest.jl

@@ -10,7 +10,7 @@
   i = rand(0.01:0.01:π/6)
   T = 2π*√(a^3/μs["Earth"])
   prop_time = 2T
-  n = 50
+  n = 30
 
   # A simple orbit raising
   start = oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"])
@@ -19,8 +19,8 @@
   new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"])
 
   # This should be close enough to 0.6
-  x0 = repeat([0.59, 0., 0.], n)
+  x0 = repeat([0.55, 0., 0.], n)
-  result = Thesis.single_shoot(start, final, sc, μs["Earth"], 0.0, prop_time, x0)
+  result = nlp_solve(start, final, sc, μs["Earth"], 0.0, prop_time, x0)
 
   # Test and plot
   @test converged(result)
@@ -29,7 +29,7 @@
   path3 = prop(zeros((100,3)), path2[end,:], sc, μs["Earth"], new_T)[1]
   path4 = prop(zeros((100,3)), final, sc, μs["Earth"], new_T)[1]
   savefig(plot_orbits([path1, path2, path3, path4],
-                      labels=["inital", "transit", "after transit", "final"],
+                      labels=["initial", "transit", "after transit", "final"],
                       colors=["#FFFFFF","#FF4444","#44FF44","#4444FF"]),
           "../plots/find_closest_test.html")
   if converged(result)

julia/test/monotonic_basin_hopping.jl

@@ -0,0 +1,29 @@
+@testset "Monotonic Basin Hopping" begin
+
+  using Thesis: mbh
+
+  # Initial Setup
+  sc = Sc("test")
+  a = rand(15000:1.:40000)
+  e = rand(0.01:0.01:0.5)
+  i = rand(0.01:0.01:π/6)
+  T = 2π*√(a^3/μs["Earth"])
+  prop_time = 2T
+  n = 25
+
+  # A simple orbit raising
+  start = oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"])
+  ΔVs = repeat([0.6, 0., 0.]', outer=(n,1))
+  final = prop(ΔVs, start, sc, μs["Earth"], prop_time)[1][end,:]
+  new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"])
+
+  # This should be close enough to 0.6
+  best, archive = mbh(start, final, sc, μs["Earth"], 0.0, prop_time, n)
+
+  # Test and plot
+  @test converged(best)
+  for path in archive
+   @test converged(path)
+  end
+
+end

julia/test/runtests.jl

@@ -11,6 +11,7 @@ using Thesis
   include("propagator.jl")
   include("plotting.jl")
   include("find_closest.jl")
+  include("monotonic_basin_hopping.jl")
 end
 
 print()