Small improvements

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
 

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)

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
+    try
-      phase_params = x[ 5 + (3n+7) * i : 5 + (3n+7) * (i+1) - 1 ]
+      for flyby in flybys
-      v∞_in = phase_params[1:3]
+        phase_params = x[ 5 + (3n+7) * i : 5 + (3n+7) * (i+1) - 1 ]
-      v∞_out = phase_params[4:6]
+        v∞_in = phase_params[1:3]
-      tof = phase_params[7]
+        v∞_out = phase_params[4:6]
-      thrusts = reshape(phase_params[8:end], (n,3))
+        tof = phase_params[7]
-      current_planet = flyby
+        thrusts = reshape(phase_params[8:end], (n,3))
-      time += tof
+        current_planet = flyby
-      goal = state(current_planet, time, v∞_in)
+        time += tof
-      final = prop(reshape(thrusts, (n,3)), start, guess.sc, tof)[2]
+        goal = state(current_planet, time, v∞_in)
-      if flyby != flybys[end]
+        final = prop(reshape(thrusts, (n,3)), start, guess.sc, tof)[2]
-        g[1+8i:6+8i] .= (final[1:6] .- goal[1:6]) .* [1., 1., 1., 10_000., 10_000., 10_000.]
+        if flyby != flybys[end]
-        g[7+8i] = (norm(v∞_out) - norm(v∞_in)) * 10_000.
+          g[1+8i:6+8i] .= (final[1:6] .- goal[1:6]) .* [1., 1., 1., 10_000., 10_000., 10_000.]
-        δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
+          g[7+8i] = (norm(v∞_out) - norm(v∞_in)) * 10_000.
-        g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) - current_planet.r
+          δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
-      else
+          g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) - current_planet.r
-        g[8*(length(flybys)-1)+1:8*(length(flybys)-1)+3] .= final[1:3] .- goal[1:3]
+        else
-        g[8*(length(flybys)-1)+4] = final[7] - guess.sc.dry_mass
+          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_iter" => 100_000,
-                       "max_cpu_time" => 300.,
+                       "max_cpu_time" => 60.,
-                       "print_level" => 0)
+                       "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
 

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

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.)

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),
+  dates = [DateTime(2023,5,23),
-           DateTime(2023,10,30),
+           DateTime(2023,10,21),
-           DateTime(2024,9,2),
+           DateTime(2024,8,24),
-           DateTime(2025,2,10),
+           DateTime(2025,2,13),
-           DateTime(2026,11,26),
+           DateTime(2026,11,22),
-           DateTime(2029,10,16)]
+           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")