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)
Created
May 7, 2020 03:31
-
-
Save mforets/fec21129ac7675963f65c6686bc94e17 to your computer and use it in GitHub Desktop.
EMBrake_scriptverif.md
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment