From 27702ba3f8207dc593edb4cac0a4e935106e4ef1 Mon Sep 17 00:00:00 2001 From: Connor Date: Sat, 25 Sep 2021 17:32:11 -0600 Subject: [PATCH] Added a sun test case and interpolation --- .gitlab-ci.yml | 3 ++- julia/src/errors.jl | 3 +++ julia/src/inner_loop/phase.jl | 3 ++- julia/src/inner_loop/propagator.jl | 36 +++++++++---------------- julia/test/inner_loop/phase.jl | 42 +++++++++++++++++++++++++----- notes.md | 6 +++++ 6 files changed, 62 insertions(+), 31 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 8bad5ca..9c18e0d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -14,5 +14,6 @@ unit-test-job: paths: - julia/plots/plot_test_earth.html - julia/plots/plot_test_sun.html - - julia/plots/nlp_test.html + - julia/plots/nlp_test_earth.html + - julia/plots/nlp_test_sun.html expire_in: 10 weeks diff --git a/julia/src/errors.jl b/julia/src/errors.jl index fd41b51..78aa617 100644 --- a/julia/src/errors.jl +++ b/julia/src/errors.jl @@ -2,6 +2,9 @@ struct LaGuerreConway_Error <: Exception end struct ΔVsize_Error <: Exception end +struct Convergence_Error <: Exception end +Base.showerror(io::IO, e::Convergence_Error) = print(io, "NLP solver didn't converge...") + struct GenOrbit_Error <: Exception end Base.showerror(io::IO, e::GenOrbit_Error) = print(io, "Infinite Loop trying to generate the init orbit") diff --git a/julia/src/inner_loop/phase.jl b/julia/src/inner_loop/phase.jl index a8c364f..656c34d 100644 --- a/julia/src/inner_loop/phase.jl +++ b/julia/src/inner_loop/phase.jl @@ -22,7 +22,7 @@ function solve_phase( start::Vector{Float64}, F[1:6, 1] .= prop(tanh.(x), start, craft, tof, primary)[2][1:6] .- final[1:6] catch e # If the error is due to something natural, just imply a penalty - if isa(Mass_Error, e) || isa(PropOne_Error, e) + if isa(e, Mass_Error) || isa(e, PropOne_Error) F .= 10000000.0 else rethrow() @@ -31,6 +31,7 @@ function solve_phase( start::Vector{Float64}, end result = nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters) + converged(result) || throw(Convergence_Error()) result.zero = tanh.(result.zero) return result diff --git a/julia/src/inner_loop/propagator.jl b/julia/src/inner_loop/propagator.jl index 763a173..b0b2d9e 100644 --- a/julia/src/inner_loop/propagator.jl +++ b/julia/src/inner_loop/propagator.jl @@ -38,40 +38,30 @@ function prop(ΔVs::Matrix{T}, state::Vector{Float64}, craft::Sc, time::Float64, - primary::Body=Sun) where T <: Real + primary::Body=Sun; + interpolate::Bool=false) where T <: Real size(ΔVs)[2] == 3 || throw(ΔVsize_Error()) n = size(ΔVs)[1] - x_states = Vector{T}() - y_states = Vector{T}() - z_states = Vector{T}() - dx_states = Vector{T}() - dy_states = Vector{T}() - dz_states = Vector{T}() - masses = Vector{T}() + states = [ Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}() ] - push!(x_states, state[1]) - push!(y_states, state[2]) - push!(z_states, state[3]) - push!(dx_states, state[4]) - push!(dy_states, state[5]) - push!(dz_states, state[6]) - push!(masses, state[7]) + for i in 1:7 push!(states[i], state[i]) end for i in 1:n + if interpolate + interpolated_state = copy(state) + for j in 1:49 + interpolated_state = [laguerre_conway(interpolated_state, time/50n, primary); 0.0] + for k in 1:7 push!(states[k], interpolated_state[k]) end + end + end state = prop_one(ΔVs[i,:], state, craft, time/n, primary) - push!(x_states, state[1]) - push!(y_states, state[2]) - push!(z_states, state[3]) - push!(dx_states, state[4]) - push!(dy_states, state[5]) - push!(dz_states, state[6]) - push!(masses, state[7]) + for j in 1:7 push!(states[j], state[j]) end state[7] >= craft.dry_mass || throw(Mass_Error(state[7])) end - return [x_states, y_states, z_states, dx_states, dy_states, dz_states, masses], state + return states, state end diff --git a/julia/test/inner_loop/phase.jl b/julia/test/inner_loop/phase.jl index e422f9c..103f378 100644 --- a/julia/test/inner_loop/phase.jl +++ b/julia/test/inner_loop/phase.jl @@ -4,11 +4,12 @@ using NLsolve, PlotlyJS + # First we'll test an Earth case # Initial Setup - T = rand( 2hour : second : 12hour) - revs = 5 - n = 10revs - start_mass = 12_000. + T = rand( 8hour : second : day ) + revs = 3 + n = 10 * revs + start_mass = 10_000. # A simple orbit raising start = gen_orbit(T, start_mass, Earth) @@ -19,7 +20,7 @@ # This should be close enough to converge thrust_guess = spiral(0.88, n, start, test_sc, revs*T, Earth) result = solve_phase(start, final, test_sc, revs*T, thrust_guess, Earth) - calc_path, calc_final = prop(result.zero, start, test_sc, revs*T, Earth) + calc_path, calc_final = prop(result.zero, start, test_sc, revs*T, Earth, interpolate=true) # Test @test converged(result) @@ -34,6 +35,35 @@ fig = plot_orbits(paths, Earth, labels=["init", "transit", "post-transit", "final"], colors=["#FFF","#F44","#4F4","#44F"]) - savefig(fig, "../plots/nlp_test.html") + savefig(fig, "../plots/nlp_test_earth.html") + + # Now the sun case + T = rand( 0.5year : hour : 1.5year ) + 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, revs*T) + final = prop(thrust, start, bepi, revs*T)[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, revs*T) + result = solve_phase(start, final, bepi, revs*T, thrust_guess) + calc_path, calc_final = prop(result.zero, start, bepi, revs*T, interpolate=true) + + # Test + @test converged(result) + @test norm(calc_final[1:6] - final[1:6]) < 1e-5 + + # 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 diff --git a/notes.md b/notes.md index 644f203..fcb0d04 100644 --- a/notes.md +++ b/notes.md @@ -28,3 +28,9 @@ I need to determine the list of decision variables for each step: - turning angle - Decided by NLP - thrust profile for each phase + +Todo: + +1. ~~Add interpolate function~~ +2. Add NLP solver sun test +3. Start on MBH