From b8bc61359622f98e4cd7ecf14b11a3ec6daefcce Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 16 Nov 2021 14:53:43 +1300 Subject: [PATCH 01/11] [docs] add background on infeasibility certificates Add tests for dual infeasibility cert --- docs/make.jl | 1 + .../background/infeasibility_certificates.md | 153 ++++++++++++++++++ src/Test/test_infeasibility_certificates.jl | 115 +++++++++++++ 3 files changed, 269 insertions(+) create mode 100644 docs/src/background/infeasibility_certificates.md create mode 100644 src/Test/test_infeasibility_certificates.jl 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..1e9bd842c8 --- /dev/null +++ b/docs/src/background/infeasibility_certificates.md @@ -0,0 +1,153 @@ +```@meta +CurrentModule = MathOptInterface +DocTestSetup = quote + using MathOptInterface + const MOI = MathOptInterface +end +DocTestFilters = [r"MathOptInterface|MOI"] +``` + +# Infeasibility certificates + +When given an infeasible or unbounded model, some (conic) solvers can produce a +certificate or infeasibility. This page explains what a certificate of +infeasibility is for unbounded and infeasible models. + +## Conic duality + +MathOptInterface uses conic duality, which we use to define infeasibility +certificates. A full explanation is given in the section [Duality](@ref), but +here is a brief overview. + +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. + +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 models + +If a model is unbounded, two things are true: + 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` and then 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`` + * [`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 models + +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. + +Therefore, 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``. Note that conic duality is defined such that +this inequality is independent of whether the objective sense is minimization or +maximization. + +If the solver has found a certificate of 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``. + +!!! 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 dual +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 dual certificate of the variable bounds can be computed using the +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 +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``. 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/Test/test_infeasibility_certificates.jl b/src/Test/test_infeasibility_certificates.jl new file mode 100644 index 0000000000..517883ff1b --- /dev/null +++ b/src/Test/test_infeasibility_certificates.jl @@ -0,0 +1,115 @@ +for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) + o = iszero(offset) ? "" : "_offset" + f_name = Symbol("test_unbounded_$(Symbol(sense))$(o)") + @eval begin + function $(f_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) + 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) + return + end + function setup_test( + ::typeof($(f_name)), + mock::MOIU.MockOptimizer, + config::Config, + ) + x = [$sense == MOI.MIN_SENSE ? -1.0 : 1.0] + flag = mock.eval_objective_value + mock.eval_objective_value = false + MOIU.set_mock_optimize!( + mock, + (mock::MOIU.MockOptimizer) -> begin + MOI.set(mock, MOI.ObjectiveValue(), 2.2 * x[1]) + MOIU.mock_optimize!( + mock, + MOI.DUAL_INFEASIBLE, + (MOI.INFEASIBILITY_CERTIFICATE, x), + MOI.NO_SOLUTION, + ) + end, + ) + return () -> mock.eval_objective_value = flag + end + end +end + +for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 0.2) + o = iszero(offset) ? "" : "_offset" + f_name = Symbol("test_infeasible_$(Symbol(sense))$(o)") + @eval begin + function $(f_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, 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()) + @test dl > config.atol + @test du < -config.atol + @test obj < -config.atol + @test isapprox(dl + du, 0.0, config) + @test isapprox(1.4 * dl + 2.5 * du, obj, config) + return + end + function setup_test( + ::typeof($(f_name)), + mock::MOIU.MockOptimizer, + config::Config, + ) + flag_obj = mock.eval_dual_objective_value + flag_var = mock.eval_variable_constraint_dual + mock.eval_dual_objective_value = false + mock.eval_variable_constraint_dual = false + MOIU.set_mock_optimize!( + mock, + (mock::MOIU.MockOptimizer) -> begin + MOI.set(mock, MOI.DualObjectiveValue(), -1.1) + MOIU.mock_optimize!( + mock, + MOI.INFEASIBLE, + MOI.NO_SOLUTION, + MOI.INFEASIBILITY_CERTIFICATE, + (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [1.0], + (MOI.VariableIndex, MOI.LessThan{Float64}) => [-1.0], + ) + end, + ) + return () -> begin + mock.eval_dual_objective_value = flag_obj + mock.eval_variable_constraint_dual = flag_var + end + end + end +end From e800a943943452eaf2f5526065bf71bfbdfa3e8b Mon Sep 17 00:00:00 2001 From: odow Date: Wed, 17 Nov 2021 12:17:12 +1300 Subject: [PATCH 02/11] More updates --- src/Test/test_infeasibility_certificates.jl | 31 ++++++++++++--------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/Test/test_infeasibility_certificates.jl b/src/Test/test_infeasibility_certificates.jl index 517883ff1b..7ed537721d 100644 --- a/src/Test/test_infeasibility_certificates.jl +++ b/src/Test/test_infeasibility_certificates.jl @@ -1,8 +1,12 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) o = iszero(offset) ? "" : "_offset" - f_name = Symbol("test_unbounded_$(Symbol(sense))$(o)") + f_unbd_name = Symbol("test_unbounded_$(Symbol(sense))$(o)") + f_infeas_name = Symbol("test_infeasible_$(Symbol(sense))$(o)") @eval begin - function $(f_name)(model::MOI.ModelLike, config::Config) + """ + 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) @@ -28,8 +32,9 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) @test isapprox(2.2 * d, obj, config) return end + function setup_test( - ::typeof($(f_name)), + ::typeof($(f_unbd_name)), mock::MOIU.MockOptimizer, config::Config, ) @@ -50,14 +55,11 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) ) return () -> mock.eval_objective_value = flag end - end -end -for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 0.2) - o = iszero(offset) ? "" : "_offset" - f_name = Symbol("test_infeasible_$(Symbol(sense))$(o)") - @eval begin - function $(f_name)(model::MOI.ModelLike, config::Config) + """ + 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) @@ -83,8 +85,9 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 0.2) @test isapprox(1.4 * dl + 2.5 * du, obj, config) return end + function setup_test( - ::typeof($(f_name)), + ::typeof($(f_infeas_name)), mock::MOIU.MockOptimizer, config::Config, ) @@ -101,8 +104,10 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 0.2) MOI.INFEASIBLE, MOI.NO_SOLUTION, MOI.INFEASIBILITY_CERTIFICATE, - (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [1.0], - (MOI.VariableIndex, MOI.LessThan{Float64}) => [-1.0], + (MOI.VariableIndex, MOI.GreaterThan{Float64}) => + [1.0], + (MOI.VariableIndex, MOI.LessThan{Float64}) => + [-1.0], ) end, ) From 9913b19e2fb3b390b50c2a487f2d8dd83e72a08e Mon Sep 17 00:00:00 2001 From: odow Date: Wed, 17 Nov 2021 12:58:27 +1300 Subject: [PATCH 03/11] More updates --- .../background/infeasibility_certificates.md | 17 +++++------ src/Bridges/Constraint/scalarize.jl | 1 - src/Test/test_infeasibility_certificates.jl | 28 +++++++++---------- src/Utilities/results.jl | 9 ++++-- 4 files changed, 29 insertions(+), 26 deletions(-) diff --git a/docs/src/background/infeasibility_certificates.md b/docs/src/background/infeasibility_certificates.md index 1e9bd842c8..52d0c982bb 100644 --- a/docs/src/background/infeasibility_certificates.md +++ b/docs/src/background/infeasibility_certificates.md @@ -102,22 +102,22 @@ 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. -Therefore, a dual improving ray is some vector ``d`` such that for all -``\eta > 0``: +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 \\ +\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, +-\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``. Note that conic duality is defined such that -this inequality is independent of whether the objective sense is minimization or -maximization. +``-\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 infeasibility: @@ -125,7 +125,8 @@ If the solver has found a certificate of infeasibility: * [`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``. + ``-\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 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 index 7ed537721d..9989bea2e0 100644 --- a/src/Test/test_infeasibility_certificates.jl +++ b/src/Test/test_infeasibility_certificates.jl @@ -38,22 +38,21 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) mock::MOIU.MockOptimizer, config::Config, ) - x = [$sense == MOI.MIN_SENSE ? -1.0 : 1.0] - flag = mock.eval_objective_value - mock.eval_objective_value = false MOIU.set_mock_optimize!( mock, (mock::MOIU.MockOptimizer) -> begin - MOI.set(mock, MOI.ObjectiveValue(), 2.2 * x[1]) MOIU.mock_optimize!( mock, MOI.DUAL_INFEASIBLE, - (MOI.INFEASIBILITY_CERTIFICATE, x), + ( + MOI.INFEASIBILITY_CERTIFICATE, + [$sense == MOI.MIN_SENSE ? -1.0 : 1.0], + ), MOI.NO_SOLUTION, ) end, ) - return () -> mock.eval_objective_value = flag + return end """ @@ -78,11 +77,16 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) 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 obj < -config.atol @test isapprox(dl + du, 0.0, config) - @test isapprox(1.4 * dl + 2.5 * du, obj, config) return end @@ -91,14 +95,11 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) mock::MOIU.MockOptimizer, config::Config, ) - flag_obj = mock.eval_dual_objective_value flag_var = mock.eval_variable_constraint_dual - mock.eval_dual_objective_value = false mock.eval_variable_constraint_dual = false MOIU.set_mock_optimize!( mock, (mock::MOIU.MockOptimizer) -> begin - MOI.set(mock, MOI.DualObjectiveValue(), -1.1) MOIU.mock_optimize!( mock, MOI.INFEASIBLE, @@ -111,10 +112,7 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) ) end, ) - return () -> begin - mock.eval_dual_objective_value = flag_obj - mock.eval_variable_constraint_dual = flag_var - end + return () -> mock.eval_variable_constraint_dual = flag_var end end end diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl index 8aa68e4173..e758c85835 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 MOI.get(model, MOI.PrimalStatus()) == MOI.INFEASIBILITY_CERTIFICATE + # Dual infeasibiltiy certificates do not include the primal objective + # constant. + obj -= MOI.constant(f, Int) + end + return obj end function constraint_constant( From ac48fdec1165c6ff22fd6f9781e4897e626c7a64 Mon Sep 17 00:00:00 2001 From: odow Date: Wed, 17 Nov 2021 14:06:09 +1300 Subject: [PATCH 04/11] Update and more tests --- .../background/infeasibility_certificates.md | 35 +++++----- src/Test/test_infeasibility_certificates.jl | 65 +++++++++++++++++++ 2 files changed, 85 insertions(+), 15 deletions(-) diff --git a/docs/src/background/infeasibility_certificates.md b/docs/src/background/infeasibility_certificates.md index 52d0c982bb..4f869d3fdb 100644 --- a/docs/src/background/infeasibility_certificates.md +++ b/docs/src/background/infeasibility_certificates.md @@ -9,9 +9,9 @@ DocTestFilters = [r"MathOptInterface|MOI"] # Infeasibility certificates -When given an infeasible or unbounded model, some (conic) solvers can produce a -certificate or infeasibility. This page explains what a certificate of -infeasibility is for unbounded and infeasible models. +When given a conic problem that is infeasible or unbounded, some solvers can +produce a certificate or infeasibility. This page explains what a certificate of +infeasibility is, and the related conventions that MathOptInterface adopts. ## Conic duality @@ -19,6 +19,8 @@ MathOptInterface uses conic duality, which we use 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} @@ -40,6 +42,8 @@ and the dual is a maximization problem in standard conic form: 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} @@ -59,14 +63,14 @@ and the dual is a minimization problem in standard conic form: \end{align} ``` -## Unbounded models +## Unbounded problems -If a model is unbounded, two things are true: +If 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` and then optimizing. Therefore, +[`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` @@ -96,7 +100,7 @@ If the solver has found a certificate of dual infeasibility: The choice of whether to scale the ray ``d`` to have magnitude `1` is left to the solver. -## Infeasible models +## 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 @@ -119,7 +123,7 @@ for any feasible dual solution ``y``. The latter simplifies to 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 infeasibility: +If the solver has found a certificate of primal infeasibility: * [`TerminationStatus`](@ref) must be `INFEASIBLE` * [`DualStatus`](@ref) must be `INFEASIBILITY_CERTIFICATE` @@ -143,12 +147,13 @@ l_A \le A x \le u_A \\ l_x \le x \le u_x, \end{align} ``` -the dual certificate of the variable bounds can be computed using the -certificate associated with the affine constraints ``d``. (Note that ``d`` will +the dual certificate of the variable bounds can be computed using the +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 -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`.) +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``. 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\}``. +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/Test/test_infeasibility_certificates.jl b/src/Test/test_infeasibility_certificates.jl index 9989bea2e0..0e9ed4c1d9 100644 --- a/src/Test/test_infeasibility_certificates.jl +++ b/src/Test/test_infeasibility_certificates.jl @@ -2,6 +2,7 @@ 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. @@ -114,5 +115,69 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) ) return () -> mock.eval_variable_constraint_dual = flag_var 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) + g = 1.0 * x + cl = MOI.add_constraint(model, g, MOI.GreaterThan(2.5)) + cu = MOI.add_constraint(model, g, 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.ScalarAffineFunction{Float64}, + MOI.LessThan{Float64}, + ) => [-1.0], + ) + end, + ) + return + end end end From daece17c07a13eeed805c9fb43b86de7ed55ba84 Mon Sep 17 00:00:00 2001 From: odow Date: Wed, 17 Nov 2021 14:53:55 +1300 Subject: [PATCH 05/11] Updates --- src/Test/test_infeasibility_certificates.jl | 23 +++++++++------------ 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/Test/test_infeasibility_certificates.jl b/src/Test/test_infeasibility_certificates.jl index 0e9ed4c1d9..5a6364f195 100644 --- a/src/Test/test_infeasibility_certificates.jl +++ b/src/Test/test_infeasibility_certificates.jl @@ -66,7 +66,7 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) 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, x, MOI.LessThan(1.4)) + cu = MOI.add_constraint(model, 1.0 * x, MOI.LessThan(1.4)) MOI.optimize!(model) @requires( MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE, @@ -96,8 +96,6 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) mock::MOIU.MockOptimizer, config::Config, ) - flag_var = mock.eval_variable_constraint_dual - mock.eval_variable_constraint_dual = false MOIU.set_mock_optimize!( mock, (mock::MOIU.MockOptimizer) -> begin @@ -108,12 +106,14 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) MOI.INFEASIBILITY_CERTIFICATE, (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [1.0], - (MOI.VariableIndex, MOI.LessThan{Float64}) => - [-1.0], + ( + MOI.ScalarAffineFunction{Float64}, + MOI.LessThan{Float64}, + ) => [-1.0], ) end, ) - return () -> mock.eval_variable_constraint_dual = flag_var + return end """ @@ -126,9 +126,8 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) MOI.set(model, MOI.ObjectiveSense(), $sense) f = 2.2 * x + $offset MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) - g = 1.0 * x - cl = MOI.add_constraint(model, g, MOI.GreaterThan(2.5)) - cu = MOI.add_constraint(model, g, MOI.LessThan(1.4)) + 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, @@ -170,10 +169,8 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, ) => [1.0], - ( - MOI.ScalarAffineFunction{Float64}, - MOI.LessThan{Float64}, - ) => [-1.0], + (MOI.VariableIndex, MOI.LessThan{Float64}) => + [-1.0], ) end, ) From 2f368b7fa5ba0eca1a1b3a533d219a4faa255fb8 Mon Sep 17 00:00:00 2001 From: odow Date: Thu, 18 Nov 2021 08:47:08 +1300 Subject: [PATCH 06/11] Add ConstraintPrimal --- docs/src/background/infeasibility_certificates.md | 1 + src/Test/test_infeasibility_certificates.jl | 7 +++++++ src/Utilities/results.jl | 14 ++++++++------ 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/docs/src/background/infeasibility_certificates.md b/docs/src/background/infeasibility_certificates.md index 4f869d3fdb..c13e07fc91 100644 --- a/docs/src/background/infeasibility_certificates.md +++ b/docs/src/background/infeasibility_certificates.md @@ -93,6 +93,7 @@ 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`. diff --git a/src/Test/test_infeasibility_certificates.jl b/src/Test/test_infeasibility_certificates.jl index 5a6364f195..9eb1f81d90 100644 --- a/src/Test/test_infeasibility_certificates.jl +++ b/src/Test/test_infeasibility_certificates.jl @@ -13,6 +13,11 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) 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, @@ -31,6 +36,8 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) @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 diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl index e758c85835..b124d205c3 100644 --- a/src/Utilities/results.jl +++ b/src/Utilities/results.jl @@ -23,7 +23,7 @@ function get_fallback(model::MOI.ModelLike, attr::MOI.ObjectiveValue) vi -> MOI.get(model, MOI.VariablePrimal(attr.result_index), vi), f, ) - if MOI.get(model, MOI.PrimalStatus()) == MOI.INFEASIBILITY_CERTIFICATE + if is_ray(MOI.get(model, MOI.PrimalStatus())) # Dual infeasibiltiy certificates do not include the primal objective # constant. obj -= MOI.constant(f, Int) @@ -178,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, Int) + end + return c end ################ Constraint Dual for Variable-wise constraints ################# From 42bfdff6cfeacad9382b74cb9711449c453a536d Mon Sep 17 00:00:00 2001 From: odow Date: Thu, 18 Nov 2021 08:50:00 +1300 Subject: [PATCH 07/11] Fix formatting --- src/Test/test_infeasibility_certificates.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Test/test_infeasibility_certificates.jl b/src/Test/test_infeasibility_certificates.jl index 9eb1f81d90..560dbb8c99 100644 --- a/src/Test/test_infeasibility_certificates.jl +++ b/src/Test/test_infeasibility_certificates.jl @@ -13,7 +13,7 @@ for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0.0, 1.2) 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 + 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)) From 631cf0c0b8c2fa62d60224644eb1f8a387d66cdc Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Thu, 18 Nov 2021 11:34:07 +1300 Subject: [PATCH 08/11] Update infeasibility_certificates.md --- docs/src/background/infeasibility_certificates.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/background/infeasibility_certificates.md b/docs/src/background/infeasibility_certificates.md index c13e07fc91..4ba2edbac1 100644 --- a/docs/src/background/infeasibility_certificates.md +++ b/docs/src/background/infeasibility_certificates.md @@ -65,7 +65,7 @@ and the dual is a minimization problem in standard conic form: ## Unbounded problems -If a problem is unbounded, if and only if: +A problem is unbounded if and only if: 1. there exists a feasible primal solution 2. the dual is infeasible. @@ -139,7 +139,7 @@ If the solver has found a certificate of primal infeasibility: ### Infeasibility certificates of variable bounds -Many linear solvers (e.g., Gurobi) do not provide explicit access to the dual +Many linear solvers (e.g., Gurobi) do not provide explicit access to the infeasibility certificate of a variable bound. However, given a set of linear constraints: ```math @@ -148,7 +148,7 @@ l_A \le A x \le u_A \\ l_x \le x \le u_x, \end{align} ``` -the dual certificate of the variable bounds can be computed using the +the certificate of the variable bounds can be computed using the 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 From ec71d9732716d353db8fc84f6aeac8b201ea5fcc Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 19 Nov 2021 10:54:39 +1300 Subject: [PATCH 09/11] Update src/Utilities/results.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Benoît Legat --- src/Utilities/results.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl index b124d205c3..3a81e133a5 100644 --- a/src/Utilities/results.jl +++ b/src/Utilities/results.jl @@ -182,7 +182,7 @@ function get_fallback( return MOI.get(model, MOI.VariablePrimal(attr.result_index), vi) end if is_ray(MOI.get(model, MOI.PrimalStatus())) - c -= MOI.constant(f, Int) + c -= MOI.constant(f, typeof(c)) end return c end From 7b8c0aa0105ca2c52338a78fb5436570ca8ff597 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 19 Nov 2021 10:55:09 +1300 Subject: [PATCH 10/11] Update src/Utilities/results.jl --- src/Utilities/results.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl index 3a81e133a5..8098e8a262 100644 --- a/src/Utilities/results.jl +++ b/src/Utilities/results.jl @@ -26,7 +26,7 @@ function get_fallback(model::MOI.ModelLike, attr::MOI.ObjectiveValue) if is_ray(MOI.get(model, MOI.PrimalStatus())) # Dual infeasibiltiy certificates do not include the primal objective # constant. - obj -= MOI.constant(f, Int) + obj -= MOI.constant(f, typeof(obj)) end return obj end From 25c132542c9093bc354863af07bebb00f64984df Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 19 Nov 2021 12:29:13 +1300 Subject: [PATCH 11/11] Update infeasibility_certificates.md --- .../src/background/infeasibility_certificates.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/src/background/infeasibility_certificates.md b/docs/src/background/infeasibility_certificates.md index 4ba2edbac1..1219666f05 100644 --- a/docs/src/background/infeasibility_certificates.md +++ b/docs/src/background/infeasibility_certificates.md @@ -10,14 +10,14 @@ DocTestFilters = [r"MathOptInterface|MOI"] # Infeasibility certificates When given a conic problem that is infeasible or unbounded, some solvers can -produce a certificate or infeasibility. This page explains what a certificate of +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, which we use to define infeasibility -certificates. A full explanation is given in the section [Duality](@ref), but -here is a brief overview. +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 @@ -69,7 +69,7 @@ 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 +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. @@ -111,7 +111,7 @@ 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 \\ +-\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} ``` @@ -139,7 +139,7 @@ If the solver has found a certificate of primal infeasibility: ### Infeasibility certificates of variable bounds -Many linear solvers (e.g., Gurobi) do not provide explicit access to the +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 @@ -148,7 +148,7 @@ l_A \le A x \le u_A \\ l_x \le x \le u_x, \end{align} ``` -the certificate of the variable bounds can be computed using the +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