# Intersection of `SciMLBase.EnsembleSolution`s

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…

1. Given two `EnsembleSolution`s, find the intersection between them
2. In reality, these `EnsembleSolution`s 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`.

Why not just use `saveat` so that the intersection is the whole thing?

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)

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)
``````

No, you can tell the solver to save at pre-defined times and intervals.

1 Like

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.

1 Like

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.

The region in magenta encloses the intersection between the manifolds (in the x-y plane). For example, the intersection area is obtained as

``````julia> sum(area, Xint)
0.0005889003496004155
``````
1 Like

Wait, intersection of `sol1.t` and `sol2.t`?

Wait, intersection of `sol1.t` and `sol2.t` ?

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!

1 Like