From 58ad3c529372fbf312d3b0382779af33c87b5f30 Mon Sep 17 00:00:00 2001 From: Connor Date: Wed, 1 Sep 2021 11:08:43 -0600 Subject: [PATCH] Still not working. Up to the loop though. --- julia/src/conversions.jl | 3 + julia/src/find_closest.jl | 129 +++++++++++++------------------------ julia/test/find_closest.jl | 21 +++--- 3 files changed, 59 insertions(+), 94 deletions(-) diff --git a/julia/src/conversions.jl b/julia/src/conversions.jl index cae9f5f..8e0d6f3 100644 --- a/julia/src/conversions.jl +++ b/julia/src/conversions.jl @@ -69,6 +69,9 @@ function xyz_to_oe(cart_vec::Vector,μ::Real) end +""" +Converts a series of thrust vectors from R,Θ,H frame to cartesian +""" function conv_T(Tr::Vector{Float64}, TΘ::Vector{Float64}, Th::Vector{Float64}, diff --git a/julia/src/find_closest.jl b/julia/src/find_closest.jl index 4af1899..dc4b448 100644 --- a/julia/src/find_closest.jl +++ b/julia/src/find_closest.jl @@ -2,39 +2,6 @@ using JuMP, Ipopt export nlp_solve -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}, final::Vector{Float64}, craft::Sc, @@ -43,7 +10,8 @@ function nlp_solve(start::Vector{Float64}, tf::Float64, Txi::Vector{Float64}, Tyi::Vector{Float64}, - Tzi::Vector{Float64}, + Tzi::Vector{Float64}; + solver_options=(), tol=1e-6) n::Int = length(Txi) @@ -53,69 +21,58 @@ function nlp_solve(start::Vector{Float64}, throw("Bad number of Tzs") end - model = Model(Ipopt.Optimizer) - set_optimizer_attribute(model, "max_cpu_time", 30.) + # First propagate the initial guess + path, masses, final_state = prop(hcat(Txi, Tyi, Tzi), + start, + craft, + μ, + tf-t0) + @assert final ≈ final_state + model = Model(optimizer_with_attributes(Ipopt.Optimizer, solver_options...)) @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] + x[i=1:n], (start=path[1][i]) + y[i=1:n], (start=path[2][i]) + z[i=1:n], (start=path[3][i]) + dx[i=1:n], (start=path[4][i]) + dy[i=1:n], (start=path[5][i]) + dz[i=1:n], (start=path[6][i]) + mass[i=1:n], (start=masses[i]) 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) + # Fix initial conditions + fix(x[1], start[1], force = true) + fix(y[1], start[2], force = true) + fix(z[1], start[3], force = true) + fix(dx[1], start[4], force = true) + fix(dy[1], start[5], force = true) + fix(dz[1], start[6], force = true) + fix(mass[1], craft.mass, force = true) - register(model, :prop_one_simple, 11, prop_one_simple; autodiff = true) - @NLexpression(model, ) + # Fix final conditions + fix(x[n], final[1], force = true) + fix(y[n], final[2], force = true) + fix(z[n], final[3], force = true) + fix(dx[n], final[4], force = true) + fix(dy[n], final[5], force = true) + fix(dz[n], final[6], force = true) -# @NLexpression() + @NLexpression(model, r_mag[i = 1:n], sqrt(x[i]^2 + y[i]^2 + z[i]^2)) + @NLexpression(model, v_mag[i = 1:n], sqrt(dx[i]^2 + dy[i]^2 + dz[i]^2)) + @NLexpression(model, σ0[i = 1:n], (x[i]*dx[i] + y[i]*dy[i] + z[i]*dz[i])/sqrt(μ)) + @NLexpression(model, a[i = 1:n], 1/(2/r_mag[i] - v_mag[i]^2/μ)) + + # Elliptical Specific + @NLexpression(model, ΔM_ell[i = 1:n], sqrt(μ) / sqrt(a[i]^3) * (tf-t0)/(2n)) + + # Parabolic/Hyperbolic Specific + @NLexpression(model, ΔN_hyp[i = 1:n], sqrt(μ) / sqrt(-a[i]^3) * (tf-t0)/(2n)) 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, x[i] == a[i-1]) @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) diff --git a/julia/test/find_closest.jl b/julia/test/find_closest.jl index 9ccc622..f5aa331 100644 --- a/julia/test/find_closest.jl +++ b/julia/test/find_closest.jl @@ -1,7 +1,6 @@ @testset "Find Closest" begin using JuMP - using Thesis: treat_inputs # Initial Setup sc = Sc("test") @@ -14,12 +13,17 @@ # A simple orbit raising start = oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"]) - ΔVs = repeat([0.6, 0., 0.]', outer=(n,1)) - final = prop(ΔVs, start, sc, μs["Earth"], prop_time)[3] + Tx, Ty, Tz = conv_T(repeat([0.6], n), repeat([0.], n), repeat([0.], n), + start, + sc.mass, + sc, + prop_time, + μs["Earth"]) + final = prop(hcat(Tx, Ty, Tz), start, sc, μs["Earth"], prop_time)[3] new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"]) # This should be close enough to 0.6 - Tr, TΘ, Th = conv_T(repeat([0.6], n), repeat([0.], n), repeat([0.], n), + Tx, Ty, Tz = conv_T(repeat([0.6], n), repeat([0.], n), repeat([0.], n), start, sc.mass, sc, @@ -31,9 +35,10 @@ μs["Earth"], 0.0, prop_time, - Tr, - TΘ, - Th) + Tx, + Ty, + Tz) + # solver_options=("max_cpu_time" => 30.)) # Test and plot @test JuMP.termination_status(result) == MOI.OPTIMAL @@ -44,7 +49,7 @@ savefig(plot_orbits([path1, path2, path3, path4], labels=["initial", "transit", "after transit", "final"], colors=["#FFFFFF","#FF4444","#44FF44","#4444FF"]), - "../plots/find_closest_test.html") + "../plots/find_closest_test.html") # if termination_status(result) == :OPTIMAL # @test norm(calc_final - final) < 1e-4 # end