Skip to content

Instantly share code, notes, and snippets.

@mforets

mforets/EMBrake_new.md

Last active May 7, 2020 19:59
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mforets/60ed38d4f48a311597f678159bf307ba to your computer and use it in GitHub Desktop.
Save mforets/60ed38d4f48a311597f678159bf307ba to your computer and use it in GitHub Desktop.
Download ZIP
EMBrake_new
Raw
EMBrake_new.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_ρ([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 = not# BETTER function run(prob, k, δ; ngens=4, max_order=1) @time sol = solve(prob, alg=GLGM06(δ=δ, max_order=max_order, static=true, dim=4, ngens=ngens), max_j

NEW RESULTS

-----------------------------------

NO PV

-----------------------------------

GLGM06, delta 1e-7, no jitter

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

julia> sol = run(prob, 1_000, 1e-7);
  0.186884 seconds (212.45 k allocations: 223.639 MiB, 4.15% gc time)

julia> analyze(sol)
final diameter of current I = 13.707464080123568
final diameter of position x = 0.0007351923316104431
tspan(vsol[first_in]) = [0.0887895, 0.0887897]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04850572743425879
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0492794261847055
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.0492794261847055
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.048000002027475806
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.048000002027475806
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04854423385309506
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0007994567734929507
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0007994567734929507
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0005874379912341613
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0006327526764030105
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0003450783028747313
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.0003450783028747313

GLGM06, delta 1e-8, no jitter

julia> sol = run(prob, 1_000, 1e-8);
  0.962470 seconds (192.44 k allocations: 1.681 GiB, 21.09% gc time)

julia> analyze(sol)
final diameter of current I = 1.3690297645113247
final diameter of position x = 7.34268190453885e-5
tspan(vsol[first_in]) = [0.0858076, 0.0858077]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04804576890888626
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.048945954868741794
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.048945954868741794
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.04800000000863762
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.04800000000863762
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.048872528049696405
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.000810869838680885
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.000810869838680885
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.00047936971342896125
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0007957831211492866
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0004551640963182902
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.0004551640963182902

GLGM06, delta 1e-9, no jitter

julia> sol = run(prob, 1_000, 1e-9);
 11.436564 seconds (182.44 k allocations: 16.441 GiB, 26.82% gc time)

for the analysis it is sufficient to compute the first 500 jumps..

julia> sol = run(prob, 500, 1e-9);
 11.005672 seconds (20.75 M allocations: 9.260 GiB, 11.05% gc time)

julia> analyze(sol, ε=0.01)
final diameter of current I = 0.025543376754228575
final diameter of position x = 1.3703165056694333e-6
tspan(vsol[first_in]) = [0.044504, 0.0445041]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04000112239745922
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.04193973300622029
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.04193973300622029
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.04000000034658299
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.04000000034658299
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04193836268971462
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.003844712032707981
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.003844712032707981
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0031099220509360377
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.003844342161346994
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0031094704221888373
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.0031094704221888373

julia> analyze(sol, ε=0.02)
final diameter of current I = 0.025543376754228575
final diameter of position x = 1.3703165056694333e-6
tspan(vsol[first_in]) = [0.0264032, 0.0264033]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.030000540301651234
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.04193973300622029
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.04193973300622029
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.030000000287769595
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.030000000287769595
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04193836268971462
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.007635381215482037
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.007635381215482037
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0031099220509360377
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.007635202833590987
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0031094704221888373
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.0031094704221888373

GLGM06, delta 1e-7, with jitter

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

julia> sol = run(prob, 1_000, 1e-7);
  1.546368 seconds (781.20 k allocations: 252.965 MiB, 2.69% gc time)

<NO VERIFICA PARA eps = 0.002 >

julia> analyze(sol, ε=0.005)
final diameter of current I = 54.714099517638516
final diameter of position x = 0.002934566941062061
tspan(vsol[first_in]) = [0.0647861, 0.0647863]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04590597079607615
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.05038266992049454
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.05038266992049454
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.04500000591771298
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.04500000591771298
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04744810297943248
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0019268337748910233
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0019268337748910233
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0009485891432081041
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0016282491695393281
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = -1.8802930953156145e-5
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = -1.8802930953156145e-5

GLGM06, delta1e-8, with jitter

julia> sol = run(prob, 1_000, 1e-8);
  1.011069 seconds (271.44 k allocations: 1.691 GiB, 19.08% gc time)

julia> analyze(sol, ε=0.002)
final diameter of current I = 17.746604024816392
final diameter of position x = 0.0009518262957221718
tspan(vsol[first_in]) = [0.0900594, 0.0900595]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04868287562433759
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.04938757641401575
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.04938757641401575
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.04800000018332063
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.04800000018332063
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.048435750118293576
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0007950689857855407
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0007950689857855407
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0006232115341202552
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0005699644062819035
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.00030943639703057016
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.00030943639703057016

GLGM06, delta1e-9, with jitter

julia> sol = run_split(prob, 500, 1e-9);
  4.486814 seconds (131.43 k allocations: 8.260 GiB, 15.82% gc time)
  3.543590 seconds (491.75 k allocations: 8.278 GiB, 23.98% gc time)

julia> sol = run(prob, 500, 1e-9);
  5.333848 seconds (621.16 k allocations: 8.284 GiB, 29.06% gc time)

julia> analyze(sol, ε=0.01)
final diameter of current I = 3.07825752613752
final diameter of position x = 0.0001651368400567288
tspan(vsol[first_in]) = [0.0446576, 0.0446577]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.040135932226251995
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.042029402282405356
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.042029402282405356
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.040000000265738404
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.040000000265738404
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04186426544234863
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.003841292610311951
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.003841292610311951
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0031339388579762107
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0037964848221430296
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0030795126322087064
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.0030795126322087064


julia> analyze(sol, ε=0.02)
final diameter of current I = 3.07825752613752
final diameter of position x = 0.0001651368400567288
tspan(vsol[first_in]) = [0.0264323, 0.0264324]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0300650756225516
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.042029402282405356
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.042029402282405356
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.03000000026506613
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.03000000026506613
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04186426544234863
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.007633674947251064
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.007633674947251064
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0031339388579762107
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.007612179929403613
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0030795126322087064
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.0030795126322087064

------------------------------------

PV1, delta = 1e-8 in all cases

------------------------------------

julia> prob = embrake_pv_1(ζ=0.0, Tsample=1e-4);

pv1, no jitter, order 1

julia> sol = run_pv(prob, 1000, 1e-8);
  8.761450 seconds (110.43 M allocations: 10.233 GiB, 14.67% gc time)

julia> analyze(sol, ε=0.005)
final diameter of current I = 137.24871950836993
final diameter of position x = 0.007305231540715523
tspan(vsol[first_in]) = [0.0705282, 0.0705283]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04772041880648274
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.052561856370724205
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.052561856370724205
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.045000000277171784
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.045000000277171784
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04525662483000868
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0018929469979134374
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0018929469979134374
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.001680604592986933
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0009743555252359067
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = -0.0007460701510317649
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = -0.0007460701510317649

pv1, no jitter, order 2

julia> sol = run_pv(prob, 1_000, 1e-8, max_order=2, ngens=8);
 38.064782 seconds (408.67 M allocations: 38.310 GiB, 14.08% gc time)

julia> analyze(sol, ε=0.002, ngens=8)
final diameter of current I = 4.247976812392487
final diameter of position x = 0.00018632095848780672
tspan(vsol[first_in]) = [0.0869794, 0.0869795]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.048229899529071606
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.04900240107961035
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.04900240107961035
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.048000000540886344
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.048000000540886344
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04881608012112254
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0008156348818995737
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0008156348818995737
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0005048211525579121
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0007211955880443044
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.000429713289397256
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.000429713289397256

julia> analyze(sol, ε=0.005, ngens=8)
final diameter of current I = 4.247976812392487
final diameter of position x = 0.00018632095848780672
tspan(vsol[first_in]) = [0.0632176, 0.0632177]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04530511435248417
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.04900240107961035
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.04900240107961035
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.04500000088282174
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.04500000088282174
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04881608012112254
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0019591495774039388
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0019591495774039388
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0005048211525579121
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0018237933999296018
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.000429713289397256
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.000429713289397256

pv1, no jitter, order3

julia> sol = run_pv(prob, 1000, 1e-8, max_order=3, ngens=12);
 41.275445 seconds (400.72 M allocations: 45.447 GiB, 21.44% gc time)

julia> analyze(sol, ε=0.002, ngens=12)
final diameter of current I = 2.94144391861186
final diameter of position x = 0.00012295781449649101
tspan(vsol[first_in]) = [0.0865408, 0.0865409]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.048161937383422726
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.04897071950761469
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.04897071950761469
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.04800000005867319
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.04800000005867319
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.0488477616931182
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0008165062855391979
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0008165062855391979
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0004932708430653205
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.000746079933035146
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0004412635988898476
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.0004412635988898476


julia> analyze(sol, ε=0.005, ngens=12)
final diameter of current I = 2.94144391861186
final diameter of position x = 0.00012295781449649101
tspan(vsol[first_in]) = [0.0630525, 0.0630526]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04524247328829584
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.04897071950761469
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.04897071950761469
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0450000011747487
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0450000011747487
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.0488477616931182
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0019599372628811428
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0019599372628811428
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0004932708430653205
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0018467174391622986
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.0004412635988898476
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = 0.0004412635988898476

pv1, with jitter, order 1

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


julia> sol = run_pv(prob, 1000, 1e-8, max_order=1, ngens=4);
  8.719063 seconds (110.52 M allocations: 10.243 GiB, 15.37% gc time)

julia> analyze(sol, ε=0.006, ngens=4)
final diameter of current I = 154.207482098115
final diameter of position x = 0.008210036200221565
tspan(vsol[first_in]) = [0.0634428, 0.0634429]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.04640244452887622
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.05301736140115831
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.05301736140115831
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.044000001406632294
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.044000001406632294
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.04480732520093674
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0022822198430200577
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0022822198430200577
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.0018293238731492435
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.0014637541655901564
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = -0.0008971962975791285
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = -0.0008971962975791285

------------------------------------

PV2

------------------------------------

julia> prob = embrake_pv_2(ζ=0.0, Tsample=1e-4);

pv2, no jitter, order 1

julia> sol = run_pv(prob, 1000, 1e-8, max_order=1, ngens=4);
  9.868345 seconds (116.07 M allocations: 10.508 GiB, 13.49% gc time)

NO VERIFICA

pv2, no jitter, order 2

julia> sol = run_pv(prob, 1000, 1e-8, max_order=2, ngens=8);
 37.802983 seconds (408.67 M allocations: 38.310 GiB, 14.24% gc time)

julia> analyze(sol, ε=0.02, ngens=8)
final diameter of current I = 883.8561589564646
final diameter of position x = 0.03551083334722416
tspan(vsol[first_in]) = [0.084553, 0.0845531]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.06601599126763294
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0667406691547255
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.06672200378026918
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.030000000089024127
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.030000000089024127
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.031211170433045027
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.008802205420732711
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.008802205420732711
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.008258793407075833
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = -0.007184654970621208
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = -0.00736879504116029
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = -0.007368539408617422


julia> analyze(sol, ε=0.05, ngens=8)
julia> sol = nothing; GC.gc(); sol = run_pv(prob, 1000, 1e-8, max_order=3, ngens=12);
 41.129208 seconds (411.24 M allocations: 45.861 GiB, 19.64% gc time)

julia> analyze(sol, ε=0.02, ngens=12)
final diameter of current I = 574.0077157938906
final diameter of position x = 0.022593942314352897
tspan(vsol[first_in]) = [0.0580367, 0.0580368]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.05847325209977274
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.060613313128198155
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.060263558263833554
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.030000001614109116
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.030000001614109116
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.03766961594948066
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.008901796674425007
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.008901796674425007
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.005519600669287567
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = -0.004423141273818343
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = -0.00493014619883479
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = -0.004629346670829155
final diameter of current I = 883.8561589564646
final diameter of position x = 0.03551083334722416
tspan(vsol[first_in]) = [0, 1.00001e-08]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.0667406691547255
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.06672200378026918
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = -0.0
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = -1.1191039075817737e-14
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.031211170433045027
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.0
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.019055876291717594
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.008258793407075833
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = -0.0
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = -0.00736879504116029
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = -0.007368539408617422

pv2, no jiter, order 3

julia> sol = nothing; GC.gc(); sol = run_pv(prob, 1000, 1e-8, max_order=3, ngens=12);
 41.129208 seconds (411.24 M allocations: 45.861 GiB, 19.64% gc time)

julia> analyze(sol, ε=0.02, ngens=12)
final diameter of current I = 574.0077157938906
final diameter of position x = 0.022593942314352897
tspan(vsol[first_in]) = [0.0580367, 0.0580368]
max values of x
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.05847325209977274
ρ([1.0, 0.0], Flowpipe(vsol.Xk[first_in:end])) = 0.060613313128198155
ρ([1.0, 0.0], Flowpipe(vsol.Xk[end:end])) = 0.060263558263833554
min values of x
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:first_in]))) = 0.030000001614109116
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[first_in:end]))) = 0.030000001614109116
-(ρ([-1.0, 0.0], Flowpipe(vsol.Xk[end:end]))) = 0.03766961594948066
max values of v
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:first_in])) = 0.008901796674425007
ρ([0.0, 1.0], Flowpipe(vsol.Xk[first_in:end])) = 0.008901796674425007
ρ([0.0, 1.0], Flowpipe(vsol.Xk[end:end])) = 0.005519600669287567
min values of v
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:first_in]))) = -0.004423141273818343
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[first_in:end]))) = -0.00493014619883479
-(ρ([0.0, -1.0], Flowpipe(vsol.Xk[end:end]))) = -0.004629346670829155

pv2, with jitter, order 1

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

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

julia> sol = run_pv(prob, 1000, 1e-8, max_order=1, ngens=4);
  9.154629 seconds (110.67 M allocations: 10.250 GiB, 14.45% gc time)

NO VERIFICA
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment