connor/thesis
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Seems like a lot changed...

Browse Source
This commit is contained in:
Connor
2022-02-19 12:30:54 -07:00
parent 5d5c79c909
commit 329a2990c5
16 changed files with 532 additions and 273 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
julia/src/Thesis.jl
View File

@@ -7,7 +7,7 @@ module Thesis
using PlotlyJS
using Distributed
using SPICE
using Dates: DateTime, Millisecond, Dates, Second, format, datetime2unix, unix2datetime
using Dates: DateTime, Millisecond, Dates, Second, format, datetime2julian, julian2datetime
try
furnsh("../../spice_files/naif0012.tls")

37
julia/src/nlp_solver.jl
View File

@@ -2,30 +2,27 @@ using SNOW
export solve_mission
const pos_scale = ones(3)
# const vel_scale = 100. * ones(3)
const vel_scale = ones(3)
# const v∞_scale = norm(vel_scale)
const v∞_scale = 1.
const state_scale = [pos_scale; vel_scale ]
function constraint_bounds(guess::Mission_Guess)
low_constraint = Vector{Float64}()
high_constraint = Vector{Float64}()
for phase in guess.phases
state_constraint = [100., 100., 100., 0.005, 0.005, 0.005] .* state_scale
# Constraints on the final state
state_constraint = [10_000., 10_000., 10_000., 0.05, 0.05, 0.05]
push!(low_constraint, (-1 * state_constraint)... )
push!(high_constraint, state_constraint... )
if phase != guess.phases[end]
push!(low_constraint, -0.005*v∞_scale, 100.)
push!(high_constraint, 0.005*v∞_scale, 1e7)
# Constraints on the v∞ and turning angle
push!(low_constraint, -0.05, 100.)
push!(high_constraint, 0.05, 1e7)
else
# Constraints on mass
push!(low_constraint, 0.)
push!(high_constraint, guess.start_mass - guess.sc.dry_mass)
end
end
# Constraints on C3 and V∞
push!(low_constraint, 0., 0.)
push!(high_constraint, 1e3, 1e3)
@@ -57,8 +54,8 @@ function solve_mission( guess::Mission_Guess,
# Establish initial conditions
v∞_out = x[2:4]
current_planet = Earth
launch_date = unix2datetime(x[1])
time = utc2et(format(launch_date,"yyyy-mm-ddTHH:MM:SS"))
launch_date = julian2datetime(x[1])
time = Dates.datetime2julian(launch_date)
start = state(current_planet, time, v∞_out, guess.start_mass)
final = zeros(7)
@@ -80,16 +77,16 @@ function solve_mission( guess::Mission_Guess,
thrusts = reshape(phase_params[8:end], (n,3))
# Propagate
final = prop(thrusts, start, guess.sc, tof)[2]
final = prop(thrusts, start, guess.sc, tof)[end,:]
current_planet = flyby
time += tof
time += tof/86400
goal = state(current_planet, time, v∞_in)[1:6]
start = state(current_planet, time, v∞_out, final[7])
# Do Checks
g[1+8i:6+8i] .= (final[1:6] .- goal[1:6]) .* state_scale
g[1+8i:6+8i] .= (final[1:6] .- goal[1:6])
if flyby != flybys[end]
g[7+8i] = (norm(v∞_out) - norm(v∞_in)) * v∞_scale
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 ) - current_planet.r
else
@@ -100,7 +97,7 @@ function solve_mission( guess::Mission_Guess,
i += 1
end
return 1.0
return -final[end]
catch e
@@ -114,7 +111,7 @@ function solve_mission( guess::Mission_Guess,
end
end
max_time = datetime2unix(latest_arrival) - datetime2unix(launch_window[1])
max_time = (datetime2julian(latest_arrival) - datetime2julian(launch_window[1]))*86400
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) + 9
@@ -123,7 +120,7 @@ function solve_mission( guess::Mission_Guess,
"acceptable_constr_viol_tol" => 100tol,
"bound_relax_factor" => 0.,
"max_iter" => 100_000,
"max_cpu_time" => 5. * length(guess.phases),
"max_cpu_time" => 10. * length(guess.phases),
"print_level" => print_level)
options = Options(solver=IPOPT(ipopt_options), derivatives=ForwardFD())

100
julia/src/types/bodies.jl
View File

@@ -1,22 +1,106 @@
export state, period
struct Body
μ::Float64
r::Float64 # radius
μ::Real
r::Real # radius
color::String
id::Int # SPICE id
a::Float64 # semi-major axis
a::Real # semi-major axis
line_color::AbstractString
name::AbstractString
end
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(format(t,"yyyy-mm-ddTHH:MM:SS"))
[ spkssb(p.id, time, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
function state(p::Body, t::DateTime, add_v∞::AbstractVector=[0., 0., 0.], add_mass::Float64=1e10)
time = Dates.datetime2julian(t)
[ basic_ephemeris(p,time); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
end
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 ]
function state(p::Body, t::Real, add_v∞::AbstractVector=[0., 0., 0.], add_mass::Real=1e10)
[ basic_ephemeris(p,t); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
end
"""
Provides a really basic, but fast, ephemeris for most planets.
"""
function basic_ephemeris(planet::Body,date::Real)
century = (date-2451545)/36525
if planet.name == "Mercury"
table = [ 252.250906 149472.6746358 0.00000535 0.000000002 ;
0.387098310 0 0 0 ;
0.20563175 0.000020406 0.0000000284 0.00000000017 ;
7.004986 -0.0059516 0.00000081 0.000000041 ;
48.330893 -0.1254229 -0.00008833 -0.000000196 ;
77.456119 0.1588643 -0.00001343 0.000000039]
elseif planet.name == "Venus"
table = [ 181.979801 58517.8156760 0.00000165 -0.000000002 ;
0.72332982 0 0 0 ;
0.00677188 -0.000047766 0.0000000975 0.00000000044 ;
3.394662 -0.0008568 -0.00003244 0.000000010 ;
76.679920 -0.2780080 -0.00014256 -0.000000198 ;
131.563707 0.0048646 -0.00138232 -0.000005332]
elseif planet.name == "Earth"
table = [ 100.466449 35999.3728519 -0.00000568 0 ;
1.000001018 0 0 0 ;
0.01670862 -0.000042037 -0.0000001236 0.00000000004 ;
0 0.0130546 -0.00000931 -0.000000034 ;
174.873174 -0.2410908 0.00004067 -0.000001327 ;
102.937348 0.3225557 0.00015026 0.000000478]
elseif planet.name == "Mars"
table = [ 355.433275 19140.2993313 0.00000261 -0.000000003 ;
1.523679342 0 0 0 ;
0.09340062 0.000090483 -0.0000000806 -0.00000000035 ;
1.849726 -0.0081479 -0.00002255 -0.000000027 ;
49.558093 -0.2949846 -0.00063993 -0.000002143 ;
336.060234 0.4438898 -0.00017321 0.000000300]
elseif planet.name == "Jupiter"
table = [ 34.351484 3034.9056746 -0.00008501 0.000000004 ;
5.202603191 0.0000001913 0 0 ;
0.04849485 0.000163244 -0.0000004719 -0.00000000197 ;
1.303270 -0.0019872 0.00003318 0.000000092 ;
100.464441 0.1766828 0.00090387 -0.000007032 ;
14.331309 0.2155525 0.00072252 -0.000004590]
elseif planet.name == "Saturn"
table = [ 50.077471 1222.1137943 0.00021004 -0.000000019 ;
9.554909596 -0.0000021389 0 0 ;
0.05550862 -0.000346818 -0.0000006456 0.00000000338 ;
2.488878 0.0025515 -0.00004903 0.000000018 ;
113.665524 -0.2566649 -0.00018345 0.000000357 ;
93.056787 0.5665496 0.00052809 0.000004882]
elseif planet.name == "Uranus"
table = [ 314.055005 429.8640561 0.00030434 0.000000026 ;
19.218446062 -0.0000000372 0.00000000098 0 ;
0.04629590 -0.000027337 0.0000000790 0.00000000025 ;
0.773196 0.0007744 0.00003749 -0.000000092 ;
74.005947 0.5211258 0.00133982 0.000018516 ;
173.005159 1.4863784 0.0021450 0.000000433]
elseif planet.name == "Neptune"
table = [ 304.348665 219.8833092 0.00030926 0.000000018 ;
30.110386869 -0.0000001663 0.00000000069 0 ;
0.00898809 0.000006408 -0.0000000008 -0.00000000005 ;
1.769952 -0.0093082 -0.00000708 0.000000028 ;
131.784057 1.1022057 0.00026006 -0.000000636 ;
48.123691 1.4262677 0.00037918 -0.000000003]
elseif planet.name == "Pluto"
table = [ 238.92903833 145.20780515 0 0 ;
39.48211675 -0.00031596 0 0 ;
0.24882730 0.00005170 0 0 ;
17.14001206 0.00004818 0 0 ;
110.30393684 -0.01183482 0 0 ;
224.06891629 -0.04062942 0 0]
else
return "No such planet"
end
poly = [1,century,century^2,century^3]
L,a,e,i,Ω,P = table * poly
L,i,Ω,P = [L,i,Ω,P] * (π/180)
ω = P-Ω
M = L-P
θ = M + (2e - 0.25e^3 + (5/96)*e^5)*sin(M) +
(1.25e^2 - (11/24)*e^4)*sin(2M) + ((13/12)*e^3 - (43/64)*e^5)*sin(3M) +
(103/96)*e^4*sin(4M) + (1097/960)*e^5*sin(5M)
oe = [a*AU,e,i,Ω,ω,θ]
return oe_to_xyz(oe, Sun.μ)
end

16
julia/src/types/mission.jl
View File

@@ -49,7 +49,7 @@ This is the opposite of the function below. It takes a vector and any other nece
"""
function Mission_Guess(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body})
# Variable mission params
launch_date = unix2datetime(x[1])
launch_date = julian2datetime(x[1])
launch_v∞ = x[2:4]
# Try to intelligently determine n
@@ -78,7 +78,7 @@ information from the guess though.
"""
function Base.Vector(g::Mission_Guess)
result = Vector{Float64}()
push!(result, datetime2unix(g.launch_date))
push!(result, datetime2julian(g.launch_date))
push!(result, g.launch_v∞...)
for phase in g.phases
push!(result,phase.v∞_in...)
@@ -91,7 +91,7 @@ end
function lowest_mission_vector(launch_window::Tuple{DateTime,DateTime}, num_phases::Int, n::Int)
result = Vector{Float64}()
push!(result, datetime2unix(launch_window[1]))
push!(result, datetime2julian(launch_window[1]))
push!(result, -10*ones(3)...)
for i in 1:num_phases
push!(result, -10*ones(3)...)
@@ -104,7 +104,7 @@ end
function highest_mission_vector(launch_window::Tuple{DateTime,DateTime}, mission_length::Float64, num_phases::Int, n::Int)
result = Vector{Float64}()
push!(result, datetime2unix(launch_window[2]))
push!(result, datetime2julian(launch_window[2]))
push!(result, 10*ones(3)...)
for i in 1:num_phases
push!(result, 10*ones(3)...)
@@ -133,14 +133,14 @@ function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, p
force=false)
# First do some checks to make sure that it's valid
if !force
time = utc2et(format(date,"yyyy-mm-ddTHH:MM:SS"))
time = datetime2julian(date)
current_planet = Earth
start = state(current_planet, time, v∞, mass)
for phase in phases
#Propagate
final = prop(phase.thrust_profile, start, sc, phase.tof)[2]
final = prop(phase.thrust_profile, start, sc, phase.tof)[end,:]
current_planet = phase.planet
time += phase.tof
time += phase.tof/86400
goal = state(current_planet, time, phase.v∞_in)
start = state(current_planet, time, phase.v∞_out, final[7])
@@ -178,7 +178,7 @@ end
function Mission(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body})
# Variable mission params
launch_date = unix2datetime(x[1])
launch_date = julian2datetime(x[1])
launch_v∞ = x[2:4]
# Try to intelligently determine n

2
julia/src/utilities/laguerre-conway.jl
View File

@@ -1,4 +1,4 @@
function laguerre_conway(state::Vector{<:Real}, time::Float64, primary::Body=Sun)
function laguerre_conway(state::AbstractArray, time::Real, primary::Body=Sun)
μ = primary.μ
n = 5 # Choose LaGuerre-Conway "n"

26
julia/src/utilities/plotting.jl
View File

@@ -59,27 +59,27 @@ end
"""
A fundamental function to generate a plot for a single orbit path
"""
function gen_plot(path::Vector{Vector{Float64}};
function gen_plot(path::AbstractArray;
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]),
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]),
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]]),
push!(traces, scatter3d(;x=([path[1,1]]),y=([path[1,2]]),z=([path[1,3]]),
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]]),
push!(traces, scatter3d(;x=([path[end,1]]),y=([path[end,2]]),z=([path[end,3]]),
mode="markers", showlegend=false,
marker=attr(color=color, size=3, symbol="diamond")))
end
@@ -222,7 +222,7 @@ end
"""
Generic plotting function for a bunch of paths. Useful for debugging.
"""
function plot(paths::Vector{Vector{Vector{Float64}}};
function plot(paths::Vector{Matrix{Real}};
labels::Union{Vector{String},Nothing}=nothing,
colors::Union{Vector{String},Nothing}=nothing,
markers::Bool=true,
@@ -243,7 +243,8 @@ function plot(paths::Vector{Vector{Vector{Float64}}};
push!(traces, gen_plot(Sun)[1]...)
# Determine layout details
layout = standard_layout(limit, title)
# layout = standard_layout(limit, title)
layout = standard_layout(1e9, title)
# Plot
PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
@@ -253,19 +254,20 @@ 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
time = datetime2julian(m.launch_date)
mass = m.start_mass
current_planet = Earth
earth_trace, limit = gen_plot(current_planet, time, year)
earth_trace, limit = gen_plot(current_planet, julian2datetime(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)
path = prop(phase.thrust_profile, start, m.sc, phase.tof)
final = path[end,:]
mass = final[7]
time += Second(floor(phase.tof))
time += phase.tof/86400
current_planet = phase.planet
start = state(current_planet, time, phase.v∞_out, mass)
path_trace, new_limit = gen_plot(path, label="Phase "*string(i))
@@ -273,7 +275,7 @@ function plot(m::Union{Mission, Mission_Guess}; title::String="Mision Plot")
limit = max(limit, new_limit)
# Then the planet
planet_trace, new_limit = gen_plot(current_planet, time, -phase.tof)
planet_trace, new_limit = gen_plot(current_planet, julian2datetime(time), -phase.tof)
push!(traces, planet_trace...)
limit = max(limit, new_limit)

83
julia/src/utilities/propagator.jl
View File

@@ -3,83 +3,52 @@ export prop
"""
Maximum ΔV that a spacecraft can impulse for a given single time step
"""
function max_ΔV(duty_cycle::Float64,
function max_ΔV(duty_cycle::Real,
num_thrusters::Int,
max_thrust::Float64,
tf::Float64,
t0::Float64,
mass::T) where T <: Real
return duty_cycle*num_thrusters*max_thrust*(tf-t0)/mass
end
"""
A convenience function for using spacecraft. Note that this function outputs a sc instead of a mass
"""
function prop_one(ΔV_unit::Vector{<:Real},
state::Vector{<:Real},
craft::Sc,
time::Float64,
primary::Body=Sun)
ΔV = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time, 0., state[7]) * ΔV_unit
halfway = laguerre_conway(state, time/2, primary) + [zeros(3); ΔV]
final = laguerre_conway(halfway, time/2, primary)
return [final; state[7] - mfr(craft)*norm(ΔV_unit)*time]
max_thrust::Real,
tf::Real,
t0::Real)
return duty_cycle*num_thrusters*max_thrust*(tf-t0)
end
"""
The propagator function
"""
function prop(ΔVs::Matrix{T},
state::Vector{Float64},
function prop(ΔVs::AbstractArray,
state::AbstractArray,
craft::Sc,
time::Float64,
primary::Body=Sun;
interpolate::Bool=false) where T <: Real
time::Real;
primary::Body=Sun)
size(ΔVs)[2] == 3 || throw(ΔVsize_Error())
size(ΔVs)[2] == 3 || throw(ΔVsize_Error(size(ΔVs)))
n = size(ΔVs)[1]
states = [ Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}() ]
for i in 1:7 push!(states[i], state[i]) end
output = Array{Real,2}(undef,n,7)
current = copy(state)
full_thrust = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time/n, 0.)
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
halfway = laguerre_conway(current[1:6], time/(2n), primary)
if time > 0
halfway += [zeros(3); full_thrust * ΔVs[i,:]] / current[7]
current[7] -= mfr(craft) * norm(ΔVs[i,:]) * time/n
else
current[7] -= mfr(craft) * norm(ΔVs[i,:]) * time/n
halfway += [zeros(3); full_thrust * ΔVs[i,:]] / current[7]
end
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
current[1:6] = laguerre_conway(halfway, time/(2n), primary)
output[i,:] = current
end
state[7] > craft.dry_mass || throw(Mass_Error(state[7]))
current[7] > craft.dry_mass || throw(Mass_Error(current[7]))
return states, state
return output
end
"""
Convenience function for propagating a state with no thrust
"""
prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1.], no_thrust, t, p)[1]
prop(x::AbstractArray, t::Real, p::Body=Sun) = prop(zeros(1000,3), x, no_thrust, t, p)
"""
This is solely for the purposes of getting the final state of a mission or guess
@@ -91,7 +60,7 @@ function prop(m::Union{Mission, Mission_Guess})
final = zeros(7)
for phase in m.phases
final = prop(phase.thrust_profile, start, m.sc, phase.tof)[2]
final = prop(phase.thrust_profile, start, m.sc, phase.tof)[end,:]
mass = final[7]
current_planet = phase.planet
time += Second(floor(phase.tof))