search

JuliaReach / ReachabilityAnalysis.jl

Computing reachable states of dynamical systems in Julia
https://juliareach.github.io/ReachabilityAnalysis.jl/
MIT License
187 stars 17 forks source link

Seemingly unnecessary overapproximation given zonotope starting set #464

Open goretkin opened 3 years ago

goretkin commented 3 years ago

Describe the bug

using an initial set of type Zonotope causes an axis-aligned overapproximation when using GLM06

To Reproduce

Steps to reproduce the behavior:

  1. Version of the package used.

    (ReachabilityAnalysisZonotopeMWE) pkg> st
  Status `~/projects/ReachabilityAnalysisZonotopeMWE/Project.toml`
[b4f0291d] LazySets v1.44.2
[91a5bcdd] Plots v1.12.0
[1e97bd63] ReachabilityAnalysis v0.12.2 `dev/ReachabilityAnalysis`

  2. Minimum data to reproduce the error.

    
using LazySets: UnionSetArray, LineSegment, BallInf, CartesianProduct
using ReachabilityAnalysis
using Plots

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=:concrete)

) tspan = (0.0, 1.0) sol = solve(prob; tspan, alg)

project out time dimension

sol_p = LazySets.project(sol, [1, 2])

axis-aligned zonotope, but why?

@show sol_p.Xk[1].X

reach_set_all = set(sol_p) plot(reach_set_all)


**Expected behavior**

A clear and concise description of what you expected to happen.

**Screenshots**

`source_set` is
![image](https://user-images.githubusercontent.com/581043/116165315-9bf39000-a6c9-11eb-9d7a-42816e85f88d.png)

`reach_set_all` is
![image](https://user-images.githubusercontent.com/581043/116165351-b0d02380-a6c9-11eb-9915-e77a11d25f3c.png)

**Additional context**

the `solve` call does many things, including:
```julia
prob_norm = ReachabilityAnalysis._normalize(prob)
prob_discrete = ReachabilityAnalysis.discretize(prob_norm, alg.δ, alg.approx_model)

which has quite a complicated x0:

julia> ReachabilityAnalysis.initial_state(prob_discrete)
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}}}}(Zonotope{Float64, Vector{Float64}, Matrix{Float64}}([0.0, 0.0], [0.5 0.5; 0.0 0.5]), 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}}}(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}}}}(LinearMap{Float64, Zonotope{Float64, Vector{Float64}, Matrix{Float64}}, Float64, Matrix{Float64}}([1.0 0.0; 0.0 1.0], Zonotope{Float64, Vector{Float64}, Matrix{Float64}}([0.0, 0.0], [0.5 0.5; 0.0 0.5])), 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]))), Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}([0.0, 0.0], [0.0, 0.0])))

It concretizes into:

julia> x0_c = LazySets.concretize(prob_discrete.x0)
VPolygon{Float64, Vector{Float64}}([[-1.0, -0.5], [0.0, -0.5], [1.0, 0.5], [0.0, 0.5]])

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:

https://github.com/JuliaReach/ReachabilityAnalysis.jl/blob/dbfacbd03d4b65deb1b73d8d8089472ba6702fcd/src/Discretization/Forward.jl#L130

It gets simplified quite well:

julia> LazySets.concretize(prob_discrete.x0.X)
Zonotope{Float64, Vector{Float64}, Matrix{Float64}}([0.0, 0.0], [0.5 0.5; 0.0 0.5])

julia> LazySets.concretize(prob_discrete.x0.Y)
Zonotope{Float64, Vector{Float64}, Matrix{Float64}}([0.0, 0.0], [0.5 0.5; 0.0 0.5])

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:

using LazySets: ConvexHull, concretize, Zonotope
concretize_skip(s::LazySets.ConvexHull) = ConvexHull(concretize(s.X), concretize(s.Y))
julia> overapproximate(concretize_skip(prob_discrete.x0), Zonotope)
Zonotope{Float64, Vector{Float64}, Matrix{Float64}}([0.0, 0.0], [0.5 0.5; 0.0 0.5])
goretkin commented 3 years ago

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)
mforets commented 3 years ago

The reason that GLGM06 returns a box is that the algorithm (here) is not smart enough to realize that the set after the conservative time discretization is indeed a zonotope.

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 _convert_or_overapproximate, one of them could be to try concretize ....

mforets commented 3 years ago

Let me mention that if the affine term is constant you can add it as a variable ("homogeneization") and use the algorithm VREP, which propagates the sets using the vertex representation (VPolygon in 2D and VPolytope in higher dimension). In particular, with VREP the discretized set is propagated exactly (unlike with eg. zonotopes or hyperrectangles).

...

# 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 homogeneize to make the transformation automatically, but several methods are still missing to be added (there's a start in https://github.com/JuliaReach/ReachabilityAnalysis.jl/pull/372).

mforets commented 3 years ago

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 setops, which applies the rule that you proposed above. Maybe also to apply the overapprox rule for the CH with a zonotope.

goretkin commented 3 years ago

The reason that GLGM06 returns a box is that the algorithm (here) is not smart enough to realize that the set after the conservative time discretization is indeed a zonotope.

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 MinkowskiSum, perhaps a MinkowskiSumArray would give the compiler an easier time. But I suppose that's not the cause for this issue. The two children of the top-level ConvexHull do simplify to Zonotopes.

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

Indeed. I mentioned this in the original post :-)

In particular, with VREP the discretized set is propagated exactly (unlike with eg. zonotopes or hyperrectangles).

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.

# the following set is what we call the initial set after discretization: and it represents the reach-set
# in the time *interval* [0, δ]

Hm, so that is the purpose of https://github.com/JuliaReach/ReachabilityAnalysis.jl/blob/dbfacbd03d4b65deb1b73d8d8089472ba6702fcd/src/Discretization/Forward.jl#L130 ?

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 [0, δ], as opposed to needing to do this at the "zeroth" step.

mforets commented 3 years ago

I might be misunderstanding what discretizing the affine ODE entails,

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 NoBloating() discretization. In some sense, this result is analogue to traditional numerical integration, but instead of obtaining a sequence of vectors in R^n, you obtain a sequence of sets in R^n. (If you start with a singleton, you obtain a sequence of singletons, try the ORBIT algorithm).

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 Forward() or the other methods. Once this initial reach-set has been computed, we propagate it according to the system's dynamics. Perhaps that should clarify why the sol[1] of flowpipes is not X0: if you do tspan(sol[1]), it is [0, delta].

Is there a benefit to modeling the conservative time discretization initial set differently?

In the case of Forward(), the default return type (setops=:lazy) is good because in that way we don't make any ad-hoc overapproximation: the algorithm such as GLGM06 or other algorithm can do whatever it wants with the set, maybe it needs to overapproximate it to propagate it using the representation of that algorithm. Let me remark that there are methods that can work with the lazy set as given, i mentioned VREP above (the problem is that it doesn't scale so it is only relevant in low dimensions); the other method is LGG09 which does scale and it uses support functions to propagate the set.