EMBrake_scriptverif.md

EMBrake_scriptverif.md

using Plots, BenchmarkTools, StaticArrays
using Revise, ReachabilityAnalysis
const RA = ReachabilityAnalysis
include("/home/mforets/.julia/dev/ReachabilityAnalysis/test/models/hybrid/embrake.jl")
LazySets.deactivate_assertions()

function LazySets.diameter(R::ReachabilityAnalysis.AbstractLazyReachSet; vars::Int)
    overapproximate(project(R, vars=vars), Interval) |> set |> diameter
end

function final_diams(sol)
    println("final diameter of current I = ",  diameter(sol[end][end], vars=(1)))
    println("final diameter of position x = ",  diameter(sol[end][end], vars=(2)))
end

function velocity2(sol; α=1.7680855258330557e-5, ngens=4)
    # [I, x, xe, xc]
    # we project along (x, v)
    M = SMatrix{2, 4}([0 1 0 0;
                       α 0 0 0])

    ZT = Zonotope{Float64, SVector{2, Float64}, SMatrix{2, ngens, Float64, 2*ngens}}
    numsets = length(sol[1]) * length(sol)
    RT = ReachSet{Float64, ZT}
    out = Vector{RT}()
    sizehint!(out, numsets)

    for Fi in sol.F
        for (i, Ri) in enumerate(Fi)
            πR = linear_map(M, set(Ri))
            Δt = tspan(Fi, i) # porque Fi es un ShiftedFlowpipe
            push!(out, ReachSet(πR, Δt))
        end
    end
    return Flowpipe(out)
end



function safety_property(sol; x0=0.05, ε=0.002, ngens=4)
    # project the solution along variables (x, v)
    vsol = velocity2(sol, ngens=ngens)

    # define the safe set P
    Hlow = HalfSpace([-1.0, 0.0], -(x0 - ε)) # x >= x0 - ε
    Hhigh = HalfSpace([1.0, 0.0], x0 + ε);   # x <= x0 + ε
    P = HPolyhedron([Hlow, Hhigh]);

    # find the first set that is inside P
    @show first_in = findfirst(X -> set(X) ⊆ P, vsol.Xk)

    # check that all states starting from first_in are contained in P
    @show all(X -> set(X) ⊆ P, vsol.Xk[first_in:end])

    # check the time span where the flowpipe enters P for the first time
    @show tspan(vsol[first_in]) # 86ms

    # check maximum values of x at different parts of the flowpipe
    println("max values of x")
    @show ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))
    @show ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))
    @show ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end]))

    # check minimum values of x at different parts of the flowpipe
    println("min values of x")
    @show -ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))
    @show -ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))
    @show -ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))

    # check maximum values of v at different parts of the flowpipe
    println("max values of v")
    @show ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in]))
    @show ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end]))
    @show ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end]))

    # check minimum values of v at different parts of the flowpipe
    println("min values of v")
    @show -ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))
    @show -ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))
    @show -ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))
end



#sol = nothing
sol2 = nothing
sol = nothing
GC.gc()



prob = embrake_no_pv(ζ=0.0, Tsample=1e-4);



@time sol = solve(prob, alg=GLGM06(δ=1e-8, max_order=1, static=true, dim=4, ngens=4), max_jumps=2);
final_diams(sol)

safety_property(sol)



function run_once(prob)
    @time sol = solve(prob, alg=GLGM06(δ=1e-9, max_order=1, static=true, dim=4, ngens=4), max_jumps=250);
    @show final_diams(sol)
    return sol
end


sol1 = run_once(prob);

prob2 = IVP(prob.s, set(sol1[end][end]));
sol2 = run_once(prob2);

prob2 = IVP(prob.s, set(sol1[end][end]));
sol2 = run_once(prob2);

prob = embrake_no_pv(ζ=[-1e-8, 1e-7], Tsample=1e-4);

# GLGM06, order 1, static
@time sol = solve(prob, alg=GLGM06(δ=1e-8, max_order=1, static=true, dim=4, ngens=4), max_jumps=1000);
final_diams(sol)

safety_property(sol, ngens=4)