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-24
-------------------------------------------------------------------
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.878303e+00 -4.434982e+00  1.936009e-01      1400   1
        20  -4.262885e+00 -4.399265e+00  3.260260e-01      2800   1
        30  -3.075162e+00 -4.382527e+00  4.524379e-01      4200   1
        40  -3.761147e+00 -4.369587e+00  5.851030e-01      5600   1
        50  -4.323162e+00 -4.362199e+00  7.249870e-01      7000   1
        60  -3.654943e+00 -4.358401e+00  9.212289e-01      8400   1
        70  -4.010883e+00 -4.357368e+00  1.066253e+00      9800   1
        80  -4.314412e+00 -4.355714e+00  1.215074e+00     11200   1
        90  -4.542422e+00 -4.353708e+00  1.369564e+00     12600   1
       100  -4.178952e+00 -4.351685e+00  1.517777e+00     14000   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 1.517777e+00
total solves   : 14000
best bound     : -4.351685e+00
simulation ci  : -4.246786e+00 ± 8.703997e-02
numeric issues : 0
-------------------------------------------------------------------