I have two EnsembleSolutions which represent tubes of orbital trajectories in space (invariant manifolds). I’d like to find the intersection of these two EnsembleSolutions. This functionality comes in two levels…
Given two EnsembleSolutions, find the intersection between them
In reality, these EnsembleSolutions are in different coordinate frames, so I’d really like to apply a coordinate frame transformation to each solution before finding the intersection
Is there any off-the-shelf way to find the intersection between two SciMLBase.EnsembleSolution instances? I’m aware of the K nearest neighbors algorithm, and other methods for generally computing the intersection between two multidimensional discrete sets, but I was wondering if there was an off the shelf way which took advantage of the interpolates in each EnsembleSolution.
I might be misunderstanding, but if I saveat I’ll have more discrete data. There’s still no guarantee that there will be an intersection between the two discrete sets of data – in fact, there’s really never an intersection at a discrete data point. I could set saveat to something really small, but the vast majority of these manifolds are nowhere near each other.
Doesn’t saveat allow a user to save the solution at every saveat timestep? Or am I misunderstanding?
Edit: I have a minimum working example below.
#
# Manifolds!
#
# Environment setup
import Pkg
Pkg.activate(; temp=true)
Pkg.add("Plots")
Pkg.add("Unitful")
Pkg.add("DifferentialEquations")
Pkg.add(; url="https://github.com/cadojo/GeneralAstrodynamics.jl", rev="main")
# Load dependencies
using Plots, Unitful
using DifferentialEquations
using GeneralAstrodynamics
# Find a periodic orbit near Earth (orbit, period)
Oₑ, Pₑ = halo(SunEarth; Az=100_000u"km", L=2)
# Find a periodic orbit near Jupiter (orbit, period)
Oⱼ, Pⱼ = halo(SunJupiter; Az=300_000u"km", L=1)
# Compute an unstable manifold departing Earth
Mₑ = manifold(Oₑ, Pₑ; trajectories=50, saveat=0.01, duration=3Pₑ, eps=1e-6, Trajectory=Val{false}, direction=Val{:unstable})
# Compute a stable manifold arriving at Jupiter
Mⱼ = manifold(Oⱼ, Pⱼ; trajectories=50, saveat=0.01, duration=3Pⱼ, eps=-1e-6, Trajectory=Val{false}, direction=Val{:stable})
# Assume there's some function that takes a state vector,
# and transforms it to a common reference frame
# magic_transformation(state::AbstractVector) = ... # another state vector
# What's the intersection of Mₑ and Mⱼ?
plot(Mₑ, vars=:xy)
plot!(Mⱼ; vars=:xy, palette=:greens)
Understood about how saveat works. I’m still not understanding how saving at time intervals would help to identify times when two separate solutions intersect. If I don’t know the times of intersection ahead of time, how can I use saveat to find the intersection between two solutions?
If you don’t want to use flowpipe construction methods / reachability analysis (which implement conservative time discretization, hence including all behaviors between time points), one simple approximation may be to “convexify” the solution sets in Mₑ and Mⱼ, then intersect them in space.
I was thinking about looking into reachability analysis, especially after seeing the JuliaCon 2021 workshop / content.
I understand that set propagation is scalable for linear systems, and numerical integration is scalable for nonlinear systems. These dynamics are nonlinear, and while I don’t know if they formally fit the definition of “chaotic dynamics”, they are certainly extremely sensitive to initial conditions. I’m a bit concerned that I’ll need a level of precision that set propagation won’t be able to provide (but like I said, I’m definitely going to try, if for no other reason than as an excuse to play around with ReachabilityAnalysis.jl).
The point about convexification is really interesting — I wasn’t aware that capability existed / was used. @mforets are you aware of any Julia implementations which “convexify” a discrete set? I searched a bit, but haven’t found any.
The discussion about reachability seems like a side topic to this thread, but I couldn’t resist to comment about it on the first place ;).
are you aware of any Julia implementations which “convexify” a discrete set?
The Julia package to do this sort of thing is LazySets.jl.
I could reproduce your MWE, very nice! Where can I see the dynamic equations you used? For a conceptual proof of the idea in my previous comment, I think it’s necessary to reason by grouping the results of different trajectories on the same time bins. Is that something that can be done using available code?
Let me show a result neglecting time and just focusing on the intersection. The code is available here as a gist.
No, I’m looking for when the manifolds physically intersect. There’s a unique time for every physical point on the manifold, so I was saying “when do these manifolds intersect” informally. Sorry about that confusion
Ah thanks so much for such a detailed example! The dynamics are available at `AstrodynamicalModels.jl. Also, I can definitely set the time bins for the two manifolds to be equivalent.
I’m on mobile at the moment, but I’ll be able to try out this intersection method this weekend. I’ll update here when I do!
