StochDynamicProgramming: the multistock problem

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 StochDynamicProgramming.jl.

using SDDP, HiGHS, Test

function test_multistock_example()
    model = SDDP.LinearPolicyGraph(
        stages = 5,
        lower_bound = -5.0,
        optimizer = HiGHS.Optimizer,
    ) do subproblem, stage
        @variable(
            subproblem,
            0 <= stock[i = 1:3] <= 1,
            SDDP.State,
            initial_value = 0.5
        )
        @variables(subproblem, begin
            0 <= control[i = 1:3] <= 0.5
            ξ[i = 1:3]  # Dummy for RHS noise.
        end)
        @constraints(
            subproblem,
            begin
                sum(control) - 0.5 * 3 <= 0
                [i = 1:3], stock[i].out == stock[i].in + control[i] - ξ[i]
            end
        )
        Ξ = collect(
            Base.product((0.0, 0.15, 0.3), (0.0, 0.15, 0.3), (0.0, 0.15, 0.3)),
        )[:]
        SDDP.parameterize(subproblem, Ξ) do ω
            return JuMP.fix.(ξ, ω)
        end
        @stageobjective(subproblem, (sin(3 * stage) - 1) * sum(control))
    end
    SDDP.train(
        model;
        iteration_limit = 100,
        cut_type = SDDP.SINGLE_CUT,
        log_frequency = 10,
    )
    @test SDDP.calculate_bound(model) ≈ -4.349 atol = 0.01

    simulation_results = SDDP.simulate(model, 5000)
    @test length(simulation_results) == 5000
    μ = SDDP.Statistics.mean(
        sum(data[:stage_objective] for data in simulation) for
        simulation in simulation_results
    )
    @test μ ≈ -4.349 atol = 0.1
    return
end

test_multistock_example()

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
problem
  nodes           : 5
  state variables : 3
  scenarios       : 1.43489e+07
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [13, 13]
  AffExpr in MOI.EqualTo{Float64}         : [3, 3]
  AffExpr in MOI.LessThan{Float64}        : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [7, 7]
  VariableRef in MOI.LessThan{Float64}    : [6, 7]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [3e-01, 2e+00]
  bounds range     [5e-01, 5e+00]
  rhs range        [2e+00, 2e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10  -3.950563e+00 -4.425944e+00  1.837828e-01      1400   1
        20  -4.274542e+00 -4.397454e+00  2.969389e-01      2800   1
        30  -3.068570e+00 -4.380102e+00  4.216728e-01      4200   1
        40  -3.786070e+00 -4.369668e+00  5.492370e-01      5600   1
        50  -4.259563e+00 -4.362339e+00  6.916220e-01      7000   1
        60  -3.647455e+00 -4.359239e+00  8.365018e-01      8400   1
        70  -4.012512e+00 -4.357262e+00  9.828730e-01      9800   1
        80  -4.293203e+00 -4.354708e+00  1.125953e+00     11200   1
        90  -4.485010e+00 -4.352846e+00  1.272223e+00     12600   1
       100  -4.178952e+00 -4.350527e+00  1.419191e+00     14000   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 1.419191e+00
total solves   : 14000
best bound     : -4.350527e+00
simulation ci  : -4.259588e+00 ± 8.809614e-02
numeric issues : 0
-------------------------------------------------------------------