StructDualDynProg: Problem 5.2, 3 stages

This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.

This example comes from StochasticDualDynamicProgramming.jl.

using SDDP, HiGHS, Test

function test_prob52_3stages()
    model = SDDP.LinearPolicyGraph(
        stages = 3,
        lower_bound = 0.0,
        optimizer = HiGHS.Optimizer,
    ) do sp, t
        n = 4
        m = 3
        i_c = [16, 5, 32, 2]
        C = [25, 80, 6.5, 160]
        T = [8760, 7000, 1500] / 8760
        D2 = [diff([0, 3919, 7329, 10315]) diff([0, 7086, 9004, 11169])]
        p2 = [0.9, 0.1]
        @variable(sp, x[i = 1:n] >= 0, SDDP.State, initial_value = 0.0)
        @variables(sp, begin
            y[1:n, 1:m] >= 0
            v[1:n] >= 0
            penalty >= 0
            ξ[j = 1:m]
        end)
        @constraints(sp, begin
            [i = 1:n], x[i].out == x[i].in + v[i]
            [i = 1:n], sum(y[i, :]) <= x[i].in
            [j = 1:m], sum(y[:, j]) + penalty >= ξ[j]
        end)
        @stageobjective(sp, i_c'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
                    JuMP.fix(ξ[j], D2[j, ω])
                end
            end
        end
        if t == 3
            @constraint(sp, sum(v) == 0)
        end
    end

    det = SDDP.deterministic_equivalent(model, HiGHS.Optimizer)
    set_silent(det)
    JuMP.optimize!(det)
    @test JuMP.objective_value(det) ≈ 406712.49 atol = 0.1

    SDDP.train(model; log_frequency = 10)
    @test SDDP.calculate_bound(model) ≈ 406712.49 atol = 0.1
    return
end

test_prob52_3stages()

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
problem
  nodes           : 3
  state variables : 4
  scenarios       : 4.00000e+00
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [29, 29]
  AffExpr in MOI.EqualTo{Float64}         : [4, 5]
  AffExpr in MOI.GreaterThan{Float64}     : [3, 3]
  AffExpr in MOI.LessThan{Float64}        : [4, 4]
  VariableRef in MOI.EqualTo{Float64}     : [3, 3]
  VariableRef in MOI.GreaterThan{Float64} : [22, 22]
  VariableRef in MOI.LessThan{Float64}    : [1, 1]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [1e+00, 1e+05]
  bounds range     [0e+00, 0e+00]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10   4.403329e+05  3.509666e+05  1.354313e-02        92   1
        20   4.055335e+05  4.054833e+05  2.418017e-02       172   1
        30   3.959476e+05  4.067125e+05  3.695416e-02       264   1
        40   3.959476e+05  4.067125e+05  4.958916e-02       344   1
        47   3.959476e+05  4.067125e+05  5.932808e-02       400   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 5.932808e-02
total solves   : 400
best bound     :  4.067125e+05
simulation ci  :  2.695623e+07 ± 3.645336e+07
numeric issues : 0
-------------------------------------------------------------------