Skip to content

Commit

Permalink
update generated tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@schillic
schillic committed Feb 27, 2024
1 parent ac330ab commit 1acd3d7
Show file tree
Hide file tree
Showing 18 changed files with 521 additions and 634 deletions.
6 changes: 3 additions & 3 deletions examples/ISS/ISS.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,8 +144,8 @@ dim(πsol_ISSF01)
using Plots, Plots.PlotMeasures, LaTeXStrings #!jl
#!jl import DisplayAs #hide

old_ztol = LazySets._ztol(Float64)
LazySets.set_ztol(Float64, 1e-8); # use higher precision for the plots
old_ztol = LazySets._ztol(Float64) #!jl
LazySets.set_ztol(Float64, 1e-8); # use higher precision for the plots #!jl

#-

Expand Down Expand Up @@ -178,7 +178,7 @@ fig = Plots.plot(πsol_ISSC01; vars=(0, 1), linecolor=:blue, color=:blue, alpha=

#-

LazySets.set_ztol(Float64, old_ztol); # reset precision
LazySets.set_ztol(Float64, old_ztol); # reset precision #!jl

# ## References

Expand Down
20 changes: 10 additions & 10 deletions test/models/generated/Brusselator.jl
View file
Original file line number Diff line number Diff line change
@@ -1,16 +1,16 @@
const A = 1.0
const B = 1.5
const B1 = B + 1

@taylorize function brusselator!(du, u, p, t)
local A = 1.0
local B = 1.5
local B1 = B + 1

x, y = u
x² = x * x
aux = x² * y
du[1] = A + aux - B1 * x
du[2] = B * x - aux

x²y = x^2 * y
du[1] = A + x²y - B1 * x
du[2] = B * x - x²y
return du
end

U0(r) = Singleton([1.0, 1.0]) ⊕ BallInf(zeros(2), r)
U0(r) = BallInf([1.0, 1.0], r);

bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim:2)
bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim:2);
113 changes: 47 additions & 66 deletions test/models/generated/Building.jl
View file
Original file line number Diff line number Diff line change
@@ -1,112 +1,93 @@
using ReachabilityAnalysis, SparseArrays, JLD2
using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath
using ReachabilityBase.Arrays: SingleEntryVector

LazySets.set_ztol(Float64, 1e-14)
const x25 = SingleEntryVector(25, 48, 1.0)
const x25e = SingleEntryVector(25, 49, 1.0)

const x25 = [zeros(24); 1.0; zeros(23)]
const x25e = vcat(x25, 0.0)
path = @current_path("Building", "building.jld2")

examples_dir = normpath(@__DIR__, "..", "..", "..", "examples")
building_path = joinpath(examples_dir, "Building", "building.jld2")
function building_common()
@load path A B

c = [fill(0.000225, 10); fill(0.0, 38)]
r = [fill(0.000025, 10); fill(0.0, 14); 0.0001; fill(0.0, 23)]
X0 = Hyperrectangle(c, r)

function building_BLDF01()
@load building_path A B
n = size(A, 1)
U = Interval(0.8, 1.0)
S = @system(x' = Ax + Bu, u ∈ U, x ∈ Universe(n))

# initial states
center_X0 = [fill(0.000225, 10); fill(0.0, 38)]
radius_X0 = [fill(0.000025, 10); fill(0.0, 14); 0.0001; fill(0.0, 23)]
X0 = Hyperrectangle(center_X0, radius_X0)
return A, B, X0, U
end

return prob_BLDF01 = InitialValueProblem(S, X0)
function building_BLDF01()
A, B, X0, U = building_common()
n = size(A, 1)
S = @system(x' = A * x + B * u, u ∈ U, x ∈ Universe(n))
prob_BLDF01 = InitialValueProblem(S, X0)
return prob_BLDF01
end

using ReachabilityAnalysis: add_dimension

function building_BLDC01()
@load building_path A B
A, B, X0, U = building_common()
n = size(A, 1)
U = Interval(0.8, 1.0)

# initial states
center_X0 = [fill(0.000225, 10); fill(0.0, 38)]
radius_X0 = [fill(0.000025, 10); fill(0.0, 14); 0.0001; fill(0.0, 23)]
X0 = Hyperrectangle(center_X0, radius_X0)

Ae = add_dimension(A)
Ae[1:n, end] = B
return prob_BLDC01 = @ivp(x' = Ae * x, x(0) ∈ X0 × U)
end

using Plots
prob_BLDC01 = @ivp(x' = Ae * x, x(0) ∈ X0 × U)
return prob_BLDC01
end;

prob_BLDF01 = building_BLDF01()
prob_BLDF01 = building_BLDF01();

sol_BLDF01_dense = solve(prob_BLDF01; T=20.0, alg=LGG09(; δ=0.004, vars=(25), n=48));
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,
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)

ρ(x25, sol_BLDF01_dense) <= 5.1e-3 # BLDF01 - BDS01

ρ(x25, sol_BLDF01_dense) <= 4e-3 # BLDF01 - BDU01

@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)) <= -0.78e-3 # BLDF01 - BDU02

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,
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)

ρ(x25, sol_BLDF01_discrete) <= 5.1e-3 # BLDF01 - BDS01

ρ(x25, sol_BLDF01_discrete)

ρ(x25, sol_BLDF01_discrete) <= 4e-3 # BLDF01 - BDU01

@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)) <= -0.78e-3 # BLDF01 - BDU02

prob_BLDC01 = building_BLDC01()
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,
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)

ρ(x25e, sol_BLDC01_dense) <= 5.1e-3 # BLDC01 - BDS01

ρ(x25e, sol_BLDC01_dense) <= 4e-3 # BLDC01 - BDU01

@assert !(ρ(x25e, sol_BLDC01_dense(20.0)) <= -0.78e-3) "the property should not be proven" # BLDC01 - BDU02
ρ(x25, sol_BLDF01_discrete(20.0))

ρ(x25e, sol_BLDC01_dense(20.0)) <= -0.78e-3 # BLDC01 - BDU02

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")

ρ(x25e, sol_BLDC01_discrete)

ρ(x25e, sol_BLDC01_discrete) <= 5.1e-3 # BLDC01 - BDS01
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)

ρ(x25e, sol_BLDC01_discrete) <= 4e-3 # BLDC01 - BDU01

@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)) <= -0.78e-3 # BLDC01 - BDU02
14 changes: 8 additions & 6 deletions test/models/generated/DuffingOscillator.jl
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using ReachabilityAnalysis, Plots
const ω = 1.2

@taylorize function duffing!(du, u, p, t)
local α = -1.0
Expand All @@ -7,16 +7,18 @@ using ReachabilityAnalysis, Plots
local γ = 0.37

x, v = u

f = γ * cos(ω * t)

du[1] = u[2]
return du[2] = -α * x - δ * v - β * x^3 + f
du[2] = -α * x - δ * v - β * x^3 + f
return du
end

ω = 1.2
T = 2 * pi / ω
X0 = Singleton([1.0, 0.0]) ⊕ BallInf(zeros(2), 0.1)
prob = @ivp(x' = duffing!(x), x(0) ∈ X0, dim = 2);
X0 = Hyperrectangle([1.0, 0.0], [0.1, 0.1])
prob = @ivp(x' = duffing!(x), x(0) ∈ X0, dim:2)

T = 2 * pi / ω;

sol = solve(prob; tspan=(0.0, 20 * T), alg=TMJets21a());

Expand Down
82 changes: 21 additions & 61 deletions test/models/generated/EMBrake.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,19 @@
ReachabilityBase.Comparison.set_rtol(Float64, 1e-12)
ReachabilityBase.Comparison.set_ztol(Float64, 1e-13)
function embrake_common(; A, Tsample, ζ, x0)
# continuous system
EMbrake = @system(x' = A * x)

# Fixed parameters
# initial condition
X₀ = Singleton([0.0, 0, 0, 0])

# reset map
Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1.0, 1.0, -1.0, -Tsample, 1.0], 4, 4)
br = sparsevec([3, 4], [x0, Tsample * x0], 4)
reset_map(X) = Ar * X + br

# hybrid system with clocked affine dynamics
ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
return IVP(ha, X₀)
end;

function embrake_no_pv(; Tsample=1.E-4, ζ=1e-6, x0=0.05)
# model's constants
Expand All @@ -19,26 +31,8 @@ function embrake_no_pv(; Tsample=1.E-4, ζ=1e-6, x0=0.05)
0 0 0 0;
0 0 0 0])

EMbrake = @system(x' = Ax)

# initial conditions
I₀ = Singleton([0.0])
x₀ = Singleton([0.0])
xe₀ = Singleton([0.0])
xc₀ = Singleton([0.0])
X₀ = I₀ × x₀ × xe₀ × xc₀

# reset map
Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1.0, 1.0, -1.0, -Tsample, 1.0], 4, 4)
br = sparsevec([3, 4], [x0, Tsample * x0], 4)
reset_map(X) = Ar * X + br

# hybrid system with clocked linear dynamics
ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
return IVP(ha, X₀)
end

# Parameter variation
return embrake_common(; A=A, Tsample=Tsample, ζ=ζ, x0=x0)
end;

function embrake_pv_1(; Tsample=1.E-4, ζ=1e-6, Δ=3.0, x0=0.05)
# model's constants
Expand All @@ -57,26 +51,8 @@ function embrake_pv_1(; Tsample=1.E-4, ζ=1e-6, Δ=3.0, x0=0.05)
0 0 0 0;
0 0 0 0])

EMbrake = @system(x' = Ax)

# initial conditions
I₀ = Singleton([0.0])
x₀ = Singleton([0.0])
xe₀ = Singleton([0.0])
xc₀ = Singleton([0.0])
X₀ = I₀ × x₀ × xe₀ × xc₀

# reset map
Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1.0, 1.0, -1.0, -Tsample, 1.0], 4, 4)
br = sparsevec([3, 4], [x0, Tsample * x0], 4)
reset_map(X) = Ar * X + br

# hybrid system with clocked linear dynamics
ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
return IVP(ha, X₀)
end

# Extended parameter variation
return embrake_common(; A=A, Tsample=Tsample, ζ=ζ, x0=x0)
end;

function embrake_pv_2(; Tsample=1.E-4, ζ=1e-6, x0=0.05, χ=5.0)
# model's constants
Expand All @@ -95,21 +71,5 @@ function embrake_pv_2(; Tsample=1.E-4, ζ=1e-6, x0=0.05, χ=5.0)
0 0 0 0;
0 0 0 0])

EMbrake = @system(x' = Ax)

# initial conditions
I₀ = Singleton([0.0])
x₀ = Singleton([0.0])
xe₀ = Singleton([0.0])
xc₀ = Singleton([0.0])
X₀ = I₀ × x₀ × xe₀ × xc₀

# reset map
Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1.0, 1.0, -1.0, -Tsample, 1.0], 4, 4)
br = sparsevec([3, 4], [x0, Tsample * x0], 4)
reset_map(X) = Ar * X + br

# hybrid system with clocked linear dynamics
ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
return IVP(ha, X₀)
end
return embrake_common(; A=A, Tsample=Tsample, ζ=ζ, x0=x0)
end;
Loading

0 comments on commit 1acd3d7

Please sign in to comment.