diff --git a/Project.toml b/Project.toml index 29c01c562..a22221352 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SDDP" uuid = "f4570300-c277-11e8-125c-4912f86ce65d" -authors = ["Oscar Dowson "] version = "1.13.0" +authors = ["Oscar Dowson "] [deps] Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -9,6 +9,7 @@ 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" diff --git a/src/SDDP.jl b/src/SDDP.jl index 1fefa6840..da9db5350 100644 --- a/src/SDDP.jl +++ b/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 diff --git a/src/algorithm.jl b/src/algorithm.jl index 68543a3ac..ec6664455 100644 --- a/src/algorithm.jl +++ b/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; diff --git a/test/algorithm.jl b/test/algorithm.jl index f16967629..e52f7d803 100644 --- a/test/algorithm.jl +++ b/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()