From b29afbce04414b56eb4827edb4832f57d6498284 Mon Sep 17 00:00:00 2001 From: Connor Date: Sun, 10 Oct 2021 20:53:14 -0600 Subject: [PATCH] NLP seems to be doing pretty well --- .gitlab-ci.yml | 13 +- julia/src/Thesis.jl | 8 +- julia/src/bodies.jl | 11 +- julia/src/constants.jl | 21 +- julia/src/errors.jl | 15 +- julia/src/inner_loop/nlp_solver.jl | 99 +++++---- julia/src/inner_loop/propagator.jl | 17 +- julia/src/mission.jl | 84 ++++++-- julia/src/plotting.jl | 315 +++++++++++++++++++++++----- julia/src/spacecraft.jl | 2 +- julia/test/inner_loop/nlp_solver.jl | 64 +++--- julia/test/plotting.jl | 67 ++++-- julia/test/runtests.jl | 6 +- 13 files changed, 524 insertions(+), 198 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 9c18e0d..61997c1 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -12,8 +12,13 @@ unit-test-job: - julia --project=julia/Project.toml -E 'using Pkg; Pkg.test()' artifacts: paths: - - julia/plots/plot_test_earth.html - - julia/plots/plot_test_sun.html - - julia/plots/nlp_test_earth.html - - julia/plots/nlp_test_sun.html + - julia/plots/nlp_test_1_phase.html + - julia/plots/nlp_test_2_phase.html + - julia/plots/plot_test_mission.html + - julia/plots/plot_test_mission_guess.html + - julia/plots/plot_test_path.html + - julia/plots/plot_test_paths.html + - julia/plots/plot_test_planet_move.html + - julia/plots/plot_test_planet_no_move.html + - julia/plots/plot_test_planet_solar_system.html expire_in: 10 weeks diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl index 7f51b72..631ef88 100644 --- a/julia/src/Thesis.jl +++ b/julia/src/Thesis.jl @@ -19,14 +19,14 @@ module Thesis include("./bodies.jl") include("./constants.jl") include("./spacecraft.jl") - include("./mission.jl") - include("./conversions.jl") - include("./lamberts.jl") - include("./plotting.jl") include("./inner_loop/laguerre-conway.jl") include("./inner_loop/propagator.jl") + include("./mission.jl") + include("./conversions.jl") + include("./plotting.jl") include("./inner_loop/nlp_solver.jl") # include("./inner_loop/monotonic_basin_hopping.jl") # include("./outer_loop.jl") + include("./lamberts.jl") end diff --git a/julia/src/bodies.jl b/julia/src/bodies.jl index 6cf6f48..4d55379 100644 --- a/julia/src/bodies.jl +++ b/julia/src/bodies.jl @@ -1,17 +1,22 @@ -export state +export state, period struct Body μ::Float64 r::Float64 # radius color::String id::Int # SPICE id + a::Float64 # semi-major axis + line_color::AbstractString + name::AbstractString end -function state(p::Body, t::DateTime, add_v∞::Vector{Float64}=[0., 0., 0.], add_mass::Float64=1e10) +period(b::Body) = 2π * √(b.a^3 / Sun.μ) + +function state(p::Body, t::DateTime, add_v∞::Vector{T}=[0., 0., 0.], add_mass::Float64=1e10) where T <: Real time = utc2et(Dates.format(t,"yyyy-mm-ddTHH:MM:SS")) [ spkssb(p.id, time, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ] end -function state(p::Body, t::Float64, add_v∞::Vector{Float64}=[0., 0., 0.], add_mass::Float64=1e10) +function state(p::Body, t::S, add_v∞::Vector{T}=[0., 0., 0.], add_mass::Float64=1e10) where {T,S <: Real} [ spkssb(p.id, t, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ] end diff --git a/julia/src/constants.jl b/julia/src/constants.jl index 7f07e85..c65e2a7 100644 --- a/julia/src/constants.jl +++ b/julia/src/constants.jl @@ -7,17 +7,16 @@ export Jupiter, Saturn, Uranus, Neptune, Pluto export G, AU, init_STM, hour, day, year, second export Pathlist -const Sun = Body(1.32712440018e11, 696000., "Electric", 10) -const Mercury = Body(2.2032e4, 2439., "heat", 1) -const Venus = Body(3.257e5, 6052., "turbid", 2) -const Earth = Body(3.986004415e5, 6378.1363, "Blues", 399) -const Moon = Body(4.902799e3, 1738., "Greys", 301) -const Mars = Body(4.305e4, 3397.2, "Reds", 4) -const Jupiter = Body(1.266865361e8, 71492., "solar", 5) -const Saturn = Body(3.794e7, 60268., "turbid", 6) -const Uranus = Body(5.794e6, 25559., "haline", 7) -const Neptune = Body(6.809e6, 24764., "ice", 8) -const Pluto = Body(9e2, 1151., "matter", 9) +const Sun = Body(1.32712440018e11, 696000., "YlOrRd", 10, 0., "", "Sun") +const Mercury = Body(2.2032e4, 2439., "YlOrRd", 1, 57909083., "#FF2299", "Mercury") +const Venus = Body(3.257e5, 6052., "Portland", 2, 108208601., "#FFCC00", "Venus") +const Earth = Body(3.986004415e5, 6378.1363, "Earth", 399, 149598023., "#0055FF", "Earth") +const Mars = Body(4.305e4, 3397.2, "Reds", 4, 227939186., "#FF0000", "Mars") +const Jupiter = Body(1.266865361e8, 71492., "YlOrRd", 5, 778298361., "#FF7700", "Jupiter") +const Saturn = Body(3.794e7, 60268., "Portland", 6, 1429394133., "#FFFF00", "Saturn") +const Uranus = Body(5.794e6, 25559., "YlGnBu", 7, 2875038615., "#0000FF", "Uranus") +const Neptune = Body(6.809e6, 24764., "Blues", 8, 4504449769., "#6666FF", "Neptune") +const Pluto = Body(9e2, 1151., "Picnic", 9, 5915799000., "#FFCCCC", "Pluto") const G = 6.67430e-20 #universal gravity parameter const AU = 149597870.691 #km diff --git a/julia/src/errors.jl b/julia/src/errors.jl index 0b4f9ad..a3561a2 100644 --- a/julia/src/errors.jl +++ b/julia/src/errors.jl @@ -15,7 +15,10 @@ struct V∞_Error <: Exception v∞_in::AbstractVector v∞_out::AbstractVector end -Base.showerror(io::IO, e::V∞_Error) = print(io, "v∞s weren't the same: ", e.v∞_in, " and ", e.v∞_out) +Base.showerror(io::IO, e::V∞_Error) = print(io, """v∞s weren't the same: + v1 = $(e.v∞_in) and + v2 = $(e.v∞_out) + difference: $(norm(e.v∞_in) - norm(e.v∞_out))""") struct PropOne_Error <: Exception ΔV_unit::AbstractVector @@ -31,3 +34,13 @@ 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") + +struct Planet_Match_Error{T} <: Exception where T <: Real + state::Vector{T} + planet::Vector{T} +end +Base.showerror(io::IO, e::Planet_Match_Error) = print(io, """" + Spacecraft didn't make it to planet + spacecraft_pos = $(e.state) + planet_pos = $(e.planet) + difference = $(abs(norm(e.state - e.planet)))""") diff --git a/julia/src/inner_loop/nlp_solver.jl b/julia/src/inner_loop/nlp_solver.jl index eace95e..cad6fd1 100644 --- a/julia/src/inner_loop/nlp_solver.jl +++ b/julia/src/inner_loop/nlp_solver.jl @@ -5,19 +5,36 @@ export solve_mission struct Result converged::Bool info::Symbol - sol::Matrix{Float64} + sol::Mission_Guess end +function constraint_bounds(guess::Mission_Guess) + low_constraint = Vector{Float64}() + high_constraint = Vector{Float64}() + + for phase in guess.phases + if phase != guess.phases[end] + push!(low_constraint, -100., -100., -100., -0.005, -0.005, -0.005, -30., 100.) + push!(high_constraint, 100., 100., 100., 0.005, 0.005, 0.005, 30., 100_000.) + else + push!(low_constraint, -100., -100., -100., 0.) + push!(high_constraint, 100., 100., 100., guess.start_mass - guess.sc.dry_mass) + end + end + return low_constraint, high_constraint + +end """ This function will take an initial guess to a mission profile and leverage Ipopt to solve the mission + +NOTE: This function will output a proper mission if it succeeded and just re-return the original +guess otherwise """ -function solve_mission( guess::Mission_Guess; - tol=1e-8, +function solve_mission( guess::Mission_Guess, + launch_window::Vector{DateTime}, + latest_arrival::DateTime; + tol=1e-10, 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) @@ -52,59 +69,39 @@ function solve_mission( guess::Mission_Guess; 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) + g[1+8i:6+8i] .= (final[1:6] .- goal[1:6]) .* [1., 1., 1., 10_000., 10_000., 10_000.] + g[7+8i] = (norm(v∞_out) - norm(v∞_in)) * 10_000. δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) ) - g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) + g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) - current_planet.r else - g[8length(flybys)+1:8length(flybys)+3] .= final[1:3] .- goal[1:3] - g[8length(flybys)+4] = final[7] + g[8*(length(flybys)-1)+1:8*(length(flybys)-1)+3] .= final[1:3] .- goal[1:3] + g[8*(length(flybys)-1)+4] = final[7] - guess.sc.dry_mass end start = state(current_planet, time, v∞_out, final[7]) i += 1 - return 1.0 end + return 1.0 end - g = zeros(8length(flybys)+4) - optimizer!(g,x0) - return g + max_time = Dates.datetime2unix(latest_arrival) - Dates.datetime2unix(launch_window[1]) + lower_x = lowest_mission_vector(launch_window, length(guess.phases), n) + upper_x = highest_mission_vector(launch_window, max_time, length(guess.phases), n) + num_constraints = 8*(length(guess.phases)-1) + 4 + g_low, g_high = constraint_bounds(guess) + ipopt_options = Dict("constr_viol_tol" => tol, + "acceptable_constr_viol_tol" => 100tol, + "max_iter" => 10_000, + "max_cpu_time" => 300., + "print_level" => 5) + options = Options(solver=IPOPT(ipopt_options), derivatives=ForwardFD()) + x, _, info = minimize(optimizer!, x0, num_constraints, lower_x, upper_x, g_low, g_high, options) + if verbose println(info) end - # 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 + if info in [:Solve_Succeeded, :Solved_To_Acceptable_Level] + return Mission(x, guess.sc, guess.start_mass, flybys) + else + return guess + end end diff --git a/julia/src/inner_loop/propagator.jl b/julia/src/inner_loop/propagator.jl index 1acc502..0e9d282 100644 --- a/julia/src/inner_loop/propagator.jl +++ b/julia/src/inner_loop/propagator.jl @@ -56,7 +56,20 @@ function prop(ΔVs::Matrix{T}, for k in 1:7 push!(states[k], interpolated_state[k]) end end end - state = prop_one(ΔVs[i,:], state, craft, time/n, primary) + try + state = prop_one(ΔVs[i,:], state, craft, time/n, primary) + catch e + if isa(e,PropOne_Error) + for val in e.ΔV_unit + # If this isn't true, then we just let it slide + if abs(val) > 1.0001 + rethrow() + end + end + else + rethrow() + end + end for j in 1:7 push!(states[j], state[j]) end state[7] >= craft.dry_mass || throw(Mass_Error(state[7])) end @@ -68,4 +81,4 @@ end """ Convenience function for propagating a state with no thrust """ -prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1e10], no_thrust, t, p)[1] +prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1.], no_thrust, t, p)[1] diff --git a/julia/src/mission.jl b/julia/src/mission.jl index 65ab7c5..ee607c0 100644 --- a/julia/src/mission.jl +++ b/julia/src/mission.jl @@ -1,6 +1,6 @@ export Phase, Mission_Guess, Mission, Bad_Mission export test_phase1, test_phase2 -export test_mg, test_mission +export test_mg export Vector struct Phase @@ -11,7 +11,7 @@ struct Phase thrust_profile::Matrix{Float64} end -const test_phase1 = Phase(Venus, [9., 10., 0.], [10., 9, 0.], 1.30464e7, zeros(20,3)) +const test_phase1 = Phase(Venus, [√20., -√80., 0.], [√50., -√50., 0.], 1.30464e7, zeros(20,3)) const test_phase2 = Phase(Jupiter, [0.3, 7.1, 0.2], [0.3, 7.1, 0.2], 3.9year, zeros(20,3)) const test_phases = [test_phase1, test_phase2] @@ -46,12 +46,6 @@ function Mission_Guess(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float mass_used += mass_consumption(sc, phase) mass - mass_used > sc.dry_mass || throw(Mass_Error(mass - mass_used)) v∞_in, v∞_out = phase.v∞_in, phase.v∞_out - if phase != phases[end] - norm(v∞_in) ≈ norm(v∞_out) || throw(V∞_Error(v∞_in, v∞_out)) - δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) ) - periapsis = (phase.planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) - periapsis > 1.1phase.planet.r || throw(HitPlanet_Error()) - end end Mission_Guess(sc, mass, date, v∞, phases, false) end @@ -101,6 +95,32 @@ function Base.Vector(g::Mission_Guess) return result end +function lowest_mission_vector(launch_window::Vector{DateTime}, num_phases::Int, n::Int) + result = Vector{Float64}() + push!(result, Dates.datetime2unix(launch_window[1])) + push!(result, -10*ones(3)...) + for i in 1:num_phases + push!(result, -10*ones(3)...) + push!(result, -10*ones(3)...) + push!(result, 86_400.0) + push!(result, -1 * ones(n,3)...) + end + return result +end + +function highest_mission_vector(launch_window::Vector{DateTime}, mission_length::Float64, num_phases::Int, n::Int) + result = Vector{Float64}() + push!(result, Dates.datetime2unix(launch_window[2])) + push!(result, 10*ones(3)...) + for i in 1:num_phases + push!(result, 10*ones(3)...) + push!(result, 10*ones(3)...) + push!(result, mission_length) + push!(result, ones(n,3)...) + 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 @@ -113,27 +133,47 @@ struct Mission end """ -Constructor for a mission. Generally mission guesses are converged +Constructor for a mission. Generally mission guesses are converged, so we check everything """ -function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, phases::Vector{Phase}) +function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, phases::Vector{Phase}; + force=false) # First do some checks to make sure that it's valid - mass_used = 0 - for phase in phases - mass_used += mass_consumption(sc, phase) - mass - mass_used > sc.dry_mass || throw(Mass_Error(mass - mass_used)) - v∞_in, v∞_out = phase.v∞_in, phase.v∞_out - if phase != phases[end] - norm(v∞_in) ≈ norm(v∞_out) || throw(V∞_Error(v∞_in, v∞_out)) - δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) ) - periapsis = (phase.planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) - periapsis > 1.1phase.planet.r || throw(HitPlanet_Error()) + if !force + time = date + current_planet = Earth + start = state(current_planet, time, v∞, mass) + for phase in phases + final = prop(phase.thrust_profile, start, sc, phase.tof)[2] + mass = final[7] + mass > sc.dry_mass || throw(Mass_Error(mass - mass_used)) + current_planet = phase.planet + time += Dates.Second(floor(phase.tof)) + p_pos = state(current_planet, time)[1:3] + abs(norm(final[1:3] - p_pos)) < √30_000. || throw(Planet_Match_Error(final[1:3], p_pos)) + v∞_in, v∞_out = phase.v∞_in, phase.v∞_out + if phase != phases[end] + abs(norm(v∞_in) - norm(v∞_out)) < 0.005 || throw(V∞_Error(v∞_in, v∞_out)) + δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) ) + periapsis = (phase.planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) + periapsis > 1.1phase.planet.r || throw(HitPlanet_Error()) + end end end Mission(sc, mass, date, v∞, phases, true) end -const test_mission = Mission(bepi, 12_000., DateTime(1992,11,19), [-3.4,1.2,0.1], test_phases) - +""" +BE CAREFUL!! This just makes a guess converged, whether true or not +""" +function Mission(g::Mission_Guess) + return Mission(g.sc, g.start_mass, g.launch_date, g.launch_v∞, g.phases, force=true) +end + +function Mission(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body}) + guess = Mission_Guess(x, sc, mass, flybys::Vector{Body}) + return Mission(guess) +end + struct Bad_Mission message::String converged::Bool diff --git a/julia/src/plotting.jl b/julia/src/plotting.jl index 18b5856..d16b606 100644 --- a/julia/src/plotting.jl +++ b/julia/src/plotting.jl @@ -3,7 +3,7 @@ using Random, PlotlyJS, Base.Iterators -export plot_orbits +export plot function random_color() num1 = rand(0:255) @@ -12,60 +12,281 @@ function random_color() return "#"*string(num1, base=16, pad=2)*string(num2, base=16, pad=2)*string(num3, base=16, pad=2) end -function get_true_max(mat::Vector{Array{Float64,2}}) - return maximum(abs.(flatten(mat))) +# Some convenience functions for maximums +""" +Gets the maximum absolute value of a path +""" +function Base.maximum(path::Vector{Vector{Float64}}) + return maximum( abs.( reduce( vcat,path ) ) ) end -function plot_orbits(paths::Vector{Vector{Vector{Float64}}}, - primary::Body=Sun; - labels::Vector{String}=Vector{String}(), - title::String="Spacecraft Position", - colors::Vector{String}=Vector{String}()) +""" +Gets the maximum absolute value of a list of paths +""" +function Base.maximum(paths::Vector{Vector{Vector{Float64}}}) + return maximum( abs.( reduce( vcat,reduce( vcat,paths ) ) ) ) +end +""" +Gets the maximum absolute value of a surface +""" +function Base.maximum(xs::Matrix{Float64},ys::Matrix{Float64},zs::Matrix{Float64}) + return maximum([ vec(xs), vec(ys), vec(zs) ]) +end + +""" +Basically rewriting Julia at this point +""" +Base.getindex(x::Nothing, i::Int) = nothing + +""" +Generates a layout that works for most plots. Requires a title and a limit. +""" +function standard_layout(limit::Float64, title::AbstractString) + limit *= 1.1 + Layout(title=title, + # width=1400, + # height=800, + paper_bgcolor="rgba(5,10,40,1.0)", + plot_bgcolor="rgba(100,100,100,0.01)", + scene = attr(xaxis = attr(autorange = false,range=[-limit,limit]), + yaxis = attr(autorange = false,range=[-limit,limit]), + zaxis = attr(autorange = false,range=[-limit,limit]), + aspectratio=attr(x=1,y=1,z=1), + aspectmode="manual")) +end + +""" +A fundamental function to generate a plot for a single orbit path +""" +function gen_plot(path::Vector{Vector{Float64}}; + label::Union{AbstractString, Nothing} = nothing, + color::AbstractString = random_color(), + markers=true) + + # Plot the orbit + if label === nothing + path_trace = scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]), + mode="lines", showlegend=false, line_color=color, line_width=3) + else + path_trace = scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]), + mode="lines", name=label, line_color=color, line_width=3) + end + traces = [ path_trace ] + + # Plot the markers + if markers + push!(traces, scatter3d(;x=([path[1][1]]),y=([path[2][1]]),z=([path[3][1]]), + mode="markers", showlegend=false, + marker=attr(color=color, size=3, symbol="circle"))) + push!(traces, scatter3d(;x=([path[1][end]]),y=([path[2][end]]),z=([path[3][end]]), + mode="markers", showlegend=false, + marker=attr(color=color, size=3, symbol="diamond"))) + end + + return traces, maximum(path[1:6]) + +end + +""" +A fundamental function to plot a single orbit path +""" +function plot(path::Vector{Vector{Float64}}; + color::AbstractString = random_color(), + title::AbstractString = "Single Orbit Path", + markers=true) + + # Generate the plot + p, limit = gen_plot(path, color=color, markers=markers) + + # Now use the standard layout + layout = standard_layout(limit, title) + + # Generate the plot + PlotlyJS.plot( PlotlyJS.Plot( p, layout ) ) + +end + +""" +This will generate the plot of a planetary body in it's correct colors and position in a +sun-centered frame at a certain time + +If the body is the sun, no time is necessary + +An optional prop time is allowed, which will also provide an orbit +""" +function gen_plot(b::Body, + time::Union{DateTime, Nothing}=nothing, + prop_time::Union{Float64, Nothing}=nothing; + scale::Int=100) + + # Gotta be a little weird with Sun Scale + if b == Sun && scale > 25 scale = 25 end + + # This is how many points to poll to generate the sphere N = 32 θ = collect(range(0,length=N,stop=2π)) ϕ = collect(range(0,length=N,stop=π)) - x = cos.(θ) * sin.(ϕ)' - y = sin.(θ) * sin.(ϕ)' - z = repeat(cos.(ϕ)',outer=[N, 1]) - ps = primary.r .* (x,y,z) - x_p,y_p,z_p = ps - t1 = [] - for i = 1:length(paths) - path = paths[i] - label = labels != [] ? labels[i] : "orbit" - color = colors != [] ? colors[i] : random_color() - push!(t1, scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]), - mode="lines", name=label, line_color=color, line_width=3)) - push!(t1, scatter3d(;x=([path[1][1]]),y=([path[2][1]]),z=([path[3][1]]), - mode="markers", showlegend=false, - marker=attr(color=color, size=3, symbol="circle"))) - push!(t1, scatter3d(;x=([path[1][end]]),y=([path[2][end]]),z=([path[3][end]]), - mode="markers", showlegend=false, - marker=attr(color=color, size=3, symbol="diamond"))) + # Convert polar -> cartesian + unscaled_xs = cos.(θ) * sin.(ϕ)' + unscaled_ys = sin.(θ) * sin.(ϕ)' + unscaled_zs = repeat(cos.(ϕ)',outer=[N, 1]) + + # Scale by the size of the body (and a scaling factor for visualization) + xs,ys,zs = scale * b.r .* (unscaled_xs,unscaled_ys,unscaled_zs) + + # Add an offset if we're not plotting the sun + body_state = zeros(6) + if b != Sun body_state = state(b, time) end + + max_value = maximum(xs .+ body_state[1], ys .+ body_state[2], zs .+ body_state[3]) + + # Generate the surface trace + body_surface = surface(;x = xs .+ body_state[1], + y = ys .+ body_state[2], + z = zs .+ body_state[3], + showscale = false, + colorscale = b.color) + outputs = [ body_surface ] + + if prop_time !== nothing + p, maybe_max = gen_plot(prop(body_state, prop_time), label=b.name, color=b.line_color, markers=false) + push!(outputs, p...) + if maybe_max > max_value max_value = maybe_max end end - limit = max(maximum(abs.(flatten(flatten(paths)))), - maximum(abs.(flatten(ps)))) * 1.1 - t2 = surface(;x=(x_p), - y=(y_p), - z=(z_p), - showscale=false, - colorscale = primary.color) - - layout = Layout(title=title, - width=1000, - height=600, - paper_bgcolor="rgba(5,10,40,1.0)", - plot_bgcolor="rgba(100,100,100,0.01)", - scene = attr(xaxis = attr(autorange = false,range=[-limit,limit]), - yaxis = attr(autorange = false,range=[-limit,limit]), - zaxis = attr(autorange = false,range=[-limit,limit]), - aspectratio=attr(x=1,y=1,z=1), - aspectmode="manual")) - - p = Plot([t1...,t2],layout) - plot(p) + return outputs, max_value + +end + +""" +This will plot a planetary body in it's correct colors and position in a +sun-centered frame at a certain time + +If the body is the sun, no time is necessary + +An optional prop time is allowed, which will also provide an orbit +""" +function plot(b::Body, + time::Union{DateTime, Nothing}=nothing, + prop_time::Union{Float64, Nothing}=nothing; + title::AbstractString="Single Planet Plot", + scale::Int=100) + + # Generate the plot + p, limit = gen_plot(b,time,prop_time,scale=scale) + + # Now use the standard layout + layout = standard_layout(limit, title) + + # Generate the plot + PlotlyJS.plot( PlotlyJS.Plot( p, layout ) ) + +end + +""" +This will plot a vector of planetary bodies at a certain time + +The Sun is plotted automatically + +An optional orbits bool is allowed, which shows the orbits +""" +function plot(bodies::Vector{Body}, + time::Union{DateTime, Nothing}; + orbits::Bool=true, + title::AbstractString="Multiple Planet Plot", + scale::Int=100) + + # Generate the plots + trace, limit = gen_plot(Sun, scale=scale) + traces = [ trace... ] + for b in bodies + if orbits + trace, new_limit = gen_plot(b,time,period(b),scale=scale) + else + trace, new_limit = gen_plot(b,time,scale=scale) + end + push!(traces, trace...) + limit = max(limit, new_limit) + end + + # Now use the standard layout + layout = standard_layout(limit, title) + + # Generate the plot + PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) ) + +end + +""" +Generic plotting function for a bunch of paths. Useful for debugging. +""" +function plot(paths::Vector{Vector{Vector{Float64}}}; + labels::Union{Vector{String},Nothing}=nothing, + colors::Union{Vector{String},Nothing}=nothing, + markers::Bool=true, + title::String="Spacecraft Position") + + # Generate the path traces + color = colors === nothing ? random_color() : colors[1] + trace, limit = gen_plot(paths[1],label=labels[1],color=color,markers=markers) + traces = [ trace... ] + for i = 2:length(paths) + color = colors === nothing ? random_color() : colors[i] + trace, new_limit = gen_plot(paths[i],label=labels[i],color=color,markers=markers) + push!(traces, trace...) + limit = max(limit, new_limit) + end + + # Generate the sun trace + push!(traces, gen_plot(Sun)[1]...) + + # Determine layout details + layout = standard_layout(limit, title) + + # Plot + PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) ) + +end + +function plot(m::Union{Mission, Mission_Guess}; title::String="Mision Plot") + # First plot the earth + # Then plot phase plus planet phase + time = m.launch_date + mass = m.start_mass + current_planet = Earth + earth_trace, limit = gen_plot(current_planet, time, year) + start = state(current_planet, time, m.launch_v∞, mass) + traces = [ earth_trace... ] + + i = 1 + for phase in m.phases + # First do the path + path, final = prop(phase.thrust_profile, start, m.sc, phase.tof, interpolate=true) + mass = final[7] + time += Dates.Second(floor(phase.tof)) + current_planet = phase.planet + start = state(current_planet, time, phase.v∞_out, mass) + path_trace, new_limit = gen_plot(path, label="Phase "*string(i)) + push!(traces, path_trace...) + limit = max(limit, new_limit) + + # Then the planet + planet_trace, new_limit = gen_plot(current_planet, time, -phase.tof) + push!(traces, planet_trace...) + limit = max(limit, new_limit) + + i += 1 + end + + # Generate the sun trace + push!(traces, gen_plot(Sun)[1]...) + + # Determine layout details + layout = standard_layout(limit, title) + + # Plot + PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) ) end diff --git a/julia/src/spacecraft.jl b/julia/src/spacecraft.jl index 25620d3..4e8d502 100644 --- a/julia/src/spacecraft.jl +++ b/julia/src/spacecraft.jl @@ -10,4 +10,4 @@ end const test_sc = Sc(8000., 0.00025/(2000*0.00981), 0.00025, 50, 0.9) const bepi = Sc(8000., 2*0.00025/(2000*0.00981), 0.00025, 2, 0.9) -const no_thrust = Sc(8000., 0.01, 0., 0, 0.) +const no_thrust = Sc(0., 0.01, 0., 0, 0.) diff --git a/julia/test/inner_loop/nlp_solver.jl b/julia/test/inner_loop/nlp_solver.jl index 49787f1..242c7c9 100644 --- a/julia/test/inner_loop/nlp_solver.jl +++ b/julia/test/inner_loop/nlp_solver.jl @@ -1,36 +1,48 @@ @testset "Phase" begin + using PlotlyJS: savefig + 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" + # Test the optimizer for a one-phase mission + # The lambert's solver said this should be pretty valid + launch_window = [DateTime(1992,11,1), DateTime(1992,12,1)] + latest_arrival = DateTime(1993,6,1) leave, arrive = DateTime(1992,11,19), DateTime(1993,4,1) + test_leave = DateTime(1992,11,12) 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) + # We can get the thrust profile and tof pretty wrong and still be ok + phase = Phase(Venus, 1.1v∞_in, v∞_in, 0.9*tof, 0.1*ones(20,3)) + guess = Mission_Guess(bepi, 12_000., test_leave, 0.9*v∞_out, [phase]) + m = solve_mission(guess, launch_window, latest_arrival, verbose=true) + @test typeof(m) == Mission + @test m.converged == true + + # Now we can plot the results to check visually + p = plot(m, title="NLP Test Solution") + savefig(p,"../plots/nlp_test_1_phase.html") + + # Now we can look at a more complicated trajectory + flybys = [Earth, Venus, Mars] + launch_window = [DateTime(2021,10,1), DateTime(2021,12,1)] + latest_arrival = DateTime(2023,1,1) + dates = [DateTime(2021,11,1), DateTime(2022,3,27), DateTime(2022,8,28)] + phases = Vector{Phase}() + launch_v∞, _, tof1 = Thesis.lamberts(flybys[1], flybys[2], dates[1], dates[2]) + println(launch_v∞) + for i in 1:length(dates)-2 + v∞_out1, v∞_in1, tof1 = Thesis.lamberts(flybys[i], flybys[i+1], dates[i], dates[i+1]) + v∞_out2, v∞_in2, tof2 = Thesis.lamberts(flybys[i+1], flybys[i+2], dates[i+1], dates[i+2]) + push!(phases, Phase(flybys[i+1], v∞_in1, v∞_out2, tof1, 0.01*ones(20,3))) + end + v∞_out, v∞_in, tof = Thesis.lamberts(flybys[end-1], flybys[end], dates[end-1], dates[end]) + push!(phases, Phase(flybys[end], v∞_in, v∞_in, tof, 0.01*ones(20,3))) + # This isn't quite right. V∞ discrepancy is too big + guess = Mission_Guess(bepi, 12_000., dates[1], launch_v∞, phases) + m = solve_mission(guess, launch_window, latest_arrival, verbose=true) + p = plot(m, title="NLP Test Solution (2 Phases)") + savefig(p,"../plots/nlp_test_2_phase.html") end diff --git a/julia/test/plotting.jl b/julia/test/plotting.jl index 739a27e..da305ff 100644 --- a/julia/test/plotting.jl +++ b/julia/test/plotting.jl @@ -2,34 +2,55 @@ println("Testing plotting features") - using PlotlyJS + using PlotlyJS: savefig, SyncPlot - # Plot an earth plot - T = rand(2hour : 1 : 4hour) - revs = 30 - n = revs*100 - - start = gen_orbit(T, 12_000., Earth) - thrust = spiral(0.9, n, start, test_sc, revs*T, Earth) - path = prop(thrust, start, test_sc, revs*T, Earth)[1] - - p = plot_orbits([path], Earth) - savefig(p,"../plots/plot_test_earth.html") - @test typeof(p) == PlotlyJS.SyncPlot - - # Now change a little bit and plot around the Sun - # This also checks that the spacecraft are configured right: - # They really shouldn't run out of fuel in 4 years + # Generate a couple of useful paths T = rand(year : hour : 4year) tof = 4year + n = 4_000 start = gen_orbit(T, 12_000.) thrust = spiral(0.9, n, start, bepi, tof) - sun_paths = Vector{Vector{Vector{Float64}}}() - push!(sun_paths, prop(zeros(100,3), start, bepi, tof)[1]) - push!(sun_paths, prop(thrust, start, bepi, tof)[1]) - p = plot_orbits(sun_paths) - savefig(p,"../plots/plot_test_sun.html") - @test typeof(p) == PlotlyJS.SyncPlot + spiral_path = prop(thrust, start, bepi, tof)[1] + no_thrust_path = prop(zeros(n,3), start, bepi, tof)[1] + paths = [ spiral_path, no_thrust_path ] + # First let's plot a single basic path + p = plot(spiral_path) + savefig(p,"../plots/plot_test_path.html") + @test typeof(p) == SyncPlot + + # Next I just want to plot a planet both not moving... + p = plot(Jupiter, Dates.now(), title="Jupiter") + savefig(p,"../plots/plot_test_planet_nomove.html") + @test typeof(p) == SyncPlot + + # ...and moving + p = plot(Venus, Dates.now(), 0.95period(Venus), title="Venus") + savefig(p,"../plots/plot_test_planet_move.html") + @test typeof(p) == SyncPlot + + # Next I want to plot the entire solar system + p = plot([Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto], + Dates.now(), + title="The Solar System") + savefig(p,"../plots/plot_test_solar_system.html") + @test typeof(p) == SyncPlot + + # Next I want to plot a vector of states + p = plot(paths, + labels = ["No Thrust Orbit","Thrust Spiral"], + title = "Vector of Paths Plot Test") + savefig(p,"../plots/plot_test_paths.html") + @test typeof(p) == SyncPlot + + # And finally, I'd like to plot a mission guess... + p = plot(test_mg) + savefig(p,"../plots/plot_test_mission_guess.html") + @test typeof(p) == SyncPlot + + # ...and a mission + p = plot(Mission(test_mg)) + savefig(p,"../plots/plot_test_mission.html") + @test typeof(p) == SyncPlot end diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 841f503..fba9d59 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -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("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