Hash risk-set cuts in multi-cut

docs/src/tutorial/16_performance.md

@@ -37,16 +37,21 @@ speed-ups.
   algorithm (e.g., [`Method=1` in Gurobi](https://www.gurobi.com/documentation/8.1/refman/method.html)
   ) is usually a performance win.
 
-## Average cut vs multi-cut
+### Single-cut vs. multi-cut
 
-There are two competing ways that cuts can be created in SDDP: _average_ cut and
-_multi_-cut. These can be specified as follows
-```julia
-SDDP.train(model; cut_type = SDDP.AVERAGE_CUT)
-SDDP.train(model; cut_type = SDDP.MULTI_CUT)
-```
+There are two competing ways that cuts can be created in SDDP: _single_-cut and
+_multi_-cut. By default, `SDDP.jl` uses the _single-cut_ version of SDDP.
 
 The performance of each method is problem-dependent. We recommend that you try
-both in order to see which one performs better. In general, the _average_ cut
+both in order to see which one performs better. In general, the _single_-cut
 method works better when the number of realizations of the stagewise-independent
-random variable is large.
+random variable is large, whereas the multi-cut method works better on small
+problems. However, the multi-cut method can cause numerical stability problems,
+particularly if used in conjunction with objective or belief state variables.
+
+You can switch between the methods by passing the relevant flag to `cut_type` in
+[`SDDP.train`](@ref).
+```julia
+SDDP.train(model; cut_type = SDDP.SINGLE_CUT)
+SDDP.train(model; cut_type = SDDP.MULTI_CUT)
+```

examples/FAST_production_management.jl

@@ -39,5 +39,5 @@ function fast_production_management(;cut_type)
     @test SDDP.calculate_bound(model) ≈ -23.96 atol = 1e-2
 end
 
-fast_production_management(cut_type = SDDP.AVERAGE_CUT)
+fast_production_management(cut_type = SDDP.SINGLE_CUT)
 fast_production_management(cut_type = SDDP.MULTI_CUT)

examples/StochDynamicProgramming.jl_multistock.jl

@@ -21,10 +21,10 @@
 using SDDP, GLPK, Test
 
 function test_multistock_example()
-    model = SDDP.PolicyGraph(SDDP.LinearGraph(5),
-            optimizer = with_optimizer(GLPK.Optimizer),
-            bellman_function = SDDP.BellmanFunction(lower_bound = -5)
-                                    ) do subproblem, stage
+    model = SDDP.LinearPolicyGraph(
+            stages = 5,
+            lower_bound = -5.0,
+            optimizer = with_optimizer(GLPK.Optimizer)) do subproblem, stage
         @variable(subproblem,
             0 <= stock[i=1:3] <= 1, SDDP.State, initial_value = 0.5)
         @variables(subproblem, begin
@@ -47,7 +47,8 @@ function test_multistock_example()
             (sin(3 * stage) - 1) * sum(control)
         )
     end
-    SDDP.train(model, iteration_limit = 100, print_level = 0)
+    SDDP.train(model, iteration_limit = 100, print_level = 0,
+        cut_type = SDDP.SINGLE_CUT)
     @test SDDP.calculate_bound(model) ≈ -4.349 atol = 0.01
 
     simulation_results = SDDP.simulate(model, 5000)

examples/StructDualDynProg.jl_prob5.2_3stages.jl

@@ -11,10 +11,10 @@
 using SDDP, GLPK, Test
 
 function test_prob52_3stages()
-    model = SDDP.PolicyGraph(SDDP.LinearGraph(3),
-                bellman_function = SDDP.BellmanFunction(lower_bound = 0.0),
-                optimizer = with_optimizer(GLPK.Optimizer)
-                ) do sp, t
+    model = SDDP.LinearPolicyGraph(
+                stages = 3,
+                lower_bound = 0.0,
+                optimizer = with_optimizer(GLPK.Optimizer)) do sp, t
         n = 4
         m = 3
         ic = [16, 5, 32, 2]
@@ -34,7 +34,7 @@ function test_prob52_3stages()
             [i=1:n], sum(y[i, :]) <= x[i].in
             [j=1:m], sum(y[:, j]) + penalty >= ξ[j]
         end)
-        @stageobjective(sp, ic'v + C' * y * T + 1e6 * penalty)
+        @stageobjective(sp, ic'v + C' * y * T + 1e5 * penalty)
         if t != 1 # no uncertainty in first stage
             SDDP.parameterize(sp, 1:size(D2, 2), p2) do ω
                 for j in 1:m
@@ -48,9 +48,6 @@ function test_prob52_3stages()
     end
     SDDP.train(model, iteration_limit = 30, print_level = 0)
     @test SDDP.calculate_bound(model) ≈ 406712.49 atol = 0.1
-    # sim = simulate(mod, 1, [:x, :penalty])
-    # @test length(sim) == 1
-    # @test isapprox(sim[1][:x][1], [2986,0,7329,854])
 end
 
 test_prob52_3stages()

examples/asset_management_stagewise.jl

@@ -60,5 +60,5 @@ function asset_management_stagewise(; cut_type)
     @test SDDP.calculate_bound(model) ≈ 1.278 atol=1e-3
 end
 
-asset_management_stagewise(cut_type = SDDP.AVERAGE_CUT)
+asset_management_stagewise(cut_type = SDDP.SINGLE_CUT)
 asset_management_stagewise(cut_type = SDDP.MULTI_CUT)

examples/belief.jl

@@ -47,7 +47,7 @@ function inventory_management_problem()
 
     # Train the policy.
     Random.seed!(123)
-    SDDP.train(model; iteration_limit = 100, print_level = 0)
+    SDDP.train(model; iteration_limit = 100, print_level = 1, cut_type = SDDP.SINGLE_CUT)
     results = SDDP.simulate(model, 500)
     objectives = [
         sum(s[:stage_objective] for s in simulation) for simulation in results

examples/infinite_horizon_hydro_thermal.jl

@@ -56,5 +56,5 @@ function infinite_hydro_thermal(; cut_type)
     @test sample_mean - sample_ci <= 119.167 <= sample_mean + sample_ci
 end
 
-infinite_hydro_thermal(cut_type = SDDP.AVERAGE_CUT)
+infinite_hydro_thermal(cut_type = SDDP.SINGLE_CUT)
 infinite_hydro_thermal(cut_type = SDDP.MULTI_CUT)

examples/objective_state_newsvendor.jl

@@ -68,5 +68,5 @@ function newsvendor_example(;cut_type)
     @test round(Statistics.mean(objectives); digits = 2) ≈ 4.04 atol=0.1
 end
 
-newsvendor_example(cut_type = SDDP.AVERAGE_CUT)
+newsvendor_example(cut_type = SDDP.SINGLE_CUT)
 newsvendor_example(cut_type = SDDP.MULTI_CUT)

src/algorithm.jl

@@ -202,7 +202,9 @@ Write the subproblem contained in `node` to the file `filename`.
 
 `format` should be one of `:mps`, `:lp`, or `:both`.
 """
-function write_subproblem_to_file(node::Node, filename::String; format::Symbol=:both)
+function write_subproblem_to_file(
+        node::Node, filename::String;
+        format::Symbol=:both, throw_error::Bool = false)
     if format ∉ (:mps, :lp, :both)
         error("Invalid `format=$(format)`. Must be `:mps`, `:lp`, or `:both`.")
     end
@@ -216,6 +218,16 @@ function write_subproblem_to_file(node::Node, filename::String; format::Symbol=:
         MOI.copy_to(lp, JuMP.backend(node.subproblem))
         MOI.write_to_file(lp, filename * ".lp")
     end
+    if throw_error
+        error("Unable to retrieve dual solution from ", node.index, ".",
+              "\n  Termination status: ", JuMP.termination_status(node.subproblem),
+              "\n  Primal status:      ", JuMP.primal_status(node.subproblem),
+              "\n  Dual status:        ", JuMP.dual_status(node.subproblem),
+              ".\n An MPS file was written to `subproblem.mps` and an LP file ",
+              "written to `subproblem.lp`. See ",
+              "https://odow.github.io/SDDP.jl/latest/tutorial/06_warnings/#Numerical-stability-1",
+              " for more information.")
+    end
 end
 
 """
@@ -245,29 +257,16 @@ function solve_subproblem(model::PolicyGraph{T},
     parameterize(node, noise)
     JuMP.optimize!(node.subproblem)
     # Test for primal feasibility.
-    primal_status = JuMP.primal_status(node.subproblem)
-    if primal_status != JuMP.MOI.FEASIBLE_POINT
-        write_subproblem_to_file(node, "subproblem")
-        error("Unable to retrieve primal solution from ", node.index, ".",
-              "\n  Termination status: ", JuMP.termination_status(node.subproblem),
-              "\n  Primal status:      ", primal_status,
-              "\n  Dual status:        ", JuMP.dual_status(node.subproblem),
-              ".\n An MPS file was written to `subproblem.mps` and an LP file ",
-              "written to `subproblem.lp`.")
+    if JuMP.primal_status(node.subproblem) != JuMP.MOI.FEASIBLE_POINT
+        write_subproblem_to_file(node, "subproblem", throw_error = true)
     end
     # If require_duals = true, check for dual feasibility and return a dict with
     # the dual on the fixed constraint associated with each incoming state
     # variable. If require_duals=false, return an empty dictionary for
     # type-stability.
     dual_values = if require_duals
-        dual_status = JuMP.dual_status(node.subproblem)
-        if dual_status != JuMP.MOI.FEASIBLE_POINT
-            write_subproblem_to_file(node, "subproblem.mps")
-            error("Unable to retrieve dual solution from ", node.index, ".",
-                  "\n  Termination status: ", JuMP.termination_status(node.subproblem),
-                  "\n  Primal status:      ", primal_status,
-                  "\n  Dual status:        ", dual_status,
-                  ".\n An MPS file was written to `subproblem.mps`")
+        if JuMP.dual_status(node.subproblem) != JuMP.MOI.FEASIBLE_POINT
+            write_subproblem_to_file(node, "subproblem", throw_error = true)
         end
         get_dual_variables(node)
     else
@@ -729,8 +728,8 @@ Train the policy for `model`. Keyword arguments:
 
  - `log_file`: filepath at which to write a log of the training progress
 
- - run_numerical_stability_report: generate a numerical stability report prior
-   to solve
+ - `run_numerical_stability_report`: generate a numerical stability report prior
+    to solve
 
  - `refine_at_similar_nodes::Bool`: if SDDP can detect that two nodes have the
     same children, it can cheaply add a cut discovered at one to the other. In
@@ -744,7 +743,9 @@ Train the policy for `model`. Keyword arguments:
  - `risk_measure`: the risk measure to use at each node.
 
  - `sampling_scheme`: a sampling scheme to use on the forward pass of the
-   algorithm. Defaults to InSampleMonteCarlo().
+    algorithm. Defaults to InSampleMonteCarlo().
+
+ - `cut_type`: choose between `SINGLE_CUT` and `MULTI_CUT` versions of SDDP.
 
 There is also a special option for infinite horizon problems
 
@@ -761,7 +762,7 @@ function train(model::PolicyGraph;
                stopping_rules = AbstractStoppingRule[],
                risk_measure = SDDP.Expectation(),
                sampling_scheme = SDDP.InSampleMonteCarlo(),
-               cut_type = SDDP.AVERAGE_CUT,
+               cut_type = SDDP.SINGLE_CUT,
                cycle_discretization_delta = 0.0,
                refine_at_similar_nodes = true,
                cut_deletion_minimum = 1
@@ -810,7 +811,7 @@ function train(model::PolicyGraph;
         cycle_discretization_delta,
         refine_at_similar_nodes
     )
-    # Update the nodes with the selected cut type (AVERAGE_CUT or MULTI_CUT)
+    # Update the nodes with the selected cut type (SINGLE_CUT or MULTI_CUT)
     # and the cut deletion minimum.
     if cut_deletion_minimum < 0
         cut_deletion_minimum = typemax(Int)