diff --git a/src/MOI/MOI_wrapper.jl b/src/MOI/MOI_wrapper.jl index 135d73f6..9ad4a2c0 100644 --- a/src/MOI/MOI_wrapper.jl +++ b/src/MOI/MOI_wrapper.jl @@ -2787,15 +2787,63 @@ function MOI.get( return _dual_multiplier(model) * p[] end -# function MOI.get( -# 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) -# pi = model.cached_solution.quadratic_dual[_info(model, c).row] -# return _dual_multiplier(model) * pi -# end +function MOI.get( + model::Optimizer, + attr::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, <:Any}, +) + # For more information on QCP duals, see + # 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) + # 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 + # checked to numeric tolerances. We use `cone_top_tol`. + cone_top, cone_top_tol = true, 1e-6 + x = zeros(length(model.variable_info)) + ret = CPXgetx(model.env, model.lp, x, 0, length(x) - 1) + _check_ret(model, ret) + ∇f = zeros(length(x)) + a_i, a_v, qrow, qcol, qval = _CPXgetqconstr(model, c) + for (i, j, v) in zip(qrow, qcol, qval) + ∇f[i + 1] += v * x[j + 1] + ∇f[j + 1] += v * x[i + 1] + if abs(x[i + 1]) > cone_top_tol || abs(x[j + 1]) > cone_top_tol + cone_top = false + end + end + for (i, v) in zip(a_i, a_v) + ∇f[i + 1] += v + if abs(x[i + 1]) > cone_top_tol + cone_top = false + end + end + # TODO(odow): if at top of cone (x = 0) dual multiplier is ill-formed. + if cone_top + return NaN + end + qind = Cint(_info(model, c).row - 1) + nz_p, surplus_p = Ref{Cint}(), Ref{Cint}() + CPXgetqconstrdslack( + model.env, model.lp, qind, nz_p, C_NULL, C_NULL, 0, surplus_p + ) + ind = Vector{Cint}(undef, -surplus_p[]) + val = Vector{Cdouble}(undef, -surplus_p[]) + ret = CPXgetqconstrdslack( + model.env, model.lp, qind, nz_p, ind, val, -surplus_p[], surplus_p + ) + _check_ret(model, ret) + ∇f_max, ∇f_i = findmax(abs.(∇f)) + if ∇f_max > cone_top_tol + for (i, v) in zip(ind, val) + if i + 1 == ∇f_i + return _dual_multiplier(model) * v / ∇f[∇f_i] + end + end + end + return 0.0 +end function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue) _throw_if_optimize_in_progress(model, attr) @@ -3464,3 +3512,27 @@ function MOI.set( return end +function MOI.get( + model::Optimizer, + ::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.SecondOrderCone}, +) + f = MOI.get(model, MOI.ConstraintFunction(), c) + qind = Cint(_info(model, c).row - 1) + surplus_p = Ref{Cint}() + CPXgetqconstrdslack( + model.env, model.lp, qind, C_NULL, C_NULL, C_NULL, 0, surplus_p + ) + ind = Vector{Cint}(undef, -surplus_p[]) + val = Vector{Cdouble}(undef, -surplus_p[]) + ret = CPXgetqconstrdslack( + model.env, model.lp, qind, C_NULL, ind, val, -surplus_p[], surplus_p + ) + _check_ret(model, ret) + slack = zeros(length(model.variable_info)) + for (i, v) in zip(ind, val) + slack[i + 1] += v + end + z = _dual_multiplier(model) + return [z * slack[_info(model, v).column] for v in f.variables] +end diff --git a/test/MathOptInterface/MOI_wrapper.jl b/test/MathOptInterface/MOI_wrapper.jl index b15416ba..3f8c5635 100644 --- a/test/MathOptInterface/MOI_wrapper.jl +++ b/test/MathOptInterface/MOI_wrapper.jl @@ -12,14 +12,11 @@ const CONFIG = MOIT.TestConfig() const OPTIMIZER = CPLEX.Optimizer() MOI.set(OPTIMIZER, MOI.Silent(), true) -const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(OPTIMIZER, Float64) +# Turn off presolve reductions so CPLEX will generate infeasibility +# certificates. +MOI.set(OPTIMIZER, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) -const CERTIFICATE_OPTIMIZER = CPLEX.Optimizer() -MOI.set(CERTIFICATE_OPTIMIZER, MOI.Silent(), true) -MOI.set(CERTIFICATE_OPTIMIZER, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) -MOI.set(CERTIFICATE_OPTIMIZER, MOI.RawParameter("CPX_PARAM_PRELINEAR"), 0) -const BRIDGED_CERTIFICATE_OPTIMIZER = - MOI.Bridges.full_bridge_optimizer(CERTIFICATE_OPTIMIZER, Float64) +const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(OPTIMIZER, Float64) function test_basic_constraint_tests() MOIT.basic_constraint_tests(BRIDGED_OPTIMIZER, CONFIG; exclude = [ @@ -62,12 +59,7 @@ function test_unittest() # TODO(odow): bug! We can't delete a vector of variables if one is in # a second order cone. "delete_soc_variables", - - # CPLEX returns INFEASIBLE_OR_UNBOUNDED without extra parameters. - # See below for the test. - "solve_unbounded_model", ]) - MOIT.solve_unbounded_model(CERTIFICATE_OPTIMIZER, CONFIG) end function test_modificationtest() @@ -76,18 +68,11 @@ end function test_contlineartest() MOIT.contlineartest(BRIDGED_OPTIMIZER, CONFIG, [ - # These tests require extra parameters to be set. - "linear8a", "linear8b", "linear8c", - # TODO(odow): This test requests the infeasibility certificate of a # variable bound. "linear12" ]) - MOIT.linear8atest(CERTIFICATE_OPTIMIZER, CONFIG) - MOIT.linear8btest(CERTIFICATE_OPTIMIZER, CONFIG) - MOIT.linear8ctest(CERTIFICATE_OPTIMIZER, CONFIG) - MOIT.linear12test(OPTIMIZER, MOIT.TestConfig(infeas_certificates=false)) end @@ -101,35 +86,34 @@ end function test_contquadratictest() MOIT.contquadratictest( - BRIDGED_CERTIFICATE_OPTIMIZER, - # TODO(odow): duals for quadratic problems. - MOIT.TestConfig(duals = false, atol = 1e-3, rtol = 1e-3), + BRIDGED_OPTIMIZER, + MOIT.TestConfig(atol = 1e-3, rtol = 1e-3), ["ncqcp"], # CPLEX doesn't support non-convex problems ) end function test_contconic() - MOIT.lintest(BRIDGED_OPTIMIZER, CONFIG, [ - # These tests require extra parameters to be set. - "lin3", "lin4" - ]) - - MOIT.lin3test(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG) - MOIT.lin4test(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG) - - # TODO(odow): duals for SOC constraints. - soc_config = MOIT.TestConfig(duals = false, atol=5e-3) + MOIT.lintest(BRIDGED_OPTIMIZER, CONFIG) - MOIT.soctest(BRIDGED_OPTIMIZER, soc_config, [ - "soc3" - ]) + soc_config = MOIT.TestConfig(atol=5e-3) + # TODO(odow): investigate why infeasibility certificates not generated for + # SOC. + MOIT.soctest(BRIDGED_OPTIMIZER, soc_config, ["soc3"]) MOIT.soc3test( BRIDGED_OPTIMIZER, - MOIT.TestConfig(duals = false, infeas_certificates = false, atol = 1e-3) + MOIT.TestConfig(atol = 1e-3, infeas_certificates = false) ) - MOIT.rsoctest(BRIDGED_OPTIMIZER, soc_config) + MOIT.rsoctest(BRIDGED_OPTIMIZER, soc_config, ["rotatedsoc2"]) + MOIT.rotatedsoc2test( + BRIDGED_OPTIMIZER, + # Need for `duals = false` fixed by https://github.com/jump-dev/MathOptInterface.jl/pull/1171 + # Remove in a future MOI 0.9.18+ release. + MOIT.TestConfig( + atol = 1e-3, infeas_certificates = false, duals = false + ), + ) MOIT.geomeantest(BRIDGED_OPTIMIZER, soc_config) end