Support function of lazy intersection via line search is not conservative

[schillic]

The intersection is flat (here in 1D: a point), but isempty detects that the set is non-empty.

using LazySets, LazySets.Approximations, Optim

julia> X = Interval(4., 5.);

julia> Y = HalfSpace([1.], 5.);

julia> Z = HalfSpace([-1.], -5.);

julia> isempty(X ∩ Y ∩ Z)
false

Now let us box-approximate in different orders. First step by step.

julia> overapproximate(X ∩ overapproximate(Y ∩ Z, Hyperrectangle), Hyperrectangle)
Hyperrectangle{Float64}([5.0], [0.0])  # good

julia> overapproximate(overapproximate(X ∩ Y, Hyperrectangle) ∩ Z, Hyperrectangle)
Hyperrectangle{Float64}([5.0], [0.0])  # good

Now a nested lazy intersection. The order matters!

julia> overapproximate(X ∩ (Y ∩ Z), Hyperrectangle)
Hyperrectangle{Float64}([5.0], [0.0])  # good

julia> overapproximate((X ∩ Y) ∩ Z, Hyperrectangle)
EmptySet{Float64}()  # bad

[schillic]

The reason for the difference is that we use different implementations of the support function.

julia> ρ([1.], (X ∩ (Y ∩ Z)))
5.0

X is an AbstractPolytope. Hence we use the implementation

return minimum([ρ(d, (Y ∩ Z) ∩ Hi) for Hi in constraints_list(X)])

On the other hand, Y and Z are HalfSpaces, and so we use line search in the second case. It now happens to be the case that the line search is not conservative.

julia> ρ([1.], ((X ∩ Y) ∩ Z))
4.999999999534339  # < 5.0 (bad)

julia> ρ([-1.], ((X ∩ Y) ∩ Z))
-5.000000000465661  # < -5.0 (good)

[mforets]

This issue is probably related to https://github.com/JuliaReach/LazySets.jl/issues/1348

[schillic]

The original issue in the OP is resolved, but the problem in the second post is still present. Below is another (simpler) example where the support-function approximation is strongly underestimating (-0.05 instead of 0). (In this example, the intersection is actually empty; but the sets are very close, so it is fine that we do not detect this.)

julia> B = Hyperrectangle([0.009231278330571413, 0.0], [0.009231221669425305, 1.0]);

julia> H = HalfSpace([1.0, 0.0], 0.0);

julia> vertices_list(intersection(B, H))
1-element Vector{Vector{Float64}}:
 [0.0, 1.0]

julia> ρ([1.0, 0.0], intersection(B, H))
0.0

# ---

julia> ρ([1.0, 0.0], B ∩ H)
-0.05666108859434239