From 83886188db5dffd441577770e8fce5fce828f191 Mon Sep 17 00:00:00 2001 From: Connor Date: Thu, 2 Sep 2021 17:21:40 -0600 Subject: [PATCH] Got MBH working! It can find better-than-basic-spiral trajectories --- julia/Manifest.toml | 183 ++++++++--------------- julia/Project.toml | 3 +- julia/src/conversions.jl | 33 +++-- julia/src/find_closest.jl | 86 ++--------- julia/src/monotonic_basin_hopping.jl | 97 +++++++----- julia/src/propagator.jl | 204 +++++++++++++------------- julia/test/find_closest.jl | 28 ++-- julia/test/monotonic_basin_hopping.jl | 65 ++++++-- 8 files changed, 318 insertions(+), 381 deletions(-) diff --git a/julia/Manifest.toml b/julia/Manifest.toml index 0cefb85..0caedc8 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 = "85d03b60274807181bae7549bb22b2204b6e5a0e" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "3.1.30" + [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -21,60 +21,24 @@ 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" @@ -108,12 +72,6 @@ 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" @@ -139,6 +97,12 @@ 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" @@ -156,6 +120,12 @@ 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" @@ -186,6 +156,11 @@ 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" @@ -196,18 +171,6 @@ 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" @@ -236,18 +199,6 @@ 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" @@ -284,6 +235,12 @@ 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" @@ -297,18 +254,6 @@ 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" @@ -319,18 +264,6 @@ 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" @@ -353,18 +286,24 @@ 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" @@ -378,12 +317,6 @@ 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" @@ -395,6 +328,12 @@ 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" @@ -475,6 +414,12 @@ git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" version = "1.6.1" +[[Static]] +deps = ["IfElse"] +git-tree-sha1 = "854b024a4a81b05c0792a4b45293b85db228bd27" +uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +version = "0.3.1" + [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb" @@ -485,6 +430,11 @@ 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" @@ -509,12 +459,6 @@ 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" @@ -530,6 +474,11 @@ 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" @@ -551,12 +500,6 @@ 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 8c25aa9..01d2154 100644 --- a/julia/Project.toml +++ b/julia/Project.toml @@ -5,9 +5,8 @@ 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/conversions.jl b/julia/src/conversions.jl index 8e0d6f3..cb66a64 100644 --- a/julia/src/conversions.jl +++ b/julia/src/conversions.jl @@ -72,9 +72,9 @@ 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}, +function conv_T(Tm::Vector{Float64}, + Ta::Vector{Float64}, + Tb::Vector{Float64}, init_state::Vector{Float64}, m::Float64, craft::Sc, @@ -84,29 +84,32 @@ function conv_T(Tr::Vector{Float64}, Txs = Float64[] Tys = Float64[] Tzs = Float64[] - state = init_state - for i in 1:length(Tr) - mag, α, β = Tr[i], TΘ[i], Th[i] + n::Int = length(Tm) + for i in 1:n + mag, α, β = Tm[i], Ta[i], Tb[i] if mag > 1 || mag < 0 - throw(ErrorException("ΔV input is too high: $mag")) + @error "ΔV input is too high: $mag" elseif α > π || α < -π - throw(ErrorException("α angle is incorrect: $α")) + @error "α angle is incorrect: $α" elseif β > π/2 || β < -π/2 - throw(ErrorException("β angle is incorrect: $β")) + @error "β angle is incorrect: $β" end + end + state = init_state + for i in 1:n + mag, α, β = Tm[i], Ta[i], Tb[i] 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), μ) + sΩ*cθ+cΩ*ci*sθ -sΩ*sθ+cΩ*ci*cθ -cΩ*si; + si*sθ si*cθ ci ] + Tx, Ty, Tz = DCM*thrust_rθh + + state = prop_one([Tx, Ty, Tz], state, craft, μ, time/n)[1] push!(Txs, Tx) push!(Tys, Ty) push!(Tzs, Tz) diff --git a/julia/src/find_closest.jl b/julia/src/find_closest.jl index dc4b448..518c95b 100644 --- a/julia/src/find_closest.jl +++ b/julia/src/find_closest.jl @@ -1,6 +1,15 @@ -using JuMP, Ipopt +using NLsolve -export nlp_solve +export nlp_solve, mass_est + +function mass_est(T) + ans = 0 + n = Int(length(T)/3) + for i in 1:n + ans += norm(T[i,:]) + end + return ans/n +end function nlp_solve(start::Vector{Float64}, final::Vector{Float64}, @@ -8,77 +17,14 @@ function nlp_solve(start::Vector{Float64}, μ::Float64, t0::Float64, tf::Float64, - Txi::Vector{Float64}, - Tyi::Vector{Float64}, - Tzi::Vector{Float64}; - solver_options=(), + x0::Matrix{Float64}; tol=1e-6) - n::Int = length(Txi) - if length(Tyi) != n - throw("Bad number of Tys") - elseif length(Tzi) != n - throw("Bad number of Tzs") + function f!(F,x) + F .= 0.0 + F[1:6, 1] .= prop_nlsolve(tanh.(x), start, craft, μ, tf-t0) .- final end - # 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[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) - - # 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) - - # 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(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] == 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) - end - - optimize!(model) - return model + return nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=1_000) end \ No newline at end of file diff --git a/julia/src/monotonic_basin_hopping.jl b/julia/src/monotonic_basin_hopping.jl index cca4c08..dbb3c41 100644 --- a/julia/src/monotonic_basin_hopping.jl +++ b/julia/src/monotonic_basin_hopping.jl @@ -1,19 +1,20 @@ -function perturb(x::AbstractVector, n::Int) - perturb_vector = 0.02 * rand(Float64, (3n)) .- 0.01 - return x + perturb_vector +function pareto(α::Float64, n::Int) + s = rand((-1,1), (n,3)) + r = rand(Float64, (n,3)) + ϵ = 1e-10 + return (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) end -function mass_better(x_star::AbstractVector, - x_current::AbstractVector, - start::AbstractVector, - final::AbstractVector, - craft::Sc, - μ::AbstractFloat, - t0::AbstractFloat, - tf::AbstractFloat) - mass_star = prop(treat_inputs(x_star), start, craft, μ, tf-t0)[2] - mass_current = prop(treat_inputs(x_current), start, craft, μ, tf-t0)[2] - return mass_star > mass_current +function perturb(x::AbstractMatrix, n::Int) + ans = x + pareto(1.01, n) + map!(elem -> elem > 1.0 ? 1.0 : elem, ans, ans) + map!(elem -> elem < -1.0 ? -1.0 : elem, ans, ans) + return ans +end + + +function new_x(n::Int) + 2.0 * rand(Float64, (n,3)) .- 1. end function mbh(start::AbstractVector, @@ -22,35 +23,51 @@ function mbh(start::AbstractVector, μ::AbstractFloat, t0::AbstractFloat, tf::AbstractFloat, - n::Int, - num_iters::Int=10, - tol=1e-6) - i::Int = 0 + n::Int; + num_iters=50, + patience_level::Int=400, + tol=1e-6, + verbose=false) + archive = [] - x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol) - while converged(x_star) == false - x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol) - end - x_current = x_star - push!(archive, x_current) - - while i < num_iters - x_star = nlp_solve(start, final, craft, μ, t0, tf, perturb(x_current.zero,n), tol=tol) - if converged(x_star) && mass_better(x_star.zero, x_current.zero, start, final, craft, μ, t0, tf) - x_current = x_star - push!(archive, x_star) - else - while converged(x_star) == false - x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol) - end - if mass_better(x_star.zero, x_current.zero, start, final, craft, μ, t0, tf) - x_current = x_star - push!(archive, x_star) - end - end + i = 0 + if verbose println("Current Iteration") end + while true i += 1 + if verbose print("\r",i) end + impatience = 0 + x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol) + while converged(x_star) == false + x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol) + end + if converged(x_star) + x_current = x_star + while impatience < patience_level + x_star = nlp_solve(start, final, craft, μ, t0, tf, perturb(tanh.(x_current.zero),n), tol=tol) + if converged(x_star) && mass_est(tanh.(x_star.zero)) < mass_est(tanh.(x_current.zero)) + x_current = x_star + impatience = 0 + else + impatience += 1 + end + end + push!(archive, x_current) + end + if i >= num_iters + break + end end - return x_current, archive + current_best_mass = 1e8 + best = archive[1] + for candidate in archive + if mass_est(tanh.(candidate.zero)) < current_best_mass + current_best_mass = mass_est(tanh.(candidate.zero)) + best = candidate + end + end + + return best, archive + end \ No newline at end of file diff --git a/julia/src/propagator.jl b/julia/src/propagator.jl index 5294285..fbf0c85 100644 --- a/julia/src/propagator.jl +++ b/julia/src/propagator.jl @@ -16,38 +16,19 @@ end This function propagates the spacecraft forward in time 1 Sim-Flanagan step (of variable length of time), applying a thrust in the center. """ -function prop_one(ΔV::Vector{T}, - state::Vector{T}, +function prop_one(thrust_unit::Vector{<:Real}, + state::Vector{<:Real}, duty_cycle::Float64, num_thrusters::Int, max_thrust::Float64, mass::T, mass_flow_rate::Float64, μ::Float64, - time::Float64) where T + time::Float64) where T<:Real - mag, α, β = ΔV - - # 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(β)] - a,e,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 - - thrust = max_ΔV(duty_cycle, num_thrusters, max_thrust, time, 0., mass) * ΔV - halfway = laguerre_conway(state, μ, time/2) + [0., 0., 0., thrust[1], thrust[2], thrust[3]] - return laguerre_conway(halfway, μ, time/2), mass - mass_flow_rate*norm(ΔV)*time + ΔV = max_ΔV(duty_cycle, num_thrusters, max_thrust, time, 0., mass) * thrust_unit + halfway = laguerre_conway(state, μ, time/2) + [0., 0., 0., ΔV...] + return laguerre_conway(halfway, μ, time/2), mass - mass_flow_rate*norm(thrust_unit)*time end @@ -64,6 +45,99 @@ 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 + +""" +This propagates over a given time period, with a certain number of intermediate steps +""" +function prop(ΔVs::Matrix{T}, + state::Vector{Float64}, + duty_cycle::Float64, + num_thrusters::Int, + max_thrust::Float64, + mass::Float64, + mass_flow_rate::Float64, + μ::Float64, + time::Float64) where T <: Real + + if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end + n = size(ΔVs)[i] + + for i in 1:n + state, mass = prop_one(ΔVs[i,:], state, duty_cycle, num_thrusters, max_thrust, mass, + mass_flow_rate, μ, time/n) + end + + return state, mass + +end + +""" +The same function, using Scs +""" +function prop(ΔVs::Matrix{T}, + state::Vector{S}, + craft::Sc, + μ::Float64, + time::Float64) where {T <: Real, S <: Real} + + if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end + n = size(ΔVs)[1] + + 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) + 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 [x_states, y_states, z_states, dx_states, dy_states, dz_states], masses, state + +end + +function prop_nlsolve(ΔVs::Matrix{T}, + state::Vector{S}, + craft::Sc, + μ::Float64, + time::Float64) where {T <: Real, S <: Real} + + n = size(ΔVs)[1] + for i in 1:n + state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n) + end + + return 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 + function prop_one_simple(Tx, Ty, Tz, x, y, z, dx, dy, dz, t, μ) # perform laguerre_conway once @@ -154,83 +228,5 @@ function prop_one_simple(Tx, Ty, Tz, x, y, z, dx, dy, dz, t, μ) 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 -""" -function prop(ΔVs::Matrix{T}, - state::Vector{Float64}, - duty_cycle::Float64, - num_thrusters::Int, - max_thrust::Float64, - mass::Float64, - mass_flow_rate::Float64, - μ::Float64, - time::Float64) where T <: Real - - if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end - n = size(ΔVs)[i] - - for i in 1:n - state, mass = prop_one(ΔVs[i,:], state, duty_cycle, num_thrusters, max_thrust, mass, - mass_flow_rate, μ, time/n) - end - - return state, mass - -end - -""" -The same function, using Scs -""" -function prop(Δ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] - - 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) - 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 [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 +end \ No newline at end of file diff --git a/julia/test/find_closest.jl b/julia/test/find_closest.jl index f5aa331..f053130 100644 --- a/julia/test/find_closest.jl +++ b/julia/test/find_closest.jl @@ -1,6 +1,6 @@ @testset "Find Closest" begin - using JuMP + using NLsolve, PlotlyJS # Initial Setup sc = Sc("test") @@ -9,10 +9,11 @@ i = rand(0.01:0.01:π/6) T = 2π*√(a^3/μs["Earth"]) prop_time = 2T - n = 30 + n = 20 # A simple orbit raising start = oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"]) +# T_craft = hcat(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, @@ -23,35 +24,26 @@ new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"]) # This should be close enough to 0.6 - Tx, Ty, Tz = conv_T(repeat([0.6], n), repeat([0.], n), repeat([0.], n), + Tx, Ty, Tz = conv_T(repeat([0.59], n), repeat([0.01], 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, - Tx, - Ty, - Tz) - # solver_options=("max_cpu_time" => 30.)) + result = nlp_solve(start, final, sc, μs["Earth"], 0.0, prop_time, hcat(Tx, Ty, Tz)) # Test and plot - @test JuMP.termination_status(result) == MOI.OPTIMAL + @test converged(result) path1 = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1] - path2, mass, calc_final = prop(treat_inputs(JuMP.value.(solution)), start, sc, μs["Earth"], prop_time) + path2, mass, calc_final = prop(tanh.(result.zero), 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 termination_status(result) == :OPTIMAL - # @test norm(calc_final - final) < 1e-4 - # end + if converged(result) + @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 419abe4..452053d 100644 --- a/julia/test/monotonic_basin_hopping.jl +++ b/julia/test/monotonic_basin_hopping.jl @@ -8,22 +8,63 @@ e = rand(0.01:0.01:0.5) i = rand(0.01:0.01:π/6) T = 2π*√(a^3/μs["Earth"]) - prop_time = 2T - n = 25 + prop_time = 0.75T + n = 10 # 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"]) +# T_craft = hcat(repeat([0.6], n), repeat([0.], n), repeat([0.], n)) + Tx, Ty, Tz = conv_T(repeat([0.8], n), repeat([0.], n), repeat([0.], n), + start, + sc.mass, + sc, + prop_time, + μs["Earth"]) + nominal_path, normal_mass, final = prop(hcat(Tx, Ty, Tz), start, sc, μs["Earth"], prop_time) + 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) + # Find the best solution + best, archive = mbh(start, + final, + sc, + μs["Earth"], + 0.0, + prop_time, + n, + num_iters=5, + patience_level=50, + verbose=true) # Test and plot - @test_skip converged(best) - #for path in archive - # @test_skip converged(path) - #end + @test converged(best) + transit, best_masses, calc_final = prop(tanh.(best.zero), start, sc, μs["Earth"], prop_time) + initial_path = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1] + after_transit = prop(zeros((100,3)), calc_final, sc, μs["Earth"], new_T)[1] + final_path = prop(zeros((100,3)), final, sc, μs["Earth"], new_T)[1] + savefig(plot_orbits([initial_path, nominal_path, final_path], + labels=["initial", "nominal transit", "final"], + colors=["#FF4444","#44FF44","#4444FF"]), + "../plots/mbh_nominal.html") + savefig(plot_orbits([initial_path, transit, after_transit, final_path], + labels=["initial", "transit", "after transit", "final"], + colors=["#FFFFFF", "#FF4444","#44FF44","#4444FF"]), + "../plots/mbh_best.html") + i = 0 + best_mass = best_masses[end] + nominal_mass = normal_mass[end] + masses = [] + for candidate in archive + @test converged(candidate) + path2, cand_ms, calc_final = prop(tanh.(candidate.zero), start, sc, μs["Earth"], prop_time) + push!(masses, cand_ms[end]) + @test norm(calc_final - final) < 1e-4 + end + @test best_mass == maximum(masses) + + # This won't always work since the test is reduced in fidelity, + # but hopefully will usually work: + @test (sc.mass - best_mass) < 1.1 * (sc.mass - nominal_mass) + @show best_mass + @show nominal_mass end