Merge pull request #721 from JuliaReach/mforets/720

#720 - add upper bounds

docs/src/lib/approximations.md

@@ -66,3 +66,9 @@ BoxDiagDirections
 ```
 
 See also `overapproximate(X::LazySet, dir::AbstractDirections)::HPolytope`.
+
+## Upper bounds
+
+```@docs
+ρ_upper_bound
+```

src/Approximations/Approximations.jl

@@ -37,5 +37,6 @@ include("box_approximations.jl")
 include("template_directions.jl")
 include("overapproximate.jl")
 include("decompositions.jl")
+include("upper_bounds.jl")
 
 end # module

src/Approximations/upper_bounds.jl

@@ -0,0 +1,113 @@
+export ρ_upper_bound
+
+"""
+    ρ_upper_bound(d::AbstractVector{N},
+                  X::LazySet{N};
+                  kwargs...) where {N<:Real}
+
+Return an upper bound of the support function of a given set.
+
+### Input
+
+- `d` -- direction
+- `X` -- set
+
+### Output
+
+An upper bound of the support function of the given set.
+
+### Algorithm
+
+The default implementation of `ρ_upper_bound` is the exact `ρ(d, X)`.
+"""
+function ρ_upper_bound(d::AbstractVector{N},
+                       X::LazySet{N};
+                       kwargs...) where {N<:Real}
+    return ρ(d, X)
+end
+
+"""
+    ρ_upper_bound(d::AbstractVector{N},
+                  cap::Intersection{N}; kwargs...) where {N<:Real}
+
+Return an upper bound of the support function of the intersection of two sets.
+
+### Input
+
+- `d`   -- direction
+- `cap` -- intersection
+
+### Output
+
+An upper bound of the support function of the given intersection.
+
+### Algorithm
+
+The support function of an intersection of ``X`` and ``Y`` is upper bounded by
+the minimum of the support functions of ``X`` and ``Y``.
+"""
+function ρ_upper_bound(d::AbstractVector{N},
+                       cap::Intersection{N}; kwargs...) where {N<:Real}
+    return min(ρ_upper_bound(d, cap.X; kwargs...), ρ_upper_bound(d, cap.Y; kwargs...))
+end
+
+"""
+    ρ_upper_bound(d::AbstractVector{N},
+                  cap::Intersection{N, <:LazySet{N}, <:AbstractPolytope{N}};
+                  kwargs...) where {N<:Real}
+
+Return an upper bound of the intersection between a compact set and a
+polytope along a given direction.
+
+### Input
+
+- `d`      -- direction
+- `cap`    -- intersection of a compact set and a polytope
+- `kwargs` -- additional arguments that are passed to the support function algorithm
+
+### Output
+
+An upper bound of the support function of the given intersection.
+
+### Algorithm
+
+The idea is to solve the univariate optimization problem `ρ(di, X ∩ Hi)` for each
+half-space in the set `P` and then take the minimum. This gives an overapproximation
+of the exact support function.
+
+This algorithm is inspired from [G. Frehse, R. Ray. Flowpipe-Guard Intersection
+for Reachability Computations with Support
+Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809).
+
+### Notes
+
+This method relies on having available the `constraints_list` of the polytope
+`P`.
+
+This method of overapproximation can return a non-empty set even if the original
+intersection is empty.
+"""
+function ρ_upper_bound(d::AbstractVector{N},
+                       cap::Intersection{N, <:LazySet{N}, <:AbstractPolytope{N}};
+                       kwargs...) where {N<:Real}
+
+    X = cap.X    # compact set
+    P = cap.Y    # polytope
+    return minimum([ρ(d, X ∩ Hi, kwargs...) for Hi in constraints_list(P)])
+end
+
+# disambiguation
+function ρ_upper_bound(d::AbstractVector{N},
+                       cap::Intersection{N, <:AbstractPolytope{N}, <:AbstractPolytope{N}};
+                       kwargs...) where {N<:Real}
+    X = cap.X    # compact set
+    P = cap.Y    # polytope
+    return minimum([ρ(d, X ∩ Hi, kwargs...) for Hi in constraints_list(P)])
+end
+
+# symmetric function
+function ρ_upper_bound(d::AbstractVector{N},
+                       cap::Intersection{N, <:AbstractPolytope{N}, <:LazySet{N}};
+                       kwargs...) where {N<:Real}
+    return ρ_upper_bound(d, cap.Y ∩ cap.X; kwargs...)
+end

test/runtests.jl

@@ -95,6 +95,7 @@ if test_suite_basic
     @time @testset "LazySets.Approximations.ballinf_approximation" begin include("unit_ballinf_approximation.jl") end
     @time @testset "LazySets.Approximations.radiusdiameter" begin include("unit_radiusdiameter.jl") end
     @time @testset "LazySets.Approximations.decompose" begin include("unit_decompose.jl") end
+    @time @testset "LazySets.Approximations.upper_bounds" begin include("unit_upper_bounds.jl") end
 
     # ========================
     # Testing method ambiguity

test/unit_upper_bounds.jl

@@ -0,0 +1,10 @@
+import LazySets.Approximations.ρ_upper_bound
+
+X = Hyperrectangle(low=[0.0, 0.0], high=[1.0, 1.0])
+d = [1.0, 2.0]
+@test ρ(d, X) == ρ_upper_bound(d, X)
+
+using Optim
+@test ρ(d, X) == ρ_upper_bound(d, X ∩ X)
+M = [1.0 1.0; 1.0 1.0]
+@test isapprox(ρ_upper_bound(d, X ∩ (M * X)), 3.0, atol=1e-6)