Add regularization for deterministic first-stage (#624)

docs/src/apireference.md

@@ -79,6 +79,7 @@ SDDP.RevisitingForwardPass
 SDDP.RiskAdjustedForwardPass
 SDDP.AlternativeForwardPass
 SDDP.AlternativePostIterationCallback
+SDDP.RegularizedForwardPass
 ```
 
 ### Risk Measures

src/plugins/forward_passes.jl

@@ -290,3 +290,78 @@ function forward_pass(
     )
     return pass
 end
+
+"""
+    RegularizedForwardPass(;
+        rho::Float64 = 0.05,
+        forward_pass::AbstractForwardPass = DefaultForwardPass(),
+    )
+
+A forward pass that regularizes the outgoing first-stage state variables with an
+L-infty trust-region constraint about the previous iteration's solution.
+Specifically, the bounds of the outgoing state variable `x` are updated from
+`(l, u)` to `max(l, x^k - rho * (u - l)) <= x <= min(u, x^k + rho * (u - l))`,
+where `x^k` is the optimal solution of `x` in the previous iteration. On the
+first iteration, the value of the state at the root node is used.
+
+By default, `rho` is set to 5%, which seems to work well empirically.
+
+Pass a different `forward_pass` to control the forward pass within the
+regularized forward pass.
+
+This forward pass is largely intended to be used for investment problems in
+which the first stage makes a series of capacity decisions that then influence
+the rest of the graph. An error is thrown if the first stage problem is not
+deterministic, and states are silently skipped if they do not have finite
+bounds.
+"""
+mutable struct RegularizedForwardPass{T<:AbstractForwardPass} <:
+               AbstractForwardPass
+    forward_pass::T
+    trial_centre::Dict{Symbol,Float64}
+    ρ::Float64
+
+    function RegularizedForwardPass(;
+        rho::Float64 = 0.05,
+        forward_pass::AbstractForwardPass = DefaultForwardPass(),
+    )
+        centre = Dict{Symbol,Float64}()
+        return new{typeof(forward_pass)}(forward_pass, centre, rho)
+    end
+end
+
+function forward_pass(
+    model::PolicyGraph,
+    options::Options,
+    fp::RegularizedForwardPass,
+)
+    if length(model.root_children) != 1
+        error(
+            "RegularizedForwardPass cannot be applied because first-stage is " *
+            "not deterministic",
+        )
+    end
+    node = model[model.root_children[1].term]
+    if length(node.noise_terms) > 1
+        error(
+            "RegularizedForwardPass cannot be applied because first-stage is " *
+            "not deterministic",
+        )
+    end
+    old_bounds = Dict{Symbol,Tuple{Float64,Float64}}()
+    for (k, v) in node.states
+        if has_lower_bound(v.out) && has_upper_bound(v.out)
+            old_bounds[k] = (l, u) = (lower_bound(v.out), upper_bound(v.out))
+            x = get(fp.trial_centre, k, model.initial_root_state[k])
+            set_lower_bound(v.out, max(l, x - fp.ρ * (u - l)))
+            set_upper_bound(v.out, min(u, x + fp.ρ * (u - l)))
+        end
+    end
+    pass = forward_pass(model, options, fp.forward_pass)
+    for (k, (l, u)) in old_bounds
+        fp.trial_centre[k] = pass.sampled_states[1][k]
+        set_lower_bound(node.states[k].out, l)
+        set_upper_bound(node.states[k].out, u)
+    end
+    return pass
+end

test/plugins/forward_passes.jl

@@ -8,6 +8,7 @@ module TestForwardPasses
 using SDDP
 using Test
 import HiGHS
+import Random
 
 function runtests()
     for name in names(@__MODULE__, all = true)
@@ -222,6 +223,58 @@ function test_DefaultForwardPass_acyclic_include_last_node()
     return
 end
 
+function test_RegularizedForwardPass()
+    function main(capacity_cost, forward_pass, hint)
+        Random.seed!(1245)
+        graph = SDDP.LinearGraph(2)
+        SDDP.add_edge(graph, 2 => 2, 0.95)
+        model = SDDP.PolicyGraph(
+            graph;
+            sense = :Min,
+            lower_bound = 0.0,
+            optimizer = HiGHS.Optimizer,
+        ) do sp, node
+            @variable(sp, 0 <= x <= 400, SDDP.State, initial_value = hint)
+            @variable(sp, 0 <= y, SDDP.State, initial_value = 0)
+            if node == 1
+                @stageobjective(sp, capacity_cost * x.out)
+                @constraint(sp, y.out == y.in)
+            else
+                @variable(sp, 0 <= u_prod <= 200)
+                @variable(sp, u_overtime >= 0)
+                @stageobjective(sp, 100u_prod + 300u_overtime + 50y.out)
+                @constraint(sp, x.out == x.in)
+                @constraint(sp, y.out <= x.in)
+                @constraint(sp, c_bal, y.out == y.in + u_prod + u_overtime)
+                SDDP.parameterize(sp, [100, 300]) do ω
+                    set_normalized_rhs(c_bal, -ω)
+                    return
+                end
+            end
+            return
+        end
+        SDDP.train(
+            model;
+            print_level = 0,
+            forward_pass = forward_pass,
+            iteration_limit = 10,
+        )
+        return SDDP.calculate_bound(model)
+    end
+    for (cost, hint) in [(0, 400), (200, 100), (400, 0)]
+        fp = SDDP.RegularizedForwardPass()
+        reg_bound = main(cost, fp, hint)
+        bound = main(cost, SDDP.DefaultForwardPass(), hint)
+        @test reg_bound >= bound - 1e-6
+    end
+    # Test that initializingn with a bad guess performs poorly
+    fp = SDDP.RegularizedForwardPass()
+    reg_bound = main(400, fp, 400)
+    bound = main(400, SDDP.DefaultForwardPass(), 0)
+    @test reg_bound < bound
+    return
+end
+
 end  # module
 
 TestForwardPasses.runtests()