diff --git a/docs/src/apireference.md b/docs/src/apireference.md index 248e80e94..d1d032543 100644 --- a/docs/src/apireference.md +++ b/docs/src/apireference.md @@ -34,6 +34,7 @@ SDDP.termination_status SDDP.write_cuts_to_file SDDP.read_cuts_from_file SDDP.write_log_to_csv +SDDP.set_numerical_difficulty_callback ``` ### [Stopping rules](@id api_stopping_rules) diff --git a/src/algorithm.jl b/src/algorithm.jl index d6fdc3825..1b84bc355 100644 --- a/src/algorithm.jl +++ b/src/algorithm.jl @@ -292,19 +292,77 @@ function _has_primal_solution(node::Node) return status in (JuMP.FEASIBLE_POINT, JuMP.NEARLY_FEASIBLE_POINT) end -function attempt_numerical_recovery(model::PolicyGraph, node::Node) - if JuMP.mode(node.subproblem) == JuMP.DIRECT - @warn( - "Unable to recover in direct mode! Remove `direct = true` when " * - "creating the policy graph." - ) - else - model.ext[:numerical_issue] = true - MOI.Utilities.reset_optimizer(node.subproblem) - optimize!(node.subproblem) +function _has_dual_solution(node::Node) + status = JuMP.dual_status(node.subproblem) + return status in (JuMP.FEASIBLE_POINT, JuMP.NEARLY_FEASIBLE_POINT) +end + +""" + set_numerical_difficulty_callback( + model::PolicyGraph, + callback::Function, + ) + +Set a callback function `callback(::PolicyGraph, ::Node; require_dual::Bool)` +that is run when the optimizer terminates without finding a primal solution (and +dual solution if `require_dual` is `true`). + +## Default callback + +The default callback is a small variation of: +```julia +function callback(::PolicyGraph, node::Node; require_dual::Bool) + MOI.Utilities.reset_optimizer(node.subproblem) + optimize!(node.subproblem) + return +end +``` +This callback is the default because a common issue is solvers declaring the +infeasible because of numerical issues related to the large number of cutting +planes. Resetting the subproblem---and therefore starting from a fresh problem +instead of warm-starting from the previous solution---is often enough to fix the +problem and allow more iterations. + +## Other callbacks + +In cases where the problem is truely infeasible (not because of numerical issues +), it may be helpful to write out the irreducible infeasible subsystem (IIS) for +debugging. For this use-case, use a callback as follows: +```julia +function callback(::PolicyGraph, node::Node; require_dual::Bool) + JuMP.compute_conflict!(node.suprobblem) + status = JuMP.get_attribute(node.subproblem, MOI.ConflictStatus()) + if status == MOI.CONFLICT_FOUND + iis_model, _ = JuMP.copy_conflict(node.subproblem) + print(iis_model) end - if !_has_primal_solution(node) - model.ext[:numerical_issue] = true + return +end +SDDP.set_numerical_difficulty_callback(model, callback) +``` +""" +function set_numerical_difficulty_callback( + model::PolicyGraph, + callback::Function, +) + model.ext[:numerical_difficulty_callback] = callback + return +end + +function attempt_numerical_recovery( + model::PolicyGraph, + node::Node; + require_dual::Bool = false, +) + model.ext[:numerical_issue] = true + callback = get( + model.ext, + :numerical_difficulty_callback, + default_numerical_difficulty_callback, + ) + callback(model, node; require_dual) + missing_dual_solution = require_dual && !_has_dual_solution(node) + if !_has_primal_solution(node) || missing_dual_solution # We use the `node.index` in the filename because two threads could both # try to write the cuts to file at the same time. If, after writing this # file, a second thread finds an infeasibility of the same node, it @@ -321,6 +379,23 @@ function attempt_numerical_recovery(model::PolicyGraph, node::Node) return end +function default_numerical_difficulty_callback( + model::PolicyGraph, + node::Node; + kwargs..., +) + if JuMP.mode(node.subproblem) == JuMP.DIRECT + @warn( + "Unable to recover in direct mode! Remove `direct = true` when " * + "creating the policy graph." + ) + return + end + MOI.Utilities.reset_optimizer(node.subproblem) + optimize!(node.subproblem) + return +end + """ _initialize_solver(node::Node; throw_error::Bool) diff --git a/src/plugins/duality_handlers.jl b/src/plugins/duality_handlers.jl index b1fbc86a1..dc3391d3c 100644 --- a/src/plugins/duality_handlers.jl +++ b/src/plugins/duality_handlers.jl @@ -91,25 +91,10 @@ min Cᵢ(x̄, u, w) + θᵢ """ struct ContinuousConicDuality <: AbstractDualityHandler end -function _has_dual_solution(node::Node) - status = JuMP.dual_status(node.subproblem) - return status in (JuMP.FEASIBLE_POINT, JuMP.NEARLY_FEASIBLE_POINT) -end - function get_dual_solution(node::Node, ::ContinuousConicDuality) if !_has_dual_solution(node) - # Attempt to recover by resetting the optimizer and re-solving. - if JuMP.mode(node.subproblem) != JuMP.DIRECT - MOI.Utilities.reset_optimizer(node.subproblem) - optimize!(node.subproblem) - end - end - if !_has_dual_solution(node) - write_subproblem_to_file( - node, - "subproblem.mof.json"; - throw_error = true, - ) + model = node.subproblem.ext[:sddp_policy_graph] + attempt_numerical_recovery(model, node; require_dual = true) end # Note: due to JuMP's dual convention, we need to flip the sign for # maximization problems. diff --git a/test/algorithm.jl b/test/algorithm.jl index 756f45990..8e555affb 100644 --- a/test/algorithm.jl +++ b/test/algorithm.jl @@ -359,6 +359,34 @@ function test_log_frequency_argument_error() return end +function test_numerical_difficulty_callback() + model = SDDP.LinearPolicyGraph(; + stages = 2, + lower_bound = 0.0, + optimizer = HiGHS.Optimizer, + ) do node, stage + @variable(node, x >= 0, SDDP.State, initial_value = 0.0) + if stage == 2 + @constraint(node, c_infeasible, x.in >= 2) + end + @stageobjective(node, x.out) + end + callback_called_from = Int[] + function my_callback(model, node; require_dual) + push!(callback_called_from, node.index) + JuMP.delete(node.subproblem, node.subproblem[:c_infeasible]) + JuMP.optimize!(node.subproblem) + return + end + SDDP.set_numerical_difficulty_callback(model, my_callback) + SDDP.train(model; iteration_limit = 3) + @test callback_called_from == [2] + @test model.most_recent_training_results.status == :iteration_limit + log = model.most_recent_training_results.log + @test map(l -> l.serious_numerical_issue, log) == [1, 0, 0] + return +end + end # module TestAlgorithm.runtests()