Still not working. Up to the loop though.

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},

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)
+  # First propagate the initial guess 
-  set_optimizer_attribute(model, "max_cpu_time", 30.)
+  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]
+    x[i=1:n], (start=path[1][i])
-    y[1:n]
+    y[i=1:n], (start=path[2][i])
-    z[1:n]
+    z[i=1:n], (start=path[3][i])
-    dx[1:n]
+    dx[i=1:n], (start=path[4][i])
-    dy[1:n]
+    dy[i=1:n], (start=path[5][i])
-    dz[1:n]
+    dz[i=1:n], (start=path[6][i])
-    mass[1:n]
+    mass[i=1:n], (start=masses[i])
   end)
 
-  @constraints(model,begin
+  # Fix initial conditions
-    x[1] == start[1]
+  fix(x[1], start[1], force = true)
-    y[1] == start[2]
+  fix(y[1], start[2], force = true)
-    z[1] == start[3]
+  fix(z[1], start[3], force = true)
-    dx[1] == start[4]
+  fix(dx[1], start[4], force = true)
-    dy[1] == start[5]
+  fix(dy[1], start[5], force = true)
-    dz[1] == start[6]
+  fix(dz[1], start[6], force = true)
-    mass[1] == craft.mass
+  fix(mass[1], craft.mass, force = true)
-    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)
+  # Fix final conditions
-  @NLexpression(model, )
+  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] == a[i-1])
-   #  @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)

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))
+  Tx, Ty, Tz = conv_T(repeat([0.6], n), repeat([0.], n), repeat([0.], n),
-  final = prop(ΔVs, start, sc, μs["Earth"], prop_time)[3]
+                      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,
+                               Tx,
-                               TΘ,
+                               Ty,
-                               Th)
+                               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