From 504d83b3901cf19001aa267b1529c133c85e0bad Mon Sep 17 00:00:00 2001 From: Connor Date: Sun, 10 Oct 2021 21:53:58 -0600 Subject: [PATCH] Small improvements --- julia/src/errors.jl | 1 + julia/src/inner_loop/laguerre-conway.jl | 4 +- julia/src/inner_loop/nlp_solver.jl | 64 +++++++++++++++---------- julia/src/inner_loop/propagator.jl | 1 - julia/src/spacecraft.jl | 2 +- julia/test/inner_loop/nlp_solver.jl | 23 ++++----- 6 files changed, 55 insertions(+), 40 deletions(-) diff --git a/julia/src/errors.jl b/julia/src/errors.jl index a3561a2..ea80e8d 100644 --- a/julia/src/errors.jl +++ b/julia/src/errors.jl @@ -1,4 +1,5 @@ struct LaGuerreConway_Error <: Exception end +Base.showerror(io::IO, e::LaGuerreConway_Error) = print(io, "LaGuerre-Conway didn't converge") struct ΔVsize_Error <: Exception end diff --git a/julia/src/inner_loop/laguerre-conway.jl b/julia/src/inner_loop/laguerre-conway.jl index b14fdce..12472a6 100644 --- a/julia/src/inner_loop/laguerre-conway.jl +++ b/julia/src/inner_loop/laguerre-conway.jl @@ -23,7 +23,7 @@ function laguerre_conway(state::Vector{<:Real}, time::Float64, primary::Body=Sun 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(LaGuerreConway_Error("LaGuerre-Conway did not converge!")) end + if i > 100 throw(LaGuerreConway_Error()) end end F = 1 - a/r0_mag * (1-cos(ΔE)) G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE) @@ -42,7 +42,7 @@ function laguerre_conway(state::Vector{<:Real}, time::Float64, primary::Body=Sun 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(LaGuerreConway_Error("LaGuerre-Conway did not converge!")) end + if i > 100 throw(LaGuerreConway_Error()) end end F = 1 - a/r0_mag * (1-cos(ΔH)) G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH) diff --git a/julia/src/inner_loop/nlp_solver.jl b/julia/src/inner_loop/nlp_solver.jl index 83735aa..82de693 100644 --- a/julia/src/inner_loop/nlp_solver.jl +++ b/julia/src/inner_loop/nlp_solver.jl @@ -34,7 +34,8 @@ function solve_mission( guess::Mission_Guess, launch_window::Vector{DateTime}, latest_arrival::DateTime; tol=1e-10, - verbose::Bool=false ) + verbose::Bool=false, + print_level=0 ) # First we define our starting point x0 = Vector(guess) @@ -58,29 +59,39 @@ function solve_mission( guess::Mission_Guess, # - 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]) .* [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 ) - current_planet.r - else - 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 + try + 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]) .* [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 ) - current_planet.r + else + 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]) + mass = final[7] + i += 1 + end + return 1.0 + catch e + if isa(e,LaGuerreConway_Error) + g[1:8*(length(flybys)-1)+4] .= 1e10 + return 1e10 + else + rethrow() end - start = state(current_planet, time, v∞_out, final[7]) - i += 1 end - return 1.0 end max_time = Dates.datetime2unix(latest_arrival) - Dates.datetime2unix(launch_window[1]) @@ -90,9 +101,9 @@ function solve_mission( guess::Mission_Guess, 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" => 0) + "max_iter" => 100_000, + "max_cpu_time" => 60., + "print_level" => print_level) options = Options(solver=IPOPT(ipopt_options), derivatives=ForwardFD()) x, _, info = minimize(optimizer!, x0, num_constraints, lower_x, upper_x, g_low, g_high, options) @@ -101,6 +112,9 @@ function solve_mission( guess::Mission_Guess, if info in [:Solve_Succeeded, :Solved_To_Acceptable_Level] return Mission(x, guess.sc, guess.start_mass, flybys) else + g = zeros(num_constraints) + optimizer!(g,x) + println(g) return guess end diff --git a/julia/src/inner_loop/propagator.jl b/julia/src/inner_loop/propagator.jl index 0e9d282..f179b72 100644 --- a/julia/src/inner_loop/propagator.jl +++ b/julia/src/inner_loop/propagator.jl @@ -71,7 +71,6 @@ function prop(ΔVs::Matrix{T}, end end for j in 1:7 push!(states[j], state[j]) end - state[7] >= craft.dry_mass || throw(Mass_Error(state[7])) end return states, state diff --git a/julia/src/spacecraft.jl b/julia/src/spacecraft.jl index 4e8d502..c647680 100644 --- a/julia/src/spacecraft.jl +++ b/julia/src/spacecraft.jl @@ -9,5 +9,5 @@ mutable struct Sc 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 bepi = Sc(2000., 2*0.00025/(2800*0.00981), 0.00025, 2, 0.9) 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 b4f98d8..765b9d7 100644 --- a/julia/test/inner_loop/nlp_solver.jl +++ b/julia/test/inner_loop/nlp_solver.jl @@ -15,10 +15,9 @@ v∞_out, v∞_in, tof = Thesis.lamberts(Earth, Venus, leave, arrive) # 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]) + guess = Mission_Guess(bepi, 3_600., 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") @@ -38,8 +37,9 @@ 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))) - guess = Mission_Guess(bepi, 12_000., dates[1], launch_v∞, phases) + guess = Mission_Guess(bepi, 3_600., dates[1], launch_v∞, phases) m = solve_mission(guess, launch_window, latest_arrival, verbose=true) + @test typeof(m) == Mission p = plot(m, title="NLP Test Solution (2 Phases)") savefig(p,"../plots/nlp_test_2_phase.html") @@ -47,23 +47,24 @@ flybys = [Earth, Venus, Earth, Mars, Earth, Jupiter] launch_window = [DateTime(2023,1,1), DateTime(2024,1,1)] latest_arrival = DateTime(2031,1,1) - dates = [DateTime(2023,6,28), - DateTime(2023,10,30), - DateTime(2024,9,2), - DateTime(2025,2,10), - DateTime(2026,11,26), - DateTime(2029,10,16)] + dates = [DateTime(2023,5,23), + DateTime(2023,10,21), + DateTime(2024,8,24), + DateTime(2025,2,13), + DateTime(2026,11,22), + DateTime(2032,1,1)] phases = Vector{Phase}() launch_v∞, _, tof1 = Thesis.lamberts(flybys[1], flybys[2], dates[1], dates[2]) 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))) + push!(phases, Phase(flybys[i+1], 1.02v∞_in1, 0.98v∞_out2, 1.02tof1, 0.02*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))) - guess = Mission_Guess(bepi, 12_000., dates[1], launch_v∞, phases) + guess = Mission_Guess(bepi, 3_600., dates[1], launch_v∞, phases) m = solve_mission(guess, launch_window, latest_arrival, verbose=true) + @test typeof(m) == Mission p = plot(m, title="NLP Test Solution (5 Phases)") savefig(p,"../plots/nlp_test_5_phase.html")