diff --git a/src/infeasibility.jl b/src/infeasibility.jl index dcb13ad..a93edf6 100644 --- a/src/infeasibility.jl +++ b/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,9 +452,17 @@ 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) @@ -459,10 +470,27 @@ function iis_elastic_filter(original_model::MOI.ModelLike, optimizer) 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 diff --git a/test/infeasibility.jl b/test/infeasibility.jl index 4d5b7c0..d9d7dee 100644 --- a/test/infeasibility.jl +++ b/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)