From 8e6d5dd36585dba085db497493175dbae3f42f5a Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 1 Mar 2024 08:12:45 +0100 Subject: [PATCH 1/3] better script selection (#jl) in examples --- examples/Brusselator/Brusselator.jl | 16 +++---- examples/Building/Building.jl | 35 +++++++------- .../DuffingOscillator/DuffingOscillator.jl | 8 ++-- examples/EMBrake/EMBrake.jl | 2 +- examples/Heat3D/Heat3D.jl | 2 +- examples/ISS/ISS.jl | 30 ++++++------ examples/LaubLoomis/LaubLoomis.jl | 37 +++++++-------- examples/Lorenz/Lorenz.jl | 20 ++++---- examples/LotkaVolterra/LotkaVolterra.jl | 16 +++---- examples/OpAmp/OpAmp.jl | 44 +++++++++--------- examples/Platoon/Platoon.jl | 26 +++++------ .../ProductionDestruction.jl | 14 +++--- examples/Quadrotor/Quadrotor.jl | 17 +++---- examples/SEIR/SEIR.jl | 12 ++--- examples/Spacecraft/Spacecraft.jl | 46 +++++++++---------- .../SquareWaveOscillator.jl | 22 ++++----- examples/TransmissionLine/TransmissionLine.jl | 25 +++++----- examples/VanDerPol/VanDerPol.jl | 8 ++-- 18 files changed, 185 insertions(+), 195 deletions(-) diff --git a/examples/Brusselator/Brusselator.jl b/examples/Brusselator/Brusselator.jl index 8700ec371..50aceea06 100644 --- a/examples/Brusselator/Brusselator.jl +++ b/examples/Brusselator/Brusselator.jl @@ -25,7 +25,7 @@ # The numerical values for the model's constants (in their respective units) are # ``A = 1, B = 1.5``. -using ReachabilityAnalysis #!jl +using ReachabilityAnalysis @taylorize function brusselator!(du, u, p, t) local A = 1.0 @@ -49,10 +49,10 @@ end # The initial set is ``x ∈ [0.8, 1]``, ``y ∈ [0, 0.2]``. The time horizon is # ``18``. These settings are taken from [^CAS13]. -U₀ = (0.8 .. 1.0) × (0.0 .. 0.2) #!jl -prob = @ivp(u' = brusselator!(u), u(0) ∈ U₀, dim:2) #!jl +U₀ = (0.8 .. 1.0) × (0.0 .. 0.2) +prob = @ivp(u' = brusselator!(u), u(0) ∈ U₀, dim:2) -T = 18.0; #!jl +T = 18.0; # ## Analysis @@ -60,8 +60,8 @@ T = 18.0; #!jl # second-order expansion in the spatial variables. The version `TMJets20` is # more efficient here. -alg = TMJets20(; orderT=6, orderQ=2) #!jl -sol = solve(prob; T=T, alg=alg); #!jl +alg = TMJets20(; orderT=6, orderQ=2) +sol = solve(prob; T=T, alg=alg); # ## Results @@ -97,7 +97,7 @@ bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim:2); # First we solve for ``r = 0.01``: -sol_01 = solve(bruss(0.01); T=30.0, alg=alg); #!jl +sol_01 = solve(bruss(0.01); T=30.0, alg=alg); #- @@ -115,7 +115,7 @@ plot!(fig, U0(0.01); color=:orange, lab="Uo", xlims=(0.6, 1.3)) #!jl # flowpipe zoomed to the last portion and compare ``r = 0.01`` with a larger set # of initial states, ``r = 0.1``. -sol_1 = solve(bruss(0.1); T=30.0, alg=alg); #!jl +sol_1 = solve(bruss(0.1); T=30.0, alg=alg); #- diff --git a/examples/Building/Building.jl b/examples/Building/Building.jl index 03076239d..8ecc6b0f1 100644 --- a/examples/Building/Building.jl +++ b/examples/Building/Building.jl @@ -45,6 +45,7 @@ using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath using ReachabilityBase.Arrays: SingleEntryVector +using ReachabilityAnalysis: add_dimension const x25 = SingleEntryVector(25, 48, 1.0) const x25e = SingleEntryVector(25, 49, 1.0) @@ -71,8 +72,6 @@ function building_BLDF01() return prob_BLDF01 end -using ReachabilityAnalysis: add_dimension - function building_BLDC01() A, B, X0, U = building_common() n = size(A, 1) @@ -127,8 +126,8 @@ prob_BLDF01 = building_BLDF01(); sol_BLDF01_dense = solve(prob_BLDF01; T=20.0, alg=LGG09(; δ=0.004, vars=(25), n=48)) -fig = plot(sol_BLDF01_dense; vars=(0, 25), linecolor=:blue, color=:blue, - alpha=0.8, lw=1.0, xlab="t", ylab="x25") +fig = plot(sol_BLDF01_dense; vars=(0, 25), linecolor=:blue, color=:blue, #!jl + alpha=0.8, lw=1.0, xlab="t", ylab="x25") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -136,20 +135,20 @@ fig = plot(sol_BLDF01_dense; vars=(0, 25), linecolor=:blue, color=:blue, @assert ρ(x25, sol_BLDF01_dense) <= 5.1e-3 "the property should be proven" # BLDF01 - BDS01 @assert !(ρ(x25, sol_BLDF01_dense) <= 4e-3) "the property should not be proven" # BLDF01 - BDU01 -ρ(x25, sol_BLDF01_dense) +ρ(x25, sol_BLDF01_dense) #!jl #- @assert !(ρ(x25, sol_BLDF01_dense(20.0)) <= -0.78e-3) "the property should not be proven" # BLDF01 - BDU02 -ρ(x25, sol_BLDF01_dense(20.0)) +ρ(x25, sol_BLDF01_dense(20.0)) #!jl # #### Discrete time sol_BLDF01_discrete = solve(prob_BLDF01; T=20.0, alg=LGG09(; δ=0.01, vars=(25), n=48, approx_model=NoBloating())); -fig = plot(sol_BLDF01_discrete; vars=(0, 25), linecolor=:blue, color=:blue, - alpha=0.8, lw=1.0, xlab="t", ylab="x25") +fig = plot(sol_BLDF01_discrete; vars=(0, 25), linecolor=:blue, color=:blue, #!jl + alpha=0.8, lw=1.0, xlab="t", ylab="x25") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -157,12 +156,12 @@ fig = plot(sol_BLDF01_discrete; vars=(0, 25), linecolor=:blue, color=:blue, @assert ρ(x25, sol_BLDF01_discrete) <= 5.1e-3 "the property should be proven" # BLDF01 - BDS01 @assert !(ρ(x25, sol_BLDF01_discrete) <= 4e-3) "the property should not be proven" # BLDF01 - BDU01 -ρ(x25, sol_BLDF01_discrete) +ρ(x25, sol_BLDF01_discrete) #!jl #- @assert !(ρ(x25, sol_BLDF01_discrete(20.0)) <= -0.78e-3) "the property should not be proven" # BLDF01 - BDU02 -ρ(x25, sol_BLDF01_discrete(20.0)) +ρ(x25, sol_BLDF01_discrete(20.0)) #!jl # ### BLDC01 @@ -172,8 +171,8 @@ prob_BLDC01 = building_BLDC01(); sol_BLDC01_dense = solve(prob_BLDC01; T=20.0, alg=LGG09(; δ=0.005, vars=(25), n=49)) -fig = plot(sol_BLDC01_dense; vars=(0, 25), linecolor=:blue, color=:blue, - alpha=0.8, lw=1.0, xlab="t", ylab="x25") +fig = plot(sol_BLDC01_dense; vars=(0, 25), linecolor=:blue, color=:blue, #!jl + alpha=0.8, lw=1.0, xlab="t", ylab="x25") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -181,20 +180,20 @@ fig = plot(sol_BLDC01_dense; vars=(0, 25), linecolor=:blue, color=:blue, @assert ρ(x25e, sol_BLDC01_dense) <= 5.1e-3 "the property should be proven" # BLDC01 - BDS01 @assert !(ρ(x25e, sol_BLDC01_dense) <= 4e-3) "the property should not be proven" # BLDC01 - BDU01 -ρ(x25e, sol_BLDC01_dense) +ρ(x25e, sol_BLDC01_dense) #!jl #- @assert !(ρ(x25e, sol_BLDC01_dense(20.0)) <= -0.78e-3) "the property should not be proven" # BLDC01 - BDU02 -ρ(x25, sol_BLDF01_discrete(20.0)) +ρ(x25, sol_BLDF01_discrete(20.0)) #!jl # #### Discrete time sol_BLDC01_discrete = solve(prob_BLDC01; T=20.0, alg=LGG09(; δ=0.01, vars=(25), n=49, approx_model=NoBloating())) -fig = plot(sol_BLDC01_discrete; vars=(0, 25), linecolor=:blue, color=:blue, - alpha=0.8, lw=1.0, xlab="t", ylab="x25") +fig = plot(sol_BLDC01_discrete; vars=(0, 25), linecolor=:blue, color=:blue, #!jl + alpha=0.8, lw=1.0, xlab="t", ylab="x25") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -202,12 +201,12 @@ fig = plot(sol_BLDC01_discrete; vars=(0, 25), linecolor=:blue, color=:blue, @assert ρ(x25e, sol_BLDC01_discrete) <= 5.1e-3 "the property should be proven" # BLDC01 - BDS01 @assert !(ρ(x25e, sol_BLDC01_discrete) <= 4e-3) "the property should not be proven" # BLDC01 - BDU01 -ρ(x25e, sol_BLDC01_discrete) +ρ(x25e, sol_BLDC01_discrete) #!jl #- @assert !(ρ(x25e, sol_BLDC01_discrete(20.0)) <= -0.78e-3) "the property should not be proven" # BLDC01 - BDU02 -ρ(x25e, sol_BLDC01_discrete(20.0)) +ρ(x25e, sol_BLDC01_discrete(20.0)) #!jl # ## References diff --git a/examples/DuffingOscillator/DuffingOscillator.jl b/examples/DuffingOscillator/DuffingOscillator.jl index 28458a448..f39cb8efc 100644 --- a/examples/DuffingOscillator/DuffingOscillator.jl +++ b/examples/DuffingOscillator/DuffingOscillator.jl @@ -7,7 +7,7 @@ # ## Model description -using ReachabilityAnalysis #!jl +using ReachabilityAnalysis const ω = 1.2 @@ -42,18 +42,18 @@ sol = solve(prob; tspan=(0.0, 20 * T), alg=TMJets21a()); using Plots #!jl #!jl import DisplayAs #hide -fig = plot(sol; vars=(0, 1), xlab="t", ylab="x", lw=0.2, color=:blue) +fig = plot(sol; vars=(0, 1), xlab="t", ylab="x", lw=0.2, color=:blue) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide #- -fig = plot(sol; vars=(0, 2), xlab="t", ylab="v", lw=0.2, color=:blue) +fig = plot(sol; vars=(0, 2), xlab="t", ylab="v", lw=0.2, color=:blue) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide #- -fig = plot(sol; vars=(1, 2), xlab="x", ylab="v", lw=0.5, color=:red) +fig = plot(sol; vars=(1, 2), xlab="x", ylab="v", lw=0.5, color=:red) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/EMBrake/EMBrake.jl b/examples/EMBrake/EMBrake.jl index 4b4cb17f2..6ed447470 100644 --- a/examples/EMBrake/EMBrake.jl +++ b/examples/EMBrake/EMBrake.jl @@ -33,7 +33,7 @@ # - **Constant inputs**: The inputs are only uncertain in the initial value, and # constant over time: ``u(0) ∈ \mathcal{U}``, ``\dot{u}(t) = 0.`` -using ReachabilityAnalysis, SparseArrays #!jl +using ReachabilityAnalysis, SparseArrays # ## Common functionality between model variants diff --git a/examples/Heat3D/Heat3D.jl b/examples/Heat3D/Heat3D.jl index 0265942c6..8094ce663 100644 --- a/examples/Heat3D/Heat3D.jl +++ b/examples/Heat3D/Heat3D.jl @@ -24,6 +24,6 @@ # - **Constant inputs**: The inputs are only uncertain in the initial value, and # constant over time: ``u(0) ∈ \mathcal{U}``, ``\dot{u}(t) = 0.`` -using ReachabilityAnalysis, SparseArrays, JLD2, ReachabilityBase.CurrentPath #!jl +using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath path = @current_path("Heat3D", "HEAT01.jld2") diff --git a/examples/ISS/ISS.jl b/examples/ISS/ISS.jl index 7dd6f9095..5edd556c9 100644 --- a/examples/ISS/ISS.jl +++ b/examples/ISS/ISS.jl @@ -35,7 +35,7 @@ # - **ISSC01** (constant inputs): The inputs are only uncertain in their initial # value, and constant over time: ``u(0) ∈ \mathcal{U}``, ``\dot{u}(t) = 0``. -using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath #!jl +using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath using ReachabilityAnalysis: add_dimension path = @current_path("ISS", "ISS.jld2") @@ -121,17 +121,17 @@ sol_ISSF01 = solve(prob_ISSF01; T=20.0, # The solution `sol_ISSF01` is a 270-dimensional set that (only) contains # template reach sets for the linear combination `C_3 x(t)`. -dim(sol_ISSF01) +dim(sol_ISSF01) #!jl # For visualization, it is necessary to specify that we want to plot "time" vs. # ``y_3(t)``. We can transform the flowpipe on the output ``y_3(t)`` by # "flattening" the flowpipe along directions `dirs`. -πsol_ISSF01 = flatten(sol_ISSF01); +πsol_ISSF01 = flatten(sol_ISSF01); #!jl # Now `πsol_ISSF01` is a one-dimensional flowpipe. -dim(πsol_ISSF01) +dim(πsol_ISSF01) #!jl #md # !!! tip "Performance tip" #md # Note that projecting the solution along direction ``C_3`` corresponds @@ -149,11 +149,11 @@ LazySets.set_ztol(Float64, 1e-8); # use higher precision for the plots #!jl #- -fig = Plots.plot(πsol_ISSF01[1:10:end]; vars=(0, 1), linecolor=:blue, color=:blue, - alpha=0.8, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0], - ytick=[-0.00075, -0.0005, -0.00025, 0, 0.00025, 0.0005, 0.00075], - xlims=(0.0, 20.0), ylims=(-0.00075, 0.00075), - bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) +fig = plot(πsol_ISSF01[1:10:end]; vars=(0, 1), linecolor=:blue, color=:blue, #!jl + alpha=0.8, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0], #!jl + ytick=[-0.00075, -0.0005, -0.00025, 0, 0.00025, 0.0005, 0.00075], #!jl + xlims=(0.0, 20.0), ylims=(-0.00075, 0.00075), #!jl + bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide # ### ISSC01 @@ -165,15 +165,15 @@ sol_ISSC01 = solve(prob_ISSC01; T=20.0, # We can flatten the flowpipe to the output ``y_3(t)`` as before: -πsol_ISSC01 = flatten(sol_ISSC01); +πsol_ISSC01 = flatten(sol_ISSC01); #!jl #- -fig = Plots.plot(πsol_ISSC01; vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8, - lw=1.0, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0], - ytick=[-0.0002, -0.0001, 0.0, 0.0001, 0.0002], - xlims=(0.0, 20.0), ylims=(-0.0002, 0.0002), - bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) +fig = plot(πsol_ISSC01; vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8, #!jl + lw=1.0, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0], #!jl + ytick=[-0.0002, -0.0001, 0.0, 0.0001, 0.0002], #!jl + xlims=(0.0, 20.0), ylims=(-0.0002, 0.0002), #!jl + bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide #- diff --git a/examples/LaubLoomis/LaubLoomis.jl b/examples/LaubLoomis/LaubLoomis.jl index b5d933257..872f148f0 100644 --- a/examples/LaubLoomis/LaubLoomis.jl +++ b/examples/LaubLoomis/LaubLoomis.jl @@ -11,7 +11,7 @@ # enzymatic activities. The dynamics can be defined by the following # seven-dimensional system of differential equations. -using ReachabilityAnalysis #!jl +using ReachabilityAnalysis @taylorize function laubloomis!(dx, x, p, t) dx[1] = 1.4 * x[3] - 0.9 * x[1] @@ -68,11 +68,11 @@ sol_1z = overapproximate(sol_1, Zonotope); # We verify that the specification holds: @assert ρ(e4, sol_1z) < 4.5 "the property should be proven" -ρ(e4, sol_1z) +ρ(e4, sol_1z) #!jl # To compute the width of the final box, we use the support function: -ρ(e4, sol_1z[end]) + ρ(-e4, sol_1z[end]) +ρ(e4, sol_1z[end]) + ρ(-e4, sol_1z[end]) #!jl # ### Case 2: Intermediate initial set @@ -87,11 +87,11 @@ sol_2z = overapproximate(sol_2, Zonotope); # We verify that the specification holds: @assert ρ(e4, sol_2z) < 4.5 "the property should be proven" -ρ(e4, sol_2z) +ρ(e4, sol_2z) #!jl # Width of final box: -ρ(e4, sol_2z[end]) + ρ(-e4, sol_2z[end]) +ρ(e4, sol_2z[end]) + ρ(-e4, sol_2z[end]) #!jl # ### Case 3: Large initial set @@ -106,30 +106,27 @@ sol_3z = overapproximate(sol_3, Zonotope); # We verify that the specification holds: @assert ρ(e4, sol_3z) < 5.0 "the property should be proven" -ρ(e4, sol_3z) +ρ(e4, sol_3z) #!jl # Width of final box: -ρ(e4, sol_3z[end]) + ρ(-e4, sol_3z[end]) +ρ(e4, sol_3z[end]) + ρ(-e4, sol_3z[end]) #!jl # ## Results using Plots, Plots.PlotMeasures, LaTeXStrings #!jl #!jl import DisplayAs #hide -fig = plot(sol_3z; vars=(0, 4), linecolor="green", color=:green, alpha=0.8, lab="W = 0.1") -plot!(fig, sol_2z; vars=(0, 4), linecolor="blue", color=:blue, alpha=0.8, lab="W = 0.05") -plot!(fig, sol_1z; vars=(0, 4), linecolor="yellow", color=:yellow, alpha=0.8, lab="W = 0.01", - tickfont=font(10, "Times"), guidefontsize=15, - xlab=L"t", - ylab=L"x_4", - xtick=[0.0, 5.0, 10.0, 15.0, 20.0], ytick=[1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0], - xlims=(0.0, 20.0), ylims=(1.5, 5.02), - bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm, - size=(600, 600)) - -plot!(fig, x -> x, x -> 4.5, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="") -plot!(fig, x -> x, x -> 5.0, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="") +fig = plot(sol_3z; vars=(0, 4), linecolor="green", color=:green, alpha=0.8, lab="W = 0.1") #!jl +plot!(fig, sol_2z; vars=(0, 4), linecolor="blue", color=:blue, alpha=0.8, lab="W = 0.05") #!jl +plot!(fig, sol_1z; vars=(0, 4), linecolor="yellow", color=:yellow, alpha=0.8, lab="W = 0.01", #!jl + tickfont=font(10, "Times"), guidefontsize=15, xlab=L"t", ylab=L"x_4", #!jl + xtick=[0.0, 5.0, 10.0, 15.0, 20.0], ytick=[1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0], #!jl + xlims=(0.0, 20.0), ylims=(1.5, 5.02), #!jl + bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm, #!jl + size=(600, 600)) #!jl +plot!(fig, x -> x, x -> 4.5, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="") #!jl +plot!(fig, x -> x, x -> 5.0, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/Lorenz/Lorenz.jl b/examples/Lorenz/Lorenz.jl index bd110fc98..f46b1bdc6 100644 --- a/examples/Lorenz/Lorenz.jl +++ b/examples/Lorenz/Lorenz.jl @@ -26,7 +26,7 @@ # ``σ``, ``ρ``, and ``β`` are proportional to the Prandtl number, Rayleigh # number, and certain physical dimensions of the layer itself. -using ReachabilityAnalysis #!jl +using ReachabilityAnalysis @taylorize function lorenz!(du, u, p, t) local σ = 10.0 @@ -54,9 +54,7 @@ prob = @ivp(x' = lorenz!(x), dim:3, x(0) ∈ X0); # ``n_Q=2``. alg = TMJets(; abstol=1e-15, orderT=10, orderQ=2, maxsteps=50_000) - sol = solve(prob; T=10.0, alg=alg) - solz = overapproximate(sol, Zonotope); # ## Results @@ -64,15 +62,15 @@ solz = overapproximate(sol, Zonotope); using Plots #!jl #!jl import DisplayAs #hide -fig = plot(solz; vars=(0, 1), xlab="t", ylab="x") +fig = plot(solz; vars=(0, 1), xlab="t", ylab="x") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide # It is apparent by inspection that variable ``x(t)`` does not exceed ``20`` in # the computed time span: -fig = plot(solz(0.0 .. 1.5); vars=(0, 1), xlab="t", ylab="x", lw=0) -plot!(fig, x -> 20.0; c=:red, xlims=(0.0, 1.5), lab="") +fig = plot(solz(0.0 .. 1.5); vars=(0, 1), xlab="t", ylab="x", lw=0) #!jl +plot!(fig, x -> 20.0; c=:red, xlims=(0.0, 1.5), lab="") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -83,11 +81,11 @@ plot!(fig, x -> 20.0; c=:red, xlims=(0.0, 1.5), lab="") #- -ρ([1.0, 0, 0], solz(0 .. 1.5)) +ρ([1.0, 0, 0], solz(0 .. 1.5)) #!jl # In a similar fashion, we can compute extremal values of variable ``y(t)``: -fig = plot(solz; vars=(0, 2), xlab="t", ylab="y") +fig = plot(solz; vars=(0, 2), xlab="t", ylab="y") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -95,11 +93,11 @@ fig = plot(solz; vars=(0, 2), xlab="t", ylab="y") # quantities are a lower bound on the exact minimum (resp. an upper bound on the # exact maximum): --ρ([0.0, -1.0, 0.0], solz), ρ([0.0, 1.0, 0.0], solz) +-ρ([0.0, -1.0, 0.0], solz), ρ([0.0, 1.0, 0.0], solz) #!jl #- -fig = plot(solz; vars=(0, 3), xlab="t", ylab="z") +fig = plot(solz; vars=(0, 3), xlab="t", ylab="z") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -107,6 +105,6 @@ fig = plot(solz; vars=(0, 3), xlab="t", ylab="z") # Below we plot the flowpipe projected on the ``x``/``z``-plane. -fig = plot(solz; vars=(1, 3), xlab="x", ylab="z") +fig = plot(solz; vars=(1, 3), xlab="x", ylab="z") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/LotkaVolterra/LotkaVolterra.jl b/examples/LotkaVolterra/LotkaVolterra.jl index 611e34d05..d1b2ba4e1 100644 --- a/examples/LotkaVolterra/LotkaVolterra.jl +++ b/examples/LotkaVolterra/LotkaVolterra.jl @@ -27,7 +27,7 @@ # # We set these parameters to ``α = 1.5``, ``β = 1``, ``γ = 3``, and ``δ = 1``. -using ReachabilityAnalysis #!jl +using ReachabilityAnalysis @taylorize function lotkavolterra!(du, u, p, t) local α, β, γ, δ = 1.5, 1.0, 3.0, 1.0 @@ -56,9 +56,9 @@ solz = overapproximate(sol, Zonotope); using Plots #!jl #!jl import DisplayAs #hide -fig = plot(solz; vars=(1, 2), alpha=0.3, lw=0.0, xlab="x", ylab="y", - lab="Flowpipe", legend=:bottomright) -plot!(fig, X0; label="X(0)") +fig = plot(solz; vars=(1, 2), alpha=0.3, lw=0.0, xlab="x", ylab="y", #!jl + lab="Flowpipe", legend=:bottomright) #!jl +plot!(fig, X0; label="X(0)") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -95,8 +95,8 @@ solz = overapproximate(sol, Zonotope); # ## Results -fig = plot(solz; vars=(1, 2), lw=0.3, title="Uncertain parameters", - lab="abstol = 1e-15", xlab="x", ylab="y") +fig = plot(solz; vars=(1, 2), lw=0.3, title="Uncertain parameters", #!jl + lab="abstol = 1e-15", xlab="x", ylab="y") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -123,8 +123,8 @@ solz = overapproximate(sol, Zonotope); # ### Results -fig = plot(solz; vars=(1, 2), color=:orange, lw=0.3, lab="ϵ = 0.05", - title="Uncertain u0 and uncertain parameters", xlab="x", ylab="y") +fig = plot(solz; vars=(1, 2), color=:orange, lw=0.3, lab="ϵ = 0.05", #!jl + title="Uncertain u0 and uncertain parameters", xlab="x", ylab="y") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/OpAmp/OpAmp.jl b/examples/OpAmp/OpAmp.jl index dc64a695b..22f574f7a 100644 --- a/examples/OpAmp/OpAmp.jl +++ b/examples/OpAmp/OpAmp.jl @@ -73,7 +73,7 @@ # It is convenient to define constants ``α`` and ``β`` such that # ``\dfrac{d e_{out}(t)}{dt} = α e_{out}(t) + β e_{in}(t)``. -using ReachabilityAnalysis #!jl +using ReachabilityAnalysis function opamp_nondet(; X0=Singleton([0.0]), R₁=2.0, R₂=6.0, C=1.e-3, @@ -101,8 +101,8 @@ sol_nondet = solve(opamp_nondet(); T=0.1, alg=INT(; δ=1e-4)); using Plots, LaTeXStrings #!jl #!jl import DisplayAs #hide -fig = plot(sol_nondet; vars=(0, 1), xlab=L"t", ylab=L"e_{out}", - title="Solution for nondeterministic input", lw=0.2) +fig = plot(sol_nondet; vars=(0, 1), xlab=L"t", ylab=L"e_{out}", #!jl + title="Solution for nondeterministic input", lw=0.2) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -110,18 +110,18 @@ fig = plot(sol_nondet; vars=(0, 1), xlab=L"t", ylab=L"e_{out}", # we solve for three different initial conditions of increasing width, for a # shorter time horizon. -Δt = 0 .. 0.04 +Δt = 0 .. 0.04 #!jl -fig = plot(; xlab=L"t", ylab=L"e_{out}(t)", title="Solution for nondeterministic input") +fig = plot(; xlab=L"t", ylab=L"e_{out}(t)", title="Solution for nondeterministic input") #!jl sol_nondet = solve(opamp_nondet(; X0=Interval(-1.0, 1.0)); T=0.1, alg=INT(; δ=1e-4)) -plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = -1 .. 1", lw=0.2) +plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = -1 .. 1", lw=0.2) #!jl sol_nondet = solve(opamp_nondet(; X0=Interval(-0.5, 0.5)); T=0.1, alg=INT(; δ=1e-4)) -plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = -0.5 .. 0.5", lw=0.2) +plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = -0.5 .. 0.5", lw=0.2) #!jl sol_nondet = solve(opamp_nondet(); T=0.1, alg=INT(; δ=1e-4)) -plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = 0", lw=0.2) +plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = 0", lw=0.2) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -160,7 +160,7 @@ plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = 0", lw=0.2) # ``e_{in}(0) ∈ E_{in, 0}``, ``e_{in}(t)`` constant and ``γ = δ = 0``. We will # consider the other cases in the next section. -using Symbolics #!jl +using Symbolics function opamp_with_saturation(; X0=Singleton(zeros(2)), R₁=2.0, R₂=6.0, C=1.e-3, @@ -210,16 +210,16 @@ sol_const = solve(prob_const; T=0.1, alg=BOX(; δ=1e-4)); # ### Results -fig = plot(sol_const; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", - title="Constant input signal", lw=0.2) -plot!(fig, x -> x, x -> 1.5, 0.0, 0.1; color="red", lw=2, ls=:dash, lab="") +fig = plot(sol_const; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", #!jl + title="Constant input signal", lw=0.2) #!jl +plot!(fig, x -> x, x -> 1.5, 0.0, 0.1; color="red", lw=2, ls=:dash, lab="") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide #- -fig = plot(sol_const; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", - title="Solution for constant input", lw=0.2) +fig = plot(sol_const; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", #!jl + title="Solution for constant input", lw=0.2) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -251,29 +251,29 @@ sol_exp = solve(prob_exp; T=0.1, alg=BOX(; δ=1e-4)); # ## Results -fig = plot(sol_lin; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", - title="Linear input with saturation", lw=0.2) +fig = plot(sol_lin; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", #!jl + title="Linear input with saturation", lw=0.2) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide #- -fig = plot(sol_exp; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", - title="Exponential input", lw=0.2) +fig = plot(sol_exp; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", #!jl + title="Exponential input", lw=0.2) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide #- -fig = plot(sol_lin; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", - title="Solution for linear input", lw=0.2) +fig = plot(sol_lin; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", #!jl + title="Solution for linear input", lw=0.2) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide #- -fig = plot(sol_exp; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", - title="Solution for exponential input", lw=0.2) +fig = plot(sol_exp; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", #!jl + title="Solution for exponential input", lw=0.2) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/Platoon/Platoon.jl b/examples/Platoon/Platoon.jl index 32cf82db7..3e10d0281 100644 --- a/examples/Platoon/Platoon.jl +++ b/examples/Platoon/Platoon.jl @@ -34,7 +34,7 @@ # communication in zero time. We will consider PLAN01: ``t_b = 10``, # ``t_c = 20``, and ``t_r = 20``. -using ReachabilityAnalysis, SparseArrays, Symbolics #!jl +using ReachabilityAnalysis, SparseArrays, Symbolics const var = @variables x[1:9] t; @@ -195,7 +195,7 @@ prob_PLAD01 = platoon(); # coordinate: boxdirs = BoxDirections(10) -length(boxdirs) +length(boxdirs) #!jl #- @@ -217,27 +217,27 @@ sol_PLAD01_BND42 = solve(prob_PLAD01; # Minimum of ``x_1(t)``: --ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND42) +-ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND42) #!jl #- # Minimum of ``x_4(t)``: --ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND42) +-ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND42) #!jl #- # Minimum of ``x_7(t)``: --ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND42) +-ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND42) #!jl # Next, we plot variable ``x_1`` over time. using Plots, LaTeXStrings #!jl #!jl import DisplayAs #hide -fig = plot(sol_PLAD01_BND42; vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND42", lw=0.1) -plot!(x -> x, x -> -42.0, 0.0, 20.0; linewidth=2, color="red", ls=:dash, leg=nothing) +fig = plot(sol_PLAD01_BND42; vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND42", lw=0.1) #!jl +plot!(x -> x, x -> -42.0, 0.0, 20.0; linewidth=2, color="red", ls=:dash, leg=nothing) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -252,7 +252,7 @@ plot!(x -> x, x -> -42.0, 0.0, 20.0; linewidth=2, color="red", ls=:dash, leg=not # are 200 such directions: octdirs = OctDirections(10) -length(octdirs) +length(octdirs) #!jl #md # !!! tip "Performance tip" #md # The increase in the number of directions implies an increase in @@ -281,23 +281,23 @@ sol_PLAD01_BND30 = solve(prob_PLAD01; # Minimum of ``x_1(t)``: --ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND30) +-ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND30) #!jl #- # Minimum of ``x_4(t)``: --ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND30) +-ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND30) #!jl #- # Minimum of ``x_7(t)``: --ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND30) +-ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND30) #!jl # Finally, we plot variable ``x_1`` over time again. -fig = plot(sol_PLAD01_BND30; vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND30", lw=0.1) -plot!(x -> x, x -> -30.0, 0.0, 20.0; linewidth=2, color="red", ls=:dash, leg=nothing) +fig = plot(sol_PLAD01_BND30; vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND30", lw=0.1) #!jl +plot!(x -> x, x -> -30.0, 0.0, 20.0; linewidth=2, color="red", ls=:dash, leg=nothing) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/ProductionDestruction/ProductionDestruction.jl b/examples/ProductionDestruction/ProductionDestruction.jl index e51b23ce8..3054c9229 100644 --- a/examples/ProductionDestruction/ProductionDestruction.jl +++ b/examples/ProductionDestruction/ProductionDestruction.jl @@ -61,7 +61,7 @@ # beginning by defining the positive orthant in three-dimensional space, # `positive_orthant`, as an unbounded polyhedron in constraint representation. -using ReachabilityAnalysis, Symbolics #!jl +using ReachabilityAnalysis, Symbolics @variables x y z positive_orthant = HPolyhedron([x >= 0, y >= 0, z >= 0], [x, y, z]); @@ -141,7 +141,7 @@ sol = solve(prob; T=100.0, alg=TMJets(; abstol=1e-12, orderT=6, orderQ=1)); property, vol = prod_dest_verif(sol) @assert property "the property should be proven" -vol +vol #!jl # #### Results @@ -150,7 +150,7 @@ vol using Plots #!jl #!jl import DisplayAs #hide -fig = plot(sol; vars=(0, 3), lc=:orange, c=:orange, alpha=0.3, lab="I", xlab="t", ylab="z") +fig = plot(sol; vars=(0, 3), lc=:orange, c=:orange, alpha=0.3, lab="I", xlab="t", ylab="z") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -183,13 +183,13 @@ sol = solve(prob; T=100.0, alg=TMJets(; abstol=9e-13, orderT=6, orderQ=1)); property, vol = prod_dest_verif(sol) @assert property "the property should be proven" -vol +vol #!jl # #### Results # We plot ``z`` over time. -fig = plot(sol; vars=(0, 3), lc=:blue, c=:blue, alpha=0.3, lab="P", xlab="t", ylab="z") +fig = plot(sol; vars=(0, 3), lc=:blue, c=:blue, alpha=0.3, lab="P", xlab="t", ylab="z") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -211,13 +211,13 @@ sol = solve(prob; T=100.0, alg=TMJets(; abstol=1e-12, orderT=6, orderQ=1)); property, vol = prod_dest_verif(sol) @assert property "the property should be proven" -vol +vol #!jl # #### Results # We plot ``z`` over time. -fig = plot(sol; vars=(0, 3), lc=:red, c=:red, alpha=0.3, lab="I & P", xlab="t", ylab="z") +fig = plot(sol; vars=(0, 3), lc=:red, c=:red, alpha=0.3, lab="I & P", xlab="t", ylab="z") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/Quadrotor/Quadrotor.jl b/examples/Quadrotor/Quadrotor.jl index 20a7f88f2..5eba8c522 100644 --- a/examples/Quadrotor/Quadrotor.jl +++ b/examples/Quadrotor/Quadrotor.jl @@ -67,6 +67,7 @@ # We leave the heading uncontrolled by setting ``τ_ψ = 0``. using ReachabilityAnalysis +using ReachabilityBase.Arrays: SingleEntryVector ## parameters of the model const g = 9.81 # gravity constant in m/s^2 @@ -156,8 +157,6 @@ end; # settings the same. No goals are specified for these cases: the objective # instead is to understand the scalability of the tool. -using ReachabilityBase.Arrays: SingleEntryVector - const T = 5.0 const v3 = SingleEntryVector(3, 12, 1.0) @@ -192,8 +191,6 @@ end; # ## Analysis -cases = ["Δ=0.1", "Δ=0.4", "Δ=0.8"]; - # ### Case 1: smaller uncertainty Wpos = 0.1 @@ -241,12 +238,12 @@ solz3 = overapproximate(sol, Zonotope); using Plots #!jl #!jl import DisplayAs #hide -fig = plot(solz3; vars=(0, 3), linecolor="green", color=:green, alpha=0.8) -plot!(solz2; vars=(0, 3), linecolor="blue", color=:blue, alpha=0.8) -plot!(solz1; vars=(0, 3), linecolor="yellow", color=:yellow, alpha=0.8, - xlab="t", ylab="x3", - xtick=[0.0, 1.0, 2.0, 3.0, 4.0, 5.0], ytick=[-1.0, -0.5, 0.0, 0.5, 1.0, 1.5], - xlims=(0.0, 5.0), ylims=(-1.0, 1.5)) +fig = plot(solz3; vars=(0, 3), linecolor="green", color=:green, alpha=0.8) #!jl +plot!(solz2; vars=(0, 3), linecolor="blue", color=:blue, alpha=0.8) #!jl +plot!(solz1; vars=(0, 3), linecolor="yellow", color=:yellow, alpha=0.8, #!jl + xlab="t", ylab="x3", #!jl + xtick=[0.0, 1.0, 2.0, 3.0, 4.0, 5.0], ytick=[-1.0, -0.5, 0.0, 0.5, 1.0, 1.5], #!jl + xlims=(0.0, 5.0), ylims=(-1.0, 1.5)) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/SEIR/SEIR.jl b/examples/SEIR/SEIR.jl index 8e6c4e457..cb44c429e 100644 --- a/examples/SEIR/SEIR.jl +++ b/examples/SEIR/SEIR.jl @@ -25,7 +25,7 @@ # exposed population, ``I`` is the stock of infected, ``R`` is the stock of # removed population (either by death or recovery), with ``S + E + I + R = N``. -using ReachabilityAnalysis #!jl +using ReachabilityAnalysis @taylorize function seir!(dx, x, p, t) S, E, I, R, α, β, γ = x @@ -71,10 +71,10 @@ solz = overapproximate(sol, Zonotope); using Plots #!jl #!jl import DisplayAs #hide -fig = plot(; legend=:outerright) -plot!(fig, solz; vars=(0, 1), color=:blue, lw=0.0, lab="S") -plot!(fig, solz; vars=(0, 2), color=:green, lw=0.0, lab="E") -plot!(fig, solz; vars=(0, 3), color=:red, lw=0.0, lab="I") -plot!(fig, solz; vars=(0, 4), color=:grey, lw=0.0, lab="R") +fig = plot(; legend=:outerright) #!jl +plot!(fig, solz; vars=(0, 1), color=:blue, lw=0.0, lab="S") #!jl +plot!(fig, solz; vars=(0, 2), color=:green, lw=0.0, lab="E") #!jl +plot!(fig, solz; vars=(0, 3), color=:red, lw=0.0, lab="I") #!jl +plot!(fig, solz; vars=(0, 4), color=:grey, lw=0.0, lab="R") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/Spacecraft/Spacecraft.jl b/examples/Spacecraft/Spacecraft.jl index bdfb648e7..b67c21939 100644 --- a/examples/Spacecraft/Spacecraft.jl +++ b/examples/Spacecraft/Spacecraft.jl @@ -187,7 +187,7 @@ end # To model the system as a hybrid automaton, it is useful to work with symbolic # state variables. -using Symbolics #!jl +using Symbolics const var = @variables x y vx vy t @@ -350,7 +350,7 @@ solz = overapproximate(sol, Zonotope); using Plots, Plots.PlotMeasures, LaTeXStrings #!jl #!jl import DisplayAs #hide -fig = plot(solz; vars=(1, 2), xlab="x", ylab="y") +fig = plot(solz; vars=(1, 2), xlab="x", ylab="y") #!jl #!jl fig = DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -360,32 +360,32 @@ fig = plot(solz; vars=(1, 2), xlab="x", ylab="y") # using different colors. The plotting code is more sophisticated than before, # and is an example of finer control of the options used to plot the flowpipes. -idx_approaching = findall(x -> x == 1, location.(solz)) -idx_attempt = findall(x -> x == 2, location.(solz)) -idx_aborting = findall(x -> x == 3, location.(solz)) +idx_approaching = findall(x -> x == 1, location.(solz)) #!jl +idx_attempt = findall(x -> x == 2, location.(solz)) #!jl +idx_aborting = findall(x -> x == 3, location.(solz)) #!jl -fig = plot(; legend=:bottomright, tickfont=font(10, "Times"), guidefontsize=15, - xlab=L"x", ylab=L"y", lw=0.0, xtick=[-750, -500, -250, 0, 250.0], - ytick=[-400, -300, -200, -100, 0.0], xlims=(-1000.0, 400.0), ylims=(-450.0, 0.0), - size=(600, 600), bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) +fig = plot(; legend=:bottomright, tickfont=font(10, "Times"), guidefontsize=15, #!jl + xlab=L"x", ylab=L"y", lw=0.0, xtick=[-750, -500, -250, 0, 250.0], #!jl + ytick=[-400, -300, -200, -100, 0.0], xlims=(-1000.0, 400.0), ylims=(-450.0, 0.0), #!jl + size=(600, 600), bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) #!jl -plot!(fig, solz[idx_approaching[1]]; vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1, - lab="Approaching") -for k in idx_approaching[2:end] - plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1) -end +plot!(fig, solz[idx_approaching[1]]; vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1, #!jl + lab="Approaching") #!jl +for k in idx_approaching[2:end] #!jl + plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1) #!jl +end #!jl -plot!(fig, solz[idx_attempt[1]]; vars=(1, 2), lw=0.0, color=:red, alpha=1, lab="Attempt") -for k in idx_attempt[2:end] - plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:red, alpha=1) -end +plot!(fig, solz[idx_attempt[1]]; vars=(1, 2), lw=0.0, color=:red, alpha=1, lab="Attempt") #!jl +for k in idx_attempt[2:end] #!jl + plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:red, alpha=1) #!jl +end #!jl -plot!(fig, solz[idx_aborting[1]]; vars=(1, 2), lw=0.0, color=:cyan, alpha=1, lab="Aborting") -for k in idx_aborting[2:end] - plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:cyan, alpha=1) -end +plot!(fig, solz[idx_aborting[1]]; vars=(1, 2), lw=0.0, color=:cyan, alpha=1, lab="Aborting") #!jl +for k in idx_aborting[2:end] #!jl + plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:cyan, alpha=1) #!jl +end #!jl -fig +fig #!jl #!jl fig = DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/SquareWaveOscillator/SquareWaveOscillator.jl b/examples/SquareWaveOscillator/SquareWaveOscillator.jl index c6b88e3c5..f0404ccfc 100644 --- a/examples/SquareWaveOscillator/SquareWaveOscillator.jl +++ b/examples/SquareWaveOscillator/SquareWaveOscillator.jl @@ -7,7 +7,7 @@ # ## Model description -using ReachabilityAnalysis, Symbolics #!jl +using ReachabilityAnalysis, Symbolics function multistable_oscillator(; X0=Interval(0.0, 0.05), V₊=+13.5, V₋=-13.5, @@ -54,7 +54,7 @@ sol = solve(prob; T=100e-4, alg=INT(; δ=1.E-6), fixpoint_check=false); # Below we print the automaton location of each flowpipe segment: -location.(sol)' +location.(sol)' #!jl # ## Results @@ -64,7 +64,7 @@ using Plots #!jl old_ztol = LazySets._ztol(Float64) #!jl LazySets.set_ztol(Float64, 1e-8); # use higher precision for the plots #!jl -fig = plot(sol; vars=(0, 1), xlab="t", ylab="v-") +fig = plot(sol; vars=(0, 1), xlab="t", ylab="v-") #!jl #!jl fig = DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -74,8 +74,8 @@ fig = plot(sol; vars=(0, 1), xlab="t", ylab="v-") # sets of the first flowpipe, we observe that only the last three actually # intersect with the guard: -fig = plot(sol[1][(end - 10):end]; vars=(0, 1), xlab="t", ylab="v-") -plot!(fig, x -> 6.75; xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red) +fig = plot(sol[1][(end - 10):end]; vars=(0, 1), xlab="t", ylab="v-") #!jl +plot!(fig, x -> 6.75; xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red) #!jl #!jl fig = DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -85,10 +85,10 @@ Xc = cluster(sol[1], [318, 319, 320], BoxClustering(1)); # Plotting this set matches with the flowpipe after the jump: -fig = plot(sol[1][(end - 10):end]; vars=(0, 1)) -plot!(fig, sol[2][1:10]; vars=(0, 1)) -plot!(fig, x -> 6.75; xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red) -plot!(fig, Xc[1]; vars=(0, 1), c=:grey) +fig = plot(sol[1][(end - 10):end]; vars=(0, 1)) #!jl +plot!(fig, sol[2][1:10]; vars=(0, 1)) #!jl +plot!(fig, x -> 6.75; xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red) #!jl +plot!(fig, Xc[1]; vars=(0, 1), c=:grey) #!jl #!jl fig = DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -100,11 +100,11 @@ plot!(fig, Xc[1]; vars=(0, 1), c=:grey) # set is contained in a previously explored starting set. sol = solve(prob; T=100e-4, alg=INT(; δ=1.E-6), fixpoint_check=true) -tspan(sol) +tspan(sol) #!jl #- -fig = plot(sol; vars=(0, 1), xlab="t", ylab="v-") +fig = plot(sol; vars=(0, 1), xlab="t", ylab="v-") #!jl #!jl fig = DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/TransmissionLine/TransmissionLine.jl b/examples/TransmissionLine/TransmissionLine.jl index e29562153..35730d397 100644 --- a/examples/TransmissionLine/TransmissionLine.jl +++ b/examples/TransmissionLine/TransmissionLine.jl @@ -114,6 +114,9 @@ # are useful constructors in the base package `LinearAlgebra` that greatly # simplify building matrices with special shape, as in our case, such as # diagonal and band matrices, with the types `Diagonal` and `Bidiagonal`. +# +# We consider the case of ``η = 20`` nodes as in [^AKS11], such that the system +# has ``n = 40`` state variables. using ReachabilityAnalysis, LinearAlgebra, SparseArrays @@ -125,7 +128,11 @@ function tline(; η=3, R=1.00, Rd=10.0, L=1e-10, C=1e-13 * 4.00) A = [A₁₁ A₁₂; A₂₁ A₂₂] B = sparse([η + 1], [1], 1 / L, 2η, 1) return A, B -end; +end + +η = 20 +n = 2η +A, B = tline(; η=η); # We can visualize the structure of the sparse coefficients matrix ``A`` for the # case ``η = 20`` with a `spy` plot: @@ -133,9 +140,8 @@ end; using Plots #!jl #!jl import DisplayAs #hide -A, _ = tline(; η=20) -fig = spy(A; legend=nothing, markersize=2.0, title="Sparsity pattern of A", - xlab="columns", ylab="rows") +fig = spy(A; legend=nothing, markersize=2.0, title="Sparsity pattern of A", #!jl + xlab="columns", ylab="rows") #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide @@ -193,13 +199,6 @@ end; # We are interested in the step response to an input voltage ``U_{in}(t)`` # arbitrarily varying over the domain ``U_{in} = [0.99, 1.01]`` for all # ``t ∈ [0, T]``. -# -# We consider the case of ``η = 20`` nodes as in [^AKS11], such that the system -# has ``n = 40`` state variables. - -η = 20 -n = 2η -A, B = tline(; η=η) Uin_ss = Interval(-0.2, 0.2) □(ϵ) = BallInf(zeros(n), ϵ) @@ -221,11 +220,11 @@ sol = solve(prob; T=0.7, alg=BOX(; δ=1e-3)); # To get the variable ``U_{out}``, we have to project onto the ``η``-th # coordinate and invert the sign of the flowpipe. -Uout_vs_t = @. (-1.0) * project(sol, η); +Uout_vs_t = @. (-1.0) * project(sol, η); #!jl # ## Results -fig = plot(Uout_vs_t; vars=(0, η), c=:blue, xlab="t", ylab="Uout", alpha=0.5, lw=0.5) +fig = plot(Uout_vs_t; vars=(0, η), c=:blue, xlab="t", ylab="Uout", alpha=0.5, lw=0.5) #!jl #!jl DisplayAs.Text(DisplayAs.PNG(fig)) #hide diff --git a/examples/VanDerPol/VanDerPol.jl b/examples/VanDerPol/VanDerPol.jl index 94c5fc284..968fdcd50 100644 --- a/examples/VanDerPol/VanDerPol.jl +++ b/examples/VanDerPol/VanDerPol.jl @@ -39,7 +39,7 @@ # The system has a stable limit cycle. This limit cycle becomes increasingly # sharp for higher values ``μ``. Here we consider the parameter ``μ = 1``. -using ReachabilityAnalysis #!jl +using ReachabilityAnalysis @taylorize function vanderpol!(du, u, p, t) x, y = u @@ -90,7 +90,7 @@ solz = overapproximate(sol, Zonotope); y_bound = ρ([0.0, 1.0], solz) @assert y_bound < 2.75 "the property should be proven" -y_bound +y_bound #!jl # This shows that the property is satisfied. Below we plot the flowpipe in the # x-y plane, together with the horizontal line ``y = 2.75``. @@ -179,11 +179,11 @@ ilast = cross_section(line, solz, [388]); # We can also calculate the length of each cross section. -lfirst = norm(ifirst.q - ifirst.p) +lfirst = norm(ifirst.q - ifirst.p) #!jl #- -llast = norm(ilast.q - ilast.p) +llast = norm(ilast.q - ilast.p) #!jl #- From f5c8433abab3e0d8e3f3e2a8bf4f165a79dac659 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 1 Mar 2024 08:19:17 +0100 Subject: [PATCH 2/3] adapt generated tests --- test/models/generated/Brusselator.jl | 14 ++++++++ test/models/generated/Building.jl | 23 +------------ test/models/generated/DuffingOscillator.jl | 8 ++--- test/models/generated/EMBrake.jl | 2 ++ test/models/generated/ISS.jl | 21 +----------- test/models/generated/LaubLoomis.jl | 25 ++------------ test/models/generated/Lorenz.jl | 19 ++--------- test/models/generated/LotkaVolterra.jl | 12 ++----- test/models/generated/OpAmp.jl | 33 +++---------------- test/models/generated/Platoon.jl | 22 ++----------- .../models/generated/ProductionDestruction.jl | 11 ++----- test/models/generated/Quadrotor.jl | 12 +------ test/models/generated/SEIR.jl | 8 ++--- test/models/generated/Spacecraft.jl | 31 ++--------------- test/models/generated/SquareWaveOscillator.jl | 17 ++-------- test/models/generated/TransmissionLine.jl | 16 +++------ test/models/generated/VanDerPol.jl | 7 ++-- 17 files changed, 47 insertions(+), 234 deletions(-) diff --git a/test/models/generated/Brusselator.jl b/test/models/generated/Brusselator.jl index 4f7cb119c..3e8052a3e 100644 --- a/test/models/generated/Brusselator.jl +++ b/test/models/generated/Brusselator.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis + @taylorize function brusselator!(du, u, p, t) local A = 1.0 local B = 1.5 @@ -11,6 +13,18 @@ return du end +U₀ = (0.8 .. 1.0) × (0.0 .. 0.2) +prob = @ivp(u' = brusselator!(u), u(0) ∈ U₀, dim:2) + +T = 18.0; + +alg = TMJets20(; orderT=6, orderQ=2) +sol = solve(prob; T=T, alg=alg); + U0(r) = BallInf([1.0, 1.0], r); bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim:2); + +sol_01 = solve(bruss(0.01); T=30.0, alg=alg); + +sol_1 = solve(bruss(0.1); T=30.0, alg=alg); diff --git a/test/models/generated/Building.jl b/test/models/generated/Building.jl index deed73bda..01392b3f2 100644 --- a/test/models/generated/Building.jl +++ b/test/models/generated/Building.jl @@ -1,5 +1,6 @@ using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath using ReachabilityBase.Arrays: SingleEntryVector +using ReachabilityAnalysis: add_dimension const x25 = SingleEntryVector(25, 48, 1.0) const x25e = SingleEntryVector(25, 49, 1.0) @@ -26,8 +27,6 @@ function building_BLDF01() return prob_BLDF01 end -using ReachabilityAnalysis: add_dimension - function building_BLDC01() A, B, X0, U = building_common() n = size(A, 1) @@ -42,52 +41,32 @@ prob_BLDF01 = building_BLDF01(); sol_BLDF01_dense = solve(prob_BLDF01; T=20.0, alg=LGG09(; δ=0.004, vars=(25), n=48)) -fig = plot(sol_BLDF01_dense; vars=(0, 25), linecolor=:blue, color=:blue, - alpha=0.8, lw=1.0, xlab="t", ylab="x25") - @assert ρ(x25, sol_BLDF01_dense) <= 5.1e-3 "the property should be proven" # BLDF01 - BDS01 @assert !(ρ(x25, sol_BLDF01_dense) <= 4e-3) "the property should not be proven" # BLDF01 - BDU01 -ρ(x25, sol_BLDF01_dense) @assert !(ρ(x25, sol_BLDF01_dense(20.0)) <= -0.78e-3) "the property should not be proven" # BLDF01 - BDU02 -ρ(x25, sol_BLDF01_dense(20.0)) sol_BLDF01_discrete = solve(prob_BLDF01; T=20.0, alg=LGG09(; δ=0.01, vars=(25), n=48, approx_model=NoBloating())); -fig = plot(sol_BLDF01_discrete; vars=(0, 25), linecolor=:blue, color=:blue, - alpha=0.8, lw=1.0, xlab="t", ylab="x25") - @assert ρ(x25, sol_BLDF01_discrete) <= 5.1e-3 "the property should be proven" # BLDF01 - BDS01 @assert !(ρ(x25, sol_BLDF01_discrete) <= 4e-3) "the property should not be proven" # BLDF01 - BDU01 -ρ(x25, sol_BLDF01_discrete) @assert !(ρ(x25, sol_BLDF01_discrete(20.0)) <= -0.78e-3) "the property should not be proven" # BLDF01 - BDU02 -ρ(x25, sol_BLDF01_discrete(20.0)) prob_BLDC01 = building_BLDC01(); sol_BLDC01_dense = solve(prob_BLDC01; T=20.0, alg=LGG09(; δ=0.005, vars=(25), n=49)) -fig = plot(sol_BLDC01_dense; vars=(0, 25), linecolor=:blue, color=:blue, - alpha=0.8, lw=1.0, xlab="t", ylab="x25") - @assert ρ(x25e, sol_BLDC01_dense) <= 5.1e-3 "the property should be proven" # BLDC01 - BDS01 @assert !(ρ(x25e, sol_BLDC01_dense) <= 4e-3) "the property should not be proven" # BLDC01 - BDU01 -ρ(x25e, sol_BLDC01_dense) @assert !(ρ(x25e, sol_BLDC01_dense(20.0)) <= -0.78e-3) "the property should not be proven" # BLDC01 - BDU02 -ρ(x25, sol_BLDF01_discrete(20.0)) sol_BLDC01_discrete = solve(prob_BLDC01; T=20.0, alg=LGG09(; δ=0.01, vars=(25), n=49, approx_model=NoBloating())) -fig = plot(sol_BLDC01_discrete; vars=(0, 25), linecolor=:blue, color=:blue, - alpha=0.8, lw=1.0, xlab="t", ylab="x25") - @assert ρ(x25e, sol_BLDC01_discrete) <= 5.1e-3 "the property should be proven" # BLDC01 - BDS01 @assert !(ρ(x25e, sol_BLDC01_discrete) <= 4e-3) "the property should not be proven" # BLDC01 - BDU01 -ρ(x25e, sol_BLDC01_discrete) @assert !(ρ(x25e, sol_BLDC01_discrete(20.0)) <= -0.78e-3) "the property should not be proven" # BLDC01 - BDU02 -ρ(x25e, sol_BLDC01_discrete(20.0)) diff --git a/test/models/generated/DuffingOscillator.jl b/test/models/generated/DuffingOscillator.jl index d064ddbcd..941062165 100644 --- a/test/models/generated/DuffingOscillator.jl +++ b/test/models/generated/DuffingOscillator.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis + const ω = 1.2 @taylorize function duffing!(du, u, p, t) @@ -21,9 +23,3 @@ prob = @ivp(x' = duffing!(x), x(0) ∈ X0, dim:2) T = 2 * pi / ω; sol = solve(prob; tspan=(0.0, 20 * T), alg=TMJets21a()); - -fig = plot(sol; vars=(0, 1), xlab="t", ylab="x", lw=0.2, color=:blue) - -fig = plot(sol; vars=(0, 2), xlab="t", ylab="v", lw=0.2, color=:blue) - -fig = plot(sol; vars=(1, 2), xlab="x", ylab="v", lw=0.5, color=:red) diff --git a/test/models/generated/EMBrake.jl b/test/models/generated/EMBrake.jl index b4f449e48..041f40f80 100644 --- a/test/models/generated/EMBrake.jl +++ b/test/models/generated/EMBrake.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis, SparseArrays + function embrake_common(; A, Tsample, ζ, x0) # continuous system EMbrake = @system(x' = A * x) diff --git a/test/models/generated/ISS.jl b/test/models/generated/ISS.jl index a179b82a6..b7f7b502d 100644 --- a/test/models/generated/ISS.jl +++ b/test/models/generated/ISS.jl @@ -1,3 +1,4 @@ +using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath using ReachabilityAnalysis: add_dimension path = @current_path("ISS", "ISS.jld2") @@ -31,27 +32,7 @@ prob_ISSF01 = ISSF01() sol_ISSF01 = solve(prob_ISSF01; T=20.0, alg=LGG09(; δ=6e-4, template=dirs, sparse=true, cache=false)); -dim(sol_ISSF01) - -πsol_ISSF01 = flatten(sol_ISSF01); - -dim(πsol_ISSF01) - -fig = Plots.plot(πsol_ISSF01[1:10:end]; vars=(0, 1), linecolor=:blue, color=:blue, - alpha=0.8, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0], - ytick=[-0.00075, -0.0005, -0.00025, 0, 0.00025, 0.0005, 0.00075], - xlims=(0.0, 20.0), ylims=(-0.00075, 0.00075), - bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) - dirs = CustomDirections([C3_ext, -C3_ext]) prob_ISSC01 = ISSC01() sol_ISSC01 = solve(prob_ISSC01; T=20.0, alg=LGG09(; δ=0.01, template=dirs, sparse=true, cache=false)); - -πsol_ISSC01 = flatten(sol_ISSC01); - -fig = Plots.plot(πsol_ISSC01; vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8, - lw=1.0, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0], - ytick=[-0.0002, -0.0001, 0.0, 0.0001, 0.0002], - xlims=(0.0, 20.0), ylims=(-0.0002, 0.0002), - bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) diff --git a/test/models/generated/LaubLoomis.jl b/test/models/generated/LaubLoomis.jl index 5ae9dbd9b..62387ae10 100644 --- a/test/models/generated/LaubLoomis.jl +++ b/test/models/generated/LaubLoomis.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis + @taylorize function laubloomis!(dx, x, p, t) dx[1] = 1.4 * x[3] - 0.9 * x[1] dx[2] = 2.5 * x[5] - 1.5 * x[2] @@ -24,9 +26,6 @@ sol_1 = solve(prob; T=20.0, alg=alg) sol_1z = overapproximate(sol_1, Zonotope); @assert ρ(e4, sol_1z) < 4.5 "the property should be proven" -ρ(e4, sol_1z) - -ρ(e4, sol_1z[end]) + ρ(-e4, sol_1z[end]) prob = laubloomis(; W=0.05) alg = TMJets(; abstol=1e-12, orderT=7, orderQ=1, adaptive=false) @@ -35,9 +34,6 @@ sol_2 = solve(prob; T=20.0, alg=alg) sol_2z = overapproximate(sol_2, Zonotope); @assert ρ(e4, sol_2z) < 4.5 "the property should be proven" -ρ(e4, sol_2z) - -ρ(e4, sol_2z[end]) + ρ(-e4, sol_2z[end]) prob = laubloomis(; W=0.1) alg = TMJets(; abstol=1e-12, orderT=7, orderQ=1, adaptive=false) @@ -46,20 +42,3 @@ sol_3 = solve(prob; T=20.0, alg=alg) sol_3z = overapproximate(sol_3, Zonotope); @assert ρ(e4, sol_3z) < 5.0 "the property should be proven" -ρ(e4, sol_3z) - -ρ(e4, sol_3z[end]) + ρ(-e4, sol_3z[end]) - -fig = plot(sol_3z; vars=(0, 4), linecolor="green", color=:green, alpha=0.8, lab="W = 0.1") -plot!(fig, sol_2z; vars=(0, 4), linecolor="blue", color=:blue, alpha=0.8, lab="W = 0.05") -plot!(fig, sol_1z; vars=(0, 4), linecolor="yellow", color=:yellow, alpha=0.8, lab="W = 0.01", - tickfont=font(10, "Times"), guidefontsize=15, - xlab=L"t", - ylab=L"x_4", - xtick=[0.0, 5.0, 10.0, 15.0, 20.0], ytick=[1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0], - xlims=(0.0, 20.0), ylims=(1.5, 5.02), - bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm, - size=(600, 600)) - -plot!(fig, x -> x, x -> 4.5, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="") -plot!(fig, x -> x, x -> 5.0, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="") diff --git a/test/models/generated/Lorenz.jl b/test/models/generated/Lorenz.jl index 420033271..a75ddcdd1 100644 --- a/test/models/generated/Lorenz.jl +++ b/test/models/generated/Lorenz.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis + @taylorize function lorenz!(du, u, p, t) local σ = 10.0 local β = 8.0 / 3.0 @@ -14,24 +16,7 @@ X0 = Hyperrectangle(; low=[0.9, 0.0, 0.0], high=[1.1, 0.0, 0.0]) prob = @ivp(x' = lorenz!(x), dim:3, x(0) ∈ X0); alg = TMJets(; abstol=1e-15, orderT=10, orderQ=2, maxsteps=50_000) - sol = solve(prob; T=10.0, alg=alg) - solz = overapproximate(sol, Zonotope); -fig = plot(solz; vars=(0, 1), xlab="t", ylab="x") - -fig = plot(solz(0.0 .. 1.5); vars=(0, 1), xlab="t", ylab="x", lw=0) -plot!(fig, x -> 20.0; c=:red, xlims=(0.0, 1.5), lab="") - @assert ρ([1.0, 0, 0], solz(0 .. 1.5)) < 20 "the property should be proven" - -ρ([1.0, 0, 0], solz(0 .. 1.5)) - -fig = plot(solz; vars=(0, 2), xlab="t", ylab="y") - --ρ([0.0, -1.0, 0.0], solz), ρ([0.0, 1.0, 0.0], solz) - -fig = plot(solz; vars=(0, 3), xlab="t", ylab="z") - -fig = plot(solz; vars=(1, 3), xlab="x", ylab="z") diff --git a/test/models/generated/LotkaVolterra.jl b/test/models/generated/LotkaVolterra.jl index d8eaa306c..bca388ddc 100644 --- a/test/models/generated/LotkaVolterra.jl +++ b/test/models/generated/LotkaVolterra.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis + @taylorize function lotkavolterra!(du, u, p, t) local α, β, γ, δ = 1.5, 1.0, 3.0, 1.0 @@ -14,10 +16,6 @@ prob = @ivp(x' = lotkavolterra!(x), dim:2, x(0) ∈ X0); sol = solve(prob; T=8.0, alg=TMJets()) solz = overapproximate(sol, Zonotope); -fig = plot(solz; vars=(1, 2), alpha=0.3, lw=0.0, xlab="x", ylab="y", - lab="Flowpipe", legend=:bottomright) -plot!(fig, X0; label="X(0)") - @taylorize function lotkavolterra_parametric!(du, u, p, t) x, y, αp, βp, γp, δp, ϵp = u xy = x * y @@ -40,9 +38,6 @@ prob = @ivp(u' = lotkavolterra_parametric!(u), dim:7, u(0) ∈ U0); sol = solve(prob; tspan=(0.0, 10.0)) solz = overapproximate(sol, Zonotope); -fig = plot(solz; vars=(1, 2), lw=0.3, title="Uncertain parameters", - lab="abstol = 1e-15", xlab="x", ylab="y") - □(ϵ) = BallInf([1.0, 1.0], ϵ) U0 = cartesian_product(□(0.05), convert(Hyperrectangle, p_int)) @@ -51,9 +46,6 @@ prob = @ivp(u' = lotkavolterra_parametric!(u), dim:7, u(0) ∈ U0); sol = solve(prob; T=10.0, alg=TMJets(; abstol=1e-10)) solz = overapproximate(sol, Zonotope); -fig = plot(solz; vars=(1, 2), color=:orange, lw=0.3, lab="ϵ = 0.05", - title="Uncertain u0 and uncertain parameters", xlab="x", ylab="y") - U0 = cartesian_product(□(0.01), convert(Hyperrectangle, p_int)) prob = @ivp(u' = lotkavolterra_parametric!(u), dim:7, u(0) ∈ U0); diff --git a/test/models/generated/OpAmp.jl b/test/models/generated/OpAmp.jl index d8558dc53..07835a4db 100644 --- a/test/models/generated/OpAmp.jl +++ b/test/models/generated/OpAmp.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis + function opamp_nondet(; X0=Singleton([0.0]), R₁=2.0, R₂=6.0, C=1.e-3, Ein=Interval(1.9, 2.1)) @@ -13,21 +15,13 @@ end; sol_nondet = solve(opamp_nondet(); T=0.1, alg=INT(; δ=1e-4)); -fig = plot(sol_nondet; vars=(0, 1), xlab=L"t", ylab=L"e_{out}", - title="Solution for nondeterministic input", lw=0.2) - -Δt = 0 .. 0.04 - -fig = plot(; xlab=L"t", ylab=L"e_{out}(t)", title="Solution for nondeterministic input") - sol_nondet = solve(opamp_nondet(; X0=Interval(-1.0, 1.0)); T=0.1, alg=INT(; δ=1e-4)) -plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = -1 .. 1", lw=0.2) sol_nondet = solve(opamp_nondet(; X0=Interval(-0.5, 0.5)); T=0.1, alg=INT(; δ=1e-4)) -plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = -0.5 .. 0.5", lw=0.2) sol_nondet = solve(opamp_nondet(); T=0.1, alg=INT(; δ=1e-4)) -plot!(fig, sol_nondet(Δt); vars=(0, 1), lab="X0 = 0", lw=0.2) + +using Symbolics function opamp_with_saturation(; X0=Singleton(zeros(2)), R₁=2.0, R₂=6.0, C=1.e-3, @@ -65,13 +59,6 @@ prob_const = opamp_with_saturation(; X0=X0, γ=0.0, δ=0.0, Es=3.0); sol_const = solve(prob_const; T=0.1, alg=BOX(; δ=1e-4)); -fig = plot(sol_const; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", - title="Constant input signal", lw=0.2) -plot!(fig, x -> x, x -> 1.5, 0.0, 0.1; color="red", lw=2, ls=:dash, lab="") - -fig = plot(sol_const; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", - title="Solution for constant input", lw=0.2) - X0 = Singleton(zeros(2)); prob_lin = opamp_with_saturation(; X0=X0, γ=0.0, δ=100.0, Es=1.0) @@ -80,15 +67,3 @@ prob_exp = opamp_with_saturation(; X0=X0, γ=-100.0, δ=100.0, Es=3.0); sol_lin = solve(prob_lin; T=0.1, alg=BOX(; δ=1e-4)) sol_exp = solve(prob_exp; T=0.1, alg=BOX(; δ=1e-4)); - -fig = plot(sol_lin; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", - title="Linear input with saturation", lw=0.2) - -fig = plot(sol_exp; vars=(0, 2), xlab=L"t", ylab=L"e_{in}(t)", - title="Exponential input", lw=0.2) - -fig = plot(sol_lin; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", - title="Solution for linear input", lw=0.2) - -fig = plot(sol_exp; vars=(0, 1), xlab=L"t", ylab=L"e_{out}(t)", - title="Solution for exponential input", lw=0.2) diff --git a/test/models/generated/Platoon.jl b/test/models/generated/Platoon.jl index 1388a5504..d9fadcded 100644 --- a/test/models/generated/Platoon.jl +++ b/test/models/generated/Platoon.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis, SparseArrays, Symbolics + const var = @variables x[1:9] t; function platoon_connected(; deterministic_switching::Bool=true, c1=5.0) @@ -117,7 +119,6 @@ end prob_PLAD01 = platoon(); boxdirs = BoxDirections(10) -length(boxdirs) alg = BOX(; δ=0.01) sol_PLAD01_BND42 = solve(prob_PLAD01; @@ -129,17 +130,7 @@ sol_PLAD01_BND42 = solve(prob_PLAD01; @assert dmin_specification(sol_PLAD01_BND42, 42) "the property should be proven" --ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND42) - --ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND42) - --ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND42) - -fig = plot(sol_PLAD01_BND42; vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND42", lw=0.1) -plot!(x -> x, x -> -42.0, 0.0, 20.0; linewidth=2, color="red", ls=:dash, leg=nothing) - octdirs = OctDirections(10) -length(octdirs) alg = LGG09(; δ=0.03, template=octdirs, approx_model=Forward(; setops=octdirs)) sol_PLAD01_BND30 = solve(prob_PLAD01; @@ -150,12 +141,3 @@ sol_PLAD01_BND30 = solve(prob_PLAD01; T=20.0); @assert dmin_specification(sol_PLAD01_BND30, 30) "the property should be proven" - --ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND30) - --ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND30) - --ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND30) - -fig = plot(sol_PLAD01_BND30; vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND30", lw=0.1) -plot!(x -> x, x -> -30.0, 0.0, 20.0; linewidth=2, color="red", ls=:dash, leg=nothing) diff --git a/test/models/generated/ProductionDestruction.jl b/test/models/generated/ProductionDestruction.jl index 415085a39..5931c3c7d 100644 --- a/test/models/generated/ProductionDestruction.jl +++ b/test/models/generated/ProductionDestruction.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis, Symbolics + @variables x y z positive_orthant = HPolyhedron([x >= 0, y >= 0, z >= 0], [x, y, z]); @@ -40,9 +42,6 @@ sol = solve(prob; T=100.0, alg=TMJets(; abstol=1e-12, orderT=6, orderQ=1)); property, vol = prod_dest_verif(sol) @assert property "the property should be proven" -vol - -fig = plot(sol; vars=(0, 3), lc=:orange, c=:orange, alpha=0.3, lab="I", xlab="t", ylab="z") @taylorize function prod_dest_IP!(du, u, p, t) x, y, z, a = u[1], u[2], u[3], u[4] @@ -64,9 +63,6 @@ sol = solve(prob; T=100.0, alg=TMJets(; abstol=9e-13, orderT=6, orderQ=1)); property, vol = prod_dest_verif(sol) @assert property "the property should be proven" -vol - -fig = plot(sol; vars=(0, 3), lc=:blue, c=:blue, alpha=0.3, lab="P", xlab="t", ylab="z") X0 = Hyperrectangle(; low=[9.5, 0.01, 0.01, 0.296], high=[10, 0.01, 0.01, 0.304]) prob = @ivp(x' = prod_dest_IP!(x), dim:4, x(0) ∈ X0); @@ -75,6 +71,3 @@ sol = solve(prob; T=100.0, alg=TMJets(; abstol=1e-12, orderT=6, orderQ=1)); property, vol = prod_dest_verif(sol) @assert property "the property should be proven" -vol - -fig = plot(sol; vars=(0, 3), lc=:red, c=:red, alpha=0.3, lab="I & P", xlab="t", ylab="z") diff --git a/test/models/generated/Quadrotor.jl b/test/models/generated/Quadrotor.jl index 484ee14df..ee4bedb14 100644 --- a/test/models/generated/Quadrotor.jl +++ b/test/models/generated/Quadrotor.jl @@ -1,4 +1,5 @@ using ReachabilityAnalysis +using ReachabilityBase.Arrays: SingleEntryVector # parameters of the model const g = 9.81 # gravity constant in m/s^2 @@ -75,8 +76,6 @@ const u₃ = 0.0 return dx end; -using ReachabilityBase.Arrays: SingleEntryVector - const T = 5.0 const v3 = SingleEntryVector(3, 12, 1.0) @@ -109,8 +108,6 @@ function quadrotor(; Wpos, Wvel) return prob end; -cases = ["Δ=0.1", "Δ=0.4", "Δ=0.8"]; - Wpos = 0.1 Wvel = 0.1 prob = quadrotor(; Wpos=Wpos, Wvel=Wvel) @@ -140,10 +137,3 @@ sol = solve(prob; T=T, alg=alg) solz3 = overapproximate(sol, Zonotope); @assert !quad_property(solz3) "the property should not be proven" - -fig = plot(solz3; vars=(0, 3), linecolor="green", color=:green, alpha=0.8) -plot!(solz2; vars=(0, 3), linecolor="blue", color=:blue, alpha=0.8) -plot!(solz1; vars=(0, 3), linecolor="yellow", color=:yellow, alpha=0.8, - xlab="t", ylab="x3", - xtick=[0.0, 1.0, 2.0, 3.0, 4.0, 5.0], ytick=[-1.0, -0.5, 0.0, 0.5, 1.0, 1.5], - xlims=(0.0, 5.0), ylims=(-1.0, 1.5)) diff --git a/test/models/generated/SEIR.jl b/test/models/generated/SEIR.jl index 95e1cc061..b970d52da 100644 --- a/test/models/generated/SEIR.jl +++ b/test/models/generated/SEIR.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis + @taylorize function seir!(dx, x, p, t) S, E, I, R, α, β, γ = x @@ -28,9 +30,3 @@ prob = @ivp(x' = seir!(x), dim:7, x(0) ∈ X0); sol = solve(prob; T=200.0, alg=TMJets21a(; orderT=7, orderQ=1)) solz = overapproximate(sol, Zonotope); - -fig = plot(; legend=:outerright) -plot!(fig, solz; vars=(0, 1), color=:blue, lw=0.0, lab="S") -plot!(fig, solz; vars=(0, 2), color=:green, lw=0.0, lab="E") -plot!(fig, solz; vars=(0, 3), color=:red, lw=0.0, lab="I") -plot!(fig, solz; vars=(0, 4), color=:grey, lw=0.0, lab="R") diff --git a/test/models/generated/Spacecraft.jl b/test/models/generated/Spacecraft.jl index e8b39d22f..21e843d39 100644 --- a/test/models/generated/Spacecraft.jl +++ b/test/models/generated/Spacecraft.jl @@ -118,6 +118,8 @@ end return du end +using Symbolics + const var = @variables x y vx vy t function spacecraft(; abort_time=(120.0, 150.0)) @@ -241,32 +243,3 @@ solz = overapproximate(sol, Zonotope); @assert line_of_sight(solz) "the property should be proven" @assert !collision_avoidance(solz) "the property cannot be proven with these settings" @assert velocity_constraint(solz) "the property should be proven" - -fig = plot(solz; vars=(1, 2), xlab="x", ylab="y") - -idx_approaching = findall(x -> x == 1, location.(solz)) -idx_attempt = findall(x -> x == 2, location.(solz)) -idx_aborting = findall(x -> x == 3, location.(solz)) - -fig = plot(; legend=:bottomright, tickfont=font(10, "Times"), guidefontsize=15, - xlab=L"x", ylab=L"y", lw=0.0, xtick=[-750, -500, -250, 0, 250.0], - ytick=[-400, -300, -200, -100, 0.0], xlims=(-1000.0, 400.0), ylims=(-450.0, 0.0), - size=(600, 600), bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm) - -plot!(fig, solz[idx_approaching[1]]; vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1, - lab="Approaching") -for k in idx_approaching[2:end] - plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1) -end - -plot!(fig, solz[idx_attempt[1]]; vars=(1, 2), lw=0.0, color=:red, alpha=1, lab="Attempt") -for k in idx_attempt[2:end] - plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:red, alpha=1) -end - -plot!(fig, solz[idx_aborting[1]]; vars=(1, 2), lw=0.0, color=:cyan, alpha=1, lab="Aborting") -for k in idx_aborting[2:end] - plot!(fig, solz[k]; vars=(1, 2), lw=0.0, color=:cyan, alpha=1) -end - -fig diff --git a/test/models/generated/SquareWaveOscillator.jl b/test/models/generated/SquareWaveOscillator.jl index 97f51bb0b..e95eea33c 100644 --- a/test/models/generated/SquareWaveOscillator.jl +++ b/test/models/generated/SquareWaveOscillator.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis, Symbolics + function multistable_oscillator(; X0=Interval(0.0, 0.05), V₊=+13.5, V₋=-13.5, R=20.E3, C=5.5556E-8, @@ -37,21 +39,6 @@ prob = multistable_oscillator(); sol = solve(prob; T=100e-4, alg=INT(; δ=1.E-6), fixpoint_check=false); -location.(sol)' - -fig = plot(sol; vars=(0, 1), xlab="t", ylab="v-") - -fig = plot(sol[1][(end - 10):end]; vars=(0, 1), xlab="t", ylab="v-") -plot!(fig, x -> 6.75; xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red) - Xc = cluster(sol[1], [318, 319, 320], BoxClustering(1)); -fig = plot(sol[1][(end - 10):end]; vars=(0, 1)) -plot!(fig, sol[2][1:10]; vars=(0, 1)) -plot!(fig, x -> 6.75; xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red) -plot!(fig, Xc[1]; vars=(0, 1), c=:grey) - sol = solve(prob; T=100e-4, alg=INT(; δ=1.E-6), fixpoint_check=true) -tspan(sol) - -fig = plot(sol; vars=(0, 1), xlab="t", ylab="v-") diff --git a/test/models/generated/TransmissionLine.jl b/test/models/generated/TransmissionLine.jl index e9aef886c..e15dcd05f 100644 --- a/test/models/generated/TransmissionLine.jl +++ b/test/models/generated/TransmissionLine.jl @@ -8,11 +8,11 @@ function tline(; η=3, R=1.00, Rd=10.0, L=1e-10, C=1e-13 * 4.00) A = [A₁₁ A₁₂; A₂₁ A₂₂] B = sparse([η + 1], [1], 1 / L, 2η, 1) return A, B -end; +end -A, _ = tline(; η=20) -fig = spy(A; legend=nothing, markersize=2.0, title="Sparsity pattern of A", - xlab="columns", ylab="rows") +η = 20 +n = 2η +A, B = tline(; η=η); function scale!(s, α=1.0) s.A .*= α @@ -20,10 +20,6 @@ function scale!(s, α=1.0) return s end; -η = 20 -n = 2η -A, B = tline(; η=η) - Uin_ss = Interval(-0.2, 0.2) □(ϵ) = BallInf(zeros(n), ϵ) X0 = -inv(Matrix(A)) * B * Uin_ss ⊕ □(0.001) @@ -36,7 +32,3 @@ scale!(s, α) prob = InitialValueProblem(s, X0); sol = solve(prob; T=0.7, alg=BOX(; δ=1e-3)); - -Uout_vs_t = @. (-1.0) * project(sol, η); - -fig = plot(Uout_vs_t; vars=(0, η), c=:blue, xlab="t", ylab="Uout", alpha=0.5, lw=0.5) diff --git a/test/models/generated/VanDerPol.jl b/test/models/generated/VanDerPol.jl index 6d0414ef8..77367bc41 100644 --- a/test/models/generated/VanDerPol.jl +++ b/test/models/generated/VanDerPol.jl @@ -1,3 +1,5 @@ +using ReachabilityAnalysis + @taylorize function vanderpol!(du, u, p, t) x, y = u local μ = 1.0 @@ -22,7 +24,6 @@ solz = overapproximate(sol, Zonotope); y_bound = ρ([0.0, 1.0], solz) @assert y_bound < 2.75 "the property should be proven" -y_bound line = LineSegment([1, 2.0], [2.0, 2.5]) @@ -42,8 +43,4 @@ end; ifirst = cross_section(line, solz, 1:13) ilast = cross_section(line, solz, [388]); -lfirst = norm(ifirst.q - ifirst.p) - -llast = norm(ilast.q - ilast.p) - @assert ilast ⊆ ifirst "the cross section should get smaller" From c870fa71252b6daf51c896277b36674a63b0cd1d Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 1 Mar 2024 08:44:24 +0100 Subject: [PATCH 3/3] remove unused dependencies in test/Project.toml --- test/Project.toml | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 6947eca3d..27569d99a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,17 +2,11 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" -DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" -ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" Flowstar = "a8054ddd-9dca-4d20-8ffe-ae96ec1541f1" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" -JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" -Optim = "429524aa-4258-5aef-a3af-852621145aeb" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -24,16 +18,10 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Aqua = "0.8" CDDLib = "0.5 - 0.9" DifferentialEquations = "6.17, 7" -DynamicPolynomials = "0.3 - 0.5" Expokit = "0.2" -ExponentialUtilities = "1" Flowstar = "0.2.4" IntervalArithmetic = "0.16 - 0.20, =0.20.9" # new versions require updates and are incompatible with dependencies -JuMP = "0.21 - 0.23, 1" LazySets = "2.3" -MultivariatePolynomials = "0.3 - 0.5" -Optim = "0.18 - 0.20, 1" -Plots = "0.25 - 0.29, 1" Polyhedra = "0.5 - 0.7" StaticArrays = "0.12, 1" Symbolics = "4 - 5"