diff --git a/docs/make.jl b/docs/make.jl index a0139c9511..6451813a65 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -37,6 +37,7 @@ const _PAGES = [ ], "Background" => [ "background/duality.md", + "background/infeasibility_certificates.md", "background/naming_conventions.md", ], "API Reference" => [ diff --git a/docs/src/background/infeasibility_certificates.md b/docs/src/background/infeasibility_certificates.md new file mode 100644 index 0000000000..1219666f05 --- /dev/null +++ b/docs/src/background/infeasibility_certificates.md @@ -0,0 +1,160 @@ +```@meta +CurrentModule = MathOptInterface +DocTestSetup = quote + using MathOptInterface + const MOI = MathOptInterface +end +DocTestFilters = [r"MathOptInterface|MOI"] +``` + +# Infeasibility certificates + +When given a conic problem that is infeasible or unbounded, some solvers can +produce a certificate of infeasibility. This page explains what a certificate of +infeasibility is, and the related conventions that MathOptInterface adopts. + +## Conic duality + +MathOptInterface uses conic duality to define infeasibility certificates. A full +explanation is given in the section [Duality](@ref), but here is a brief +overview. + +### Minimization problems + +For a minimization problem in geometric conic form, the primal is: +```math +\begin{align} +& \min_{x \in \mathbb{R}^n} & a_0^\top x + b_0 +\\ +& \;\;\text{s.t.} & A_i x + b_i & \in \mathcal{C}_i & i = 1 \ldots m, +\end{align} +``` +and the dual is a maximization problem in standard conic form: +```math +\begin{align} +& \max_{y_1, \ldots, y_m} & -\sum_{i=1}^m b_i^\top y_i + b_0 +\\ +& \;\;\text{s.t.} & a_0 - \sum_{i=1}^m A_i^\top y_i & = 0 +\\ +& & y_i & \in \mathcal{C}_i^* & i = 1 \ldots m, +\end{align} +``` +where each ``\mathcal{C}_i`` is a closed convex cone and ``\mathcal{C}_i^*`` is +its dual cone. + +### Maximization problems + +For a maximization problem in geometric conic form, the primal is: +```math +\begin{align} +& \max_{x \in \mathbb{R}^n} & a_0^\top x + b_0 +\\ +& \;\;\text{s.t.} & A_i x + b_i & \in \mathcal{C}_i & i = 1 \ldots m, +\end{align} +``` +and the dual is a minimization problem in standard conic form: +```math +\begin{align} +& \min_{y_1, \ldots, y_m} & \sum_{i=1}^m b_i^\top y_i + b_0 +\\ +& \;\;\text{s.t.} & a_0 + \sum_{i=1}^m A_i^\top y_i & = 0 +\\ +& & y_i & \in \mathcal{C}_i^* & i = 1 \ldots m. +\end{align} +``` + +## Unbounded problems + +A problem is unbounded if and only if: + 1. there exists a feasible primal solution + 2. the dual is infeasible. + +A feasible primal solution—if one exists—can be obtained by setting +[`ObjectiveSense`](@ref) to `FEASIBILITY_SENSE` before optimizing. Therefore, +most solvers terminate after they prove the dual is infeasible via a certificate +of dual infeasibility, but _before_ they have found a feasible primal solution. +This is also the reason that MathOptInterface defines the `DUAL_INFEASIBLE` +status instead of `UNBOUNDED`. + +A certificate of dual infeasibility is an improving ray of the primal problem. +That is, there exists some vector ``d`` such that for all ``\eta > 0``: +```math +A_i (x + \eta d) + b_i \in \mathcal{C}_i,\ \ i = 1 \ldots m, +``` +and (for minimization problems): +```math +a_0^\top (x + \eta d) + b_0 < a_0^\top x + b_0, +``` +for any feasible point ``x``. The latter simplifies to ``a_0^\top d < 0``. For +maximization problems, the inequality is reversed, so that ``a_0^\top d > 0``. + +If the solver has found a certificate of dual infeasibility: + + * [`TerminationStatus`](@ref) must be `DUAL_INFEASIBLE` + * [`PrimalStatus`](@ref) must be `INFEASIBILITY_CERTIFICATE` + * [`VariablePrimal`](@ref) must be the corresponding value of ``d`` + * [`ConstraintPrimal`](@ref) must be the corresponding value of ``A_i d`` + * [`ObjectiveValue`](@ref) must be the value ``a_0^\top d``. Note that this is + the value of the objective function at `d`, ignoring the constant `b_0`. + +!!! note + The choice of whether to scale the ray ``d`` to have magnitude `1` is left + to the solver. + +## Infeasible problems + +A certificate of primal infeasibility is an improving ray of the dual problem. +However, because infeasibility is independent of the objective function, we +first homogenize the primal problem by removing its objective. + +For a minimization problem, a dual improving ray is some vector ``d`` such that +for all ``\eta > 0``: +```math +\begin{align} +-\sum_{i=1}^m A_i^\top (y_i + \eta d_i) & = 0 \\ +(y_i + \eta d_i) & \in \mathcal{C}_i^* & i = 1 \ldots m, +\end{align} +``` +and: +```math +-\sum_{i=1}^m b_i^\top (y_i + \eta d_i) > -\sum_{i=1}^m b_i^\top y_i, +``` +for any feasible dual solution ``y``. The latter simplifies to +``-\sum_{i=1}^m b_i^\top d_i > 0``. For a maximization problem, the inequality +is ``\sum_{i=1}^m b_i^\top d_i < 0``. (Note that these are the same inequality, +modulo a `-` sign.) + +If the solver has found a certificate of primal infeasibility: + + * [`TerminationStatus`](@ref) must be `INFEASIBLE` + * [`DualStatus`](@ref) must be `INFEASIBILITY_CERTIFICATE` + * [`ConstraintDual`](@ref) must be the corresponding value of ``d`` + * [`DualObjectiveValue`](@ref) must be the value + ``-\sum_{i=1}^m b_i^\top d_i`` for minimization problems and + ``\sum_{i=1}^m b_i^\top d_i`` for maximization problems. + +!!! note + The choice of whether to scale the ray ``d`` to have magnitude `1` is left + to the solver. + +### Infeasibility certificates of variable bounds + +Many linear solvers (e.g., Gurobi) do not provide explicit access to the primal +infeasibility certificate of a variable bound. However, given a set of linear +constraints: +```math +\begin{align} +l_A \le A x \le u_A \\ +l_x \le x \le u_x, +\end{align} +``` +the primal certificate of the variable bounds can be computed using the primal +certificate associated with the affine constraints, ``d``. (Note that ``d`` will +have one element for each row of the ``A`` matrix, and that some or all of the +elements in the vectors ``l_A`` and ``u_A`` may be ``\pm \infty``. If both +``l_A`` and ``u_A`` are finite for some row, the corresponding element in ``d` + must be `0`.) + +Given ``d``, compute ``\bar{d} = d^\top A``. If the bound is finite, a +certificate for the lower variable bound of ``x_i`` is ``\max\{\bar{d}_i, 0\}``, +and a certificate for the upper variable bound is ``\min\{\bar{d}_i, 0\}``. diff --git a/src/Bridges/Constraint/scalarize.jl b/src/Bridges/Constraint/scalarize.jl index c0886af963..2982db3baf 100644 --- a/src/Bridges/Constraint/scalarize.jl +++ b/src/Bridges/Constraint/scalarize.jl @@ -175,7 +175,6 @@ function MOI.set( bridge::ScalarizeBridge, value, ) - # TODO do no add constant if the primal status is a ray like in Vectorize MOI.set.(model, attr, bridge.scalar_constraints, value .- bridge.constants) return end diff --git a/src/Test/test_infeasibility_certificates.jl b/src/Test/test_infeasibility_certificates.jl new file mode 100644 index 0000000000..560dbb8c99 --- /dev/null +++ b/src/Test/test_infeasibility_certificates.jl @@ -0,0 +1,187 @@ +for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) + o = iszero(offset) ? "" : "_offset" + f_unbd_name = Symbol("test_unbounded_$(Symbol(sense))$(o)") + f_infeas_name = Symbol("test_infeasible_$(Symbol(sense))$(o)") + f_infeas_affine_name = Symbol("test_infeasible_affine_$(Symbol(sense))$(o)") + @eval begin + """ + Test the infeasibility certificates of an unbounded model. + """ + function $(f_unbd_name)(model::MOI.ModelLike, config::Config) + @requires _supports(config, MOI.optimize!) + x = MOI.add_variable(model) + MOI.set(model, MOI.ObjectiveSense(), $sense) + f = 2.2 * x + $offset + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + c = if $sense == MOI.MIN_SENSE + MOI.add_constraint(model, 1.3 * x, MOI.LessThan(1.1)) + else + MOI.add_constraint(model, 1.3 * x, MOI.GreaterThan(1.1)) + end + MOI.optimize!(model) + @requires( + MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE, + ) + @requires( + MOI.get(model, MOI.PrimalStatus()) == + MOI.INFEASIBILITY_CERTIFICATE, + ) + d = MOI.get(model, MOI.VariablePrimal(), x) + obj = MOI.get(model, MOI.ObjectiveValue()) + if $sense == MOI.MIN_SENSE + @test d < -config.atol + @test obj < -config.atol + else + @test d > config.atol + @test obj > config.atol + end + @test isapprox(2.2 * d, obj, config) + Ad = MOI.get(model, MOI.ConstraintPrimal(), c) + @test isapprox(1.3 * d, Ad, config) + return + end + + function setup_test( + ::typeof($(f_unbd_name)), + mock::MOIU.MockOptimizer, + config::Config, + ) + MOIU.set_mock_optimize!( + mock, + (mock::MOIU.MockOptimizer) -> begin + MOIU.mock_optimize!( + mock, + MOI.DUAL_INFEASIBLE, + ( + MOI.INFEASIBILITY_CERTIFICATE, + [$sense == MOI.MIN_SENSE ? -1.0 : 1.0], + ), + MOI.NO_SOLUTION, + ) + end, + ) + return + end + + """ + Test the infeasibility certificates of an infeasible model. + """ + function $(f_infeas_name)(model::MOI.ModelLike, config::Config) + @requires _supports(config, MOI.optimize!) + x = MOI.add_variable(model) + MOI.set(model, MOI.ObjectiveSense(), $sense) + f = 2.2 * x + $offset + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + cl = MOI.add_constraint(model, x, MOI.GreaterThan(2.5)) + cu = MOI.add_constraint(model, 1.0 * x, MOI.LessThan(1.4)) + MOI.optimize!(model) + @requires( + MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE, + ) + @requires( + MOI.get(model, MOI.DualStatus()) == + MOI.INFEASIBILITY_CERTIFICATE, + ) + dl = MOI.get(model, MOI.ConstraintDual(), cl) + du = MOI.get(model, MOI.ConstraintDual(), cu) + obj = MOI.get(model, MOI.DualObjectiveValue()) + if $sense == MOI.MIN_SENSE + @test obj > config.atol + @test isapprox(-(1.4 * dl + 2.5 * du), obj, config) + else + @test obj < -config.atol + @test isapprox(1.4 * dl + 2.5 * du, obj, config) + end + @test dl > config.atol + @test du < -config.atol + @test isapprox(dl + du, 0.0, config) + return + end + + function setup_test( + ::typeof($(f_infeas_name)), + mock::MOIU.MockOptimizer, + config::Config, + ) + MOIU.set_mock_optimize!( + mock, + (mock::MOIU.MockOptimizer) -> begin + MOIU.mock_optimize!( + mock, + MOI.INFEASIBLE, + MOI.NO_SOLUTION, + MOI.INFEASIBILITY_CERTIFICATE, + (MOI.VariableIndex, MOI.GreaterThan{Float64}) => + [1.0], + ( + MOI.ScalarAffineFunction{Float64}, + MOI.LessThan{Float64}, + ) => [-1.0], + ) + end, + ) + return + end + + """ + Test the infeasibility certificates of an infeasible model with affine + constraints. + """ + function $(f_infeas_affine_name)(model::MOI.ModelLike, config::Config) + @requires _supports(config, MOI.optimize!) + x = MOI.add_variable(model) + MOI.set(model, MOI.ObjectiveSense(), $sense) + f = 2.2 * x + $offset + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + cl = MOI.add_constraint(model, 1.0 * x, MOI.GreaterThan(2.5)) + cu = MOI.add_constraint(model, x, MOI.LessThan(1.4)) + MOI.optimize!(model) + @requires( + MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE, + ) + @requires( + MOI.get(model, MOI.DualStatus()) == + MOI.INFEASIBILITY_CERTIFICATE, + ) + dl = MOI.get(model, MOI.ConstraintDual(), cl) + du = MOI.get(model, MOI.ConstraintDual(), cu) + obj = MOI.get(model, MOI.DualObjectiveValue()) + if $sense == MOI.MIN_SENSE + @test obj > config.atol + @test isapprox(-(1.4 * dl + 2.5 * du), obj, config) + else + @test obj < -config.atol + @test isapprox(1.4 * dl + 2.5 * du, obj, config) + end + @test dl > config.atol + @test du < -config.atol + @test isapprox(dl + du, 0.0, config) + return + end + + function setup_test( + ::typeof($(f_infeas_affine_name)), + mock::MOIU.MockOptimizer, + config::Config, + ) + MOIU.set_mock_optimize!( + mock, + (mock::MOIU.MockOptimizer) -> begin + MOIU.mock_optimize!( + mock, + MOI.INFEASIBLE, + MOI.NO_SOLUTION, + MOI.INFEASIBILITY_CERTIFICATE, + ( + MOI.ScalarAffineFunction{Float64}, + MOI.GreaterThan{Float64}, + ) => [1.0], + (MOI.VariableIndex, MOI.LessThan{Float64}) => + [-1.0], + ) + end, + ) + return + end + end +end diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl index 8aa68e4173..8098e8a262 100644 --- a/src/Utilities/results.jl +++ b/src/Utilities/results.jl @@ -19,11 +19,16 @@ the `ObjectiveFunction` value. function get_fallback(model::MOI.ModelLike, attr::MOI.ObjectiveValue) F = MOI.get(model, MOI.ObjectiveFunctionType()) f = MOI.get(model, MOI.ObjectiveFunction{F}()) - # TODO do not include constant if primal solution is a ray - return eval_variables( + obj = eval_variables( vi -> MOI.get(model, MOI.VariablePrimal(attr.result_index), vi), f, ) + if is_ray(MOI.get(model, MOI.PrimalStatus())) + # Dual infeasibiltiy certificates do not include the primal objective + # constant. + obj -= MOI.constant(f, typeof(obj)) + end + return obj end function constraint_constant( @@ -173,11 +178,13 @@ function get_fallback( idx::MOI.ConstraintIndex, ) f = MOI.get(model, MOI.ConstraintFunction(), idx) - # TODO do not include constant if primal solution is a ray - return eval_variables( - vi -> MOI.get(model, MOI.VariablePrimal(attr.result_index), vi), - f, - ) + c = eval_variables(f) do vi + return MOI.get(model, MOI.VariablePrimal(attr.result_index), vi) + end + if is_ray(MOI.get(model, MOI.PrimalStatus())) + c -= MOI.constant(f, typeof(c)) + end + return c end ################ Constraint Dual for Variable-wise constraints #################