From 5daf007e5be63bfceaffae332210ffcce94de896 Mon Sep 17 00:00:00 2001 From: Connor Date: Mon, 30 Aug 2021 23:28:36 -0600 Subject: [PATCH] Latest attempt at jump. Does not work yet. --- julia/Manifest.toml | 183 +++++++++++++++++--------- julia/Project.toml | 3 +- julia/src/Thesis.jl | 2 +- julia/src/conversions.jl | 44 ++++++- julia/src/find_closest.jl | 120 +++++++++++++++-- julia/src/laguerre-conway.jl | 6 +- julia/src/plotting.jl | 9 +- julia/src/propagator.jl | 143 ++++++++++++++++++-- julia/test/find_closest.jl | 33 +++-- julia/test/monotonic_basin_hopping.jl | 18 +-- julia/test/runtests.jl | 8 +- 11 files changed, 452 insertions(+), 117 deletions(-) diff --git a/julia/Manifest.toml b/julia/Manifest.toml index 1dbf6e2..0cefb85 100644 --- a/julia/Manifest.toml +++ b/julia/Manifest.toml @@ -1,14 +1,14 @@ # This file is machine-generated - editing it directly is not advised +[[ASL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "370cafc70604b2522f2c7cf9915ebcd17b4cd38b" +uuid = "ae81ac8f-d209-56e5-92de-9978fef736f9" +version = "0.1.2+0" + [[ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -[[ArrayInterface]] -deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"] -git-tree-sha1 = "baf4ef9082070477046bd98306952292bfcb0af9" -uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "3.1.25" - [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -21,24 +21,60 @@ version = "0.1.0" [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +[[BenchmarkTools]] +deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"] +git-tree-sha1 = "42ac5e523869a84eac9669eaceed9e4aa0e1587b" +uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +version = "1.1.4" + [[BinDeps]] deps = ["Libdl", "Pkg", "SHA", "URIParser", "Unicode"] git-tree-sha1 = "1289b57e8cf019aede076edab0587eb9644175bd" uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee" version = "1.0.2" +[[BinaryProvider]] +deps = ["Libdl", "Logging", "SHA"] +git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058" +uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" +version = "0.5.10" + [[Blink]] deps = ["Base64", "BinDeps", "Distributed", "JSExpr", "JSON", "Lazy", "Logging", "MacroTools", "Mustache", "Mux", "Reexport", "Sockets", "WebIO", "WebSockets"] git-tree-sha1 = "08d0b679fd7caa49e2bca9214b131289e19808c0" uuid = "ad839575-38b3-5650-b840-f874b8c74a25" version = "0.12.5" +[[Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.8+0" + +[[Calculus]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad" +uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" +version = "0.5.1" + [[ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" version = "1.3.0" +[[CodecBzip2]] +deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"] +git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7" +uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" +version = "0.7.2" + +[[CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.0" + [[ColorTypes]] deps = ["FixedPointNumbers", "Random"] git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597" @@ -72,6 +108,12 @@ git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" version = "1.7.0" +[[DataStructures]] +deps = ["Compat", "InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.18.10" + [[DataValueInterfaces]] git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464" @@ -97,12 +139,6 @@ git-tree-sha1 = "3ed8fa7178a10d1cd0f1ca524f249ba6937490c0" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" version = "1.3.0" -[[Distances]] -deps = ["LinearAlgebra", "Statistics", "StatsAPI"] -git-tree-sha1 = "abe4ad222b26af3337262b8afb28fab8d215e9f8" -uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" -version = "0.10.3" - [[Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -120,12 +156,6 @@ uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" [[FileWatching]] uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" -[[FiniteDiff]] -deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] -git-tree-sha1 = "8b3c09b56acaf3c0e581c66638b85c8650ee9dca" -uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.8.1" - [[FixedPointNumbers]] deps = ["Statistics"] git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" @@ -156,11 +186,6 @@ git-tree-sha1 = "6187bb2d5fcbb2007c39e7ac53308b0d371124bd" uuid = "9fb69e20-1954-56bb-a84f-559cc56a8ff7" version = "0.2.2" -[[IfElse]] -git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef" -uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" -version = "0.1.0" - [[IniFile]] deps = ["Test"] git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8" @@ -171,6 +196,18 @@ version = "0.5.0" deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[Ipopt]] +deps = ["BinaryProvider", "Ipopt_jll", "Libdl", "LinearAlgebra", "MathOptInterface", "MathProgBase"] +git-tree-sha1 = "380786b4929b8d18d76e909c6b2eca355b7c3bd6" +uuid = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +version = "0.7.0" + +[[Ipopt_jll]] +deps = ["ASL_jll", "Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "MUMPS_seq_jll", "OpenBLAS32_jll", "Pkg"] +git-tree-sha1 = "82124f27743f2802c23fcb05febc517d0b15d86e" +uuid = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7" +version = "3.13.4+2" + [[IrrationalConstants]] git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94" uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" @@ -199,6 +236,18 @@ git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" version = "0.21.2" +[[JSONSchema]] +deps = ["HTTP", "JSON", "ZipFile"] +git-tree-sha1 = "b84ab8139afde82c7c65ba2b792fe12e01dd7307" +uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" +version = "0.3.3" + +[[JuMP]] +deps = ["Calculus", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Printf", "Random", "SparseArrays", "SpecialFunctions", "Statistics"] +git-tree-sha1 = "4f0a771949bbe24bf70c89e8032c107ebe03f6ba" +uuid = "4076af6c-e467-56ae-b986-b466b2749572" +version = "0.21.9" + [[Kaleido_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "2ef87eeaa28713cb010f9fb0be288b6c1a4ecd53" @@ -235,12 +284,6 @@ uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" [[Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" -[[LineSearches]] -deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] -git-tree-sha1 = "f27132e551e959b3667d8c93eae90973225032dd" -uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" -version = "7.1.1" - [[LinearAlgebra]] deps = ["Libdl"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -254,6 +297,18 @@ version = "0.3.0" [[Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" +[[METIS_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "2dc1a9fc87e57e32b1fc186db78811157b30c118" +uuid = "d00139f3-1899-568f-a2f0-47f597d42d70" +version = "5.1.0+5" + +[[MUMPS_seq_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "METIS_jll", "OpenBLAS32_jll", "Pkg"] +git-tree-sha1 = "1a11a84b2af5feb5a62a820574804056cdc59c39" +uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d" +version = "5.2.1+4" + [[MacroTools]] deps = ["Markdown", "Random"] git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7" @@ -264,6 +319,18 @@ version = "0.5.7" deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" +[[MathOptInterface]] +deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"] +git-tree-sha1 = "575644e3c05b258250bb599e57cf73bbf1062901" +uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +version = "0.9.22" + +[[MathProgBase]] +deps = ["LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "9abbe463a1e9fc507f12a69e7f29346c2cdc472c" +uuid = "fdba3010-5040-5b88-9595-932c9decdf73" +version = "0.7.8" + [[MbedTLS]] deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"] git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe" @@ -286,24 +353,18 @@ git-tree-sha1 = "36995ef0d532fe08119d70b2365b7b03d4e00f48" uuid = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70" version = "1.0.10" +[[MutableArithmetics]] +deps = ["LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01" +uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +version = "0.2.20" + [[Mux]] deps = ["AssetRegistry", "Base64", "HTTP", "Hiccup", "Pkg", "Sockets", "WebSockets"] git-tree-sha1 = "82dfb2cead9895e10ee1b0ca37a01088456c4364" uuid = "a975b10e-0019-58db-a62f-e48ff68538c9" version = "0.7.6" -[[NLSolversBase]] -deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] -git-tree-sha1 = "144bab5b1443545bc4e791536c9f1eacb4eed06a" -uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" -version = "7.8.1" - -[[NLsolve]] -deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"] -git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1" -uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -version = "4.5.1" - [[NaNMath]] git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" @@ -317,6 +378,12 @@ git-tree-sha1 = "3469ef96607a6b9a1e89e54e6f23401073ed3126" uuid = "510215fc-4207-5dde-b226-833fc4488ee2" version = "0.3.3" +[[OpenBLAS32_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "ba4a8f683303c9082e84afba96f25af3c7fb2436" +uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2" +version = "0.3.12+1" + [[OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" @@ -328,12 +395,6 @@ git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" version = "1.4.1" -[[Parameters]] -deps = ["OrderedCollections", "UnPack"] -git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345" -uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" -version = "0.12.2" - [[Parsers]] deps = ["Dates"] git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779" @@ -414,12 +475,6 @@ git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" version = "1.6.1" -[[Static]] -deps = ["IfElse"] -git-tree-sha1 = "62701892d172a2fa41a1f829f66d2b0db94a9a63" -uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.3.0" - [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb" @@ -430,11 +485,6 @@ version = "1.2.12" deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -[[StatsAPI]] -git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510" -uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.0.0" - [[TOML]] deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" @@ -459,6 +509,12 @@ uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[[TranscodingStreams]] +deps = ["Random", "Test"] +git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.9.6" + [[URIParser]] deps = ["Unicode"] git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d" @@ -474,11 +530,6 @@ version = "1.3.0" deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" -[[UnPack]] -git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b" -uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" -version = "1.0.2" - [[Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" @@ -500,6 +551,12 @@ git-tree-sha1 = "eae2fbbc34a79ffd57fb4c972b08ce50b8f6a00d" uuid = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62" version = "0.6.3" +[[ZipFile]] +deps = ["Libdl", "Printf", "Zlib_jll"] +git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7" +uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" +version = "0.9.3" + [[Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" diff --git a/julia/Project.toml b/julia/Project.toml index 01d2154..8c25aa9 100644 --- a/julia/Project.toml +++ b/julia/Project.toml @@ -5,8 +5,9 @@ version = "0.1.0" [deps] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl index 468c411..69c5dcd 100644 --- a/julia/src/Thesis.jl +++ b/julia/src/Thesis.jl @@ -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") diff --git a/julia/src/conversions.jl b/julia/src/conversions.jl index 6bc40d2..cae9f5f 100644 --- a/julia/src/conversions.jl +++ b/julia/src/conversions.jl @@ -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 diff --git a/julia/src/find_closest.jl b/julia/src/find_closest.jl index 9d496c7..4af1899 100644 --- a/julia/src/find_closest.jl +++ b/julia/src/find_closest.jl @@ -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 \ No newline at end of file diff --git a/julia/src/laguerre-conway.jl b/julia/src/laguerre-conway.jl index cb2cca3..7ca4503 100644 --- a/julia/src/laguerre-conway.jl +++ b/julia/src/laguerre-conway.jl @@ -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 diff --git a/julia/src/plotting.jl b/julia/src/plotting.jl index e9fccf3..45d5fac 100644 --- a/julia/src/plotting.jl +++ b/julia/src/plotting.jl @@ -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), diff --git a/julia/src/propagator.jl b/julia/src/propagator.jl index e6ee9f6..5294285 100644 --- a/julia/src/propagator.jl +++ b/julia/src/propagator.jl @@ -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 diff --git a/julia/test/find_closest.jl b/julia/test/find_closest.jl index f7870a3..9ccc622 100644 --- a/julia/test/find_closest.jl +++ b/julia/test/find_closest.jl @@ -1,6 +1,6 @@ @testset "Find Closest" begin - using NLsolve + using JuMP using Thesis: treat_inputs # Initial Setup @@ -15,25 +15,38 @@ # 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)[1][end,:] + final = prop(ΔVs, 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 - x0 = repeat([0.55, 0., 0.], n) - result = nlp_solve(start, final, sc, μs["Earth"], 0.0, prop_time, x0) + Tr, TΘ, Th = conv_T(repeat([0.6], n), repeat([0.], n), repeat([0.], n), + start, + sc.mass, + sc, + prop_time, + μs["Earth"]) + result, solution = nlp_solve(start, + final, + sc, + μs["Earth"], + 0.0, + prop_time, + Tr, + TΘ, + Th) # Test and plot - @test converged(result) + @test JuMP.termination_status(result) == MOI.OPTIMAL path1 = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1] - path2, mass = prop(treat_inputs(result.zero), start, sc, μs["Earth"], prop_time) - path3 = prop(zeros((100,3)), path2[end,:], sc, μs["Earth"], new_T)[1] + path2, mass, calc_final = prop(treat_inputs(JuMP.value.(solution)), start, sc, μs["Earth"], prop_time) + path3 = prop(zeros((100,3)), calc_final, sc, μs["Earth"], new_T)[1] path4 = prop(zeros((100,3)), final, sc, μs["Earth"], new_T)[1] savefig(plot_orbits([path1, path2, path3, path4], labels=["initial", "transit", "after transit", "final"], colors=["#FFFFFF","#FF4444","#44FF44","#4444FF"]), "../plots/find_closest_test.html") - if converged(result) - @test norm(path2[end,:] - final) < 1e-4 - end + # if termination_status(result) == :OPTIMAL + # @test norm(calc_final - final) < 1e-4 + # end end diff --git a/julia/test/monotonic_basin_hopping.jl b/julia/test/monotonic_basin_hopping.jl index 45e8373..419abe4 100644 --- a/julia/test/monotonic_basin_hopping.jl +++ b/julia/test/monotonic_basin_hopping.jl @@ -12,18 +12,18 @@ n = 25 # 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)[1][end,:] - new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"]) + # 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)[1][end,:] + # new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"]) # This should be close enough to 0.6 - best, archive = mbh(start, final, sc, μs["Earth"], 0.0, prop_time, n) + # best, archive = mbh(start, final, sc, μs["Earth"], 0.0, prop_time, n) # Test and plot - @test converged(best) - for path in archive - @test converged(path) - end + @test_skip converged(best) + #for path in archive + # @test_skip converged(path) + #end end diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index dc3d99d..9d25283 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -6,10 +6,10 @@ using Thesis # Tests @testset "All Tests" begin - include("spacecraft.jl") - include("laguerre-conway.jl") - include("propagator.jl") - include("plotting.jl") +# include("spacecraft.jl") +# include("laguerre-conway.jl") +# include("propagator.jl") +# include("plotting.jl") include("find_closest.jl") include("monotonic_basin_hopping.jl") end