Add relax_with_penalty

Project.toml

@@ -12,7 +12,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-MathOptInterface = "1.10"
+MathOptInterface = "1.11"
 MutableArithmetics = "1"
 OrderedCollections = "1"
 julia = "1.6"

docs/make.jl

@@ -88,6 +88,7 @@ if !_FAST
     for file in [
         joinpath("getting_started", "getting_started_with_julia.md"),
         joinpath("getting_started", "getting_started_with_JuMP.md"),
+        joinpath("getting_started", "debugging.md"),
         joinpath("linear", "tips_and_tricks.md"),
     ]
         filename = joinpath(@__DIR__, "src", "tutorials", file)

docs/src/changelog.md

@@ -13,6 +13,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
  - Added a `result` keyword argument to [`solution_summary`] to allow
    summarizing models with multiple solutions (#3138)
+ - Add [`relax_with_penalty!`](@ref), which is a useful tool when debugging
+   infeasible models (#3140)
 
 ### Other
 

docs/src/reference/constraints.md

@@ -33,6 +33,7 @@ normalized_rhs
 set_normalized_rhs
 
 add_to_function_constant
+relax_with_penalty!
 ```
 
 ## Deletion

docs/src/tutorials/getting_started/debugging.jl

@@ -164,11 +164,11 @@ import HiGHS
 
 # A simple example of an infeasible model is:
 
-model = Model(HiGHS.Optimizer)
+model = Model(HiGHS.Optimizer);
 set_silent(model)
 @variable(model, x >= 0)
 @objective(model, Max, 2x + 1)
-@constraint(model, 2x - 1 <= -2)
+@constraint(model, con, 2x - 1 <= -2)
 
 # because the bound says that `x >= 0`, but we can rewrite the constraint to be
 # `x <= -1/2`. When the problem is infeasible, JuMP may return one of a number
@@ -230,6 +230,51 @@ termination_status(model)
 #     solvers such as Gurobi and CPLEX do. If the solver does support computing
 #     conflicts, read [Conflicts](@ref) for more details.
 
+# ### Penalty relaxation
+
+# Another strategy to debug sources of infeasibility is the
+# [`relax_with_penalty!`](@ref) function.
+#
+# The penalty relaxation modifies constraints of the form ``f(x) \in S`` into
+# ``f(x) + y - z \in S``, where ``y, z \ge 0``, and then it introduces a
+# penalty term into the objective of ``a \times (y + z)`` (if minimizing, else
+# ``-a``), where ``a`` is a penalty.
+
+map = relax_with_penalty!(model)
+
+# Here `map` is a dictionary which maps constraint indices to an affine
+# expression representing ``(y + z)``.
+
+# If we optimize the relaxed model, this time we get a feasible solution:
+
+optimize!(model)
+termination_status(model)
+
+# Iterate over the contents of `map` to see which constraints are violated:
+
+for (con, penalty) in map
+    violation = value(penalty)
+    if violation > 0
+        println("Constraint `$(name(con))` is violated by $violation")
+    end
+end
+
+# Once you find a violated constraint in the relaxed problem, take a look to see
+# if there is a typo or other common mistake in that particular constraint.
+
+# Consult the docstring [`relax_with_penalty!`](@ref) for information on how to
+# modify the penalty cost term `a`, either for every constraint in the model or
+# a particular subset of the constraints.
+
+# When using [`relax_with_penalty!`](@ref), you should be aware that:
+#
+#  * Variable bounds and integrality restrictions are not relaxed. If the
+#    problem is still infeasible after calling [`relax_with_penalty!`](@ref),
+#    check the variable bounds.
+#  * You cannot undo the penalty relaxation. If you need an unmodified model,
+#    rebuild the problem, or call [`copy_model`](@ref) before calling
+#    [`relax_with_penalty!`](@ref).
+
 # ## Debugging an unbounded model
 
 # A model is unbounded if there is no limit on how good the objective value can
@@ -241,7 +286,7 @@ termination_status(model)
 
 # A simple example of an unbounded model is:
 
-model = Model(HiGHS.Optimizer)
+model = Model(HiGHS.Optimizer);
 set_silent(model)
 @variable(model, x >= 0)
 @objective(model, Max, 2x + 1)
@@ -288,7 +333,7 @@ termination_status(model)
 # the variable must be less-than or equal to the expression of the objective
 # function. For example:
 
-model = Model(HiGHS.Optimizer)
+model = Model(HiGHS.Optimizer);
 set_silent(model)
 @variable(model, x >= 0)
 ## @objective(model, Max, 2x + 1)

src/constraints.jl

@@ -1344,3 +1344,107 @@ function all_constraints(
     append!(ret, all_nonlinear_constraints(model))
     return ret
 end
+
+"""
+    relax_with_penalty!(
+        model::Model,
+        [penalties::Dict{ConstraintRef,Float64}];
+        [default::Union{Nothing,Real} = nothing,]
+    )
+
+Destructively modify the model in-place to create a penalized relaxation of the
+constraints.
+
+!!! warning
+    This is a destructive routine that modifies the model in-place. If you don't
+    want to modify the original model, use [`copy_model`](@ref) to create a copy
+    before calling [`relax_with_penalty!`](@ref).
+
+## Reformulation
+
+See [`MOI.Utilities.ScalarPenaltyRelaxation`](@ref) for details of the
+reformulation.
+
+For each constraint `ci`, the penalty passed to
+[`MOI.Utilities.ScalarPenaltyRelaxation`](@ref) is `get(penalties, ci, default)`.
+If the value is `nothing`, because `ci` does not exist in `penalties` and
+`default = nothing`, then the constraint is skipped.
+
+## Return value
+
+This function returns a `Dict{ConstraintRef,AffExpr}` that maps each constraint
+index to the corresponding `y + z` as an [`AffExpr`](@ref). In an optimal
+solution, query the value of these functions to compute the violation of each
+constraint.
+
+## Relax a subset of constraints
+
+To relax a subset of constraints, pass a `penalties` dictionary and set
+`default = nothing`.
+
+## Examples
+
+```jldoctest
+julia> function new_model()
+           model = Model()
+           @variable(model, x)
+           @objective(model, Max, 2x + 1)
+           @constraint(model, c1, 2x - 1 <= -2)
+           @constraint(model, c2, 3x >= 0)
+           return model
+       end
+new_model (generic function with 1 method)
+
+julia> model_1 = new_model();
+
+julia> relax_with_penalty!(model_1; default = 2.0)
+Dict{ConstraintRef{Model, C, ScalarShape} where C, AffExpr} with 2 entries:
+  c1 : 2 x - _[3] ≤ -1.0 => _[3]
+  c2 : 3 x + _[2] ≥ 0.0  => _[2]
+
+julia> print(model_1)
+Max 2 x - 2 _[2] - 2 _[3] + 1
+Subject to
+ c2 : 3 x + _[2] ≥ 0.0
+ c1 : 2 x - _[3] ≤ -1.0
+ _[2] ≥ 0.0
+ _[3] ≥ 0.0
+
+julia> model_2 = new_model();
+
+julia> relax_with_penalty!(model_2, Dict(model_2[:c2] => 3.0))
+Dict{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.GreaterThan{Float64}}, ScalarShape}, AffExpr} with 1 entry:
+  c2 : 3 x + _[2] ≥ 0.0 => _[2]
+
+julia> print(model_2)
+Max 2 x - 3 _[2] + 1
+Subject to
+ c2 : 3 x + _[2] ≥ 0.0
+ c1 : 2 x ≤ -1.0
+ _[2] ≥ 0.0
+```
+"""
+function relax_with_penalty!(
+    model::Model,
+    penalties::Dict;
+    default::Union{Nothing,Real} = nothing,
+)
+    if default !== nothing
+        default = Float64(default)
+    end
+    moi_penalties = Dict{MOI.ConstraintIndex,Float64}(
+        index(k) => Float64(v) for (k, v) in penalties
+    )
+    map = MOI.modify(
+        backend(model),
+        MOI.Utilities.PenaltyRelaxation(moi_penalties; default = default),
+    )
+    return Dict(
+        constraint_ref_with_index(model, k) => jump_function(model, v) for
+        (k, v) in map
+    )
+end
+
+function relax_with_penalty!(model::Model; default::Real = 1.0)
+    return relax_with_penalty!(model, Dict(); default = default)
+end

test/constraint.jl

@@ -1220,6 +1220,88 @@ function test_Model_all_constraints(::Any, ::Any)
     return
 end
 
+function test_Model_relax_with_penalty!_default(::Any, ::Any)
+    model = Model()
+    @variable(model, x >= 0)
+    map = relax_with_penalty!(model)
+    @test isempty(map)
+    @constraint(model, c1, x <= 1)
+    @constraint(model, c2, x == 0)
+    map = relax_with_penalty!(model)
+    @test length(map) == 2
+    @test map[c1] isa AffExpr
+    @test map[c2] isa AffExpr
+    @test num_variables(model) == 4
+    @test objective_sense(model) == MOI.MIN_SENSE
+    @test objective_function(model) == map[c1] + map[c2]
+    return
+end
+
+function test_Model_relax_with_penalty!_max(::Any, ::Any)
+    model = Model()
+    @variable(model, x >= 0)
+    @constraint(model, c1, x <= 1)
+    @constraint(model, c2, x == 0)
+    @objective(model, Max, 1.0 * x + 2.5)
+    map = relax_with_penalty!(model)
+    @test length(map) == 2
+    @test map[c1] isa AffExpr
+    @test map[c2] isa AffExpr
+    @test num_variables(model) == 4
+    @test objective_sense(model) == MOI.MAX_SENSE
+    @test objective_function(model) == x + 2.5 - map[c1] - map[c2]
+    return
+end
+
+function test_Model_relax_with_penalty!_constant(::Any, ::Any)
+    model = Model()
+    @variable(model, x >= 0)
+    map = relax_with_penalty!(model)
+    @test isempty(map)
+    @constraint(model, c1, x <= 1)
+    @constraint(model, c2, x == 0)
+    map = relax_with_penalty!(model; default = 2)
+    @test length(map) == 2
+    @test map[c1] isa AffExpr
+    @test map[c2] isa AffExpr
+    @test num_variables(model) == 4
+    @test objective_sense(model) == MOI.MIN_SENSE
+    @test objective_function(model) == 2.0 * map[c1] + 2.0 * map[c2]
+    return
+end
+
+function test_Model_relax_with_penalty!_specific(::Any, ::Any)
+    model = Model()
+    @variable(model, x >= 0)
+    map = relax_with_penalty!(model)
+    @test isempty(map)
+    @constraint(model, c1, x <= 1)
+    @constraint(model, c2, x == 0)
+    map = relax_with_penalty!(model, Dict(c1 => 3.0))
+    @test length(map) == 1
+    @test map[c1] isa AffExpr
+    @test num_variables(model) == 2
+    @test objective_sense(model) == MOI.MIN_SENSE
+    @test objective_function(model) == 3.0 * map[c1]
+    return
+end
+
+function test_Model_relax_with_penalty!_specific_with_default(::Any, ::Any)
+    model = Model()
+    @variable(model, x >= 0)
+    map = relax_with_penalty!(model)
+    @test isempty(map)
+    @constraint(model, c1, x <= 1)
+    @constraint(model, c2, x == 0)
+    map = relax_with_penalty!(model, Dict(c1 => 3.0); default = 1)
+    @test length(map) == 2
+    @test map[c1] isa AffExpr
+    @test num_variables(model) == 4
+    @test objective_sense(model) == MOI.MIN_SENSE
+    @test objective_function(model) == 3 * map[c1] + map[c2]
+    return
+end
+
 function runtests()
     for name in names(@__MODULE__; all = true)
         if !startswith("$(name)", "test_")