Handle unboundedness in Infeasibility analysis (#37)

* Fix print of Infeasibility analysis

* format

* Handle unboundedness in Infeasibility analysis

* format

* tweaks and tests

* format

Author: joaquimg

src/infeasibility.jl

@@ -438,7 +438,10 @@ function iis_elastic_filter(original_model::MOI.ModelLike, optimizer)
 
     # handle optimize not called
     status = MOI.get(original_model, MOI.TerminationStatus())
-    if status != MOI.INFEASIBLE
+    if !(
+        status in
+        (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.ALMOST_INFEASIBLE)
+    )
         println(
             "iis resolver cannot continue because model is found to be $(status) by the solver",
         )
@@ -449,20 +452,45 @@ function iis_elastic_filter(original_model::MOI.ModelLike, optimizer)
     reference_map = MOI.copy_to(model, original_model)
     MOI.set(model, MOI.Silent(), true)
 
+    obj_sense = MOI.get(model, MOI.ObjectiveSense())
+    base_obj_type = MOI.get(model, MOI.ObjectiveFunctionType())
+    base_obj_func = MOI.get(model, MOI.ObjectiveFunction{base_obj_type}())
+
     constraint_to_affine =
         MOI.modify(model, MOI.Utilities.PenaltyRelaxation(default = 1.0))
     # might need to do something related to integers / binary
+    relaxed_obj_type = MOI.get(model, MOI.ObjectiveFunctionType())
+    relaxed_obj_func = MOI.get(model, MOI.ObjectiveFunction{relaxed_obj_type}())
+
+    pure_relaxed_obj_func = relaxed_obj_func - base_obj_func
 
     max_iterations = length(constraint_to_affine)
 
     tolerance = 1e-5
 
     de_elastisized = []
 
+    changed_obj = false
+
     for i in 1:max_iterations
         MOI.optimize!(model)
         status = MOI.get(model, MOI.TerminationStatus())
-        if status == MOI.INFEASIBLE
+        if status in ( # possibily primal unbounded
+            MOI.INFEASIBLE_OR_UNBOUNDED,
+            MOI.DUAL_INFEASIBLE,
+            MOI.ALMOST_DUAL_INFEASIBLE,
+        )
+            #try with a pure relaxation objective
+            MOI.set(
+                model,
+                MOI.ObjectiveFunction{relaxed_obj_type}(),
+                pure_relaxed_obj_func,
+            )
+            changed_obj = true
+            MOI.optimize!(model)
+        end
+        if status in
+           (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.ALMOST_INFEASIBLE)
             break
         end
         for (con, func) in constraint_to_affine
@@ -503,20 +531,58 @@ function iis_elastic_filter(original_model::MOI.ModelLike, optimizer)
         end
     end
 
+    if changed_obj
+        MOI.set(
+            model,
+            MOI.ObjectiveFunction{relaxed_obj_type}(),
+            relaxed_obj_func,
+        )
+    end
+
     # consider deleting all no iis constraints
     # be careful with intervals
 
+    obj_type = MOI.get(model, MOI.ObjectiveFunctionType())
+    obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_type}())
+    obj_sense = MOI.get(model, MOI.ObjectiveSense())
+
     # deletion filter
     cadidates = MOI.ConstraintIndex[]
     for (con, var, has_lower) in de_elastisized
         _set_bound_zero(model, var, has_lower, T)
         MOI.optimize!(model)
         status = MOI.get(model, MOI.TerminationStatus())
-        if status in (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE)
+        if status in
+           (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.ALMOST_INFEASIBLE)
             # this constraint is not in IIS
-        elseif status in (MOI.OPTIMAL, MOI.ALMOST_OPTIMAL)
+        elseif status in (
+            MOI.OPTIMAL,
+            MOI.ALMOST_OPTIMAL,
+            MOI.LOCALLY_SOLVED,
+            MOI.ALMOST_LOCALLY_SOLVED,
+        )
             push!(cadidates, con)
             _fix_to_zero(model, var, T)
+        elseif status in (
+            MOI.INFEASIBLE_OR_UNBOUNDED,
+            MOI.DUAL_INFEASIBLE,
+            MOI.ALMOST_DUAL_INFEASIBLE, # possibily primal unbounded
+        )
+            MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
+            MOI.optimize!(model)
+            primal_status = MOI.get(model, MOI.PrimalStatus())
+            if primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT)
+                # this constraint is not in IIS
+                push!(cadidates, con)
+                _fix_to_zero(model, var, T)
+                MOI.set(model, MOI.ObjectiveSense(), obj_sense)
+                MOI.set(model, MOI.ObjectiveFunction{obj_type}(), obj_func)
+            else
+                error(
+                    "IIS failed due numerical instability, got status $status,",
+                    "then, for MOI.FEASIBILITY_SENSE objective, got primal status $primal_status",
+                )
+            end
         else
             error("IIS failed due numerical instability, got status $status")
         end

test/infeasibility.jl

@@ -419,6 +419,33 @@ function test_iis()
     return
 end
 
+function test_iis_free_var()
+    model = Model(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, y)
+    @constraint(model, c1, x + y <= 1)
+    @constraint(model, c2, x + y >= 2)
+    @objective(model, Max, -2x + y)
+    optimize!(model)
+    data = ModelAnalyzer.analyze(
+        ModelAnalyzer.Infeasibility.Analyzer(),
+        model,
+        optimizer = HiGHS.Optimizer,
+    )
+    list = ModelAnalyzer.list_of_issue_types(data)
+    @test length(list) == 1
+    ret = ModelAnalyzer.list_of_issues(data, list[1])
+    @test length(ret) == 1
+    @test length(ret[].constraint) == 2
+    @test Set([ret[].constraint[1], ret[].constraint[2]]) ==
+          Set(JuMP.index.([c2, c1]))
+    iis = ModelAnalyzer.constraints(ret[], model)
+    @test length(iis) == 2
+    @test Set(iis) == Set([c2, c1])
+    return
+end
+
 function test_iis_multiple()
     model = Model(HiGHS.Optimizer)
     set_silent(model)