From 936547d7adbb55690dbcf571dc19924c9b83f74a Mon Sep 17 00:00:00 2001 From: odow Date: Wed, 4 Nov 2020 14:16:04 +1300 Subject: [PATCH] Add farkas certificate for variable bounds --- src/MOI/MOI_wrapper.jl | 77 +++++++++-- test/MathOptInterface/MOI_wrapper.jl | 198 ++++++++++++++++++++++++++- 2 files changed, 253 insertions(+), 22 deletions(-) diff --git a/src/MOI/MOI_wrapper.jl b/src/MOI/MOI_wrapper.jl index cfa4fed3..8f9284f6 100644 --- a/src/MOI/MOI_wrapper.jl +++ b/src/MOI/MOI_wrapper.jl @@ -2720,6 +2720,51 @@ 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`. + +Compute 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) + nzcnt_p, surplus_p = Ref{Cint}(), Ref{Cint}() + cmatbeg = Vector{Cint}(undef, 2) + ret = CPXgetcols( + model.env, + model.lp, + nzcnt_p, + cmatbeg, + C_NULL, + C_NULL, + Cint(0), + surplus_p, + col, + col, + ) + cmatind = Vector{Cint}(undef, -surplus_p[]) + cmatval = Vector{Cdouble}(undef, -surplus_p[]) + ret = CPXgetcols( + model.env, + model.lp, + nzcnt_p, + cmatbeg, + cmatind, + cmatval, + -surplus_p[], + surplus_p, + col, + col, + ) + _check_ret(model, ret) + return sum(v * model.certificate[i + 1] for (i, v) in zip(cmatind, cmatval)) +end + function MOI.get( model::Optimizer, attr::MOI.ConstraintDual, @@ -2728,6 +2773,10 @@ 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_dual_certificate + dual = -_farkas_variable_dual(model, col) + return min(0.0, dual) + end p = Ref{Cdouble}() ret = CPXgetdj(model.env, model.lp, p, col, col) _check_ret(model, ret) @@ -2758,6 +2807,10 @@ 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_dual_certificate + dual = -_farkas_variable_dual(model, col) + return max(0.0, dual) + end p = Ref{Cdouble}() ret = CPXgetdj(model.env, model.lp, p, col, col) _check_ret(model, ret) @@ -2783,25 +2836,16 @@ end function MOI.get( model::Optimizer, attr::MOI.ConstraintDual, - c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}, -) - _throw_if_optimize_in_progress(model, attr) - MOI.check_result_index_bounds(model, attr) - col = Cint(column(model, c) - 1) - p = Ref{Cdouble}() - ret = CPXgetdj(model.env, model.lp, p, col, col) - _check_ret(model, ret) - return _dual_multiplier(model) * p[] -end - -function MOI.get( - model::Optimizer, - attr::MOI.ConstraintDual, - c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}}, + c::MOI.ConstraintIndex{ + MOI.SingleVariable, <:Union{MOI.Interval{Float64}, MOI.EqualTo{Float64}} + }, ) _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) col = Cint(column(model, c) - 1) + if model.has_dual_certificate + return -_farkas_variable_dual(model, col) + end p = Ref{Cdouble}() ret = CPXgetdj(model.env, model.lp, p, col, col) _check_ret(model, ret) @@ -2834,6 +2878,9 @@ function MOI.get( # https://www.ibm.com/support/knowledgecenter/SSSA5P_12.10.0/ilog.odms.cplex.help/CPLEX/UsrMan/topics/cont_optim/qcp/17_QCP_duals.html _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) + if model.has_dual_certificate + error("Infeasibility certificate not available for $(c)") + end # The derivative of a quadratic f(x) = x^TQx + a^Tx + b <= 0 is # ∇f(x) = Q^Tx + Qx + a # The dual is undefined if x is at the point of the cone. This can only be diff --git a/test/MathOptInterface/MOI_wrapper.jl b/test/MathOptInterface/MOI_wrapper.jl index 12ac39e0..84fc4c6d 100644 --- a/test/MathOptInterface/MOI_wrapper.jl +++ b/test/MathOptInterface/MOI_wrapper.jl @@ -67,13 +67,7 @@ function test_modificationtest() end function test_contlineartest() - MOIT.contlineartest(BRIDGED_OPTIMIZER, CONFIG, [ - # TODO(odow): This test requests the infeasibility certificate of a - # variable bound. - "linear12" - ]) - - MOIT.linear12test(OPTIMIZER, MOIT.TestConfig(infeas_certificates=false)) + MOIT.contlineartest(BRIDGED_OPTIMIZER, CONFIG) end function test_intlineartest() @@ -542,6 +536,196 @@ function test_ZeroOne_INTERVAL() @test tmp[] == 2.0 end +function test_farkas_dual_min() + model = CPLEX.Optimizer() + MOI.set(model, MOI.Silent(), true) + MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) + 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 = CPLEX.Optimizer() + MOI.set(model, MOI.Silent(), true) + MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) + 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 = CPLEX.Optimizer() + MOI.set(model, MOI.Silent(), true) + MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) + 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 = CPLEX.Optimizer() + MOI.set(model, MOI.Silent(), true) + MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) + 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 = CPLEX.Optimizer() + MOI.set(model, MOI.Silent(), true) + MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) + 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 = CPLEX.Optimizer() + MOI.set(model, MOI.Silent(), true) + MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) + 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 + @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 + end # module TestMOIwrapper runtests(TestMOIwrapper)