From e83a705e77110d5a66934e616621759e497c0127 Mon Sep 17 00:00:00 2001 From: Connor Date: Sun, 12 Sep 2021 19:25:40 -0600 Subject: [PATCH] I think I'm pretty much at 2b now --- julia/src/inner_loop/find_closest.jl | 34 +++++++++------- julia/src/inner_loop/inner_loop.jl | 39 +++++++++++-------- julia/src/inner_loop/laguerre-conway.jl | 5 +-- .../src/inner_loop/monotonic_basin_hopping.jl | 28 ++++++------- julia/src/inner_loop/propagator.jl | 10 +++-- julia/test/inner_loop/inner_loop.jl | 10 +++-- .../inner_loop/monotonic_basin_hopping.jl | 5 ++- julia/test/runtests.jl | 10 ++--- 8 files changed, 81 insertions(+), 60 deletions(-) diff --git a/julia/src/inner_loop/find_closest.jl b/julia/src/inner_loop/find_closest.jl index 9c8faa2..8d2dd45 100644 --- a/julia/src/inner_loop/find_closest.jl +++ b/julia/src/inner_loop/find_closest.jl @@ -24,28 +24,32 @@ function nlp_solve(start::Vector{Float64}, t0::Float64, tf::Float64, x0::Matrix{Float64}; - tol=1e-6) + tol=1e-6, + num_iters=1_000) function f!(F,x) F .= 0.0 F[1:6, 1] .= prop_nlsolve(tanh.(x), start, craft, μ, tf-t0) .- final end - # return nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=1_000) - p = addprocs(1) - response = Channel(1) - @async put!(response, remotecall_fetch(nlsolve, 2, f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=1_000)) + # return nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters) - start=time() - while !isready(response) && (time() - start) < 30. - sleep(0.1) - end + return nlsolve(f!, atanh.(x0), ftol=tol, iterations=num_iters) - if isready(response) - return fetch(response) - else - rmprocs(p); - return "error" - end + # p = addprocs(1) + # response = Channel(1) + # @async put!(response, remotecall_fetch(nlsolve, 2, f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=1_000)) + + # start=time() + # while !isready(response) && (time() - start) < 30. + # sleep(0.1) + # end + + # if isready(response) + # return fetch(response) + # else + # rmprocs(p); + # return "error" + # end end \ No newline at end of file diff --git a/julia/src/inner_loop/inner_loop.jl b/julia/src/inner_loop/inner_loop.jl index fde352e..c03b4a6 100644 --- a/julia/src/inner_loop/inner_loop.jl +++ b/julia/src/inner_loop/inner_loop.jl @@ -47,23 +47,28 @@ function inner_loop(launch_date::DateTime, time = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS")) thrust_profiles = [] - for phase in phases - planet1_state = spkssb(ids[phase.from_planet], time, "J2000") - time += phase.time_of_flight - planet2_state = spkssb(ids[phase.to_planet], time, "J2000") - # TODO: Come up with improved method of calculating "n" - start = planet1_state + [0., 0., 0., phase.v∞_outgoing...] - final = planet2_state + [0., 0., 0., phase.v∞_incoming...] - if mbh_specs == nothing - best = mbh(start, final, craft, μs["Sun"], 0.0, phase.time_of_flight, 10, - verbose=verbose)[1] - else - best = mbh(start, final, craft, μs["Sun"], 0.0, phase.time_of_flight, 10, - verbose=verbose, num_iters=mbh_specs[1], patience_level=mbh_specs[2])[1] + try + for phase in phases + planet1_state = spkssb(ids[phase.from_planet], time, "J2000") + time += phase.time_of_flight + planet2_state = spkssb(ids[phase.to_planet], time, "J2000") + # TODO: Come up with improved method of calculating "n" + start = planet1_state + [0., 0., 0., phase.v∞_outgoing...] + final = planet2_state + [0., 0., 0., phase.v∞_incoming...] + if mbh_specs === nothing + best = mbh(start, final, craft, μs["Sun"], 0.0, phase.time_of_flight, 10, + verbose=verbose)[1] + else + num_iters, sil, dil = mbh_specs + best = mbh(start, final, craft, μs["Sun"], 0.0, phase.time_of_flight, 10, + verbose=verbose, num_iters=num_iters, search_patience_lim=sil, + drill_patience_lim=dil) + end + push!(thrust_profiles, best) + craft.mass = prop(tanh.(best.zero), planet1_state, sc, μs["Sun"], prop_time)[2][end] end - push!(thrust_profiles, best) - craft.mass = prop(tanh.(best.zero), planet1_state, sc, μs["Sun"], prop_time)[2][end] + return craft.mass, thrust_profiles + catch + return "One path did not converge" end - - return craft.mass, thrust_profiles end diff --git a/julia/src/inner_loop/laguerre-conway.jl b/julia/src/inner_loop/laguerre-conway.jl index 7ca4503..5f332c8 100644 --- a/julia/src/inner_loop/laguerre-conway.jl +++ b/julia/src/inner_loop/laguerre-conway.jl @@ -22,6 +22,7 @@ function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T sign = dF >= 0 ? 1 : -1 ΔE_new = ΔE - n*F / ( dF + sign * √(abs((n-1)^2*dF^2 - n*(n-1)*F*d2F ))) i += 1 + if i > 100 throw(ErrorException("LaGuerre-Conway did not converge!")) end end F = 1 - a/r0_mag * (1-cos(ΔE)) G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE) @@ -40,6 +41,7 @@ function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T sign = dF >= 0 ? 1 : -1 ΔH_new = ΔH - n*F / ( dF + sign * √(abs((n-1)^2*dF^2 - n*(n-1)*F*d2F ))) i += 1 + if i > 100 throw(ErrorException("LaGuerre-Conway did not converge!")) end end F = 1 - a/r0_mag * (1-cos(ΔH)) G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH) @@ -48,9 +50,6 @@ function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T Gt = 1 - a/r * (1-cos(ΔH)) end - if i > 100 - throw(ErrorException("LaGuerre-Conway did not converge!")) - end return [ F*state[1:3] + G*state[4:6]; Ft*state[1:3] + Gt*state[4:6]] diff --git a/julia/src/inner_loop/monotonic_basin_hopping.jl b/julia/src/inner_loop/monotonic_basin_hopping.jl index 74ba046..12c4b8e 100644 --- a/julia/src/inner_loop/monotonic_basin_hopping.jl +++ b/julia/src/inner_loop/monotonic_basin_hopping.jl @@ -27,7 +27,8 @@ function mbh(start::AbstractVector, tf::AbstractFloat, n::Int; num_iters=25, - patience_level::Int=200, + search_patience_lim::Int=200, + drill_patience_lim::Int=200, tol=1e-6, verbose=false) @@ -38,25 +39,24 @@ function mbh(start::AbstractVector, while true i += 1 if verbose print("\r",i) end - # TODO: Should this be two separate "impatience" values? - impatience = 0 - x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol) - while converged(x_star) == false && impatience < patience_level - impatience += 1 - x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol) - println(impatience) + search_impatience = 0 + drill_impatience = 0 + x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol, num_iters=100) + while converged(x_star) == false && search_impatience < search_patience_lim + search_impatience += 1 + x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol, num_iters=100) end - if impatience > patience_level break end - impatience = 0 + if drill_impatience > drill_patience_lim break end + drill_impatience = 0 if converged(x_star) x_current = x_star - while impatience < patience_level + while drill_impatience < drill_patience_lim x_star = nlp_solve(start, final, craft, μ, t0, tf, perturb(tanh.(x_current.zero),n), tol=tol) if converged(x_star) && mass_est(tanh.(x_star.zero)) < mass_est(tanh.(x_current.zero)) x_current = x_star - impatience = 0 + drill_impatience = 0 else - impatience += 1 + drill_impatience += 1 end end push!(archive, x_current) @@ -65,6 +65,8 @@ function mbh(start::AbstractVector, end if verbose println() end + if archive == [] error("there were no converged paths") end + current_best_mass = 1e8 best = archive[1] for candidate in archive diff --git a/julia/src/inner_loop/propagator.jl b/julia/src/inner_loop/propagator.jl index fbf0c85..9d9eb0b 100644 --- a/julia/src/inner_loop/propagator.jl +++ b/julia/src/inner_loop/propagator.jl @@ -113,11 +113,15 @@ function prop_nlsolve(ΔVs::Matrix{T}, time::Float64) where {T <: Real, S <: Real} n = size(ΔVs)[1] - for i in 1:n - state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n) + try + for i in 1:n + state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n) + end + return state + catch + return [0., 0., 0., 0., 0., 0.] end - return state end diff --git a/julia/test/inner_loop/inner_loop.jl b/julia/test/inner_loop/inner_loop.jl index 6a1ef8e..7913d59 100644 --- a/julia/test/inner_loop/inner_loop.jl +++ b/julia/test/inner_loop/inner_loop.jl @@ -8,7 +8,7 @@ phase1 = Phase("Earth", "Mars", - 3600*24*365*1.5, + 3600*24*365*1.85, [1., 0.3, 0.], [3., 3., 0.]) @@ -18,8 +18,12 @@ [2., 3.7416573867739413, 0.], [0.3, 1., 0.]) - inner_loop(DateTime(2024,3,5), sc, [phase1, phase2], verbose=true, mbh_specs=(5,100)) + result = inner_loop(DateTime(2024,3,5), + sc, + [phase1, phase2], + verbose=true, + mbh_specs=(5,50,100)) - @test true + @test result == "One path did not converge" end diff --git a/julia/test/inner_loop/monotonic_basin_hopping.jl b/julia/test/inner_loop/monotonic_basin_hopping.jl index bda02c1..d53233e 100644 --- a/julia/test/inner_loop/monotonic_basin_hopping.jl +++ b/julia/test/inner_loop/monotonic_basin_hopping.jl @@ -1,5 +1,7 @@ @testset "Monotonic Basin Hopping" begin + using PlotlyJS, NLsolve + println("Testing Monotonic Basin Hopper") # Initial Setup @@ -32,7 +34,8 @@ prop_time, n, num_iters=5, - patience_level=50, + search_patience_lim=50, + drill_patience_lim=100, verbose=true) # Test and plot diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 33914b4..5ca2423 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -14,11 +14,11 @@ end # Tests @testset "All Tests" begin - # include("spacecraft.jl") - # include("plotting.jl") - # include("inner_loop/laguerre-conway.jl") - # include("inner_loop/propagator.jl") - # include("inner_loop/find_closest.jl") + include("spacecraft.jl") + include("plotting.jl") + include("inner_loop/laguerre-conway.jl") + include("inner_loop/propagator.jl") + include("inner_loop/find_closest.jl") include("inner_loop/monotonic_basin_hopping.jl") include("inner_loop/inner_loop.jl") end