Update for meeting with Bosanac
This commit is contained in:
Connor
@@ -25,8 +25,8 @@ module Thesis
|
||||
include("./plotting.jl")
|
||||
include("./inner_loop/laguerre-conway.jl")
|
||||
include("./inner_loop/propagator.jl")
|
||||
include("./inner_loop/phase.jl")
|
||||
include("./inner_loop/monotonic_basin_hopping.jl")
|
||||
include("./inner_loop/nlp_solver.jl")
|
||||
# include("./inner_loop/monotonic_basin_hopping.jl")
|
||||
# include("./outer_loop.jl")
|
||||
|
||||
end
|
||||
|
||||
@@ -26,3 +26,8 @@ struct Mass_Error{T} <: Exception where T <: Real
|
||||
mass::T
|
||||
end
|
||||
Base.showerror(io::IO, e::Mass_Error) = print(io, "mass (", e.mass, ") got too low in propagation")
|
||||
|
||||
struct Mission_Creation_Error{T} <: Exception where T <: Real
|
||||
num_entries::T
|
||||
end
|
||||
Base.showerror(io::IO, e::Mission_Creation_Error) = print(io, "Tried to create a mission with $(e.num_entries) entries")
|
||||
|
||||
110
julia/src/inner_loop/nlp_solver.jl
Normal file
110
julia/src/inner_loop/nlp_solver.jl
Normal file
@@ -0,0 +1,110 @@
|
||||
using SNOW
|
||||
|
||||
export solve_mission
|
||||
|
||||
struct Result
|
||||
converged::Bool
|
||||
info::Symbol
|
||||
sol::Matrix{Float64}
|
||||
end
|
||||
|
||||
"""
|
||||
This function will take an initial guess to a mission profile and leverage Ipopt to solve the mission
|
||||
"""
|
||||
function solve_mission( guess::Mission_Guess;
|
||||
tol=1e-8,
|
||||
verbose::Bool=false )
|
||||
# First we need to define a function that propagates the mission all the way through
|
||||
# The function takes in a constraints vector (g) and an input vector (x) and spits out
|
||||
# the objective function. In my case, the objective will be constant.
|
||||
#
|
||||
|
||||
# First we define our starting point
|
||||
x0 = Vector(guess)
|
||||
n = size(guess.phases[1].thrust_profile)[1]
|
||||
flybys = [ p.planet for p in guess.phases ]
|
||||
|
||||
# And our NLP
|
||||
function optimizer!(g,x)
|
||||
# Establish initial conditions
|
||||
mass = guess.start_mass
|
||||
v∞_out = x[2:4]
|
||||
current_planet = Earth
|
||||
launch_date = Dates.unix2datetime(x[1])
|
||||
time = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS"))
|
||||
start = state(current_planet, time, v∞_out, mass)
|
||||
|
||||
# Now, for each phase we must require:
|
||||
# - That the ending state matches (unless final leg)
|
||||
# - That the v∞_out == v∞_in (unless final leg)
|
||||
# - That the minimum height is acceptable (unless final leg)
|
||||
# - That the ending position matches (if final leg)
|
||||
# - That the ending mass is acceptable (if final leg)
|
||||
i = 0
|
||||
for flyby in flybys
|
||||
phase_params = x[ 5 + (3n+7) * i : 5 + (3n+7) * (i+1) - 1 ]
|
||||
v∞_in = phase_params[1:3]
|
||||
v∞_out = phase_params[4:6]
|
||||
tof = phase_params[7]
|
||||
thrusts = reshape(phase_params[8:end], (n,3))
|
||||
current_planet = flyby
|
||||
time += tof
|
||||
goal = state(current_planet, time, v∞_in)
|
||||
final = prop(reshape(thrusts, (n,3)), start, guess.sc, tof)[2]
|
||||
if flyby != flybys[end]
|
||||
g[1+8i:6+8i] .= final[1:6] .- goal[1:6]
|
||||
g[7+8i] = norm(v∞_out) - norm(v∞_in)
|
||||
δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
|
||||
g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
|
||||
else
|
||||
g[8length(flybys)+1:8length(flybys)+3] .= final[1:3] .- goal[1:3]
|
||||
g[8length(flybys)+4] = final[7]
|
||||
end
|
||||
start = state(current_planet, time, v∞_out, final[7])
|
||||
i += 1
|
||||
return 1.0
|
||||
end
|
||||
end
|
||||
g = zeros(8length(flybys)+4)
|
||||
optimizer!(g,x0)
|
||||
return g
|
||||
|
||||
|
||||
|
||||
# n = size(x0)[1]
|
||||
|
||||
# function f!(g,x)
|
||||
# try
|
||||
# g[1:6] .= prop(reshape(x,(n,3)), start, craft, tof, primary)[2][1:6]
|
||||
# return 1.
|
||||
# catch e
|
||||
# if isa(e,LaGuerreConway_Error) || isa(e, Mass_Error) || isa(e, PropOne_Error)
|
||||
# return 10_000.
|
||||
# else
|
||||
# rethrow
|
||||
# end
|
||||
# end
|
||||
# end
|
||||
|
||||
# lower_x = -1 * ones(3n)
|
||||
# upper_x = ones(3n)
|
||||
# num_constraints = 6
|
||||
# filename = "ipopt_"*string(rand(UInt))
|
||||
# ipopt_options = Dict("constr_viol_tol" => tol,
|
||||
# "acceptable_constr_viol_tol" => 1e-4,
|
||||
# "max_iter" => 10_000,
|
||||
# "max_cpu_time" => 30.,
|
||||
# "print_level" => 0,
|
||||
# "output_file" => filename)
|
||||
# options = Options(solver=IPOPT(ipopt_options),
|
||||
# derivatives=ForwardAD())
|
||||
# x, _, info = minimize(f!, vec(x0), num_constraints, lower_x, upper_x, final[1:6], final[1:6], options)
|
||||
# rm(filename)
|
||||
# if info in [:Solve_Succeeded, :Solved_To_Acceptable_Level]
|
||||
# return Result(true, info, reshape(x,(n,3)))
|
||||
# else
|
||||
# if verbose println(info) end
|
||||
# return Result(false, info, x0)
|
||||
# end
|
||||
|
||||
end
|
||||
@@ -1,60 +0,0 @@
|
||||
using SNOW
|
||||
|
||||
export solve_phase
|
||||
|
||||
struct Result
|
||||
converged::Bool
|
||||
info::Symbol
|
||||
sol::Matrix{Float64}
|
||||
end
|
||||
|
||||
"""
|
||||
This function will take a single phase (so an initial state, and a final state) and an initial guess
|
||||
to the thrust profile and use an NLP solver to find the nearest thrust profile to that initial guess
|
||||
that satisfies the final state condition
|
||||
"""
|
||||
function solve_phase( start::Vector{Float64},
|
||||
final::Vector{Float64},
|
||||
craft::Sc,
|
||||
tof::Float64,
|
||||
x0::Matrix{Float64},
|
||||
primary::Body=Sun;
|
||||
tol=1e-8,
|
||||
verbose::Bool=false )
|
||||
n = size(x0)[1]
|
||||
|
||||
function f!(g,x)
|
||||
try
|
||||
g[1:6] .= prop(reshape(x,(n,3)), start, craft, tof, primary)[2][1:6]
|
||||
return 1.
|
||||
catch e
|
||||
if isa(e,LaGuerreConway_Error) || isa(e, Mass_Error) || isa(e, PropOne_Error)
|
||||
return 10_000.
|
||||
else
|
||||
rethrow
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
lower_x = -1 * ones(3n)
|
||||
upper_x = ones(3n)
|
||||
num_constraints = 6
|
||||
filename = "ipopt_"*string(rand(UInt))
|
||||
ipopt_options = Dict("constr_viol_tol" => tol,
|
||||
"acceptable_constr_viol_tol" => 1e-4,
|
||||
"max_iter" => 10_000,
|
||||
"max_cpu_time" => 30.,
|
||||
"print_level" => 0,
|
||||
"output_file" => filename)
|
||||
options = Options(solver=IPOPT(ipopt_options),
|
||||
derivatives=ForwardAD())
|
||||
x, _, info = minimize(f!, vec(x0), num_constraints, lower_x, upper_x, final[1:6], final[1:6], options)
|
||||
rm(filename)
|
||||
if info in [:Solve_Succeeded, :Solved_To_Acceptable_Level]
|
||||
return Result(true, info, reshape(x,(n,3)))
|
||||
else
|
||||
if verbose println(info) end
|
||||
return Result(false, info, x0)
|
||||
end
|
||||
|
||||
end
|
||||
@@ -1,6 +1,7 @@
|
||||
export Phase, Mission_Guess, Mission, Bad_Mission
|
||||
export test_phase1, test_phase2
|
||||
export test_mg, test_mission
|
||||
export Vector
|
||||
|
||||
struct Phase
|
||||
planet::Body
|
||||
@@ -55,6 +56,51 @@ function Mission_Guess(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float
|
||||
Mission_Guess(sc, mass, date, v∞, phases, false)
|
||||
end
|
||||
|
||||
"""
|
||||
This is the opposite of the function below. It takes a vector and any other necessary information and outputs a mission guess
|
||||
"""
|
||||
function Mission_Guess(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body})
|
||||
# Variable mission params
|
||||
launch_date = Dates.unix2datetime(x[1])
|
||||
launch_v∞ = x[2:4]
|
||||
|
||||
# Try to intelligently determine n
|
||||
three_n = ((length(x)-4)/length(flybys)) - 7
|
||||
three_n % 3 == 0 || throw(Mission_Creation_Error(three_n))
|
||||
n::Int = three_n/3
|
||||
|
||||
# Build the phases
|
||||
i = 0
|
||||
phases = Vector{Phase}()
|
||||
for flyby in flybys
|
||||
phase_params = x[ 5 + (3n+7) * i : 5 + (3n+7) * (i+1) - 1 ]
|
||||
v∞_in = phase_params[1:3]
|
||||
v∞_out = phase_params[4:6]
|
||||
tof = phase_params[7]
|
||||
thrusts = reshape(phase_params[8:end], (n,3))
|
||||
push!(phases, Phase(flyby, v∞_in, v∞_out, tof, thrusts))
|
||||
i += 1
|
||||
end
|
||||
return Mission_Guess(sc, mass, launch_date, launch_v∞, phases, false)
|
||||
end
|
||||
|
||||
"""
|
||||
This is a function to convert a mission guess into an nlp-ready vector. It doesn't contain all
|
||||
information from the guess though.
|
||||
"""
|
||||
function Base.Vector(g::Mission_Guess)
|
||||
result = Vector{Float64}()
|
||||
push!(result, Dates.datetime2unix(g.launch_date))
|
||||
push!(result, g.launch_v∞...)
|
||||
for phase in g.phases
|
||||
push!(result,phase.v∞_in...)
|
||||
push!(result,phase.v∞_out...)
|
||||
push!(result,phase.tof)
|
||||
push!(result,phase.thrust_profile...)
|
||||
end
|
||||
return result
|
||||
end
|
||||
|
||||
const test_mg = Mission_Guess(bepi, 12_000., DateTime(1992,11,19), [-3.4,1.2,0.1], test_phases)
|
||||
|
||||
struct Mission
|
||||
|
||||
36
julia/test/inner_loop/nlp_solver.jl
Normal file
36
julia/test/inner_loop/nlp_solver.jl
Normal file
@@ -0,0 +1,36 @@
|
||||
@testset "Phase" begin
|
||||
|
||||
println("Testing NLP solver")
|
||||
|
||||
# We'll start by testing the mission_guess -> vector function
|
||||
vec = Vector(test_mg)
|
||||
@test typeof(vec) == Vector{Float64}
|
||||
|
||||
# Now we go in the other direction
|
||||
flybys = [ p.planet for p in test_mg.phases ]
|
||||
guess = Mission_Guess(vec, test_mg.sc, test_mg.start_mass, flybys)
|
||||
@test typeof(guess) == Mission_Guess
|
||||
@test guess.sc == test_mg.sc
|
||||
@test guess.start_mass == test_mg.start_mass
|
||||
@test guess.launch_date == test_mg.launch_date
|
||||
@test guess.launch_v∞ == test_mg.launch_v∞
|
||||
@test guess.converged == test_mg.converged
|
||||
for i in 1:length(guess.phases)
|
||||
@test guess.phases[i].planet == test_mg.phases[i].planet
|
||||
@test guess.phases[i].v∞_in == test_mg.phases[i].v∞_in
|
||||
@test guess.phases[i].v∞_out == test_mg.phases[i].v∞_out
|
||||
@test guess.phases[i].tof == test_mg.phases[i].tof
|
||||
@test guess.phases[i].thrust_profile == test_mg.phases[i].thrust_profile
|
||||
end
|
||||
|
||||
# Now we test an example run of the basic "inner function"
|
||||
leave, arrive = DateTime(1992,11,19), DateTime(1993,4,1)
|
||||
earth_state = state(Earth, leave)
|
||||
venus_state = state(Venus, arrive)
|
||||
v∞_out, v∞_in, tof = Thesis.lamberts(Earth, Venus, leave, arrive)
|
||||
phase = Phase(Venus, v∞_in, v∞_in, tof, zeros(20,3))
|
||||
guess = Mission_Guess(bepi, 12_000., leave, v∞_out, [phase])
|
||||
g = solve_mission(guess)
|
||||
println(g)
|
||||
|
||||
end
|
||||
@@ -1,74 +0,0 @@
|
||||
@testset "Phase" begin
|
||||
|
||||
println("Testing NLP solver")
|
||||
|
||||
using SNOW, PlotlyJS
|
||||
|
||||
# First we'll test an Earth case
|
||||
# Initial Setup
|
||||
T = rand( 8hour : second : day )
|
||||
revs = 8
|
||||
n = 10 * revs
|
||||
start_mass = 10_000.
|
||||
|
||||
# A simple orbit raising
|
||||
start = gen_orbit(T, start_mass, Earth)
|
||||
thrust = spiral(0.9, n, start, test_sc, revs*T, Earth)
|
||||
final = prop(thrust, start, test_sc, revs*T, Earth)[2]
|
||||
new_T = 2π * √( xyz_to_oe(final, Earth.μ)[1]^3 / Earth.μ )
|
||||
|
||||
# This should be close enough to converge, but Earth orbits are picky
|
||||
thrust_guess = spiral(0.89, n, start, test_sc, revs*T, Earth)
|
||||
guess_final = prop(thrust_guess, start, test_sc, revs*T, Earth)[2]
|
||||
result = solve_phase(start, final, test_sc, revs*T, thrust_guess, Earth, verbose=true)
|
||||
calc_path, calc_final = prop(result.sol, start, test_sc, revs*T, Earth, interpolate=true)
|
||||
|
||||
# Test
|
||||
@test norm(guess_final[1:6] - final[1:6]) > 1e-4
|
||||
@test result.converged
|
||||
@test norm(calc_final[1:6] - final[1:6]) < 1e-4
|
||||
|
||||
# Plot
|
||||
paths = Pathlist()
|
||||
push!(paths, prop(start, T, Earth),
|
||||
calc_path,
|
||||
prop(calc_final, T, Earth),
|
||||
prop(final, T, Earth) )
|
||||
fig = plot_orbits(paths, Earth,
|
||||
labels=["init", "transit", "post-transit", "final"],
|
||||
colors=["#FFF","#F44","#4F4","#44F"])
|
||||
savefig(fig, "../plots/nlp_test_earth.html")
|
||||
|
||||
# Now the sun case
|
||||
T = rand( 0.5year : hour : 1.5year )
|
||||
tof = 0.7T
|
||||
n = 20
|
||||
start_mass = 12_000.
|
||||
|
||||
# A simple orbit raising again, to keep things simple
|
||||
start = gen_orbit(T, start_mass)
|
||||
thrust = spiral(0.6, n, start, bepi, tof)
|
||||
final = prop(thrust, start, bepi, tof)[2]
|
||||
new_T = 2π * √( xyz_to_oe(final, Sun.μ)[1]^3 / Sun.μ )
|
||||
|
||||
# This should be close enough to converge
|
||||
thrust_guess = spiral(0.55, n, start, bepi, tof)
|
||||
guess_final = prop(thrust_guess, start, bepi, tof)[2]
|
||||
result = solve_phase(start, final, bepi, tof, thrust_guess, verbose=true)
|
||||
calc_path, calc_final = prop(result.sol, start, bepi, tof, interpolate=true)
|
||||
|
||||
# Test
|
||||
@test norm(guess_final[1:6] - final[1:6]) > 1e-4
|
||||
@test result.converged
|
||||
@test norm(calc_final[1:6] - final[1:6]) < 1e-4
|
||||
|
||||
# Plot
|
||||
paths = Pathlist()
|
||||
push!(paths, prop(start, T), calc_path, prop(calc_final, T), prop(final, T) )
|
||||
fig = plot_orbits(paths,
|
||||
labels=["init", "transit", "post-transit", "final"],
|
||||
colors=["#FFF","#F44","#4F4","#44F"])
|
||||
savefig(fig, "../plots/nlp_test_sun.html")
|
||||
|
||||
|
||||
end
|
||||
@@ -14,9 +14,9 @@ catch
|
||||
end
|
||||
|
||||
@testset "All Tests" begin
|
||||
include("plotting.jl")
|
||||
include("inner_loop/laguerre-conway.jl")
|
||||
include("inner_loop/propagator.jl")
|
||||
include("inner_loop/phase.jl")
|
||||
include("inner_loop/monotonic_basin_hopping.jl")
|
||||
# include("plotting.jl")
|
||||
# include("inner_loop/laguerre-conway.jl")
|
||||
# include("inner_loop/propagator.jl")
|
||||
include("inner_loop/nlp_solver.jl")
|
||||
# include("inner_loop/monotonic_basin_hopping.jl")
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user