|
| 1 | +# Copyright (c) 2025: Joaquim Garcia, Oscar Dowson and contributors |
| 2 | +# |
| 3 | +# Use of this source code is governed by an MIT-style license that can be found |
| 4 | +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. |
| 5 | + |
| 6 | +struct InfeasibleBounds{T} |
| 7 | + variable::VariableRef |
| 8 | + lb::T |
| 9 | + ub::T |
| 10 | +end |
| 11 | + |
| 12 | +struct InfeasibleIntegrality{T} |
| 13 | + variable::VariableRef |
| 14 | + lb::T |
| 15 | + ub::T |
| 16 | + set::Union{MOI.Integer,MOI.ZeroOne}#, MOI.Semicontinuous{T}, MOI.Semiinteger{T}} |
| 17 | +end |
| 18 | + |
| 19 | +struct InfeasibleConstraintRange{T} |
| 20 | + constraint::ConstraintRef |
| 21 | + lb::T |
| 22 | + ub::T |
| 23 | + set::Union{MOI.EqualTo{T},MOI.LessThan{T},MOI.GreaterThan{T}} |
| 24 | +end |
| 25 | + |
| 26 | +struct IrreducibleInfeasibleSubset |
| 27 | + constraint::Vector{ConstraintRef} |
| 28 | +end |
| 29 | + |
| 30 | +Base.@kwdef mutable struct InfeasibilityData |
| 31 | + infeasible_bounds::Vector{InfeasibleBounds} = InfeasibleBounds[] |
| 32 | + infeasible_integrality::Vector{InfeasibleIntegrality} = |
| 33 | + InfeasibleIntegrality[] |
| 34 | + |
| 35 | + constraint_range::Vector{InfeasibleConstraintRange} = |
| 36 | + InfeasibleConstraintRange[] |
| 37 | + |
| 38 | + iis::Vector{IrreducibleInfeasibleSubset} = IrreducibleInfeasibleSubset[] |
| 39 | +end |
| 40 | + |
| 41 | +function infeasibility_analysis(model::Model; optimizer = nothing) |
| 42 | + T = Float64 |
| 43 | + |
| 44 | + out = InfeasibilityData() |
| 45 | + |
| 46 | + variables = Dict{VariableRef,Interval{T}}() |
| 47 | + |
| 48 | + # first layer of infeasibility analysis is bounds consistency |
| 49 | + bounds_consistent = true |
| 50 | + for var in JuMP.all_variables(model) |
| 51 | + lb = if JuMP.has_lower_bound(var) |
| 52 | + JuMP.lower_bound(var) |
| 53 | + else |
| 54 | + -Inf |
| 55 | + end |
| 56 | + ub = if JuMP.has_upper_bound(var) |
| 57 | + JuMP.upper_bound(var) |
| 58 | + else |
| 59 | + Inf |
| 60 | + end |
| 61 | + if lb > ub |
| 62 | + push!(out.infeasible_bounds, InfeasibleBounds(var, lb, ub)) |
| 63 | + bounds_consistent = false |
| 64 | + else |
| 65 | + variables[var] = Interval(lb, ub) |
| 66 | + end |
| 67 | + if JuMP.is_integer(var) |
| 68 | + if abs(ub - lb) < 1 && ceil(ub) == ceil(lb) |
| 69 | + push!( |
| 70 | + out.infeasible_integrality, |
| 71 | + InfeasibleIntegrality(var, lb, ub, MOI.Integer()), |
| 72 | + ) |
| 73 | + bounds_consistent = false |
| 74 | + end |
| 75 | + end |
| 76 | + if JuMP.is_binary(var) |
| 77 | + if lb > 0 && ub < 1 |
| 78 | + push!( |
| 79 | + out.infeasible_integrality, |
| 80 | + InfeasibleIntegrality(var, lb, ub, MOI.ZeroOne()), |
| 81 | + ) |
| 82 | + bounds_consistent = false |
| 83 | + end |
| 84 | + end |
| 85 | + end |
| 86 | + # check PSD diagonal >= 0 ? |
| 87 | + # other cones? |
| 88 | + if !bounds_consistent |
| 89 | + return out |
| 90 | + end |
| 91 | + |
| 92 | + # second layer of infeasibility analysis is constraint range analysis |
| 93 | + range_consistent = true |
| 94 | + for (F, S) in JuMP.list_of_constraint_types(model) |
| 95 | + F != JuMP.GenericAffExpr{T,JuMP.VariableRef} && continue |
| 96 | + # TODO: handle quadratics |
| 97 | + !(S in (MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T})) && |
| 98 | + continue |
| 99 | + for con in JuMP.all_constraints(model, F, S) |
| 100 | + con_obj = JuMP.constraint_object(con) |
| 101 | + interval = JuMP.value(x -> variables[x], con_obj.func) |
| 102 | + if con_obj.set isa MOI.EqualTo{T} |
| 103 | + rhs = con_obj.set.value |
| 104 | + if interval.lo > rhs || interval.hi < rhs |
| 105 | + push!( |
| 106 | + out.constraint_range, |
| 107 | + InfeasibleConstraintRange( |
| 108 | + con, |
| 109 | + interval.lo, |
| 110 | + interval.hi, |
| 111 | + con_obj.set, |
| 112 | + ), |
| 113 | + ) |
| 114 | + range_consistent = false |
| 115 | + end |
| 116 | + elseif con_obj.set isa MOI.LessThan{T} |
| 117 | + rhs = con_obj.set.upper |
| 118 | + if interval.lo > rhs |
| 119 | + push!( |
| 120 | + out.constraint_range, |
| 121 | + InfeasibleConstraintRange( |
| 122 | + con, |
| 123 | + interval.lo, |
| 124 | + interval.hi, |
| 125 | + con_obj.set, |
| 126 | + ), |
| 127 | + ) |
| 128 | + range_consistent = false |
| 129 | + end |
| 130 | + elseif con_obj.set isa MOI.GreaterThan{T} |
| 131 | + rhs = con_obj.set.lower |
| 132 | + if interval.hi < rhs |
| 133 | + push!( |
| 134 | + out.constraint_range, |
| 135 | + InfeasibleConstraintRange( |
| 136 | + con, |
| 137 | + interval.lo, |
| 138 | + interval.hi, |
| 139 | + con_obj.set, |
| 140 | + ), |
| 141 | + ) |
| 142 | + range_consistent = false |
| 143 | + end |
| 144 | + end |
| 145 | + end |
| 146 | + end |
| 147 | + |
| 148 | + if !range_consistent |
| 149 | + return out |
| 150 | + end |
| 151 | + |
| 152 | + # check if there is a optimizer |
| 153 | + # third layer is an IIS resolver |
| 154 | + if optimizer === nothing |
| 155 | + println("iis resolver cannot continue because no optimizer is provided") |
| 156 | + return out |
| 157 | + end |
| 158 | + iis = iis_elastic_filter(model, optimizer) |
| 159 | + # for now, only one iis is computed |
| 160 | + push!(out.iis, IrreducibleInfeasibleSubset(iis)) |
| 161 | + |
| 162 | + return out |
| 163 | +end |
| 164 | + |
| 165 | +function iis_elastic_filter(original_model::Model, optimizer) |
| 166 | + |
| 167 | + # if JuMP.termination_status(original_model) == MOI.OPTIMIZE_NOT_CALLED |
| 168 | + # println("iis resolver cannot continue because model is not optimized") |
| 169 | + # # JuMP.optimize!(original_model) |
| 170 | + # end |
| 171 | + |
| 172 | + status = JuMP.termination_status(original_model) |
| 173 | + if status != MOI.INFEASIBLE |
| 174 | + println( |
| 175 | + "iis resolver cannot continue because model is found to be $(status) by the solver", |
| 176 | + ) |
| 177 | + return |
| 178 | + end |
| 179 | + |
| 180 | + model, reference_map = JuMP.copy_model(original_model) |
| 181 | + JuMP.set_optimizer(model, optimizer) |
| 182 | + JuMP.set_silent(model) |
| 183 | + # TODO handle ".ext" to avoid warning |
| 184 | + |
| 185 | + constraint_to_affine = JuMP.relax_with_penalty!(model, default = 1.0) |
| 186 | + # might need to do somehting related to integers / binary |
| 187 | + |
| 188 | + JuMP.optimize!(model) |
| 189 | + |
| 190 | + max_iterations = length(constraint_to_affine) |
| 191 | + |
| 192 | + tolerance = 1e-5 |
| 193 | + |
| 194 | + for _ in 1:max_iterations |
| 195 | + if JuMP.termination_status(model) == MOI.INFEASIBLE |
| 196 | + break |
| 197 | + end |
| 198 | + for (con, func) in constraint_to_affine |
| 199 | + if length(func.terms) == 1 |
| 200 | + var = collect(keys(func.terms))[1] |
| 201 | + if value(var) > tolerance |
| 202 | + fix(var, 0.0; force = true) |
| 203 | + # or delete(model, var) |
| 204 | + delete!(constraint_to_affine, con) |
| 205 | + end |
| 206 | + elseif length(func.terms) == 2 |
| 207 | + var = collect(keys(func.terms)) |
| 208 | + coef1 = func.terms[var[1]] |
| 209 | + coef2 = func.terms[var[2]] |
| 210 | + if value(var1) > tolerance && value(var2) > tolerance |
| 211 | + error("IIS failed due numerical instability") |
| 212 | + elseif value(var[1]) > tolerance |
| 213 | + fix(var[1], 0.0; force = true) |
| 214 | + # or delete(model, var1) |
| 215 | + delete!(constraint_to_affine, con) |
| 216 | + constraint_to_affine[con] = coef2 * var[2] |
| 217 | + elseif value(var[2]) > tolerance |
| 218 | + fix(var[2], 0.0; force = true) |
| 219 | + # or delete(model, var2) |
| 220 | + delete!(constraint_to_affine, con) |
| 221 | + constraint_to_affine[con] = coef1 * var[1] |
| 222 | + end |
| 223 | + else |
| 224 | + println( |
| 225 | + "$con and relaxing function with more than two terms: $func", |
| 226 | + ) |
| 227 | + end |
| 228 | + JuMP.optimize!(model) |
| 229 | + end |
| 230 | + end |
| 231 | + |
| 232 | + pre_iis = Set(keys(constraint_to_affine)) |
| 233 | + iis = JuMP.ConstraintRef[] |
| 234 | + for con in pre_iis |
| 235 | + push!(iis, reference_map[con]) |
| 236 | + end |
| 237 | + |
| 238 | + return iis |
| 239 | +end |
0 commit comments