Skip to content

Commit 7019ca7

Browse files
odowodow
authored
Add farkas certificate for variable bounds (#337)
1 parent c04b2ad commit 7019ca7

File tree

2 files changed

+253
-22
lines changed

2 files changed

+253
-22
lines changed

‎src/MOI/MOI_wrapper.jl‎

Lines changed: 62 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -2720,6 +2720,51 @@ function _dual_multiplier(model::Optimizer)
27202720
return MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE ? 1.0 : -1.0
27212721
end
27222722

2723+
"""
2724+
_farkas_variable_dual(model::Optimizer, col::Cint)
2725+
2726+
Return a Farkas dual associated with the variable bounds of `col`.
2727+
2728+
Compute the Farkas dual as:
2729+
2730+
ā * x = λ' * A * x <= λ' * b = -β + sum(āᵢ * Uᵢ | āᵢ < 0) + sum(āᵢ * Lᵢ | āᵢ > 0)
2731+
2732+
The Farkas dual of the variable is ā, and it applies to the upper bound if ā < 0,
2733+
and it applies to the lower bound if ā > 0.
2734+
"""
2735+
function _farkas_variable_dual(model::Optimizer, col::Cint)
2736+
nzcnt_p, surplus_p = Ref{Cint}(), Ref{Cint}()
2737+
cmatbeg = Vector{Cint}(undef, 2)
2738+
ret = CPXgetcols(
2739+
model.env,
2740+
model.lp,
2741+
nzcnt_p,
2742+
cmatbeg,
2743+
C_NULL,
2744+
C_NULL,
2745+
Cint(0),
2746+
surplus_p,
2747+
col,
2748+
col,
2749+
)
2750+
cmatind = Vector{Cint}(undef, -surplus_p[])
2751+
cmatval = Vector{Cdouble}(undef, -surplus_p[])
2752+
ret = CPXgetcols(
2753+
model.env,
2754+
model.lp,
2755+
nzcnt_p,
2756+
cmatbeg,
2757+
cmatind,
2758+
cmatval,
2759+
-surplus_p[],
2760+
surplus_p,
2761+
col,
2762+
col,
2763+
)
2764+
_check_ret(model, ret)
2765+
return sum(v * model.certificate[i + 1] for (i, v) in zip(cmatind, cmatval))
2766+
end
2767+
27232768
function MOI.get(
27242769
model::Optimizer,
27252770
attr::MOI.ConstraintDual,
@@ -2728,6 +2773,10 @@ function MOI.get(
27282773
_throw_if_optimize_in_progress(model, attr)
27292774
MOI.check_result_index_bounds(model, attr)
27302775
col = Cint(column(model, c) - 1)
2776+
if model.has_dual_certificate
2777+
dual = -_farkas_variable_dual(model, col)
2778+
return min(0.0, dual)
2779+
end
27312780
p = Ref{Cdouble}()
27322781
ret = CPXgetdj(model.env, model.lp, p, col, col)
27332782
_check_ret(model, ret)
@@ -2758,6 +2807,10 @@ function MOI.get(
27582807
_throw_if_optimize_in_progress(model, attr)
27592808
MOI.check_result_index_bounds(model, attr)
27602809
col = Cint(column(model, c) - 1)
2810+
if model.has_dual_certificate
2811+
dual = -_farkas_variable_dual(model, col)
2812+
return max(0.0, dual)
2813+
end
27612814
p = Ref{Cdouble}()
27622815
ret = CPXgetdj(model.env, model.lp, p, col, col)
27632816
_check_ret(model, ret)
@@ -2783,25 +2836,16 @@ end
27832836
function MOI.get(
27842837
model::Optimizer,
27852838
attr::MOI.ConstraintDual,
2786-
c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}},
2787-
)
2788-
_throw_if_optimize_in_progress(model, attr)
2789-
MOI.check_result_index_bounds(model, attr)
2790-
col = Cint(column(model, c) - 1)
2791-
p = Ref{Cdouble}()
2792-
ret = CPXgetdj(model.env, model.lp, p, col, col)
2793-
_check_ret(model, ret)
2794-
return _dual_multiplier(model) * p[]
2795-
end
2796-
2797-
function MOI.get(
2798-
model::Optimizer,
2799-
attr::MOI.ConstraintDual,
2800-
c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}},
2839+
c::MOI.ConstraintIndex{
2840+
MOI.SingleVariable, <:Union{MOI.Interval{Float64}, MOI.EqualTo{Float64}}
2841+
},
28012842
)
28022843
_throw_if_optimize_in_progress(model, attr)
28032844
MOI.check_result_index_bounds(model, attr)
28042845
col = Cint(column(model, c) - 1)
2846+
if model.has_dual_certificate
2847+
return -_farkas_variable_dual(model, col)
2848+
end
28052849
p = Ref{Cdouble}()
28062850
ret = CPXgetdj(model.env, model.lp, p, col, col)
28072851
_check_ret(model, ret)
@@ -2834,6 +2878,9 @@ function MOI.get(
28342878
# https://www.ibm.com/support/knowledgecenter/SSSA5P_12.10.0/ilog.odms.cplex.help/CPLEX/UsrMan/topics/cont_optim/qcp/17_QCP_duals.html
28352879
_throw_if_optimize_in_progress(model, attr)
28362880
MOI.check_result_index_bounds(model, attr)
2881+
if model.has_dual_certificate
2882+
error("Infeasibility certificate not available for $(c)")
2883+
end
28372884
# The derivative of a quadratic f(x) = x^TQx + a^Tx + b <= 0 is
28382885
# ∇f(x) = Q^Tx + Qx + a
28392886
# The dual is undefined if x is at the point of the cone. This can only be

‎test/MathOptInterface/MOI_wrapper.jl‎

Lines changed: 191 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -67,13 +67,7 @@ function test_modificationtest()
6767
end
6868

6969
function test_contlineartest()
70-
MOIT.contlineartest(BRIDGED_OPTIMIZER, CONFIG, [
71-
# TODO(odow): This test requests the infeasibility certificate of a
72-
# variable bound.
73-
"linear12"
74-
])
75-
76-
MOIT.linear12test(OPTIMIZER, MOIT.TestConfig(infeas_certificates=false))
70+
MOIT.contlineartest(BRIDGED_OPTIMIZER, CONFIG)
7771
end
7872

7973
function test_intlineartest()
@@ -542,6 +536,196 @@ function test_ZeroOne_INTERVAL()
542536
@test tmp[] == 2.0
543537
end
544538

539+
function test_farkas_dual_min()
540+
model = CPLEX.Optimizer()
541+
MOI.set(model, MOI.Silent(), true)
542+
MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0)
543+
x = MOI.add_variables(model, 2)
544+
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
545+
MOI.set(
546+
model,
547+
MOI.ObjectiveFunction{MOI.SingleVariable}(),
548+
MOI.SingleVariable(x[1]),
549+
)
550+
clb = MOI.add_constraint.(
551+
model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)
552+
)
553+
c = MOI.add_constraint(
554+
model,
555+
MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
556+
MOI.LessThan(-1.0),
557+
)
558+
MOI.optimize!(model)
559+
@test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
560+
@test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
561+
clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
562+
c_dual = MOI.get(model, MOI.ConstraintDual(), c)
563+
@show clb_dual, c_dual
564+
@test clb_dual[1] > 1e-6
565+
@test clb_dual[2] > 1e-6
566+
@test c_dual[1] < -1e-6
567+
@test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
568+
@test clb_dual[2] ≈ -c_dual atol = 1e-6
569+
end
570+
571+
function test_farkas_dual_min_interval()
572+
model = CPLEX.Optimizer()
573+
MOI.set(model, MOI.Silent(), true)
574+
MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0)
575+
x = MOI.add_variables(model, 2)
576+
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
577+
MOI.set(
578+
model,
579+
MOI.ObjectiveFunction{MOI.SingleVariable}(),
580+
MOI.SingleVariable(x[1]),
581+
)
582+
clb = MOI.add_constraint.(
583+
model, MOI.SingleVariable.(x), MOI.Interval(0.0, 10.0)
584+
)
585+
c = MOI.add_constraint(
586+
model,
587+
MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
588+
MOI.LessThan(-1.0),
589+
)
590+
MOI.optimize!(model)
591+
@test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
592+
@test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
593+
clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
594+
c_dual = MOI.get(model, MOI.ConstraintDual(), c)
595+
@show clb_dual, c_dual
596+
@test clb_dual[1] > 1e-6
597+
@test clb_dual[2] > 1e-6
598+
@test c_dual[1] < -1e-6
599+
@test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
600+
@test clb_dual[2] ≈ -c_dual atol = 1e-6
601+
end
602+
603+
function test_farkas_dual_min_equalto()
604+
model = CPLEX.Optimizer()
605+
MOI.set(model, MOI.Silent(), true)
606+
MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0)
607+
x = MOI.add_variables(model, 2)
608+
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
609+
MOI.set(
610+
model,
611+
MOI.ObjectiveFunction{MOI.SingleVariable}(),
612+
MOI.SingleVariable(x[1]),
613+
)
614+
clb = MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.EqualTo(0.0))
615+
c = MOI.add_constraint(
616+
model,
617+
MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
618+
MOI.LessThan(-1.0),
619+
)
620+
MOI.optimize!(model)
621+
@test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
622+
@test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
623+
clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
624+
c_dual = MOI.get(model, MOI.ConstraintDual(), c)
625+
@show clb_dual, c_dual
626+
@test clb_dual[1] > 1e-6
627+
@test clb_dual[2] > 1e-6
628+
@test c_dual[1] < -1e-6
629+
@test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
630+
@test clb_dual[2] ≈ -c_dual atol = 1e-6
631+
end
632+
633+
function test_farkas_dual_min_ii()
634+
model = CPLEX.Optimizer()
635+
MOI.set(model, MOI.Silent(), true)
636+
MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0)
637+
x = MOI.add_variables(model, 2)
638+
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
639+
MOI.set(
640+
model,
641+
MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
642+
MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0),
643+
)
644+
clb = MOI.add_constraint.(
645+
model, MOI.SingleVariable.(x), MOI.LessThan(0.0)
646+
)
647+
c = MOI.add_constraint(
648+
model,
649+
MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
650+
MOI.LessThan(-1.0),
651+
)
652+
MOI.optimize!(model)
653+
@test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
654+
@test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
655+
clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
656+
c_dual = MOI.get(model, MOI.ConstraintDual(), c)
657+
@show clb_dual, c_dual
658+
@test clb_dual[1] < -1e-6
659+
@test clb_dual[2] < -1e-6
660+
@test c_dual[1] < -1e-6
661+
@test clb_dual[1] ≈ 2 * c_dual atol = 1e-6
662+
@test clb_dual[2] ≈ c_dual atol = 1e-6
663+
end
664+
665+
function test_farkas_dual_max()
666+
model = CPLEX.Optimizer()
667+
MOI.set(model, MOI.Silent(), true)
668+
MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0)
669+
x = MOI.add_variables(model, 2)
670+
MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
671+
MOI.set(
672+
model,
673+
MOI.ObjectiveFunction{MOI.SingleVariable}(),
674+
MOI.SingleVariable(x[1]),
675+
)
676+
clb = MOI.add_constraint.(
677+
model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)
678+
)
679+
c = MOI.add_constraint(
680+
model,
681+
MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
682+
MOI.LessThan(-1.0),
683+
)
684+
MOI.optimize!(model)
685+
@test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
686+
@test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
687+
clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
688+
c_dual = MOI.get(model, MOI.ConstraintDual(), c)
689+
@show clb_dual, c_dual
690+
@test clb_dual[1] > 1e-6
691+
@test clb_dual[2] > 1e-6
692+
@test c_dual[1] < -1e-6
693+
@test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
694+
@test clb_dual[2] ≈ -c_dual atol = 1e-6
695+
end
696+
697+
function test_farkas_dual_max_ii()
698+
model = CPLEX.Optimizer()
699+
MOI.set(model, MOI.Silent(), true)
700+
MOI.set(model, MOI.RawParameter("CPX_PARAM_REDUCE"), 0)
701+
x = MOI.add_variables(model, 2)
702+
MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
703+
MOI.set(
704+
model,
705+
MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
706+
MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0),
707+
)
708+
clb = MOI.add_constraint.(
709+
model, MOI.SingleVariable.(x), MOI.LessThan(0.0)
710+
)
711+
c = MOI.add_constraint(
712+
model,
713+
MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
714+
MOI.LessThan(-1.0),
715+
)
716+
MOI.optimize!(model)
717+
@test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
718+
@test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
719+
clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
720+
c_dual = MOI.get(model, MOI.ConstraintDual(), c)
721+
@show clb_dual, c_dual
722+
@test clb_dual[1] < -1e-6
723+
@test clb_dual[2] < -1e-6
724+
@test c_dual[1] < -1e-6
725+
@test clb_dual[1] ≈ 2 * c_dual atol = 1e-6
726+
@test clb_dual[2] ≈ c_dual atol = 1e-6
727+
end
728+
545729
end # module TestMOIwrapper
546730

547731
runtests(TestMOIwrapper)

0 commit comments

Comments
 (0)