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

Latest attempt at jump. Does not work yet.

Browse Source
This commit is contained in:
Connor
2021-08-30 23:28:36 -06:00
parent 850f05ce38
commit 5daf007e5b
11 changed files with 452 additions and 117 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
julia/src/Thesis.jl
View File

@@ -3,8 +3,8 @@ module Thesis
using LinearAlgebra, ForwardDiff, PlotlyJS
include("./constants.jl")
include("./conversions.jl")
include("./spacecraft.jl")
include("./conversions.jl")
include("./plotting.jl")
include("./laguerre-conway.jl")
include("./propagator.jl")

44
julia/src/conversions.jl
View File

@@ -1,4 +1,4 @@
export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe
export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe, conv_T
function oe_to_rθh(oe::Vector,μ::Real) :: Vector
@@ -68,3 +68,45 @@ function xyz_to_oe(cart_vec::Vector,μ::Real)
return [a,e,i,Ω,ω,ν]
end
function conv_T(Tr::Vector{Float64},
TΘ::Vector{Float64},
Th::Vector{Float64},
init_state::Vector{Float64},
m::Float64,
craft::Sc,
time::Float64,
μ::Float64)
Txs = Float64[]
Tys = Float64[]
Tzs = Float64[]
state = init_state
for i in 1:length(Tr)
mag, α, β = Tr[i], TΘ[i], Th[i]
if mag > 1 || mag < 0
throw(ErrorException("ΔV input is too high: $mag"))
elseif α > π || α < -π
throw(ErrorException("α angle is incorrect: $α"))
elseif β > π/2 || β < -π/2
throw(ErrorException("β angle is incorrect: $β"))
end
thrust_rθh = mag * [cos(β)*sin(α), cos(β)*cos(α), sin(β)]
_,_,i,Ω,ω,ν = xyz_to_oe(state, μ)
θ = ω+ν
cΩ,sΩ,ci,si,cθ,sθ = cos(Ω),sin(Ω),cos(i),sin(i),cos(θ),sin(θ)
DCM = [cΩ*cθ-sΩ*ci*sθ -cΩ*sθ-sΩ*ci*cθ sΩ*si;
sΩ*cθ+cΩ*ci*sθ -sΩ*sθ+cΩ*ci*cθ -cΩ*si;
si*sθ si*cθ ci]
ΔV = DCM*thrust_rθh
Tx, Ty, Tz = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time, 0., m) * ΔV
x, y, z, dx, dy, dz = state
state = prop_one_simple(Tx, Ty, Tz, x, y, z, dx, dy, dz, time/length(Tr), μ)
push!(Txs, Tx)
push!(Tys, Ty)
push!(Tzs, Tz)
end
return Txs, Tys, Tzs
end

120
julia/src/find_closest.jl
View File

@@ -1,10 +1,38 @@
using NLsolve
using JuMP, Ipopt
export nlp_solve
function treat_inputs(x::AbstractVector)
n::Int = length(x)/3
reshape(x,(3,n))'
function treat_inputs(x::Vector{T}) where T <: Real
thrust = T[]
α = T[]
β = T[]
n::Int = length(x)
for i in 1:n
if i % 3 == 1
push!(thrust,x[i])
elseif i % 3 == 2
push!(α,x[i])
else
push!(β,x[i])
end
end
hcat(thrust, α, β)
end
function treat_inputs(x::NTuple{90,T}) where T <: Real
thrust = T[]
α = T[]
β = T[]
for i in 1:90
if i % 3 == 1
push!(thrust,x[i])
elseif i % 3 == 2
push!(α,x[i])
else
push!(β,x[i])
end
end
hcat(thrust, α, β)
end
function nlp_solve(start::Vector{Float64},
@@ -13,15 +41,87 @@ function nlp_solve(start::Vector{Float64},
μ::Float64,
t0::Float64,
tf::Float64,
x0::AbstractVector;
Txi::Vector{Float64},
Tyi::Vector{Float64},
Tzi::Vector{Float64},
tol=1e-6)
n::Int = length(x0)/3
function f!(F,x)
F[1:6] .= prop(treat_inputs(x), start, craft, μ, tf-t0)[1][end,:] - final
F[7:3n] .= 0.
n::Int = length(Txi)
if length(Tyi) != n
throw("Bad number of Tys")
elseif length(Tzi) != n
throw("Bad number of Tzs")
end
return nlsolve(f!, x0, ftol=tol, autodiff=:forward, iterations=1_000)
model = Model(Ipopt.Optimizer)
set_optimizer_attribute(model, "max_cpu_time", 30.)
@variables(model, begin
Tx[i=1:n], (start=Txi[i])
Ty[i=1:n], (start=Tyi[i])
Tz[i=1:n], (start=Tzi[i])
x[1:n]
y[1:n]
z[1:n]
dx[1:n]
dy[1:n]
dz[1:n]
mass[1:n]
end)
@constraints(model,begin
x[1] == start[1]
y[1] == start[2]
z[1] == start[3]
dx[1] == start[4]
dy[1] == start[5]
dz[1] == start[6]
mass[1] == craft.mass
x[n] == final[1]
y[n] == final[2]
z[n] == final[3]
dx[n] == final[4]
dy[n] == final[5]
dz[n] == final[6]
end)
register(model, :prop_one_simple, 11, prop_one_simple; autodiff = true)
@NLexpression(model, )
# @NLexpression()
for i in 2:n
@NLconstraint(model, x[i] == sqrt(x[i-1]^2))
# @NLconstraint(model, x[i] == prop_one_simple(Tx[i-1], Ty[i-1], Tz[i-1],
# x[i-1], y[i-1], z[i-1],
# dx[i-1], dy[i-1], dz[i-1],
# (tf-t0)/n, μ)[1])
# @NLconstraint(model, y[i] == prop_one_simple(Tx[i-1], Ty[i-1], Tz[i-1],
# x[i-1], y[i-1], z[i-1],
# dx[i-1], dy[i-1], dz[i-1],
# (tf-t0)/n, μ)[2])
# @NLconstraint(model, z[i] == prop_one_simple(Tx[i-1], Ty[i-1], Tz[i-1],
# x[i-1], y[i-1], z[i-1],
# dx[i-1], dy[i-1], dz[i-1],
# (tf-t0)/n, μ)[3])
# @NLconstraint(model, dx[i] == prop_one_simple(Tx[i-1], Ty[i-1], Tz[i-1],
# x[i-1], y[i-1], z[i-1],
# dx[i-1], dy[i-1], dz[i-1],
# (tf-t0)/n, μ)[4])
# @NLconstraint(model, dy[i] == prop_one_simple(Tx[i-1], Ty[i-1], Tz[i-1],
# x[i-1], y[i-1], z[i-1],
# dx[i-1], dy[i-1], dz[i-1],
# (tf-t0)/n, μ)[5])
# @NLconstraint(model, dz[i] == prop_one_simple(Tx[i-1], Ty[i-1], Tz[i-1],
# x[i-1], y[i-1], z[i-1],
# dx[i-1], dy[i-1], dz[i-1],
# (tf-t0)/n, μ)[6])
@NLconstraint(model, mass[i] == mass[i-1] - craft.mass_flow_rate*√(Tx[i-1]^2 +
Ty[i-1]^2 +
Tz[i-1]^2)*(tf-t0)/n)
end
optimize!(model)
return model
end

6
julia/src/laguerre-conway.jl
View File

@@ -1,12 +1,12 @@
function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T <: Real
function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T
n = 5 # Choose LaGuerre-Conway "n"
i = 0
r0 = state[1:3] # Are we in elliptical or hyperbolic orbit?
r0_mag = norm(r0)
r0_mag = √(state[1]^2 + state[2]^2 + state[3]^2)
v0 = state[4:6]
v0_mag = norm(v0)
v0_mag = √(state[4]^2 + state[5]^2 + state[6]^2)
σ0 = (r0 ⋅ v0)/√(μ)
a = 1 / ( 2/r0_mag - v0_mag^2/μ )
coeff = 1 - r0_mag/a

9
julia/src/plotting.jl
View File

@@ -16,7 +16,7 @@ function get_true_max(mat::Vector{Array{Float64,2}})
return maximum(abs.(flatten(mat)))
end
function plot_orbits(paths::Vector{Array{Float64,2}};
function plot_orbits(paths::Vector{Vector{Vector{Float64}}};
primary::String="Earth",
plot_theme::Symbol=:juno,
labels::Vector{String}=Vector{String}(),
@@ -34,13 +34,14 @@ function plot_orbits(paths::Vector{Array{Float64,2}};
t1 = []
for i = 1:length(paths)
path = [ x for x in paths[i] ]
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]),
push!(t1, scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]),
mode="lines", name=label, line_color=color, line_width=3))
end
limit = max(maximum(abs.(flatten(paths))),maximum(abs.(flatten(ps)))) * 1.1
limit = max(maximum(abs.(flatten(flatten(paths)))),
maximum(abs.(flatten(ps)))) * 1.1
t2 = surface(;x=(x_p),
y=(y_p),

143
julia/src/propagator.jl
View File

@@ -17,14 +17,14 @@ This function propagates the spacecraft forward in time 1 Sim-Flanagan step (of
applying a thrust in the center.
"""
function prop_one(ΔV::Vector{T},
state::Vector{S},
state::Vector{T},
duty_cycle::Float64,
num_thrusters::Int,
max_thrust::Float64,
mass::S,
mass::T,
mass_flow_rate::Float64,
μ::Float64,
time::Float64) where {T <: Real, S <: Real}
time::Float64) where T
mag, α, β = ΔV
@@ -64,6 +64,100 @@ function prop_one(ΔV_unit::Vector{T},
return state, Sc(mass, craft.mass_flow_rate, craft.max_thrust, craft.num_thrusters, craft.duty_cycle)
end
function prop_one_simple(Tx, Ty, Tz, x, y, z, dx, dy, dz, t, μ)
# perform laguerre_conway once
r0_mag = √(x^2 + y^2 + z^2)
v0_mag = √(dx^2 + dy^2 + dz^2)
σ0 = ([x, y, z] ⋅ [dx, dy, dz])/√(μ)
a = 1 / ( 2/r0_mag - v0_mag^2/μ )
coeff = 1 - r0_mag/a
if a > 0 # Elliptical
ΔM = ΔE_new = √(μ) / sqrt(a^3) * t/2
ΔE = 1000
while abs(ΔE - ΔE_new) > 1e-10
ΔE = ΔE_new
F = ΔE - ΔM + σ0 / √(a) * (1-cos(ΔE)) - coeff * sin(ΔE)
dF = 1 + σ0 / √(a) * sin(ΔE) - coeff * cos(ΔE)
d2F = σ0 / √(a) * cos(ΔE) + coeff * sin(ΔE)
sign = dF >= 0 ? 1 : -1
ΔE_new = ΔE - 5*F / ( dF + sign * √(abs((5-1)^2*dF^2 - 5*(5-1)*F*d2F )))
end
F = 1 - a/r0_mag * (1-cos(ΔE))
G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE)
r = a + (r0_mag - a) * cos(ΔE) + σ0 * √(a) * sin(ΔE)
Ft = -√(a)*√(μ) / (r*r0_mag) * sin(ΔE)
Gt = 1 - a/r * (1-cos(ΔE))
else # Hyperbolic or Parabolic
ΔN = √(μ) / sqrt(-a^3) * t/2
ΔH = 0
ΔH_new = t/2 < 0 ? -1 : 1
while abs(ΔH - ΔH_new) > 1e-10
ΔH = ΔH_new
F = -ΔN - ΔH + σ0 / √(-a) * (cos(ΔH)-1) + coeff * sin(ΔH)
dF = -1 + σ0 / √(-a) * sin(ΔH) + coeff * cos(ΔH)
d2F = σ0 / √(-a) * cos(ΔH) + coeff * sin(ΔH)
sign = dF >= 0 ? 1 : -1
ΔH_new = ΔH - 5*F / ( dF + sign * √(abs((5-1)^2*dF^2 - 5*(5-1)*F*d2F )))
end
F = 1 - a/r0_mag * (1-cos(ΔH))
G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH)
r = a + (r0_mag - a) * cos(ΔH) + σ0 * √(-a) * sin(ΔH)
Ft = -√(-a)*√(μ) / (r*r0_mag) * sin(ΔH)
Gt = 1 - a/r * (1-cos(ΔH))
end
# add the thrust vector
x,y,z,dx,dy,dz = [F*[x,y,z] + G*[dx,dy,dz]; Ft*[x,y,z] + Gt*[dx,dy,dz] + [Tx, Ty, Tz]]
#perform again
r0_mag = √(x^2 + y^2 + z^2)
v0_mag = √(dx^2 + dy^2 + dz^2)
σ0 = ([x, y, z] ⋅ [dx, dy, dz])/√(μ)
a = 1 / ( 2/r0_mag - v0_mag^2/μ )
coeff = 1 - r0_mag/a
if a > 0 # Elliptical
ΔM = ΔE_new = √(μ) / sqrt(a^3) * t/2
ΔE = 1000
while abs(ΔE - ΔE_new) > 1e-10
ΔE = ΔE_new
F = ΔE - ΔM + σ0 / √(a) * (1-cos(ΔE)) - coeff * sin(ΔE)
dF = 1 + σ0 / √(a) * sin(ΔE) - coeff * cos(ΔE)
d2F = σ0 / √(a) * cos(ΔE) + coeff * sin(ΔE)
sign = dF >= 0 ? 1 : -1
ΔE_new = ΔE - 5*F / ( dF + sign * √(abs((5-1)^2*dF^2 - 5*(5-1)*F*d2F )))
end
F = 1 - a/r0_mag * (1-cos(ΔE))
G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE)
r = a + (r0_mag - a) * cos(ΔE) + σ0 * √(a) * sin(ΔE)
Ft = -√(a)*√(μ) / (r*r0_mag) * sin(ΔE)
Gt = 1 - a/r * (1-cos(ΔE))
else # Hyperbolic or Parabolic
ΔN = √(μ) / sqrt(-a^3) * t/2
ΔH = 0
ΔH_new = t/2 < 0 ? -1 : 1
while abs(ΔH - ΔH_new) > 1e-10
ΔH = ΔH_new
F = -ΔN - ΔH + σ0 / √(-a) * (cos(ΔH)-1) + coeff * sin(ΔH)
dF = -1 + σ0 / √(-a) * sin(ΔH) + coeff * cos(ΔH)
d2F = σ0 / √(-a) * cos(ΔH) + coeff * sin(ΔH)
sign = dF >= 0 ? 1 : -1
ΔH_new = ΔH - 5*F / ( dF + sign * √(abs((5-1)^2*dF^2 - 5*(5-1)*F*d2F )))
end
F = 1 - a/r0_mag * (1-cos(ΔH))
G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH)
r = a + (r0_mag - a) * cos(ΔH) + σ0 * √(-a) * sin(ΔH)
Ft = -√(-a)*√(μ) / (r*r0_mag) * sin(ΔH)
Gt = 1 - a/r * (1-cos(ΔH))
end
return [F*[x,y,z] + G*[dx,dy,dz]; Ft*[x,y,z] + Gt*[dx,dy,dz]]
return repeat([sqrt(Tx^2 + Ty^2 + Tz^2)],6)
end
"""
This propagates over a given time period, with a certain number of intermediate steps
"""
@@ -92,24 +186,51 @@ end
"""
The same function, using Scs
"""
function prop(ΔVs::AbstractArray{T},
state::Vector{Float64},
function prop(ΔVs::AbstractMatrix,
state::AbstractVector,
craft::Sc,
μ::Float64,
time::Float64) where T <: Real
time::Float64)
if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
n = size(ΔVs)[1]
states = state'
masses = craft.mass
x_states = [state[1]]
y_states = [state[2]]
z_states = [state[3]]
dx_states = [state[4]]
dy_states = [state[5]]
dz_states = [state[6]]
masses = [craft.mass]
for i in 1:n
state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n)
states = [states; state']
masses = [masses, craft.mass]
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, craft.mass)
end
return states, masses
return [x_states, y_states, z_states, dx_states, dy_states, dz_states], masses, state
end
function prop_simple(ΔVs::AbstractMatrix,
state::AbstractVector,
craft::Sc,
μ::Float64,
time::Float64)
if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
n = size(ΔVs)[1]
for i in 1:n
state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n)
end
return state
end