Add variable bounds to incoming states in LagrangianDuality (#819)

Author: odow

docs/src/guides/add_multidimensional_noise.md

@@ -26,11 +26,12 @@ julia> model = SDDP.LinearPolicyGraph(
            SDDP.parameterize(subproblem, Ω, P) do ω
                JuMP.fix(x.out, ω.value)
                @stageobjective(subproblem, ω.coefficient * x.out)
-               println("ω is: ", ω)
            end
        end;
 
-julia> SDDP.simulate(model, 1);
+julia> for s in SDDP.simulate(model, 1)[1]
+           println("ω is: ", s[:noise_term])
+       end
 ω is: (value = 1, coefficient = 4)
 ω is: (value = 1, coefficient = 3)
 ω is: (value = 2, coefficient = 4)
@@ -46,7 +47,7 @@ For sampling multidimensional random variates, it can be useful to use the
 There are several ways to go about this. If the sample space is finite, and
 small enough that it makes sense to enumerate each element, we can use
 `Base.product` and `Distributions.support` to generate the entire sample space
-`Ω` from each of the marginal distributions. 
+`Ω` from each of the marginal distributions.
 
 We can then evaluate the density function of the product distribution on each
 element of this space to get the vector of corresponding probabilities, `P`.
@@ -75,11 +76,12 @@ julia> model = SDDP.LinearPolicyGraph(
            SDDP.parameterize(subproblem, Ω, P) do ω
                JuMP.fix(x.out, ω[1])
                @stageobjective(subproblem, ω[2] * x.out + ω[3])
-               println("ω is: ", ω)
            end
        end;
 
-julia> SDDP.simulate(model, 1);
+julia> for s in SDDP.simulate(model, 1)[1]
+           println("ω is: ", s[:noise_term])
+       end
 ω is: [10, 0, 3]
 ω is: [0, 1, 6]
 ω is: [6, 0, 5]
@@ -117,11 +119,12 @@ julia> model = SDDP.LinearPolicyGraph(
            SDDP.parameterize(subproblem, Ω, P) do ω
                JuMP.fix(x.out, ω[1])
                @stageobjective(subproblem, ω[2] * x.out + ω[3])
-               println("ω is: ", ω)
            end
        end;
 
-julia> SDDP.simulate(model, 1);
+julia> for s in SDDP.simulate(model, 1)[1]
+           println("ω is: ", s[:noise_term])
+       end
 ω is: [54, 38, 19]
 ω is: [43, 3, 13]
 ω is: [43, 4, 17]

docs/src/guides/add_noise_in_the_constraint_matrix.md

@@ -15,20 +15,15 @@ julia> model = SDDP.LinearPolicyGraph(
                stages=3, lower_bound = 0, optimizer = HiGHS.Optimizer
                ) do subproblem, t
            @variable(subproblem, x, SDDP.State, initial_value = 0.0)
-           @constraint(subproblem, emissions, 1x.out <= 1)
+           @constraint(subproblem, emissions, 1 * x.out <= 1)
            SDDP.parameterize(subproblem, [0.2, 0.5, 1.0]) do ω
+               ## rewrite 1 * x.out <= 1 to ω * x.out <= 1
                JuMP.set_normalized_coefficient(emissions, x.out, ω)
-               println(emissions)
            end
            @stageobjective(subproblem, -x.out)
        end
 A policy graph with 3 nodes.
  Node indices: 1, 2, 3
-
-julia> SDDP.simulate(model, 1);
-emissions : x_out <= 1
-emissions : 0.2 x_out <= 1
-emissions : 0.5 x_out <= 1
 ```
 
 !!! note

src/plugins/duality_handlers.jl

@@ -179,6 +179,16 @@ function get_dual_solution(node::Node, lagrange::LagrangianDuality)
         x_in_value[i] = JuMP.fix_value(state.in)
         h_expr[i] = @expression(node.subproblem, state.in - x_in_value[i])
         JuMP.unfix(state.in)
+        l, u, is_integer = node.incoming_state_bounds[key]
+        if l > -Inf
+            JuMP.set_lower_bound(state.in, l)
+        end
+        if u < Inf
+            JuMP.set_upper_bound(state.in, u)
+        end
+        if is_integer
+            JuMP.set_integer(state.in)
+        end
         λ_star[i] = conic_dual[key]
     end
     # Check that the conic dual is feasible for the subproblem. Sometimes it
@@ -194,6 +204,9 @@ function get_dual_solution(node::Node, lagrange::LagrangianDuality)
             return L_k === nothing ? nothing : (s * L_k, s * h_k)
         end
     for (i, (_, state)) in enumerate(node.states)
+        if JuMP.is_integer(state.in)
+            JuMP.unset_integer(state.in)
+        end
         JuMP.fix(state.in, x_in_value[i]; force = true)
     end
     λ_solution = Dict{Symbol,Float64}(

src/user_interface.jl

@@ -655,6 +655,8 @@ mutable struct Node{T}
     ext::Dict{Symbol,Any}
     # Lock for threading
     lock::ReentrantLock
+    # (lower, upper, is_integer)
+    incoming_state_bounds::Dict{Symbol,Tuple{Float64,Float64,Bool}}
 end
 
 function Base.show(io::IO, node::Node)
@@ -987,6 +989,7 @@ function PolicyGraph(
             # The extension dictionary.
             Dict{Symbol,Any}(),
             ReentrantLock(),
+            Dict{Symbol,Tuple{Float64,Float64,Bool}}(),
         )
         subproblem.ext[:sddp_policy_graph] = policy_graph
         policy_graph.nodes[node_index] = subproblem.ext[:sddp_node] = node
@@ -1032,9 +1035,83 @@ function PolicyGraph(
     if length(graph.belief_partition) > 0
         initialize_belief_states(policy_graph, graph)
     end
+    domain = _get_incoming_domain(policy_graph)
+    for (node_name, node) in policy_graph.nodes
+        for (k, v) in domain[node_name]
+            node.incoming_state_bounds[k] = something(v, (-Inf, Inf, false))
+        end
+    end
     return policy_graph
 end
 
+function _get_incoming_domain(model::PolicyGraph{T}) where {T}
+    function _bounds(x)
+        l, u = -Inf, Inf
+        if has_lower_bound(x)
+            l = lower_bound(x)
+        end
+        if has_upper_bound(x)
+            u = upper_bound(x)
+        end
+        if is_fixed(x)
+            l = u = fix_value(x)
+        end
+        is_int = is_integer(x)
+        if is_binary(x)
+            l, u = max(something(l, 0.0), 0.0), min(something(u, 1.0), 1.0)
+            is_int = true
+        end
+        return l, u, is_int
+    end
+    outgoing_bounds = Dict{Tuple{T,Symbol},Any}(
+        (k, state_name) => nothing for (k, node) in model.nodes for
+        (state_name, _) in node.states
+    )
+    for (k, node) in model.nodes
+        for noise in node.noise_terms
+            parameterize(node, noise.term)
+            for (state_name, state) in node.states
+                domain = outgoing_bounds[(k, state_name)]
+                l_new, u_new, is_int_new = _bounds(state.out)
+                outgoing_bounds[(k, state_name)] = if domain === nothing
+                    (l_new, u_new, is_int_new)
+                else
+                    l, u, is_int = domain::Tuple{Float64,Float64,Bool}
+                    (min(l, l_new), max(u, u_new), is_int & is_int_new)
+                end
+            end
+        end
+    end
+    incoming_bounds = Dict{T,Dict{Symbol,Any}}()
+    for (k, node) in model.nodes
+        incoming_bounds[k] = Dict{Symbol,Any}(
+            state_name => nothing for (state_name, _) in node.states
+        )
+    end
+    for (parent_name, parent) in model.nodes
+        for (state_name, state) in parent.states
+            domain_new = outgoing_bounds[(parent_name, state_name)]
+            l_new, u_new, is_int_new = domain_new::Tuple{Float64,Float64,Bool}
+            for child in parent.children
+                domain = incoming_bounds[child.term][state_name]
+                incoming_bounds[child.term][state_name] = if domain === nothing
+                    (l_new, u_new, is_int_new)
+                else
+                    l, u, is_int = domain::Tuple{Float64,Float64,Bool}
+                    (min(l, l_new), max(u, u_new), is_int & is_int_new)
+                end
+            end
+        end
+    end
+    # The incoming state from the root node can be anything
+    for (state_name, value) in model.initial_root_state
+        for child in model.root_children
+            incoming_bounds[child.term][state_name] = (-Inf, Inf, false)
+        end
+    end
+    return incoming_bounds
+end
+
 # Internal function: set up ::BeliefState for each node.
 function initialize_belief_states(
     policy_graph::PolicyGraph{T},

test/plugins/sampling_schemes.jl

@@ -380,8 +380,8 @@ function test_SimulatorSamplingScheme_with_noise()
     ) do sp, node
         t, price = node
         @variable(sp, 0 <= x <= 1, SDDP.State, initial_value = 0)
-        SDDP.parameterize(sp, [(price, i) for i in 1:2]) do ω
-            return SDDP.@stageobjective(sp, price * x.out + i)
+        SDDP.parameterize(sp, [(price, i) for i in 1:2]) do (p, i)
+            return SDDP.@stageobjective(sp, p * x.out + i)
         end
     end
     sampler = SDDP.SimulatorSamplingScheme(simulator)

test/user_interface.jl

@@ -305,7 +305,8 @@ function test_parameterize()
     end
     node = model[2]
     @test length(node.noise_terms) == 3
-    @test JuMP.upper_bound(node.subproblem[:x]) == 1
+    # SDDP is allowed to call `parameterize` during construction
+    # @test JuMP.upper_bound(node.subproblem[:x]) == 1
     node.parameterize(node.noise_terms[2].term)
     @test JuMP.upper_bound(node.subproblem[:x]) == 2
     node.parameterize(3)
@@ -486,20 +487,19 @@ function test_numerical_stability_report()
 end
 
 function test_objective_state()
-    model = SDDP.LinearPolicyGraph(;
-        stages = 2,
-        lower_bound = 0,
-        optimizer = HiGHS.Optimizer,
-    ) do subproblem, t
-        @variable(subproblem, x, SDDP.State, initial_value = 0)
-        SDDP.parameterize(subproblem, [1, 2]) do ω
-            price = SDDP.objective_state(subproblem)
-            @stageobjective(subproblem, price * x.out)
-        end
-    end
     @test_throws(
         ErrorException("No objective state defined."),
-        SDDP.simulate(model, 1; parallel_scheme = SDDP.Serial()),
+        SDDP.LinearPolicyGraph(;
+            stages = 2,
+            lower_bound = 0,
+            optimizer = HiGHS.Optimizer,
+        ) do subproblem, t
+            @variable(subproblem, x, SDDP.State, initial_value = 0)
+            SDDP.parameterize(subproblem, [1, 2]) do ω
+                price = SDDP.objective_state(subproblem)
+                @stageobjective(subproblem, price * x.out)
+            end
+        end,
     )
     @test_throws(
         ErrorException("add_objective_state can only be called once."),
@@ -911,6 +911,71 @@ function test_print_problem_statistics()
     return
 end
 
+function test_incoming_state_bounds()
+    graph = SDDP.Graph(0)
+    SDDP.add_node.((graph,), 1:4)
+    SDDP.add_edge(graph, 0 => 4, 1.0)
+    SDDP.add_edge(graph, 4 => 1, 0.5)
+    SDDP.add_edge(graph, 4 => 2, 0.5)
+    SDDP.add_edge(graph, 1 => 2, 1.0)
+    SDDP.add_edge(graph, 2 => 3, 0.5)
+    SDDP.add_edge(graph, 3 => 1, 0.5)
+    SDDP.add_edge(graph, 3 => 4, 0.5)
+    model = SDDP.PolicyGraph(graph; lower_bound = 0.0) do sp, node
+        @variable(sp, x, Int, SDDP.State, initial_value = 0)
+        @variable(sp, y, SDDP.State, initial_value = -1)
+        @variable(sp, z, SDDP.State, initial_value = 1)
+        if node == 1
+            set_binary(z.out)
+            SDDP.parameterize(sp, 1:2) do w
+                JuMP.set_lower_bound(x.out, -w)
+                JuMP.set_upper_bound(y.out, w)
+                return
+            end
+        elseif node == 2
+            SDDP.parameterize(sp, 1:2) do w
+                if w == 1
+                    set_upper_bound(y.out, 1.0)
+                elseif w == 2 && has_upper_bound(y.out)
+                    delete_upper_bound(y.out)
+                end
+                return
+            end
+        elseif node == 3
+            set_lower_bound(x.out, 1)
+            set_upper_bound(x.out, 2)
+            fix(y.out, 3)
+        elseif node == 4
+            fix(x.out, 0)
+            fix(y.out, -1)
+            fix(z.out, 1)
+        end
+        @stageobjective(sp, 0)
+        return
+    end
+    @test model[1].incoming_state_bounds == Dict(
+        :x => (0.0, 2.0, true),
+        :y => (-1.0, 3.0, false),
+        :z => (-Inf, Inf, false),
+    )
+    @test model[2].incoming_state_bounds == Dict(
+        :x => (-2.0, Inf, true),
+        :y => (-Inf, 2.0, false),
+        :z => (0.0, 1.0, false),
+    )
+    @test model[3].incoming_state_bounds == Dict(
+        :x => (-Inf, Inf, true),
+        :y => (-Inf, Inf, false),
+        :z => (-Inf, Inf, false),
+    )
+    @test model[4].incoming_state_bounds == Dict(
+        :x => (-Inf, Inf, false),
+        :y => (-Inf, Inf, false),
+        :z => (-Inf, Inf, false),
+    )
+    return
+end
+
 end  # module
 
 TestUserInterface.runtests()