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-25
-------------------------------------------------------------------
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.EqualTo{Float64}     : [3, 3]
  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  -4.385228e+00 -4.442269e+00  1.901629e-01      1400   1
        20  -3.881547e+00 -4.392566e+00  3.058648e-01      2800   1
        30  -3.671354e+00 -4.376606e+00  4.332879e-01      4200   1
        40  -4.312032e+00 -4.370551e+00  6.244590e-01      5600   1
        50  -4.077684e+00 -4.362538e+00  7.617218e-01      7000   1
        60  -3.609363e+00 -4.357411e+00  9.031339e-01      8400   1
        70  -4.997551e+00 -4.355252e+00  1.049878e+00      9800   1
        80  -4.608038e+00 -4.352744e+00  1.194759e+00     11200   1
        90  -4.949477e+00 -4.350750e+00  1.339477e+00     12600   1
       100  -4.165260e+00 -4.349697e+00  1.490497e+00     14000   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 1.490497e+00
total solves   : 14000
best bound     : -4.349697e+00
simulation ci  : -4.310541e+00 ± 9.196174e-02
numeric issues : 0
-------------------------------------------------------------------