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

**URL:** <https://discourse.julialang.org/t/convexifying-a-6-dimensional-closed-trajectory-with-lazysets-jl/67764>\
**Category:** Numerics\
**Tags:** reachabilityanalysis\
**Created:** [September 6, 2021, 4:30pm UTC](https://discourse.julialang.org/t/convexifying-a-6-dimensional-closed-trajectory-with-lazysets-jl/67764 "2021-09-06T16:30:30Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [September 6, 2021, 4:30pm UTC](https://discourse.julialang.org/t/convexifying-a-6-dimensional-closed-trajectory-with-lazysets-jl/67764/1 "2021-09-06T16:30:30Z")

</div>

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](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/7), @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.

```julia

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](https://global.discourse-cdn.com/julialang/original/3X/4/d/4d8a107eb0f6cca8bbd643e3edae0d933126711a.png)

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [September 7, 2021, 1:10am UTC](https://discourse.julialang.org/t/convexifying-a-6-dimensional-closed-trajectory-with-lazysets-jl/67764/2 "2021-09-07T01:10:59Z")

</div>

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/a/d/ad12849af492998cb721da94f33007f52bd2d6c3.png)

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.

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [September 7, 2021, 11:01am UTC](https://discourse.julialang.org/t/convexifying-a-6-dimensional-closed-trajectory-with-lazysets-jl/67764/3 "2021-09-07T11:01:35Z")

</div>

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

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