diff --git a/src/bb_inits_and_defaults.jl b/src/bb_inits_and_defaults.jl
index 31c3509..43e53d6 100644
--- a/src/bb_inits_and_defaults.jl
+++ b/src/bb_inits_and_defaults.jl
@@ -220,10 +220,22 @@ function init_juniper_problem!(jp::JuniperProblem, model::MOI.AbstractOptimizer)
     jp.start_time = time()
 
     jp.nl_solver = model.options.nl_solver
+    jp.mip_model = nothing
 
     if model.options.mip_solver !== nothing
         jp.mip_solver = model.options.mip_solver
     end
+    if model.options.mip_model !== nothing
+        if MOI.get(model.options.mip_model, MOI.TerminationStatus()) ==
+           MOI.OPTIMIZE_NOT_CALLED
+            throw(
+                ErrorException(
+                    "The MIP model has not been solved, optimize it before setting `mip_model`.",
+                ),
+            )
+        end
+        jp.mip_model = model.options.mip_model.moi_backend.optimizer.model
+    end
     jp.options = model.options
     if MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
         jp.obj_sense = :Min
diff --git a/src/fpump.jl b/src/fpump.jl
index 1db7dc0..21fde53 100644
--- a/src/fpump.jl
+++ b/src/fpump.jl
@@ -337,6 +337,7 @@ function fpump(optimizer, m)
     # the tolerance can be changed => current atol
     catol = m.options.atol
     atol_counter = 0
+    mip_model = isnothing(m.mip_model) ? optimizer : m.mip_model
     while !are_type_correct(nlp_sol, m.var_type, m.disc2var_idx, catol) &&
               time() - start_fpump < tl &&
               time() - m.start_time < m.options.time_limit
@@ -350,7 +351,7 @@ function fpump(optimizer, m)
             ),
         )
             mip_status, mip_sol, mip_obj =
-                generate_mip(optimizer, m, nlp_sol, tabu_list, start_fpump)
+                generate_mip(mip_model, m, nlp_sol, tabu_list, start_fpump)
         else
             # if no linear constraints just round the discrete variables
             mip_obj = NaN
diff --git a/src/solver.jl b/src/solver.jl
index 7d688a7..00a2e17 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -39,6 +39,7 @@ function get_default_options()
     tabu_list_length = 30
     num_resolve_nlp_feasibility_pump = 1
     mip_solver = nothing
+    mip_model = nothing
     allow_almost_solved = true
     allow_almost_solved_integral = true
     registered_functions = nothing
@@ -82,6 +83,7 @@ function get_default_options()
         tabu_list_length,
         num_resolve_nlp_feasibility_pump,
         mip_solver,
+        mip_model,
         allow_almost_solved,
         allow_almost_solved_integral,
         registered_functions,
diff --git a/src/types.jl b/src/types.jl
index 197b4b1..0e69d2c 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -38,6 +38,7 @@ mutable struct SolverOptions
     tabu_list_length::Int64
     num_resolve_nlp_feasibility_pump::Int64
     mip_solver::Any
+    mip_model::Any
     allow_almost_solved::Bool
     allow_almost_solved_integral::Bool
     registered_functions::Union{Nothing,Vector{RegisteredFunction}}
@@ -97,6 +98,7 @@ mutable struct JuniperProblem
     nsolutions::Int64
 
     mip_solver::Any
+    mip_model::Union{MOI.AbstractOptimizer,Nothing}
 
     relaxation_time::Float64
     start_time::Float64
diff --git a/test/fpump.jl b/test/fpump.jl
index 3451507..5f13f52 100644
--- a/test/fpump.jl
+++ b/test/fpump.jl
@@ -292,4 +292,37 @@ include("basic/gamsworld.jl")
         end
         @test JuMP.objective_value(m) ≈ 0.0
     end
+
+    @testset "Custom linear relaxation" begin
+        optimizer = optimizer_with_attributes(
+            Juniper.Optimizer,
+            DefaultTestSolver(
+                branch_strategy = :MostInfeasible,
+                mip_solver = optimizer_with_attributes(
+                    HiGHS.Optimizer,
+                    "output_flag" => false,
+                ),
+            )...,
+        )
+        model = Model(optimizer)
+        @variable(model, a, integer = true)
+        @constraint(model, 0 <= model[:a] <= 10)
+        @NLconstraint(model, model[:a] * abs(model[:a]) >= 3)
+        @objective(model, Min, model[:a])
+
+        mip = Model(optimizer)
+        @variable(mip, a, integer = true)
+        @constraint(mip, mip[:a] * mip[:a] == 0)
+        @constraint(mip, mip[:a] <= 10)
+        @objective(mip, Min, mip[:a])
+        set_silent(mip)
+
+        set_optimizer_attribute(model, "mip_model", mip)
+        # optimize!(mip) hasn't been run so we expect an error
+        @test_throws ErrorException JuMP.optimize!(model)
+        JuMP.optimize!(mip)
+        set_optimizer_attribute(model, "mip_model", mip)
+        JuMP.optimize!(model)
+        @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
+    end
 end