Convexifying a 6-dimensional closed trajectory with `LazySets.jl`

Specific Domains Numerics
1

First off, if anyone knows how to change the tag reachability-analysi to reachability-analysis, please do so! Discourse won’t let me type reachability-analysis in the category because a previous tag with an incorrect spelling was created.

I want to convexify a closed orbit (e.g. a 6-dimensional line which closes in on itself). How can I do this? In another post, @mforets kindly offered a detailed example for convexifying a manifold. I can’t figure out how to do so with a simple “line”, though.

I can approximate the plane of the orbit, as in the image below, but I can’t figure out how to only take the edge of the orbital plane.

Some pseudo code for what I’d like to do is below. Much of this code is taken directly from @mforets’ example provided in the previously linked Discourse post.


using Plots
using LazySets
using DifferentialEquations
using ReachabilityAnalysis
using GeneralAstrodynamics

points = let
  orbit, T = halo(SunEarth; Az = 0.005, L = 2)
  trajectory = propagate(orbit, T; abstol=1e-16)
  Matrix(solution(trajectory)')
end
	
@assert points[1,:] ≈ points[end,:]
	
convexified_points = [
  Singleton(row[1:2]) for row ∈ eachrow(points)
]
	
dom = box_approximation(UnionSetArray(convexified_points))
part = split(dom, [2, 2], [4, 1])
idx_dom = []
for D in part
  idx = findall(p -> p ⊆ D, convexified_points)
  isempty(idx) && continue
  push!(idx_dom, idx)
end

myset = UnionSetArray([ConvexHullArray(convexified_points[ii]) for ii in idx_dom])

plot(myset)

img

2

If you’d like to get a set that covers the egg-shaped region:

points, times, traj = let
  orbit, T = halo(SunEarth; Az = 0.005, L = 2)
  @show trajectory = propagate(orbit, T; abstol=1e-16)
  Matrix(solution(trajectory)'), trajectory.solution.t, trajectory
end

V = VPolytope(points');
πV = project(V, 1:2)
plot(πV, c=:red)
plot!(traj, vars=(1, 2), marker=:o, c=:blue)

Screenshot from 2021-09-06 22-02-36
Screenshot from 2021-09-06 22-02-361818×1236 199 KB

This computes one polygon in vertex representation, πV, that encloses all these points. I don’t know if this is useful for your application, is it?

To clarify my comment in the other thread: it was about the case in which there is not one but a bunch of trajectories – then one may just create (“convexify”) one set per each group of points among all trajectories falling inside each time bin [(k-1)\delta, k\delta] where \delta is some predefined time step. This method produces a flowpipe-like object that can be used to make geometric intersections.

4 Likes
3

To comment on the use of a reachability solver, I took the model from AstrodynamicalModels.jl:

LazySets.set_ztol(Float64, 1e-14) # for plots

@taylorize function f!(du, u, p, t)
    local μ = 3.003480593992993e-6
    
    x, y, z, vx, vy, vz = u

    du[1] = vx
    du[2] = vy
    du[3] = vz

    aux0 = sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3
    aux01 = (sqrt(y^2 + z^2 + (x + μ)^2)^-3)*(1 - μ)

    aux1 = μ*(x + μ - 1)*aux0
    aux2 = (x + μ)*aux01
    du[4] = x + 2*vy - aux1 - aux2

    aux3 = y*(μ*aux0 + aux01)
    du[5] = y - 2*vx - aux3

    aux4 = -μ*aux0 - aux01
    du[6] = z*aux4

    return du
end

U0 = Singleton([1.005279631460777, 0.0, 0.004591549059400884, 0.0, 0.019016771516337055, -0.0])
prob = @ivp(u' = f!(u), u(0) ∈ U0, dim=6)
alg = TMJets(orderQ=2, orderT=6, abstol=1e-12)

T = 3.0293993740463603
sol = solve(prob, T=T, alg=alg);

# time span of the first ten reachable sets computed
tspan.(sol)[1:10]
10-element Vector{IntervalArithmetic.Interval{Float64}}:
 [0, 0.00440149]
    [0.00440148, 0.00880456]
    [0.00880455, 0.0132125]
    [0.0132124, 0.0176284]
    [0.0176283, 0.0220558]
    [0.0220557, 0.0264979]
    [0.0264978, 0.0309584]
    [0.0309583, 0.0354411]
    [0.035441, 0.03995]
    [0.0399499, 0.0444894]

The solution can be plotted with plot(sol, vars=(1, 2), c=:red); plot!(traj, vars=(1, 2), marker=:o, c=:blue) (they overlap – taking a larger number of orderT, gives a coarser flowpipe hence the sets can be observed in the plot more clearly).

1 Like