diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl
index 91bfd82f..0cff55a0 100644
--- a/src/MOI_wrapper.jl
+++ b/src/MOI_wrapper.jl
@@ -2512,6 +2512,33 @@ function _dual_multiplier(model::Optimizer)
     return MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE ? 1.0 : -1.0
 end
 
+"""
+    _farkas_variable_dual(model::Optimizer, col::Cint)
+
+Return a Farkas dual associated with the variable bounds of `col`.
+
+Gurobi computes the Farkas dual as:
+
+    ā * x = λ' * A * x <= λ' * b = -β + sum(āᵢ * Uᵢ | āᵢ < 0) + sum(āᵢ * Lᵢ | āᵢ > 0)
+
+The Farkas dual of the variable is ā, and it applies to the upper bound if ā < 0,
+and it applies to the lower bound if ā > 0.
+"""
+function _farkas_variable_dual(model::Optimizer, col::Cint)
+    numnzP = Ref{Cint}()
+    ret = GRBgetvars(model, numnzP, C_NULL, C_NULL, C_NULL, col, 1)
+    _check_ret(model, ret)
+    vbeg = Vector{Cint}(undef, 2)
+    vind = Vector{Cint}(undef, numnzP[])
+    vval = Vector{Cdouble}(undef, numnzP[])
+    ret = GRBgetvars(model, numnzP, vbeg, vind, vval, col, 1)
+    _check_ret(model, ret)
+    λ = Vector{Cdouble}(undef, numnzP[])
+    ret = GRBgetdblattrlist(model, "FarkasDual", length(vind), vind, λ)
+    _check_ret(model, ret)
+    return λ' * vval
+end
+
 function MOI.get(
     model::Optimizer,
     attr::MOI.ConstraintDual,
@@ -2519,10 +2546,13 @@ function MOI.get(
 )
     _throw_if_optimize_in_progress(model, attr)
     MOI.check_result_index_bounds(model, attr)
+    col = Cint(column(model, c) - 1)
+    if model.has_infeasibility_cert
+        dual = _farkas_variable_dual(model, col)
+        return min(dual, 0.0)
+    end
     reduced_cost = Ref{Cdouble}()
-    ret = GRBgetdblattrelement(
-        model, "RC", Cint(column(model, c) - 1), reduced_cost
-    )
+    ret = GRBgetdblattrelement(model, "RC", col, reduced_cost)
     _check_ret(model, ret)
     sense = MOI.get(model, MOI.ObjectiveSense())
     # The following is a heuristic for determining whether the reduced cost
@@ -2550,10 +2580,13 @@ function MOI.get(
 )
     _throw_if_optimize_in_progress(model, attr)
     MOI.check_result_index_bounds(model, attr)
+    col = Cint(column(model, c) - 1)
+    if model.has_infeasibility_cert
+        dual = _farkas_variable_dual(model, col)
+        return max(dual, 0.0)
+    end
     reduced_cost = Ref{Cdouble}()
-    ret = GRBgetdblattrelement(
-        model, "RC", Cint(column(model, c) - 1), reduced_cost
-    )
+    ret = GRBgetdblattrelement(model, "RC", col, reduced_cost)
     _check_ret(model, ret)
     sense = MOI.get(model, MOI.ObjectiveSense())
     # The following is a heuristic for determining whether the reduced cost
@@ -2581,11 +2614,12 @@ function MOI.get(
 )
     _throw_if_optimize_in_progress(model, attr)
     MOI.check_result_index_bounds(model, attr)
+    col = Cint(column(model, c) - 1)
+    if model.has_infeasibility_cert
+        return _farkas_variable_dual(model, col)
+    end
     reduced_cost = Ref{Cdouble}()
-    reduced_cost = Ref{Cdouble}()
-    ret = GRBgetdblattrelement(
-        model, "RC", Cint(column(model, c) - 1), reduced_cost
-    )
+    ret = GRBgetdblattrelement(model, "RC", col, reduced_cost)
     _check_ret(model, ret)
     return _dual_multiplier(model) * reduced_cost[]
 end
@@ -2597,41 +2631,45 @@ function MOI.get(
 )
     _throw_if_optimize_in_progress(model, attr)
     MOI.check_result_index_bounds(model, attr)
+    col = Cint(column(model, c) - 1)
+    if model.has_infeasibility_cert
+        return _farkas_variable_dual(model, col)
+    end
     reduced_cost = Ref{Cdouble}()
-    ret = GRBgetdblattrelement(
-        model, "RC", Cint(column(model, c) - 1), reduced_cost
-    )
+    ret = GRBgetdblattrelement(model, "RC", col, reduced_cost)
     _check_ret(model, ret)
     return _dual_multiplier(model) * reduced_cost[]
 end
 
 function MOI.get(
-    model::Optimizer, attr::MOI.ConstraintDual,
-    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:Any}
+    model::Optimizer,
+    attr::MOI.ConstraintDual,
+    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:Any},
 )
     _throw_if_optimize_in_progress(model, attr)
     MOI.check_result_index_bounds(model, attr)
     valueP = Ref{Cdouble}()
+    row = Cint(_info(model, c).row - 1)
     if model.has_infeasibility_cert
-        ret = GRBgetdblattrelement(
-            model, "FarkasDual", Cint(_info(model, c).row - 1), valueP
-        )
+        ret = GRBgetdblattrelement(model, "FarkasDual", row, valueP)
         _check_ret(model, ret)
-        return -_dual_multiplier(model) * valueP[]
+        return -valueP[]
     end
-    ret = GRBgetdblattrelement(
-        model, "Pi", Cint(_info(model, c).row - 1), valueP
-    )
+    ret = GRBgetdblattrelement(model, "Pi", row, valueP)
     _check_ret(model, ret)
     return _dual_multiplier(model) * valueP[]
 end
 
 function MOI.get(
-    model::Optimizer, attr::MOI.ConstraintDual,
-    c::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, <:Any}
+    model::Optimizer,
+    attr::MOI.ConstraintDual,
+    c::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, <:Any},
 )
     _throw_if_optimize_in_progress(model, attr)
     MOI.check_result_index_bounds(model, attr)
+    if model.has_infeasibility_cert
+        error("Infeasibility certificate not available for $(c).")
+    end
     valueP = Ref{Cdouble}()
     ret = GRBgetdblattrelement(
         model, "QCPi", Cint(_info(model, c).row - 1), valueP
diff --git a/test/MOI/MOI_wrapper.jl b/test/MOI/MOI_wrapper.jl
index 63169ab5..f365b5ed 100644
--- a/test/MOI/MOI_wrapper.jl
+++ b/test/MOI/MOI_wrapper.jl
@@ -70,11 +70,7 @@ function test_modificationtest()
 end
 
 function test_contlineartest()
-    MOIT.contlineartest(OPTIMIZER, MOIT.TestConfig(basis = true), [
-        # This requires an infeasiblity certificate for a variable bound.
-        "linear12"
-    ])
-    MOIT.linear12test(OPTIMIZER, MOIT.TestConfig(infeas_certificates=false))
+    MOIT.contlineartest(OPTIMIZER, MOIT.TestConfig(basis = true))
 end
 
 function test_contquadratictest()
@@ -1082,7 +1078,8 @@ function test_InterruptException()
 end
 
 function test_indicator_name()
-    MOI.empty!(OPTIMIZER)
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
     x = MOI.add_variables(model, 2)
     MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
     f = MOI.VectorAffineFunction(
@@ -1099,7 +1096,8 @@ function test_indicator_name()
 end
 
 function test_indicator_on_one()
-    MOI.empty!(OPTIMIZER)
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
     x = MOI.add_variables(model, 2)
     MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
     f = MOI.VectorAffineFunction(
@@ -1116,7 +1114,8 @@ function test_indicator_on_one()
 end
 
 function test_indicator_on_zero()
-    MOI.empty!(OPTIMIZER)
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
     x = MOI.add_variables(model, 2)
     MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
     f = MOI.VectorAffineFunction(
@@ -1133,7 +1132,8 @@ function test_indicator_on_zero()
 end
 
 function test_indicator_nonconstant_x()
-    MOI.empty!(OPTIMIZER)
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
     x = MOI.add_variables(model, 2)
     MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
     f = MOI.VectorAffineFunction(
@@ -1148,7 +1148,8 @@ function test_indicator_nonconstant_x()
 end
 
 function test_indicator_too_many_indicators()
-    MOI.empty!(OPTIMIZER)
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
     x = MOI.add_variables(model, 2)
     MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
     f = MOI.VectorAffineFunction(
@@ -1163,7 +1164,8 @@ function test_indicator_too_many_indicators()
 end
 
 function test_indicator_nonconstant()
-    MOI.empty!(OPTIMIZER)
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
     x = MOI.add_variables(model, 2)
     MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
     f = MOI.VectorAffineFunction(
@@ -1177,6 +1179,198 @@ function test_indicator_nonconstant()
     @test_throws ErrorException MOI.add_constraint(model, f, s)
 end
 
+function test_farkas_dual_min()
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
+    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
+    x = MOI.add_variables(model, 2)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(
+        model,
+        MOI.ObjectiveFunction{MOI.SingleVariable}(),
+        MOI.SingleVariable(x[1]),
+    )
+    clb = MOI.add_constraint.(
+        model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)
+    )
+    c = MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
+        MOI.LessThan(-1.0),
+    )
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
+    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
+    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
+    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
+    @show clb_dual, c_dual
+    @test clb_dual[1] > -1e-6
+    @test clb_dual[2] > -1e-6
+    @test c_dual[1] < 1e-6
+    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
+    @test clb_dual[2] ≈ -c_dual atol = 1e-6
+end
+
+function test_farkas_dual_min_interval()
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
+    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
+    x = MOI.add_variables(model, 2)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(
+        model,
+        MOI.ObjectiveFunction{MOI.SingleVariable}(),
+        MOI.SingleVariable(x[1]),
+    )
+    clb = MOI.add_constraint.(
+        model, MOI.SingleVariable.(x), MOI.Interval(0.0, 10.0)
+    )
+    c = MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
+        MOI.LessThan(-1.0),
+    )
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
+    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
+    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
+    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
+    @show clb_dual, c_dual
+    @test clb_dual[1] > -1e-6
+    @test clb_dual[2] > -1e-6
+    @test c_dual[1] < 1e-6
+    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
+    @test clb_dual[2] ≈ -c_dual atol = 1e-6
+end
+
+function test_farkas_dual_min_equalto()
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
+    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
+    x = MOI.add_variables(model, 2)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(
+        model,
+        MOI.ObjectiveFunction{MOI.SingleVariable}(),
+        MOI.SingleVariable(x[1]),
+    )
+    clb = MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.EqualTo(0.0))
+    c = MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
+        MOI.LessThan(-1.0),
+    )
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
+    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
+    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
+    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
+    @show clb_dual, c_dual
+    @test clb_dual[1] > -1e-6
+    @test clb_dual[2] > -1e-6
+    @test c_dual[1] < 1e-6
+    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
+    @test clb_dual[2] ≈ -c_dual atol = 1e-6
+end
+
+function test_farkas_dual_min_ii()
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
+    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
+    x = MOI.add_variables(model, 2)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(
+        model,
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
+        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0),
+    )
+    clb = MOI.add_constraint.(
+        model, MOI.SingleVariable.(x), MOI.LessThan(0.0)
+    )
+    c = MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
+        MOI.LessThan(-1.0),
+    )
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
+    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
+    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
+    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
+    @show clb_dual, c_dual
+    @test clb_dual[1] < 1e-6
+    @test clb_dual[2] < 1e-6
+    @test c_dual[1] < 1e-6
+    @test clb_dual[1] ≈ 2 * c_dual atol = 1e-6
+    @test clb_dual[2] ≈ c_dual atol = 1e-6
+end
+
+function test_farkas_dual_max()
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
+    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
+    x = MOI.add_variables(model, 2)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+    MOI.set(
+        model,
+        MOI.ObjectiveFunction{MOI.SingleVariable}(),
+        MOI.SingleVariable(x[1]),
+    )
+    clb = MOI.add_constraint.(
+        model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)
+    )
+    c = MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
+        MOI.LessThan(-1.0),
+    )
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
+    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
+    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
+    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
+    @show clb_dual, c_dual
+    @test clb_dual[1] > -1e-6
+    @test clb_dual[2] > -1e-6
+    @test c_dual[1] < 1e-6
+    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
+    @test clb_dual[2] ≈ -c_dual atol = 1e-6
+end
+
+function test_farkas_dual_max_ii()
+    model = Gurobi.Optimizer(GRB_ENV)
+    MOI.set(model, MOI.Silent(), true)
+    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
+    x = MOI.add_variables(model, 2)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+    MOI.set(
+        model,
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
+        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0),
+    )
+    clb = MOI.add_constraint.(
+        model, MOI.SingleVariable.(x), MOI.LessThan(0.0)
+    )
+    c = MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
+        MOI.LessThan(-1.0),
+    )
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
+    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
+    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
+    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
+    @show clb_dual, c_dual
+    @test clb_dual[1] < 1e-6
+    @test clb_dual[2] < 1e-6
+    # Confirmed bug in Gurobi <=9.1 and fixed in next release. When updating for new 
+    # release, confirm fixed and add a conditional check on `_GRB_VERSION`.
+    @test_broken c_dual[1] < 1e-6
+    @test_broken clb_dual[1] ≈ 2 * c_dual atol = 1e-6
+    @test_broken clb_dual[2] ≈ c_dual atol = 1e-6
+end
+
 end
 
 runtests(TestMOIWrapper)