diff --git a/docs/src/lib/approximations.md b/docs/src/lib/approximations.md
index 6273600771..85093af9bb 100644
--- a/docs/src/lib/approximations.md
+++ b/docs/src/lib/approximations.md
@@ -66,9 +66,3 @@ BoxDiagDirections
 ```
 
 See also `overapproximate(X::LazySet, dir::AbstractDirections)::HPolytope`.
-
-## Upper bounds
-
-```@docs
-ρ_upper_bound
-```
diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl
index 85f3fdaddd..fa0f70c51a 100644
--- a/src/Approximations/Approximations.jl
+++ b/src/Approximations/Approximations.jl
@@ -37,6 +37,5 @@ include("box_approximations.jl")
 include("template_directions.jl")
 include("overapproximate.jl")
 include("decompositions.jl")
-include("upper_bounds.jl")
 
 end # module
diff --git a/src/Approximations/box_approximations.jl b/src/Approximations/box_approximations.jl
index bf9f2bd3e3..920872c6f7 100644
--- a/src/Approximations/box_approximations.jl
+++ b/src/Approximations/box_approximations.jl
@@ -3,15 +3,13 @@
 # ===================================
 
 """
-    box_approximation(S::LazySet; upper_bound::Bool=false)::Hyperrectangle
+    box_approximation(S::LazySet)::Hyperrectangle
 
 Overapproximate a convex set by a tight hyperrectangle.
 
 ### Input
 
 - `S`           -- convex set
-- `upper_bound` -- (optional, default: `false`) use overapproximation in support
-                   function computation?
 
 ### Output
 
@@ -24,20 +22,19 @@ of the given set in the canonical directions, and the lengths of the sides can
 be recovered from the distance among support functions in the same directions.
 """
 function box_approximation(S::LazySet{N};
-                           upper_bound::Bool=false
                           )::Union{Hyperrectangle{N}, EmptySet{N}} where N<:Real
-    (c, r) = box_approximation_helper(S; upper_bound=upper_bound)
-    if upper_bound && r[1] < 0
+    (c, r) = box_approximation_helper(S)
+    if r[1] < 0
         return EmptySet{N}()
     end
     return Hyperrectangle(c, r)
 end
 
 # special case: Hyperrectangle
-box_approximation(S::Hyperrectangle; upper_bound::Bool=false) = S
+box_approximation(S::Hyperrectangle) = S
 
 # special case: other rectangle
-box_approximation(S::AbstractHyperrectangle; upper_bound::Bool=false) =
+box_approximation(S::AbstractHyperrectangle) =
     Hyperrectangle(center(S), radius_hyperrectangle(S))
 
 """
@@ -48,7 +45,7 @@ Alias for `box_approximation`.
 interval_hull = box_approximation
 
 """
-    box_approximation_symmetric(S::LazySet{N}; upper_bound::Bool=false
+    box_approximation_symmetric(S::LazySet{N}
                                )::Union{Hyperrectangle{N}, EmptySet{N}}
                                 where {N<:Real}
 
@@ -57,8 +54,6 @@ Overapproximate a convex set by a tight hyperrectangle centered in the origin.
 ### Input
 
 - `S`           -- convex set
-- `upper_bound` -- (optional, default: `false`) use overapproximation in support
-                   function computation?
 
 ### Output
 
@@ -70,11 +65,10 @@ The center of the box is the origin, and the radius is obtained by computing the
 maximum value of the support function evaluated at the canonical directions.
 """
 function box_approximation_symmetric(S::LazySet{N};
-                                     upper_bound::Bool=false
                                     )::Union{Hyperrectangle{N},
                                              EmptySet{N}} where {N<:Real}
-    (c, r) = box_approximation_helper(S; upper_bound=upper_bound)
-    if upper_bound && r[1] < 0
+    (c, r) = box_approximation_helper(S)
+    if r[1] < 0
         return EmptySet{N}()
     end
     return Hyperrectangle(zeros(N, length(c)), abs.(c) .+ r)
@@ -88,7 +82,7 @@ Alias for `box_approximation_symmetric`.
 symmetric_interval_hull = box_approximation_symmetric
 
 """
-    box_approximation_helper(S::LazySet{N}; upper_bound::Bool=false
+    box_approximation_helper(S::LazySet{N};
                             ) where {N<:Real}
 
 Common code of `box_approximation` and `box_approximation_symmetric`.
@@ -96,8 +90,6 @@ Common code of `box_approximation` and `box_approximation_symmetric`.
 ### Input
 
 - `S`           -- convex set
-- `upper_bound` -- (optional, default: `false`) use overapproximation in support
-                   function computation?
 
 ### Output
 
@@ -115,24 +107,22 @@ The lengths of the sides can be recovered from the distance among support
 functions in the same directions.
 """
 @inline function box_approximation_helper(S::LazySet{N};
-                                          upper_bound::Bool=false
                                          ) where {N<:Real}
     zero_N = zero(N)
     one_N = one(N)
-    ρ_rec = upper_bound ? ρ_upper_bound : ρ
     n = dim(S)
     c = Vector{N}(undef, n)
     r = Vector{N}(undef, n)
     d = zeros(N, n)
     @inbounds for i in 1:n
         d[i] = one_N
-        htop = ρ_rec(d, S)
+        htop = ρ(d, S)
         d[i] = -one_N
-        hbottom = -ρ_rec(d, S)
+        hbottom = -ρ(d, S)
         d[i] = zero_N
         c[i] = (htop + hbottom) / 2
         r[i] = (htop - hbottom) / 2
-        if upper_bound && r[i] < 0
+        if r[i] < 0
             # contradicting bounds => set is empty
             # terminate with first radius entry being negative
             r[1] = r[i]
@@ -143,7 +133,7 @@ functions in the same directions.
 end
 
 """
-    ballinf_approximation(S::LazySet{N}; upper_bound::Bool=false
+    ballinf_approximation(S::LazySet{N};
                          )::BallInf{N} where {N<:Real}
 
 Overapproximate a convex set by a tight ball in the infinity norm.
@@ -151,8 +141,6 @@ Overapproximate a convex set by a tight ball in the infinity norm.
 ### Input
 
 - `S`           -- convex set
-- `upper_bound` -- (optional, default: `false`) use overapproximation in support
-                   function computation?
 
 ### Output
 
@@ -164,26 +152,24 @@ The center and radius of the box are obtained by evaluating the support function
 of the given convex set along the canonical directions.
 """
 function ballinf_approximation(S::LazySet{N};
-                               upper_bound::Bool=false
                               )::Union{BallInf{N}, EmptySet{N}} where {N<:Real}
     zero_N = zero(N)
     one_N = one(N)
-    ρ_rec = upper_bound ? ρ_upper_bound : ρ
     n = dim(S)
     c = Vector{N}(undef, n)
     r = zero_N
     d = zeros(N, n)
     @inbounds for i in 1:n
         d[i] = one_N
-        htop = ρ_rec(d, S)
+        htop = ρ(d, S)
         d[i] = -one_N
-        hbottom = -ρ_rec(d, S)
+        hbottom = -ρ(d, S)
         d[i] = zero_N
         c[i] = (htop + hbottom) / 2
         rcur = (htop - hbottom) / 2
         if (rcur > r)
             r = rcur
-        elseif upper_bound && rcur < 0
+        elseif rcur < 0
             # contradicting bounds => set is empty
             return EmptySet{N}()
         end
diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl
index 83680ed173..02049a5051 100644
--- a/src/Approximations/overapproximate.jl
+++ b/src/Approximations/overapproximate.jl
@@ -30,7 +30,6 @@ Otherwise the result is an ε-close approximation as a polygon.
 - `S`           -- convex set, assumed to be two-dimensional
 - `HPolygon`    -- type for dispatch
 - `ε`           -- (optional, default: `Inf`) error bound
-- `upper_bound` -- (optional, default: `false`) currently ignored
 
 ### Output
 
@@ -38,8 +37,7 @@ A polygon in constraint representation.
 """
 function overapproximate(S::LazySet{N},
                          ::Type{<:HPolygon},
-                         ε::Real=Inf;
-                         upper_bound::Bool=false
+                         ε::Real=Inf
                         )::HPolygon where {N<:Real}
     @assert dim(S) == 2
     if ε == Inf
@@ -59,15 +57,14 @@ function overapproximate(S::LazySet{N},
 end
 
 """
-    overapproximate(S::LazySet, ε::Real; [upper_bound]::Bool=false)::HPolygon
+    overapproximate(S::LazySet, ε::Real)::HPolygon
 
-Alias for `overapproximate(S, HPolygon, ε, upper_bound=upper_bound)`.
+Alias for `overapproximate(S, HPolygon, ε)`.
 """
 function overapproximate(S::LazySet,
                          ε::Real;
-                         upper_bound::Bool=false
                         )::HPolygon
-    return overapproximate(S, HPolygon, ε; upper_bound=upper_bound)
+    return overapproximate(S, HPolygon, ε)
 end
 
 """
@@ -80,8 +77,6 @@ Return an approximation of a given set as a hyperrectangle.
 
 - `S`              -- set
 - `Hyperrectangle` -- type for dispatch
-- `upper_bound`    -- (optional, default: `false`) use overapproximation in
-                      support function computation?
 
 ### Output
 
@@ -89,9 +84,8 @@ A hyperrectangle.
 """
 function overapproximate(S::LazySet,
                          ::Type{<:Hyperrectangle};
-                         upper_bound::Bool=false
                         )::Union{Hyperrectangle, EmptySet}
-    box_approximation(S; upper_bound=upper_bound)
+    box_approximation(S)
 end
 
 """
@@ -100,14 +94,12 @@ end
 Alias for `overapproximate(S, Hyperrectangle)`.
 """
 overapproximate(S::LazySet;
-                upper_bound::Bool=false
                )::Union{Hyperrectangle, EmptySet} =
-    overapproximate(S, Hyperrectangle; upper_bound=upper_bound)
+    overapproximate(S, Hyperrectangle)
 
 """
     overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}},
-                    ::Type{<:Zonotope};
-                    [upper_bound]::Bool=false)::Zonotope where {N<:Real}
+                    ::Type{<:Zonotope})::Zonotope where {N<:Real}
 
 Overapproximate the convex hull of two zonotopes.
 
@@ -115,7 +107,6 @@ Overapproximate the convex hull of two zonotopes.
 
 - `S`           -- convex hull of two zonotopes of the same order
 - `Zonotope`    -- type for dispatch
-- `upper_bound` -- (optional, default: `false`) currently ignored
 
 ### Algorithm
 
@@ -134,8 +125,7 @@ zonotope, which is in general expensive in high dimensions. This is further inve
 in: *Zonotopes as bounding volumes, L. J. Guibas et al, Proc. of Symposium on Discrete Algorithms, pp. 803-812*.
 """
 function overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}},
-                         ::Type{<:Zonotope};
-                         upper_bound::Bool=false)::Zonotope where {N<:Real}
+                         ::Type{<:Zonotope})::Zonotope where {N<:Real}
     Z1, Z2 = S.X, S.Y
     @assert order(Z1) == order(Z2)
     center = (Z1.center+Z2.center)/2
@@ -145,8 +135,7 @@ end
 
 """
     overapproximate(X::LazySet,
-                    dir::AbstractDirections;
-                    [upper_bound]::Bool=false
+                    dir::AbstractDirections
                    )::HPolytope
 
 Overapproximating a set with template directions.
@@ -155,8 +144,6 @@ Overapproximating a set with template directions.
 
 - `X`           -- set
 - `dir`         -- direction representation
-- `upper_bound` -- (optional, default: `false`) use overapproximation in support
-                   function computation?
 
 ### Output
 
@@ -164,22 +151,19 @@ A `HPolytope` overapproximating the set `X` with the directions from `dir`.
 """
 function overapproximate(X::LazySet{N},
                          dir::AbstractDirections{N};
-                         upper_bound::Bool=false
                         )::HPolytope{N} where N
-    ρ_rec = upper_bound ? ρ_upper_bound : ρ
     halfspaces = Vector{LinearConstraint{N}}()
     sizehint!(halfspaces, length(dir))
     H = HPolytope(halfspaces)
     for d in dir
-        addconstraint!(H, LinearConstraint(d, ρ_rec(d, X)))
+        addconstraint!(H, LinearConstraint(d, ρ(d, X)))
     end
     return H
 end
 
 """
     overapproximate(S::LazySet{N},
-                    ::Type{Interval};
-                    [upper_bound]::Bool=false
+                    ::Type{Interval}
                    ) where {N<:Real}
 
 Return the overapproximation of a real unidimensional set with an interval.
@@ -188,15 +172,13 @@ Return the overapproximation of a real unidimensional set with an interval.
 
 - `S`           -- one-dimensional set
 - `Interval`    -- type for dispatch
-- `upper_bound` -- (optional, default: `false`) currently ignored
 
 ### Output
 
 An interval.
 """
 function overapproximate(S::LazySet{N},
-                         ::Type{Interval};
-                         upper_bound::Bool=false
+                         ::Type{Interval}
                         ) where {N<:Real}
     @assert dim(S) == 1
     lo = σ([-one(N)], S)[1]
@@ -207,7 +189,6 @@ end
 """
     overapproximate(cap::Intersection{N, <:LazySet, S},
                     dir::AbstractDirections{N};
-                    [upper_bound]::Bool=false,
                     kwargs...
                    ) where {N<:Real, S<:AbstractPolytope{N}}
 
@@ -218,8 +199,6 @@ polytope given a set of template directions.
 
 - `cap`         -- intersection of a compact set and a polytope
 - `dir`         -- template directions
-- `upper_bound` -- (optional, default: `false`) use overapproximation in support
-                   function computation?
 - `kwargs`      -- additional arguments that are passed to the support function
                    algorithm
 
@@ -253,7 +232,6 @@ intersection is empty.
 """
 function overapproximate(cap::Intersection{N, <:LazySet, S},
                          dir::AbstractDirections{N};
-                         upper_bound::Bool=false,
                          kwargs...
                         ) where {N<:Real, S<:AbstractPolytope{N}}
 
@@ -263,12 +241,11 @@ function overapproximate(cap::Intersection{N, <:LazySet, S},
     Hi = constraints_list(P)
     m = length(Hi)
     Q = HPolytope{N}()
-    ρ_rec = upper_bound ? ρ_upper_bound : ρ
 
     for di in dir
-        ρ_X_Hi_min = ρ_rec(di, X ∩ Hi[1], kwargs...)
+        ρ_X_Hi_min = ρ(di, X ∩ Hi[1], kwargs...)
         for i in 2:m
-            ρ_X_Hi = ρ_rec(di, X ∩ Hi[i], kwargs...)
+            ρ_X_Hi = ρ(di, X ∩ Hi[i], kwargs...)
             if ρ_X_Hi < ρ_X_Hi_min
                 ρ_X_Hi_min = ρ_X_Hi
             end
diff --git a/src/Approximations/upper_bounds.jl b/src/Approximations/upper_bounds.jl
index 9d42bbb3af..139597f9cb 100644
--- a/src/Approximations/upper_bounds.jl
+++ b/src/Approximations/upper_bounds.jl
@@ -1,243 +1,2 @@
-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},
-                  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` for lazy operations is to pass
-all keyword arguments to `ρ(d, X; kwargs...)` if that method exists.
-"""
-function ρ_upper_bound(d::AbstractVector{N},
-                       X::Union{LinearMap{N}};
-                       kwargs...) where {N<:Real}
-    return ρ(d, X, upper_bound=true, kwargs...)
-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
-
-"""
-    ρ_upper_bound(d::AbstractVector{N},
-                  cap::Intersection{N, <:LazySet{N},
-                                       <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
-                  [algorithm]::String="line_search",
-                  [check_intersection]::Bool=true,
-                  [kwargs...]) where N<:Real
-
-Return an upper bound of the support function of the intersection of a compact
-set and a half-space or a hyperplane or a line in a given direction.
-
-### Input
-
-- `d`         -- direction
-- `cap`       -- lazy intersection of a compact set and a half-space/hyperplane/
-                 line
-- `algorithm` -- (optional, default: `"line_search"`): the algorithm to
-                 calculate the support function; valid options are:
-
-    * `"line_search"` -- solve the associated univariate optimization problem
-                         using a line search method (either Brent or the
-                         Golden Section method)
-    * `"projection"`  -- only valid for intersection with a hyperplane;
-                         evaluates the support function by reducing the problem
-                         to the 2D intersection of a rank 2 linear
-                         transformation of the given compact set in the plane
-                         generated by the given direction `d` and the
-                         hyperplane's normal vector `n`
-
-- `check_intersection` -- (optional, default: `true`) if `true`, check if the
-                          intersection is empty before actually calculating the
-                          support function
-
-### Output
-
-The scalar value of the support function of the set `cap` in the given
-direction.
-
-### Notes
-
-It is assumed that the set `cap.X` is compact.
-
-The `check_intersection` flag can be useful if it is known in advance that the
-intersection is non-empty.
-
-Any additional number of arguments to the algorithm backend can be passed as
-keyword arguments.
-
-### Algorithm
-
-We fall back to the implementation of `ρ` with the same arguments.
-"""
-function ρ_upper_bound(d::AbstractVector{N},
-                       cap::Intersection{N,
-                                         <:LazySet{N},
-                                         <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
-                       algorithm::String="line_search",
-                       check_intersection::Bool=true,
-                       kwargs...
-                      ) where N<:Real
-    return ρ(d, cap; algorithm=algorithm, check_intersection=check_intersection,
-             upper_bound=true, kwargs...)
-end
-
-# symmetric case
-function ρ_upper_bound(d::AbstractVector{N},
-                       cap::Intersection{N,
-                                         <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
-                                         <:LazySet{N}};
-                       algorithm::String="line_search",
-                       check_intersection::Bool=true,
-                       kwargs...
-                      ) where N<:Real
-    return ρ(d, cap; algorithm=algorithm, check_intersection=check_intersection,
-             upper_bound=true, kwargs...)
-end
-
-# disambiguation
-function ρ_upper_bound(d::AbstractVector{N},
-                       cap::Intersection{N,
-                                         <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
-                                         <:AbstractPolytope{N}};
-                       algorithm::String="line_search",
-                       check_intersection::Bool=true,
-                       kwargs...
-                      ) where N<:Real
-    return ρ(d, cap; algorithm=algorithm, check_intersection=check_intersection,
-             upper_bound=true, kwargs...)
-end
-
-# disambiguation symmetric case
-function ρ_upper_bound(d::AbstractVector{N},
-                       cap::Intersection{N,
-                                         <:AbstractPolytope{N},
-                                         <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
-                       algorithm::String="line_search",
-                       check_intersection::Bool=true,
-                       kwargs...
-                      ) where N<:Real
-    return ρ(d, cap; algorithm=algorithm, check_intersection=check_intersection,
-             upper_bound=true, kwargs...)
-end
diff --git a/src/Intersection.jl b/src/Intersection.jl
index a29c5acc0a..3bb0d39620 100644
--- a/src/Intersection.jl
+++ b/src/Intersection.jl
@@ -110,38 +110,230 @@ direction.
 The support vector in the given direction.
 """
 function σ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}
-    # TODO document behavior if the direction has norm zero
-    # TODO error message if the intersection is empty!
-    # TODO implement
-    error("the exact support vector of an intersection is not implemented yet")
+    error("the exact support vector of an intersection is not implemented")
 end
 
 """
-    ρ(d::AbstractVector{N}, cap::Intersection{N};
-      upper_bound=false, kwargs...) where {N<:Real}
+    ρ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}
 
-Return the support function of the intersection of two convex sets in a given
-direction.
+Return an upper bound on the support function of the intersection of two convex
+sets in a given direction.
 
 ### Input
 
-- `d`           -- direction
-- `cap`         -- intersection of two convex sets
-- `upper_bound` -- (optional, default: `false`) if `false`, compute the support
-                   function exactly; otherwise use an overapproximative algorithm
+- `d`    -- direction
+- `cap`  -- intersection of two convex sets
 
 ### Output
 
-The support function in the given direction.
+An uper bound on the support function in the given direction.
+
+### 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 ρ(d::AbstractVector{N}, cap::Intersection{N};
-           upper_bound=false, kwargs...) where {N<:Real}
-    if upper_bound
-        return LazySets.Approximations.ρ_upper_bound(d, cap; kwargs...)
+function ρ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real}
+    return min(ρ(d, cap.X), ρ(d, cap.Y))
+end
+
+function ρ_helper(d::AbstractVector{N},
+                  cap::Intersection{N},
+                  algorithm::String,
+                  check_intersection::Bool;
+                  kwargs...) where N<:Real
+    X = cap.X # compact set
+    H = cap.Y # halfspace or hyperplane or line
+
+    # if the intersection is empty => stop
+    if check_intersection && is_intersection_empty(X, H, false)
+        error("the intersection is empty")
+    end
+
+    if algorithm == "line_search"
+        @assert isdefined(Main, :Optim) "the algorithm $algorithm needs " *
+                                        "the package 'Optim' to be loaded"
+        (s, _) = _line_search(d, X, H; kwargs...)
+    elseif algorithm == "projection"
+        @assert H isa Hyperplane "the algorithm $algorithm cannot be used with a
+                                  $(typeof(H)); it only works with hyperplanes"
+        s = _projection(d, X, H; kwargs...)
     else
-        error("the exact support function of an intersection is not implemented; " *
-              "try using `upper_bound=true`")
+        error("algorithm $(algorithm) unknown")
     end
+    return s
+end
+
+"""
+    ρ(d::AbstractVector{N},
+      cap::Intersection{N,
+                        <:LazySet{N},
+                        <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
+      [algorithm]::String="line_search",
+      [check_intersection]::Bool=true,
+      [kwargs...]) where N<:Real
+
+Return the support function of the intersection of a compact set and a
+half-space/hyperplane/line in a given direction.
+
+### Input
+
+- `d`         -- direction
+- `cap`       -- lazy intersection of a compact set and a half-space/hyperplane/
+                 line
+- `algorithm` -- (optional, default: `"line_search"`): the algorithm to
+                 calculate the support function; valid options are:
+
+    * `"line_search"` -- solve the associated univariate optimization problem
+                         using a line search method (either Brent or the
+                         Golden Section method)
+    * `"projection"`  -- only valid for intersection with a hyperplane;
+                         evaluates the support function by reducing the problem
+                         to the 2D intersection of a rank 2 linear
+                         transformation of the given compact set in the plane
+                         generated by the given direction `d` and the
+                         hyperplane's normal vector `n`
+
+- `check_intersection` -- (optional, default: `true`) if `true`, check if the
+                          intersection is empty before actually calculating the
+                          support function
+
+### Output
+
+The scalar value of the support function of the set `cap` in the given
+direction.
+
+### Notes
+
+It is assumed that the set `cap.X` is compact.
+
+The `check_intersection` flag can be useful if it is known in advance that the
+intersection is non-empty.
+
+Any additional number of arguments to the algorithm backend can be passed as
+keyword arguments.
+
+### Algorithm
+
+The algorithms are based on solving the associated optimization problem
+
+```math
+\\min_\\{ λ ∈ D_h \\} ρ(ℓ - λa, X) + λb.
+```
+where ``D_h = \\{ λ : λ ≥ 0 \\}`` if ``H`` is a half-space or
+``D_h = \\{ λ : λ ∈ \\mathbb{R} \\}`` if ``H`` is a hyperplane.
+
+For additional information we refer to:
+
+- [G. Frehse, R. Ray. Flowpipe-Guard Intersection for Reachability Computations with
+  Support Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809).
+- [C. Le Guernic. Reachability Analysis of Hybrid Systems with Linear Continuous
+  Dynamics, PhD thesis](https://tel.archives-ouvertes.fr/tel-00422569v2).
+- [T. Rockafellar, R. Wets.
+  Variational Analysis](https://www.springer.com/us/book/9783540627722).
+"""
+function ρ(d::AbstractVector{N},
+           cap::Intersection{N,
+                             <:LazySet{N},
+                             <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
+           algorithm::String="line_search",
+           check_intersection::Bool=true,
+           kwargs...) where N<:Real
+    return ρ_helper(d, cap, algorithm, check_intersection, kwargs...)
+end
+
+# symmetric method
+function ρ(ℓ::AbstractVector{N},
+           cap::Intersection{N,
+                             <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
+                             <:LazySet{N}};
+           algorithm::String="line_search",
+           check_intersection::Bool=true,
+           kwargs...) where N<:Real
+    return ρ_helper(d, cap.Y ∩ cap.X, algorithm, check_intersection, kwargs...)
+end
+
+"""
+    ρ(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 ρ(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
+
+# symmetric method
+function ρ(d::AbstractVector{N},
+           cap::Intersection{N, <:AbstractPolytope{N}, <:LazySet{N}};
+           kwargs...) where {N<:Real}
+    return ρ(d, cap.Y ∩ cap.X; kwargs...)
+end
+
+# disambiguation
+function ρ(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
+
+# disambiguation
+function ρ(d::AbstractVector{N},
+           cap::Intersection{N,
+                             <:AbstractPolytope{N},
+                             <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
+           algorithm::String="line_search",
+           check_intersection::Bool=true,
+           kwargs...) where N<:Real
+    return ρ_helper(d, cap, algorithm, check_intersection; kwargs...)
+end
+
+# symmetric method
+function ρ(d::AbstractVector{N},
+           cap::Intersection{N,
+                             <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
+                             <:AbstractPolytope{N}};
+           algorithm::String="line_search",
+           check_intersection::Bool=true,
+           kwargs...) where N<:Real
+    return ρ_helper(d, cap.Y ∩ cap.X, algorithm, check_intersection; kwargs...)
 end
 
 """
@@ -316,115 +508,9 @@ function ∈(x::AbstractVector{N}, ia::IntersectionArray{N})::Bool where {N<:Rea
     return true
 end
 
-"""
-    ρ(d::AbstractVector{N},
-      cap::Intersection{N,
-                        <:LazySet{N},
-                        <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
-      [algorithm]::String="line_search",
-      [check_intersection]::Bool=true,
-      [kwargs...]) where N<:Real
-
-Return the support function of the intersection of a compact set and a
-half-space/hyperplane/line in a given direction.
-
-### Input
-
-- `d`         -- direction
-- `cap`       -- lazy intersection of a compact set and a half-space/hyperplane/
-                 line
-- `algorithm` -- (optional, default: `"line_search"`): the algorithm to
-                 calculate the support function; valid options are:
-
-    * `"line_search"` -- solve the associated univariate optimization problem
-                         using a line search method (either Brent or the
-                         Golden Section method)
-    * `"projection"`  -- only valid for intersection with a hyperplane;
-                         evaluates the support function by reducing the problem
-                         to the 2D intersection of a rank 2 linear
-                         transformation of the given compact set in the plane
-                         generated by the given direction `d` and the
-                         hyperplane's normal vector `n`
-
-- `check_intersection` -- (optional, default: `true`) if `true`, check if the
-                          intersection is empty before actually calculating the
-                          support function
-
-### Output
-
-The scalar value of the support function of the set `cap` in the given
-direction.
-
-### Notes
-
-It is assumed that the set `cap.X` is compact.
-
-The `check_intersection` flag can be useful if it is known in advance that the
-intersection is non-empty.
-
-Any additional number of arguments to the algorithm backend can be passed as
-keyword arguments.
-
-### Algorithm
-
-The algorithms are based on solving the associated optimization problem
-
-```math
-\\min_\\{ λ ∈ D_h \\} ρ(ℓ - λa, X) + λb.
-```
-where ``D_h = \\{ λ : λ ≥ 0 \\}`` if ``H`` is a half-space or
-``D_h = \\{ λ : λ ∈ \\mathbb{R} \\}`` if ``H`` is a hyperplane.
-
-For additional information we refer to:
-
-- [G. Frehse, R. Ray. Flowpipe-Guard Intersection for Reachability Computations with
-  Support Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809).
-- [C. Le Guernic. Reachability Analysis of Hybrid Systems with Linear Continuous
-  Dynamics, PhD thesis](https://tel.archives-ouvertes.fr/tel-00422569v2).
-- [T. Rockafellar, R. Wets.
-  Variational Analysis](https://www.springer.com/us/book/9783540627722).
-"""
-function ρ(d::AbstractVector{N},
-           cap::Intersection{N,
-                             <:LazySet{N},
-                             <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
-           algorithm::String="line_search",
-           check_intersection::Bool=true,
-           upper_bound::Bool=false,
-           kwargs...) where N<:Real
-
-    X = cap.X    # compact set
-    H = cap.Y    # halfspace or hyperplane
-
-    # if the intersection is empty => stop
-    if check_intersection &&
-            is_intersection_empty(X, H, false; upper_bound=upper_bound)
-        error("the intersection is empty")
-    end
-    
-    if algorithm == "line_search"
-        @assert isdefined(Main, :Optim) "the algorithm $algorithm needs " *
-                                        "the package 'Optim' to be loaded"
-        (s, _) = _line_search(d, X, H; upper_bound=upper_bound, kwargs...)
-    elseif algorithm == "projection"
-        @assert H isa Hyperplane "the algorithm $algorithm cannot be used with a
-                                  $(typeof(H)); it only works with hyperplanes"
-        s = _projection(d, X, H; upper_bound=upper_bound, kwargs...)
-    else
-        error("algorithm $(algorithm) unknown")
-    end
-    return s
-end
-
-# Symmetric case
-ρ(ℓ::AbstractVector{N},
-  cap::Intersection{N,
-                    <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
-                    <:LazySet{N}};
-  algorithm::String="line_search", check_intersection::Bool=true,
-  kwargs...) where N<:Real =
-    ρ(ℓ, cap.Y ∩ cap.X; algorithm=algorithm,
-      check_intersection=check_intersection, kwargs...)
+# ==================================
+# Algorithms for lazy intersection
+# ==================================
 
 function load_optim_intersection()
 return quote
@@ -432,7 +518,7 @@ return quote
 import Optim
 
 """
-    _line_search(ℓ, X, H; [upper_bound]=false, [kwargs...])
+    _line_search(ℓ, X, H; [kwargs...])
 
 Given a compact and convex set ``X`` and a halfspace ``H = \\{x: a^T x ≤ b \\}``
 or a hyperplane ``H = \\{x: a^T x = b \\}``, calculate:
@@ -448,9 +534,6 @@ where ``D_h = \\{ λ : λ ≥ 0 \\}`` if ``H`` is a half-space or
 - `ℓ`           -- direction
 - `X`           -- set
 - `H`           -- halfspace or hyperplane
-- `upper_bound` -- (optional, default: `false`) if `false`, compute the support
-                   function exactly; otherwise use an overapproximative
-                   algorithm
 
 ### Output
 
@@ -495,14 +578,12 @@ julia> _line_search([1.0, 0.0], X, H, upper=1e3, method=GoldenSection())
 (1.0, 381.9660112501051)
 ```
 """
-function _line_search(ℓ, X, H::Union{HalfSpace, Hyperplane, Line};
-                      upper_bound::Bool=false, kwargs...)
+function _line_search(ℓ, X, H::Union{HalfSpace, Hyperplane, Line}; kwargs...)
     options = Dict(kwargs)
 
     # Initialization
     a, b = H.a, H.b
-    ρ_rec = upper_bound ? LazySets.Approximations.ρ_upper_bound : ρ
-    f(λ) = ρ_rec(ℓ - λ[1] * a, X) + λ[1] * b
+    f(λ) = ρ(ℓ - λ[1] * a, X) + λ[1] * b
 
     if haskey(options, :lower)
         lower = pop!(options, :lower)
@@ -541,7 +622,6 @@ end # load_optim
                 [lazy_linear_map]=false,
                 [lazy_2d_intersection]=true,
                 [algorithm_2d_intersection]=nothing,
-                [upper_bound]::Bool=false,
                 [kwargs...]) where {N}
 
 Given a compact and convex set ``X`` and a hyperplane ``H = \\{x: n ⋅ x = γ \\}``,
@@ -561,9 +641,6 @@ rank-2 projection ``Π_{nℓ} X`` and the line ``Lγ = \\{(x, y): x = γ \\}``.
 - `algorithm_2d_intersection` -- (optional, default: `nothing`) if given, fixes the
                                  support function algorithm used for the intersection
                                  in 2D; otherwise the default is implied
-- `upper_bound`          -- (optional, default: `false`) if `false`, compute the
-                            support function exactly; otherwise use an
-                            overapproximative algorithm
 
 ### Output
 
@@ -597,7 +674,6 @@ function _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}};
                      lazy_linear_map=false,
                      lazy_2d_intersection=true,
                      algorithm_2d_intersection=nothing,
-                     upper_bound::Bool=false,
                      kwargs...) where {N}
 
     n = H.a                  # normal vector to the hyperplane
@@ -612,10 +688,9 @@ function _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}};
     Xnℓ = lazy_linear_map ? LinearMap(Πnℓ, X) : linear_map(Πnℓ, X)
     Xnℓ⋂Lγ = lazy_2d_intersection ? Intersection(Xnℓ, Lγ) : intersection(Xnℓ, Lγ)
 
-    ρ_rec = upper_bound ? LazySets.Approximations.ρ_upper_bound : ρ
     if algorithm_2d_intersection == nothing
-        return ρ_rec(y_dir, Xnℓ⋂Lγ; kwargs...)
+        return ρ(y_dir, Xnℓ⋂Lγ; kwargs...)
     else
-        return ρ_rec(y_dir, Xnℓ⋂Lγ, algorithm=algorithm_2d_intersection; kwargs...)
+        return ρ(y_dir, Xnℓ⋂Lγ, algorithm=algorithm_2d_intersection; kwargs...)
     end
 end
diff --git a/src/is_intersection_empty.jl b/src/is_intersection_empty.jl
index fb059518f4..8436f9d4de 100644
--- a/src/is_intersection_empty.jl
+++ b/src/is_intersection_empty.jl
@@ -721,17 +721,15 @@ the half-space: ``σ(-a) ∈ hs``.
 The support vector is thus also a witness.
 
 Optional keyword arguments can be passed to the `ρ` function. In particular, if
-`X` is a lazy intersection, the option `upper_bound=true` may be useful.
+`X` is a lazy intersection, options can be passed to the line search algorithm.
 """
 function is_intersection_empty(X::LazySet{N},
                                hs::HalfSpace{N},
                                witness::Bool=false;
-                               upper_bound::Bool=false,
                                kwargs...
                                )::Union{Bool, Tuple{Bool,Vector{N}}} where N<:Real
     if !witness
-        ρ_rec = upper_bound ? LazySets.Approximations.ρ_upper_bound : ρ
-        return -ρ_rec(-hs.a, X; kwargs...) > hs.b
+        return -ρ(-hs.a, X; kwargs...) > hs.b
     end
 
     # for witness production, we compute the support vector instead
diff --git a/test/runtests.jl b/test/runtests.jl
index becffaf7c1..9e51174df5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -96,7 +96,6 @@ 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
diff --git a/test/unit_plot.jl b/test/unit_plot.jl
index bc1f848865..59c109462a 100644
--- a/test/unit_plot.jl
+++ b/test/unit_plot.jl
@@ -1,4 +1,4 @@
-using Plots
+using Plots, Optim
 
 for N in [Float64, Rational{Int}, Float32]
     p0 = zero(N)
@@ -47,8 +47,7 @@ for N in [Float64, Rational{Int}, Float32]
     # binary set operations
     ch = ConvexHull(b1, bi)
     cha = ConvexHullArray([b1, bi])
-    its = Intersection(b1, bi)
-    itsa = IntersectionArray([b1, bi])
+
     ms = MinkowskiSum(b1, bi)
     msa = MinkowskiSumArray([b1, bi])
     cms = CacheMinkowskiSum([b1, bi])
@@ -79,8 +78,6 @@ for N in [Float64, Rational{Int}, Float32]
     @test_throws ErrorException plot(l) # TODO see #576
     plot(ch)
     plot(cha)
-    @test_throws ErrorException plot(its) # TODO not implemented yet
-    @test_throws ErrorException plot(itsa) # TODO not implemented yet
     plot(sih)
     plot(lm)
     plot(ms)
@@ -109,8 +106,6 @@ for N in [Float64, Rational{Int}, Float32]
         @test_throws ErrorException plot(hs, ε) # TODO see #576/#578
         @test_throws ErrorException plot(hp, ε) # TODO see #576/#578
         @test_throws ErrorException plot(l, ε) # TODO see #576/#578
-        @test_throws ErrorException plot(its, ε) # TODO not implemented yet
-        @test_throws ErrorException plot(itsa, ε) # TODO not implemented yet
         plot(ch, ε)
         plot(cha, ε)
         plot(sih, ε)
@@ -169,3 +164,28 @@ for N in [Float64, Float32]
         @test_throws ErrorException plot(b2, ε) # TODO see #578
     end
 end
+
+# set types that only work with Float64
+for N in [Float64]
+    v0 = zeros(N, 2)
+    p1 = one(N)
+
+    b1 = Ball1(v0, p1)
+    bi = BallInf(v0, p1)
+
+    # binary set operations
+    its = Intersection(b1, bi)
+    itsa = IntersectionArray([b1, bi])
+
+    # -----
+    # plots
+    # -----
+
+    plot(its)
+    @test_throws ErrorException plot(itsa) # TODO not implemented yet
+
+    # ε-close
+    ε = N(1e-2)
+    plot(its)
+    @test_throws ErrorException plot(itsa, ε) # TODO not implemented yet
+end