Skip to content

Add farkas certificate for variable bounds #337

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
odow merged 1 commit into from
Nov 5, 2020
Merged

Add farkas certificate for variable bounds #337

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 62 additions & 15 deletions src/MOI/MOI_wrapper.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
198 changes: 191 additions & 7 deletions test/MathOptInterface/MOI_wrapper.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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)