diff --git a/docs/src/lib/binary_functions.md b/docs/src/lib/binary_functions.md
index 382cb2aeda..a9aece180e 100644
--- a/docs/src/lib/binary_functions.md
+++ b/docs/src/lib/binary_functions.md
@@ -29,6 +29,8 @@ end
 ⊆(::LineSegment{N}, ::LazySet{N}, ::Bool=false) where {N<:Real}
 ⊆(::LineSegment{N}, ::Hyperrectangle{N}, ::Bool=false) where {N<:Real}
 ⊆(::Interval, ::Interval)
+⊆(::EmptySet{N}, ::LazySet{N}, ::Bool=false) where {N<:Real}
+⊆(::LazySet{N}, ::EmptySet{N}, ::Bool=false) where {N<:Real}
 ```
 
 ## Check for emptiness of intersection
diff --git a/src/Approximations/box_approximations.jl b/src/Approximations/box_approximations.jl
index 920872c6f7..8dc90692a1 100644
--- a/src/Approximations/box_approximations.jl
+++ b/src/Approximations/box_approximations.jl
@@ -37,6 +37,9 @@ box_approximation(S::Hyperrectangle) = S
 box_approximation(S::AbstractHyperrectangle) =
     Hyperrectangle(center(S), radius_hyperrectangle(S))
 
+# special case: empty set
+box_approximation(∅::EmptySet) = ∅
+
 """
     interval_hull
 
@@ -74,6 +77,9 @@ function box_approximation_symmetric(S::LazySet{N};
     return Hyperrectangle(zeros(N, length(c)), abs.(c) .+ r)
 end
 
+# special case: empty set
+box_approximation_symmetric(∅::EmptySet) = ∅
+
 """
     symmetric_interval_hull
 
@@ -176,3 +182,6 @@ function ballinf_approximation(S::LazySet{N};
     end
     return BallInf(c, r)
 end
+
+# special case: empty set
+ballinf_approximation(∅::EmptySet) = ∅
diff --git a/src/is_subset.jl b/src/is_subset.jl
index 361af8628c..fae4a63702 100644
--- a/src/is_subset.jl
+++ b/src/is_subset.jl
@@ -482,6 +482,10 @@ function ⊆(L::LineSegment{N}, H::AbstractHyperrectangle{N}, witness::Bool=fals
     end
 end
 
+
+# --- Interval ---
+
+
 """
     ⊆(x::Interval, y::Interval)
 
@@ -499,3 +503,87 @@ Check whether an interval is contained in another interval.
 function ⊆(x::Interval, y::Interval)
     return x.dat ⊆ y.dat
 end
+
+
+# --- EmptySet ---
+
+
+"""
+    ⊆(∅::EmptySet{N}, X::LazySet{N}, [witness]::Bool=false
+     )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+
+Check whether an empty set is contained in another set.
+
+### Input
+
+- `∅`       -- empty set
+- `X`       -- another set
+- `witness` -- (optional, default: `false`) compute a witness if activated
+               (ignored, just kept for interface reasons)
+
+### Output
+
+`true`.
+"""
+function ⊆(∅::EmptySet{N}, X::LazySet{N}, witness::Bool=false
+          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    return witness ? (true, N[]) : true
+end
+
+# disambiguation
+function ⊆(∅::EmptySet{N}, H::AbstractHyperrectangle{N}, witness::Bool=false
+          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    return witness ? (true, N[]) : true
+end
+
+"""
+    ⊆(X::LazySet{N}, ∅::EmptySet{N}, [witness]::Bool=false
+     )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+
+Check whether a set is contained in an empty set.
+
+### Input
+
+- `X`       -- another set
+- `∅`       -- empty set
+- `witness` -- (optional, default: `false`) compute a witness if activated
+
+### Output
+
+`true` iff `X` is empty.
+
+### Algorithm
+
+We rely on `isempty(X)` for the emptiness check and on `an_element(X)` for
+witness production.
+"""
+function ⊆(X::LazySet{N}, ∅::EmptySet{N}, witness::Bool=false
+          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    if isempty(X)
+        return witness ? (true, N[]) : true
+    else
+        return witness ? (false, an_element(X)) : false
+    end
+end
+
+# disambiguation
+function ⊆(X::AbstractPolytope{N}, ∅::EmptySet{N}, witness::Bool=false
+          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    if isempty(X)
+        return witness ? (true, N[]) : true
+    else
+        return witness ? (false, an_element(X)) : false
+    end
+end
+function ⊆(X::AbstractSingleton{N}, ∅::EmptySet{N}, witness::Bool=false
+          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    return witness ? (false, an_element(X)) : false
+end
+function ⊆(X::LineSegment{N}, ∅::EmptySet{N}, witness::Bool=false
+          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    return witness ? (false, an_element(X)) : false
+end
+function ⊆(X::EmptySet{N}, ∅::EmptySet{N}, witness::Bool=false
+          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    return witness ? (true, N[]) : true
+end
diff --git a/test/unit_EmptySet.jl b/test/unit_EmptySet.jl
index 937d92be32..15dbba49cf 100644
--- a/test/unit_EmptySet.jl
+++ b/test/unit_EmptySet.jl
@@ -46,6 +46,12 @@ for N in [Float64, Rational{Int}, Float32]
     @test !∈(N[0], E)
     @test !∈(N[0, 0], E)
 
+    # subset
+    @test ⊆(E, B) && ⊆(E, B, true)[1]
+    subset, point = ⊆(B, E, true)
+    @test !⊆(B, E) && !subset && point ∈ B
+    @test ⊆(E, E) && ⊆(E, E, true)[1]
+
     # emptiness check
     @test isempty(E)
 
diff --git a/test/unit_ballinf_approximation.jl b/test/unit_ballinf_approximation.jl
index 1e5091507b..01c73ac16a 100644
--- a/test/unit_ballinf_approximation.jl
+++ b/test/unit_ballinf_approximation.jl
@@ -19,4 +19,8 @@ for N in [Float64, Float32, Rational{Int}]
     biexp = BallInf(N[0.5, 0], N(1.5))
     @test bi.center ≈ biexp.center
     @test bi.radius ≈ biexp.radius
+
+    # empty set
+    E = EmptySet{N}()
+    @test ballinf_approximation(E) == E
 end
diff --git a/test/unit_box_approximation.jl b/test/unit_box_approximation.jl
index a52313462d..71e1cf7b16 100644
--- a/test/unit_box_approximation.jl
+++ b/test/unit_box_approximation.jl
@@ -33,6 +33,10 @@ for N in [Float64, Rational{Int}, Float32]
     @test h.center ≈ hexp.center
     @test h.radius ≈ hexp.radius
 
+    # empty set
+    E = EmptySet{N}()
+    @test box_approximation(E) == E
+
     # ===================================================================
     # Testing box_approximation_symmetric (= symmetric interval hull)
     # ===================================================================
@@ -67,4 +71,8 @@ for N in [Float64, Rational{Int}, Float32]
 
     # Testing alias symmetric_interval_hull
     h = symmetric_interval_hull(b)
+
+    # empty set
+    E = EmptySet{N}()
+    @test box_approximation_symmetric(E) == E
 end