diff --git a/src/Initialization/init_JuMP.jl b/src/Initialization/init_JuMP.jl
index df5dc976c8..425e9345e0 100644
--- a/src/Initialization/init_JuMP.jl
+++ b/src/Initialization/init_JuMP.jl
@@ -30,25 +30,60 @@ function has_lp_infeasibility_ray(model)
     return primal_status(model) == JuMP.INFEASIBILITY_CERTIFICATE
 end
 
-# solve a linear program (in the old MathProgBase interface)
-function linprog(c, A, sense::Char, b, l::Number, u::Number, solver)
+# solve a linear program (in the MathOptInterface)
+function linprog(c, A, sense, b, l, u, solver_or_model)
     n = length(c)
-    m = length(b)
-    return linprog(c, A, fill(sense, m), b, fill(l, n), fill(u, n), solver)
-end
+    model = linprog_model(solver_or_model)
+
+    x = lin_prog_variable!(model, n, l, u)
+    lin_prog_constraints!(model, A, sense, b, x)
 
-function linprog(c, A, sense, b, l, u, solver)
-    n = length(c)
-    model = Model(solver)
-    @variable(model, l[i] <= x[i=1:n] <= u[i])
     @objective(model, Min, c' * x)
-    eq_rows, ge_rows, le_rows = sense .== '=', sense .== '>', sense .== '<'
-    @constraint(model, A[eq_rows, :] * x .== b[eq_rows])
-    @constraint(model, A[ge_rows, :] * x .>= b[ge_rows])
-    @constraint(model, A[le_rows, :] * x .<= b[le_rows])
+
     optimize!(model)
     return (status = termination_status(model),
             objval = objective_value(model),
             sol    = value.(x),
             model  = model)
 end
+
+@inline function linprog_model(model::JuMP.Model)
+    Base.empty!(model)
+    return model
+end
+@inline linprog_model(solver) = Model(solver)
+
+function lin_prog_variable!(model, n, l::Number, u::Number)
+    if isfinite(l) && isfinite(u)
+        return @variable(model, x[1:n], lower_bound=l, upper_bound=u)
+    elseif isfinite(l)
+        return @variable(model, x[1:n], lower_bound=l)
+    elseif isfinite(u)
+        return @variable(model, x[1:n], upper_bound=u)
+    else 
+        return @variable(model, x[1:n])
+    end
+end
+
+function lin_prog_variable!(model, n, l::Vector{T}, u::Vector{T}) where {T}
+    return @variable(model, l[i] <= x[i=1:n] <= u[i])
+end
+
+function lin_prog_constraints!(model, A, sense::Char, b, x)
+    if sense == '='
+        @constraint(model, A * x == b)
+    elseif sense == '>'
+        @constraint(model, A * x >= b)
+    elseif sense == '<'
+        @constraint(model, A * x <= b)
+    else
+        throw("Invalid sense ($sense) for constraints in lin_prog.")
+    end
+end
+
+function lin_prog_constraints!(model, A, sense, b, x)
+    eq_rows, ge_rows, le_rows = sense .== '=', sense .== '>', sense .== '<'
+    @constraint(model, view(A, eq_rows, :) * x == view(b, eq_rows))
+    @constraint(model, view(A, ge_rows, :) * x >= view(b, ge_rows))
+    @constraint(model, view(A, le_rows, :) * x <= view(b, le_rows))
+end
\ No newline at end of file
diff --git a/src/Utils/lp_solvers.jl b/src/Utils/lp_solvers.jl
index 5c910b5b64..7fac856c28 100644
--- a/src/Utils/lp_solvers.jl
+++ b/src/Utils/lp_solvers.jl
@@ -1,17 +1,45 @@
+const THREAD_FLOAT_LP_SOLVERs = JuMP.Model[]
+const THREAD_EXACT_LP_SOLVERs = JuMP.Model[]
+
+@inline default_lp_solver(::Type{T}) where {T} = default_lp_solver(T, Threads.threadid())
+@noinline function default_lp_solver(::Type{T}, tid::Int) where {T}
+    THREAD_LP_SOLVERs = thread_specific_lp_solvers(T)
+
+    0 < tid <= length(THREAD_LP_SOLVERs) || _rng_length_assert()
+    if @inbounds isassigned(THREAD_LP_SOLVERs, tid)
+        @inbounds LP = THREAD_LP_SOLVERs[tid]
+    else
+        LP = Model(default_lp_solver_factory(T))
+        @inbounds THREAD_LP_SOLVERs[tid] = LP
+    end
+    return LP
+end
+@noinline _rng_length_assert() =  @assert false "0 < tid <= length(THREAD_RNGs)"
+
+function init_lp_solvers()
+    # Code is heavily inspired by the global state of RNGs in Random.
+    # That code contains the following comment, although it is unclear
+    # what this pertains to (most likely empty!).
+    # "it ensures that we didn't save a bad object"
+    resize!(empty!(THREAD_FLOAT_LP_SOLVERs), Threads.nthreads())
+    resize!(empty!(THREAD_EXACT_LP_SOLVERs), Threads.nthreads())
+end
+
 # default LP solver for floating-point numbers
-function default_lp_solver(::Type{<:AbstractFloat})
+@inline function default_lp_solver_factory(::Type{<:AbstractFloat})
     return JuMP.optimizer_with_attributes(() -> GLPK.Optimizer(; method=GLPK.SIMPLEX))
 end
 
 # default LP solver for rational numbers
-function default_lp_solver(::Type{<:Rational})
+@inline function default_lp_solver_factory(::Type{<:Rational})
     return JuMP.optimizer_with_attributes(() -> GLPK.Optimizer(; method=GLPK.EXACT))
 end
 
+@inline thread_specific_lp_solvers(::Type{<:AbstractFloat}) = THREAD_FLOAT_LP_SOLVERs
+@inline thread_specific_lp_solvers(::Type{<:Rational}) = THREAD_EXACT_LP_SOLVERs
+
 # default LP solver given two possibly different numeric types
-function default_lp_solver(M::Type{<:Number}, N::Type{<:Number})
-    return default_lp_solver(promote_type(M, N))
-end
+@inline default_lp_solver(M::Type{<:Number}, N::Type{<:Number}) = default_lp_solver(promote_type(M, N))
 
 # check for Polyhedra backend (fallback method)
 function _is_polyhedra_backend(backend)
diff --git a/src/init.jl b/src/init.jl
index 6b697efc37..b3eaec5f97 100644
--- a/src/init.jl
+++ b/src/init.jl
@@ -14,6 +14,8 @@ function __init__()
     @require RangeEnclosures = "1b4d18b6-9e5d-11e9-236c-f792b01831f8" include("Initialization/init_RangeEnclosures.jl")
     @require Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" include("Initialization/init_Symbolics.jl")
     @require WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" include("Initialization/init_WriteVTK.jl")
+
+    init_lp_solvers()
 end
 
 # there is a conflict with `Revise` and `StaticArrays(Core)`, so it is loaded without `Requires`
diff --git a/test/Project.toml b/test/Project.toml
index 1bf8f43b9c..17468e7439 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,12 +4,14 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4"
 ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 # IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153"
 IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Javis = "78b212ba-a7f9-42d4-b726-60726080707e"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca"
@@ -36,7 +38,6 @@ Expokit = "0.2"
 ExponentialUtilities = "1"
 GR = "0"
 IntervalArithmetic = "0.15 - 0.20"
-# IntervalConstraintProgramming = "0.9 - 0.12"
 IntervalMatrices = "0.8"
 Ipopt = "1"
 Javis = "0.5 - 0.9"
diff --git a/test/Utils/lp_solvers.jl b/test/Utils/lp_solvers.jl
new file mode 100644
index 0000000000..480b62d8c6
--- /dev/null
+++ b/test/Utils/lp_solvers.jl
@@ -0,0 +1,47 @@
+# Check that the default model and explicit solver
+# specification works for linprog
+
+for N in [Float64]
+    p = HPolyhedron{N}()
+    c1 = LinearConstraint(N[2, 2], N(12))
+    c2 = LinearConstraint(N[-3, 3], N(6))
+    c3 = LinearConstraint(N[-1, -1], N(0))
+    c4 = LinearConstraint(N[2, -4], N(0))
+    addconstraint!(p, c3)
+    addconstraint!(p, c1)
+    addconstraint!(p, c4)
+    addconstraint!(p, c2)
+
+    P = HPolyhedron([HalfSpace(N[3 // 50, -4 // 10], N(1)),
+                     HalfSpace(N[-1 // 50, 1 // 10], N(-1))])
+
+    # With default model 
+    # support vector
+    d = N[1, 0]
+    @test σ(d, p) == N[4, 2]
+    d = N[0, 1]
+    @test σ(d, p) == N[2, 4]
+    d = N[-1, 0]
+    @test σ(d, p) == N[-1, 1]
+    d = N[0, -1]
+    @test σ(d, p) == N[0, 0]
+
+    # an_element
+    @test an_element(P) ∈ P
+
+    # With explicit solver
+    solver = optimizer_with_attributes(() -> GLPK.Optimizer(; method=GLPK.SIMPLEX))
+
+    # support vector
+    d = N[1, 0]
+    @test σ(d, p; solver=solver) == N[4, 2]
+    d = N[0, 1]
+    @test σ(d, p; solver=solver) == N[2, 4]
+    d = N[-1, 0]
+    @test σ(d, p; solver=solver) == N[-1, 1]
+    d = N[0, -1]
+    @test σ(d, p; solver=solver) == N[0, 0]
+
+    # an_element
+    @test an_element(P; solver=solver) ∈ P
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 5d6daddaa2..b1b7e52030 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 using LazySets, LazySets.Approximations, Test, LinearAlgebra, SparseArrays,
-      StaticArrays
+      StaticArrays, GLPK
+using JuMP: optimizer_with_attributes
 
 # fix random number generator seed
 using Random
@@ -287,6 +288,13 @@ if test_suite_basic
         include("Approximations/hausdorff_distance.jl")
     end
 
+    # ================================
+    # Testing shared utility functions
+    # ================================
+    @time @testset "LazySets.lp_solvers" begin
+        include("Utils/lp_solvers.jl")
+    end
+
     # ========================
     # Testing method ambiguity
     # ========================