diff --git a/src/plugins/bellman_functions.jl b/src/plugins/bellman_functions.jl index 438c8d876..4d5496c4b 100644 --- a/src/plugins/bellman_functions.jl +++ b/src/plugins/bellman_functions.jl @@ -479,6 +479,10 @@ function _refine_bellman_function_no_lock( end end +_copy_value(::Nothing) = nothing +_copy_value(state::ObjectiveState) = state.state +_copy_value(state::BeliefState) = copy(state.belief) + function _add_average_cut( node::Node, outgoing_state::Dict{Symbol,Float64}, @@ -502,10 +506,8 @@ function _add_average_cut( end # Now add the average-cut to the subproblem. We include the objective-state # component μᵀy and the belief state (if it exists). - obj_y = - node.objective_state === nothing ? nothing : node.objective_state.state - belief_y = - node.belief_state === nothing ? nothing : node.belief_state.belief + obj_y = _copy_value(node.objective_state) + belief_y = _copy_value(node.belief_state) _add_cut( node.bellman_function.global_theta, θᵏ, @@ -542,9 +544,8 @@ function _add_multi_cut( objective_realizations[i], dual_variables[i], outgoing_state, - node.objective_state === nothing ? nothing : - node.objective_state.state, - node.belief_state === nothing ? nothing : node.belief_state.belief, + _copy_value(node.objective_state), + _copy_value(node.belief_state), ) end model = JuMP.owner_model(bellman_function.global_theta.theta) diff --git a/test/plugins/bellman_functions.jl b/test/plugins/bellman_functions.jl index 75dcf0596..24448b18c 100644 --- a/test/plugins/bellman_functions.jl +++ b/test/plugins/bellman_functions.jl @@ -407,6 +407,57 @@ function test_cut_selection_flags() return end +function test_issue_892() + graph = SDDP.Graph( + :root_node, + [:Ad, :Ah, :Bd, :Bh], + [ + (:root_node => :Ad, 0.5), + (:root_node => :Bd, 0.5), + (:Ad => :Ah, 1.0), + (:Ah => :Ad, 0.9), + (:Bd => :Bh, 1.0), + (:Bh => :Bd, 0.9), + ], + ) + SDDP.add_ambiguity_set(graph, [:Ad, :Bd], 1e2) + SDDP.add_ambiguity_set(graph, [:Ah, :Bh], 1e2) + model = SDDP.PolicyGraph( + graph; + lower_bound = 0.0, + optimizer = HiGHS.Optimizer, + ) do sp, node + @variable(sp, 0 <= x <= 2, SDDP.State, initial_value = 0.0) + @variable(sp, u >= 0) + @variable(sp, w == 0) + @constraint(sp, w == x.in - x.out + u) + if node == :Ad || node == :Bd + @stageobjective(sp, u) + else + P = Dict(:Ah => [0.2, 0.8], :Bh => [0.8, 0.2]) + SDDP.parameterize(ω -> fix(w, ω), sp, 1:2, P[node]) + @stageobjective(sp, 2 * u + x.out) + end + end + SDDP.train(model; iteration_limit = 200, print_level = 0) + simulations = SDDP.simulate( + model, + 100; + sampling_scheme = SDDP.InSampleMonteCarlo(; + max_depth = 50, + terminate_on_dummy_leaf = false, + ), + ) + function calculate_objective(sim) + ρ(t) = 0.9^div(t - 1, 2) + return sum(ρ(t) * s[:stage_objective] for (t, s) in enumerate(sim)) + end + objectives = calculate_objective.(simulations) + @test minimum(objectives) < 12 + @test SDDP.calculate_bound(model) < 19 + return +end + end # module TestBellmanFunctions.runtests()