diff --git a/src/Interfaces/AbstractPolyhedron.jl b/src/Interfaces/AbstractPolyhedron.jl
index a406b60df1..273e90f638 100644
--- a/src/Interfaces/AbstractPolyhedron.jl
+++ b/src/Interfaces/AbstractPolyhedron.jl
@@ -78,7 +78,6 @@ function _linear_map_hrep_helper(M::AbstractMatrix, P::LazySet,
 end
 
 # internal functions; defined here due to optional dependencies and submodules
-function isfeasible end
 function remove_redundant_constraints end
 function remove_redundant_constraints! end
 function _ishalfspace end
diff --git a/src/Interfaces/AbstractPolyhedron_functions.jl b/src/Interfaces/AbstractPolyhedron_functions.jl
index 3bb8cce7d8..4ace2e7337 100644
--- a/src/Interfaces/AbstractPolyhedron_functions.jl
+++ b/src/Interfaces/AbstractPolyhedron_functions.jl
@@ -1,7 +1,6 @@
 export constrained_dimensions,
        remove_redundant_constraints,
        remove_redundant_constraints!,
-       isfeasible,
        addconstraint!,
        ishyperplanar
 
@@ -949,49 +948,6 @@ function _project_polyhedron(P::LazySet, block; kwargs...)
     return constraints_list(πP)
 end
 
-"""
-    isfeasible(A::AbstractMatrix, b::AbstractVector, [witness]::Bool=false;
-               [solver]=nothing)
-
-Check for feasibility of linear constraints given in matrix-vector form.
-
-### Input
-
-- `A`       -- constraints matrix
-- `b`       -- constraints vector
-- `witness` -- (optional; default: `false`) flag for witness production
-- `solver`  -- (optional; default: `nothing`) LP solver
-
-### Output
-
-If `witness` is `false`, the result is a `Bool`.
-
-If `witness` is `true`, the result is a pair `(res, w)` where `res` is a `Bool`
-and `w` is a witness point/vector.
-
-### Algorithm
-
-This implementation solves the corresponding feasibility linear program.
-"""
-function isfeasible(A::AbstractMatrix, b::AbstractVector, witness::Bool=false;
-                    solver=nothing)
-    N = promote_type(eltype(A), eltype(b))
-    # feasibility LP
-    lbounds, ubounds = -Inf, Inf
-    sense = '<'
-    obj = zeros(N, size(A, 2))
-    if isnothing(solver)
-        solver = default_lp_solver(N)
-    end
-    lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
-    if is_lp_optimal(lp.status)
-        return witness ? (true, lp.sol) : true
-    elseif is_lp_infeasible(lp.status)
-        return witness ? (false, N[]) : false
-    end
-    return error("LP returned status $(lp.status) unexpectedly")
-end
-
 # convenience function to invert the result of `isfeasible` while still including the witness result
 function _isinfeasible(constraints::AbstractVector{<:HalfSpace},
                        witness::Bool=false; solver=nothing)
diff --git a/src/Sets/HalfSpace/isfeasible.jl b/src/Sets/HalfSpace/isfeasible.jl
index d32ae49d79..a8eb3d337f 100644
--- a/src/Sets/HalfSpace/isfeasible.jl
+++ b/src/Sets/HalfSpace/isfeasible.jl
@@ -1,4 +1,3 @@
-
 """
     isfeasible(constraints::AbstractVector{<:HalfSpace}, [witness]::Bool=false;
                [solver]=nothing)
diff --git a/src/Utils/lp_solvers.jl b/src/Utils/lp_solvers.jl
index 46525b53b5..60ca0e5e02 100644
--- a/src/Utils/lp_solvers.jl
+++ b/src/Utils/lp_solvers.jl
@@ -1,3 +1,5 @@
+export isfeasible
+
 function default_lp_solver(::Type{T}) where {T}
     key = task_local_lp_solver_key(T)
     LP = get!(() -> JuMP.Model(default_lp_solver_factory(T)), task_local_storage(), key)
@@ -35,3 +37,46 @@ function default_lp_solver_polyhedra(N; kwargs...)
     require(@__MODULE__, :Polyhedra; fun_name="default_lp_solver_polyhedra")
     return error("no default solver for numeric type $N")
 end
+
+"""
+    isfeasible(A::AbstractMatrix, b::AbstractVector, [witness]::Bool=false;
+               [solver]=nothing)
+
+Check for feasibility of linear constraints given in matrix-vector form.
+
+### Input
+
+- `A`       -- constraints matrix
+- `b`       -- constraints vector
+- `witness` -- (optional; default: `false`) flag for witness production
+- `solver`  -- (optional; default: `nothing`) LP solver
+
+### Output
+
+If `witness` is `false`, the result is a `Bool`.
+
+If `witness` is `true`, the result is a pair `(res, w)` where `res` is a `Bool`
+and `w` is a witness point/vector.
+
+### Algorithm
+
+This implementation solves the corresponding feasibility linear program.
+"""
+function isfeasible(A::AbstractMatrix, b::AbstractVector, witness::Bool=false;
+                    solver=nothing)
+    N = promote_type(eltype(A), eltype(b))
+    # feasibility LP
+    lbounds, ubounds = -Inf, Inf
+    sense = '<'
+    obj = zeros(N, size(A, 2))
+    if isnothing(solver)
+        solver = default_lp_solver(N)
+    end
+    lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
+    if is_lp_optimal(lp.status)
+        return witness ? (true, lp.sol) : true
+    elseif is_lp_infeasible(lp.status)
+        return witness ? (false, N[]) : false
+    end
+    return error("LP returned status $(lp.status) unexpectedly")
+end