Write out IIS when an infeasible subproblem is encountered

Project.toml

@@ -1,14 +1,15 @@
 name = "SDDP"
 uuid = "f4570300-c277-11e8-125c-4912f86ce65d"
-authors = ["Oscar Dowson <[email protected]>"]
 version = "1.13.0"
+authors = ["Oscar Dowson <[email protected]>"]
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MathOptIIS = "8c4f8055-bd93-4160-a86b-a0c04941dbff"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -22,6 +23,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 HTTP = "1"
 JSON = "0.21, 1"
 JuMP = "1"
+MathOptIIS = "0.1.1"
 MutableArithmetics = "1"
 RecipesBase = "1"
 Reexport = "1"

src/SDDP.jl

@@ -11,6 +11,7 @@ Reexport.@reexport using JuMP
 import Distributed
 import HTTP
 import JSON
+import MathOptIIS
 import MutableArithmetics
 import Printf
 import Random

src/algorithm.jl

@@ -381,8 +381,8 @@ function attempt_numerical_recovery(
         # file, a second thread finds an infeasibility of the same node, it
         # doesn't matter if we over-write this file.
         filename = "model_infeasible_node_$(node.index).cuts.json"
-        @info "Writing cuts to the file `$filename`"
         write_cuts_to_file(model, filename)
+        write_iis_to_file(node, "subproblem_$(node.index)")
         write_subproblem_to_file(
             node,
             "subproblem_$(node.index).mof.json";
@@ -392,6 +392,58 @@ function attempt_numerical_recovery(
     return
 end
 
+function Base.getindex(map::JuMP.GenericReferenceMap, x::State)
+    return State(map[x.in], map[x.out])
+end
+
+_supports_iis(::Type{T}) where {T} = false
+_supports_iis(::Type{<:GenericVariableRef}) = true
+_supports_iis(::Type{<:GenericAffExpr}) = true
+_supports_iis(::Type{<:MOI.LessThan}) = true
+_supports_iis(::Type{<:MOI.GreaterThan}) = true
+_supports_iis(::Type{<:MOI.EqualTo}) = true
+_supports_iis(::Type{<:MOI.Interval}) = true
+_supports_iis(::Type{<:MOI.ZeroOne}) = true
+_supports_iis(::Type{<:MOI.Integer}) = true
+
+function write_iis_to_file(node, filename)
+    # For now, write out the IIS only if the subproblem is a MILP.
+    for (F, S) in list_of_constraint_types(node.subproblem)
+        if !_supports_iis(F) || !_supports_iis(S)
+            return
+        end
+    end
+    iis = MathOptIIS.Optimizer()
+    MOI.set(iis, MathOptIIS.InfeasibleModel(), backend(node.subproblem))
+    MOI.set(
+        iis,
+        MathOptIIS.InnerOptimizer(),
+        () -> MOI.instantiate(node.optimizer; with_cache_type = Float64),
+    )
+    MOI.compute_conflict!(iis)
+    @show MOI.get(iis, MOI.ConflictStatus())
+    if MOI.get(iis, MOI.ConflictStatus()) != MOI.CONFLICT_FOUND
+        return
+    end
+    # we need to cache the .ext dict, because JuMP doesn't know how to copy
+    # the elements.
+    ext_kwargs = [k => v for (k, v) in node.subproblem.ext]
+    empty!(node.subproblem.ext)
+    # Copy the model, filtering constraints which are not in the conflict
+    function filter_constraints(cref::ConstraintRef)
+        status = MOI.get(iis, MOI.ConstraintConflictStatus(), index(cref))
+        return status != MOI.NOT_IN_CONFLICT
+    end
+    new_model, _ = copy_model(node.subproblem; filter_constraints)
+    set_objective_sense(new_model, MOI.FEASIBILITY_SENSE)
+    # Restore the .ext dictionary
+    for (k, v) in ext_kwargs
+        node.subproblem.ext[k] = v
+    end
+    JuMP.write_to_file(new_model, "$(filename).iis.lp")
+    return
+end
+
 function default_numerical_difficulty_callback(
     model::PolicyGraph,
     node::Node;

test/algorithm.jl

@@ -8,6 +8,7 @@ module TestAlgorithm
 using SDDP
 using Test
 import HiGHS
+import Ipopt
 
 function runtests()
     for name in names(@__MODULE__; all = true)
@@ -472,6 +473,78 @@ function test_root_node_risk_measure()
     return
 end
 
+function test_iis_lp()
+    model = SDDP.LinearPolicyGraph(;
+        stages = 3,
+        lower_bound = 0.0,
+        optimizer = HiGHS.Optimizer,
+    ) do sp, stage
+        @variable(sp, 0 <= x <= 100, SDDP.State, initial_value = 0)
+        @variable(sp, 0 <= u_p <= 200)
+        @variable(sp, u_o >= 0)
+        @constraint(sp, x.out >= u_p)
+        @constraint(sp, u_p == 101)
+        @constraint(sp, x.out == x.in + u_p + 1.1 * u_o)
+        @stageobjective(sp, 100 * u_p + 300 * u_o + 50 * x.out)
+    end
+    root = pwd()
+    try
+        mktempdir() do path
+            cd(path)
+            try
+                SDDP.train(model)
+            catch
+            end
+            contents = read("subproblem_1.iis.lp", String)
+            @test occursin("u_p = 101", contents)
+            @test !occursin("300", contents)
+            @test !occursin("1.1", contents)
+        end
+    finally
+        cd(root)
+    end
+    return
+end
+
+function test_iis_nlp()
+    model = SDDP.LinearPolicyGraph(;
+        stages = 3,
+        lower_bound = 0.0,
+        optimizer = Ipopt.Optimizer,
+    ) do sp, stage
+        @variable(sp, 0 <= x <= 1, SDDP.State, initial_value = 0)
+        @variable(sp, 0 <= u_p <= 200)
+        @variable(sp, u_o >= 0)
+        @constraint(sp, x.out >= u_p)
+        @constraint(sp, exp(u_p) == 101)
+        @constraint(sp, x.out == x.in + u_p + 1.1 * u_o)
+        @stageobjective(sp, 100 * u_p + 300 * u_o + 50 * x.out)
+    end
+    try
+        SDDP.train(model)
+    catch
+    end
+    @test !isfile("subproblem_1.iis.lp")
+    return
+end
+
+function test_iis_unbounded()
+    model = SDDP.LinearPolicyGraph(;
+        stages = 3,
+        lower_bound = 0.0,
+        optimizer = HiGHS.Optimizer,
+    ) do sp, stage
+        @variable(sp, 0 <= x, SDDP.State, initial_value = 0)
+        @stageobjective(sp, -x.out)
+    end
+    try
+        SDDP.train(model)
+    catch
+    end
+    @test !isfile("subproblem_3.iis.lp")
+    return
+end
+
 end  # module
 
 TestAlgorithm.runtests()