Week #5: Visualizing a limit cycle

Summary of contributions

• LazySets#2257

• LazySets#2258

• LazySets#2284

• NeuralVerification#125

• NeuralVerification#132

Visualizing properties using ReachabilityAnalysis.jl

Limit cycle

A limit cycle is an isolated closed trajectory. Isolated means that neighboring trajectories are not closed, they spiral away or towards the limit cycle. A limit cycle is deemed stable if all neighboring trajectories approach it, otherwise it is unstable. The Van der Pol oscillator has a stable limit cycle, we can visualize it watching how the flowpipe evolves with time.

First let's define the dynamics of the system:

using ReachabilityAnalysis, Plots, Optim, Polyhedra

@taylorize function vanderpol!(dx, x, params, t)

    local μ = 1.0

    dx[1] = x[2]

    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]

    return dx

end

Then I have to choose a big enough set of initial conditions that will allow us to see the flowpipe shrink with time:

X0 = Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45])

prob = @ivp(x' = vanderpol!(x), dim=2, x(0) ∈ X0)

sol = solve(prob, T=7.0, alg=TMJets());

The reachability was done with the non-linear algorithm TMJets

To further analyse the flowpipe, I will overapproximate it with zonotopes.

solz = overapproximate(sol, Zonotope);

plot(solz, vars=(1, 2), lw=0.2, xlims=(-2.5, 2.5), xlab="x", ylab="y")

plot!(x -> 2.75, color=:red, lab="y = 2.75", style=:dash, legend=:bottomright)

I will highlight the list few sets and the last set such that they overlap, then I will intersect a line somewhat perpendicular to the trajectory, that will allow me to get a cross-section of the sets.

plot(solz, vars=(1, 2), lw=0.2, xlims=(0.0, 2.5), ylims=(1.6, 2.8), xlab="x", ylab="y")

plot!(X0, color=:orange, lab="X0")

plot!(solz[1:5], vars=(1, 2), color=:green, lw=1.0, alpha=0.5, lab="F[1:5]")

plot!(solz[200], vars=(1, 2), color=:red, lw=1.0, alpha=0.6, lab="F[200]")

plot!(LineSegment([1, 2.], [2., 2.5]), lw=2.0)

Here I define a function to get the cross section of the flowpipe, the function needs the flowpipe, a line segment that cuts the flowpipe and the indices of the subsets to cut

using ReachabilityAnalysis: ReachSolution

function cross_section(line::LineSegment, RS::ReachSolution, idx)

    i = reduce(convex_hull, map(x -> intersection(line, x), set.(RS[idx])))

    vl = vertices_list(i)

    return LineSegment(vl[1], vl[2])

end

Then I get the cross section of the first five sets and the last set, I will call them i1 and i2 respectively

i1 = cross_section(line, solz, 1:5)

i2 = cross_section(line, solz, [200])

Then I calculate the length of each cross section, remember that the system is 2D, so the cross section will be a line segment.

l1 = norm(i1.q - i1.p)

l2 = norm(i2.q - i2.p);

@show l1

@show l2;

plot(i1, lw=3.0, alpha=1.0, label="First subsets", legend=:bottomright)

plot!(i2, lw=5.0, alpha=1.0, label="Last subset")

i2 ⊆ i1

As we can see, the cross section of the las subset is a subset of the first few sets, thus, the cycle will continue, presumably getting smaller each revolution.