-
Notifications
You must be signed in to change notification settings - Fork 17
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Seemingly unnecessary overapproximation given zonotope starting set #464
Comments
My workaround for now: using LazySets: UnionSetArray, LineSegment, BallInf, CartesianProduct
using ReachabilityAnalysis
using Plots
# https://github.com/JuliaReach/ReachabilityAnalysis.jl/issues/464
using LazySets: ConvexHull, concretize, Zonotope
concretize_skip(s::LazySets.ConvexHull) = ConvexHull(concretize(s.X), concretize(s.Y))
concretize_skip(s) = s
ReachabilityAnalysis._apply_setops(X::LazySet, ::Val{:custom_gng}) = concretize_skip(X)
ReachabilityAnalysis._alias(setops::Val{:custom_gng}) = setops
source_set = [1 1 ; 0 1] * BallInf([0.0, 0.0], 0.5);
source_set = LazySets.concretize(source_set)
# no-op dynamics
A = zeros(2, 2)
c = zeros(2)
prob = @ivp(x' = A * x + c, x(0) ∈ source_set)
alg = GLGM06(δ=0.1, max_order=10,
approx_model=ReachabilityAnalysis.Forward(sih=:concrete, exp=:base, setops=:custom_gng)
)
tspan = (0.0, 1.0)
sol = solve(prob; tspan, alg)
# project out time dimension
sol_p = LazySets.project(sol, [1, 2])
@show sol_p.Xk[1].X
reach_set_all = set(sol_p)
plot(reach_set_all) |
The reason that You can check this by observing the following: approx_model = Forward()
δ = 0.1
probn = normalize(prob)
prob_discr = discretize(probn, δ, approx_model);
# the following set is what we call the initial set after discretization: and it represents the reach-set
# in the time *interval* [0, δ]
Ω0 = initial_state(prob_discr)
# now the type is quite complicated
@show typeof(Ω0)
typeof(Ω0) = ConvexHull{Float64, Zonotope{Float64, Vector{Float64}, Matrix{Float64}}, MinkowskiSum{Float64, MinkowskiSum{Float64, LinearMap{Float64, Zonotope{Float64, Vector{Float64}, Matrix{Float64}}, Float64, Matrix{Float64}}, MinkowskiSum{Float64, LinearMap{Float64, Singleton{Float64, Vector{Float64}}, Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}}}, Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}}}
# but if you plot it, it'll be the source set as expected
plot(Ω0)
# we can check equivalence of sets
isequivalent(Ω0, source_set) # broken, should be fixed in LazySets
....
isequivalent(concretize(Ω0), source_set) # ok
true
# the box overapproximation
ReachabilityAnalysis._convert_or_overapproximate(Ω0, Zonotope)
Zonotope{Float64, Vector{Float64}, Matrix{Float64}}([0.0, 0.0], [1.0 0.0; 0.0 0.5])
I didn't read in detail your second comment yet, but it seems in the correct direction as a way to fix this. We can add more heuristics to |
Let me mention that if the affine term is constant you can add it as a variable ("homogeneization") and use the algorithm ...
# no-op dynamics
A = zeros(2, 2)
prob = @ivp(x' = A * x, x(0) ∈ source_set)
alg = VREP(δ=0.1, approx_model=Forward(setops=:concrete))
# plots the initial source set
reach_set_all = set(sol_p)
plot(reach_set_all) There is a |
Generally speaking, it should be noted that the convex hull of two zonotopes is no longer a zonotope. But we have methods in LazySets to approximate that set by an outer zonotope (in this case, the result is probably exact). So I think we should add a new |
Is there a benefit to modeling the conservative time discretization initial set differently? For example, I notice julia> Ω0.Y.X.Y # descend down MinkowskiSum structs
MinkowskiSum{Float64, LinearMap{Float64, Singleton{Float64, Vector{Float64}}, Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}}(LinearMap{Float64, Singleton{Float64, Vector{Float64}}, Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}(
0.1 ⋅
⋅ 0.1, Singleton{Float64, Vector{Float64}}([0.0, 0.0])), Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}([0.0, 0.0], [0.0, 0.0])) Instead of nested
Indeed. I mentioned this in the original post :-)
I thought that the zonotope would also be propagated exactly under an affine transformation. I might be misunderstanding what discretizing the affine ODE entails, and perhaps the discretized dynamics are not affine.
Hm, so that is the purpose of
Perhaps it's just an edge case, but I am surprised that one needs to define the initial set differently, because I would expect the first step of propagating X0 through the discretized system to produce the reach-set in the time interval |
It depends what you want to compute. If you are interested in a sequence of reach-sets associated to time points, then the bloating factors are not necessary and you can just go with On the other hand, if you are interested in dense-time reachability, that is reach-sets that cover all possible behaviors in time intervals that are dense in R, ie. [0, delta], [delta, 2delta] and so on, then you need to use one of the conservative time discretization algorithms, such as
In the case of |
Describe the bug
using an initial set of type
Zonotope
causes an axis-aligned overapproximation when using GLM06To Reproduce
Steps to reproduce the behavior:
Expected behavior
A clear and concise description of what you expected to happen.
Screenshots
source_set
isreach_set_all
isAdditional context
the
solve
call does many things, including:which has quite a complicated
x0
:It
concretize
s into:and this gets overapproximated by the axis-aligned bounding box.
I would hope that
x0
could be a simpler set, and in this case it certainly could be, but I do not know in general.The complication is coming from:
ReachabilityAnalysis.jl/src/Discretization/Forward.jl
Line 130 in dbfacbd
It gets simplified quite well:
The convex hull of two zonotopes is not in general a zonotope, as illustrated here. In this case, the Zonotope overapproximation would be exact, however:
The text was updated successfully, but these errors were encountered: