diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index e263b3e..05a1414 100644 --- a/.gitlab-ci.yml +++ b/.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 diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl index b21e0bc..468c411 100644 --- a/julia/src/Thesis.jl +++ b/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 \ No newline at end of file diff --git a/julia/src/find_closest.jl b/julia/src/find_closest.jl index 5f12307..9d496c7 100644 --- a/julia/src/find_closest.jl +++ b/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}, - final::Vector{Float64}, - craft::Sc, - μ::Float64, - t0::Float64, - tf::Float64, - x0::AbstractVector, - tol=1e-6) +function nlp_solve(start::Vector{Float64}, + final::Vector{Float64}, + craft::Sc, + μ::Float64, + t0::Float64, + tf::Float64, + x0::AbstractVector; + 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 \ No newline at end of file diff --git a/julia/src/monotonic_basin_hopping.jl b/julia/src/monotonic_basin_hopping.jl new file mode 100644 index 0000000..cca4c08 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/julia/test/find_closest.jl b/julia/test/find_closest.jl index 0236d56..f7870a3 100644 --- a/julia/test/find_closest.jl +++ b/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) - result = Thesis.single_shoot(start, final, sc, μs["Earth"], 0.0, prop_time, x0) + x0 = repeat([0.55, 0., 0.], n) + 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) diff --git a/julia/test/monotonic_basin_hopping.jl b/julia/test/monotonic_basin_hopping.jl new file mode 100644 index 0000000..45e8373 --- /dev/null +++ b/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 diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index cd904bc..dc3d99d 100644 --- a/julia/test/runtests.jl +++ b/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()