StochDynamicProgramming: the multistock problem

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 SDDP.jl contributors, 2017-23

Problem
  Nodes           : 5
  State variables : 3
  Scenarios       : 1.43489e+07
  Existing cuts   : false
  Subproblem structure                      : (min, max)
    Variables                               : (13, 13)
    VariableRef in MOI.LessThan{Float64}    : (6, 7)
    AffExpr in MOI.LessThan{Float64}        : (1, 1)
    VariableRef in MOI.GreaterThan{Float64} : (7, 7)
    AffExpr in MOI.EqualTo{Float64}         : (3, 3)
Options
  Solver          : serial mode
  Risk measure    : SDDP.Expectation()
  Sampling scheme : SDDP.InSampleMonteCarlo

Numerical stability report
  Non-zero Matrix range     [1e+00, 1e+00]
  Non-zero Objective range  [3e-01, 2e+00]
  Non-zero Bounds range     [5e-01, 5e+00]
  Non-zero RHS range        [2e+00, 2e+00]
No problems detected

 Iteration    Simulation       Bound         Time (s)    Proc. ID   # Solves
       10   -4.487696e+00  -4.443335e+00   1.100810e-01          1       1400
       20   -4.842991e+00  -4.388941e+00   2.392249e-01          1       2800
       30   -4.193961e+00  -4.375300e+00   3.789530e-01          1       4200
       40   -4.945658e+00  -4.369811e+00   5.202959e-01          1       5600
       50   -4.805113e+00  -4.366601e+00   6.746209e-01          1       7000
       60   -4.397496e+00  -4.362569e+00   8.284979e-01          1       8400
       70   -4.350929e+00  -4.358615e+00   9.848108e-01          1       9800
       80   -4.819443e+00  -4.354454e+00   1.259897e+00          1      11200
       90   -4.519670e+00  -4.352503e+00   1.424106e+00          1      12600
      100   -3.893235e+00  -4.351248e+00   1.595077e+00          1      14000

Terminating training
  Status         : iteration_limit
  Total time (s) : 1.595077e+00
  Total solves   : 14000
  Best bound     : -4.351248e+00
  Simulation CI  : -4.352150e+00 ± 8.487769e-02
------------------------------------------------------------------------------

This page was generated using Literate.jl.