From 306efca317a88283106632e717250af3d5474aa5 Mon Sep 17 00:00:00 2001 From: odow Date: Mon, 23 Aug 2021 21:08:51 +1200 Subject: [PATCH] Fix bug in lagrangian duality --- src/algorithm.jl | 4 +++- src/plugins/duality_handlers.jl | 23 ++++++++++++----------- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/src/algorithm.jl b/src/algorithm.jl index cc5bf4531..08c3ae52b 100644 --- a/src/algorithm.jl +++ b/src/algorithm.jl @@ -382,7 +382,9 @@ function solve_subproblem( end state = get_outgoing_state(node) stage_objective = stage_objective_value(node.stage_objective) - objective, dual_values = get_dual_solution(node, duality_handler) + TimerOutputs.@timeit SDDP_TIMER "get_dual_solution" begin + objective, dual_values = get_dual_solution(node, duality_handler) + end if node.post_optimize_hook !== nothing node.post_optimize_hook(pre_optimize_ret) end diff --git a/src/plugins/duality_handlers.jl b/src/plugins/duality_handlers.jl index 2eca288e4..29c793403 100644 --- a/src/plugins/duality_handlers.jl +++ b/src/plugins/duality_handlers.jl @@ -264,15 +264,12 @@ function get_dual_solution(node::Node, lagrange::LagrangianDuality) h_expr[i] = @expression(node.subproblem, state.in - x_in_value[i]) # Relax the constraint from the problem. JuMP.unfix(state.in) - # We set new bounds here for a few reasons. First, it ensures that the - # relaxed primal problem always has a bounded optimal solution. Without - # the bounds, the primal problem could be unbounded (i.e., an infeasible - # dual solution). With bounds on x_in, it ensures that the dual problem - # is always feasible. The choice of ±1 is somewhat arbitrary. It needs - # to be large enough to avoid numerical error, but not too large to - # cause other issues. - JuMP.set_lower_bound(state.in, x_in_value[i] - 1) - JuMP.set_upper_bound(state.in, x_in_value[i] + 1) + # We need bounds to ensure that the dual problem is feasible. However, + # they can't be too tight. Let's use 1e9 as a default... + lb = has_lower_bound(state.out) ? lower_bound(state.out) : -1e9 + ub = has_upper_bound(state.out) ? upper_bound(state.out) : 1e9 + JuMP.set_lower_bound(state.in, lb) + JuMP.set_upper_bound(state.in, ub) end # Create the model for the cutting plane algorithm model = JuMP.Model(something(lagrange.optimizer, node.optimizer)) @@ -391,8 +388,12 @@ function get_dual_solution(node::Node, ::StrengthenedConicDuality) x[i] = JuMP.fix_value(state.in) h_expr[i] = @expression(node.subproblem, state.in - x[i]) JuMP.unfix(state.in) - JuMP.set_lower_bound(state.in, x[i] - 1) - JuMP.set_upper_bound(state.in, x[i] + 1) + # We need bounds to ensure that the dual problem is feasible. However, + # they can't be too tight. Let's use 1e9 as a default... + lb = has_lower_bound(state.out) ? lower_bound(state.out) : -1e9 + ub = has_upper_bound(state.out) ? upper_bound(state.out) : 1e9 + JuMP.set_lower_bound(state.in, lb) + JuMP.set_upper_bound(state.in, ub) λ_k[i] = conic_dual[key] end lagrangian_obj = _solve_primal_problem(node.subproblem, λ_k, h_expr, h_k)