I think I'm pretty much at 2b now
This commit is contained in:
@@ -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
|
||||
@@ -47,6 +47,7 @@ function inner_loop(launch_date::DateTime,
|
||||
|
||||
time = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS"))
|
||||
thrust_profiles = []
|
||||
try
|
||||
for phase in phases
|
||||
planet1_state = spkssb(ids[phase.from_planet], time, "J2000")
|
||||
time += phase.time_of_flight
|
||||
@@ -54,16 +55,20 @@ function inner_loop(launch_date::DateTime,
|
||||
# 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
|
||||
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=mbh_specs[1], patience_level=mbh_specs[2])[1]
|
||||
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
|
||||
|
||||
return craft.mass, thrust_profiles
|
||||
catch
|
||||
return "One path did not converge"
|
||||
end
|
||||
end
|
||||
|
||||
@@ -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]]
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -113,11 +113,15 @@ function prop_nlsolve(ΔVs::Matrix{T},
|
||||
time::Float64) where {T <: Real, S <: Real}
|
||||
|
||||
n = size(ΔVs)[1]
|
||||
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
|
||||
|
||||
|
||||
end
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user