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

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)


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)


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

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