SLDP: example 2

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

This example is derived from Section 4.3 of the paper: Ahmed, S., Cabral, F. G., & da Costa, B. F. P. (2019). Stochastic Lipschitz Dynamic Programming. Optimization Online. PDF

using SDDP
import HiGHS
import Test

function sldp_example_two(; first_stage_integer::Bool = true, N = 2)
    model = SDDP.LinearPolicyGraph(
        stages = 2,
        lower_bound = -100.0,
        optimizer = HiGHS.Optimizer,
    ) do sp, t
        @variable(sp, 0 <= x[1:2] <= 5, SDDP.State, initial_value = 0.0)
        if t == 1
            if first_stage_integer
                @variable(sp, 0 <= u[1:2] <= 5, Int)
                @constraint(sp, [i = 1:2], u[i] == x[i].out)
            end
            @stageobjective(sp, -1.5 * x[1].out - 4 * x[2].out)
        else
            @variable(sp, 0 <= y[1:4] <= 1, Bin)
            @variable(sp, ω[1:2])
            @stageobjective(sp, -16 * y[1] - 19 * y[2] - 23 * y[3] - 28 * y[4])
            @constraint(
                sp,
                2 * y[1] + 3 * y[2] + 4 * y[3] + 5 * y[4] <= ω[1] - x[1].in
            )
            @constraint(
                sp,
                6 * y[1] + 1 * y[2] + 3 * y[3] + 2 * y[4] <= ω[2] - x[2].in
            )
            steps = range(5, stop = 15, length = N)
            SDDP.parameterize(sp, [[i, j] for i in steps for j in steps]) do φ
                return JuMP.fix.(ω, φ)
            end
        end
    end
    if get(ARGS, 1, "") == "--write"
        # Run `$ julia sldp_example_two.jl --write` to update the benchmark
        # model directory
        model_dir = joinpath(@__DIR__, "..", "..", "..", "benchmarks", "models")
        SDDP.write_to_file(
            model,
            joinpath(model_dir, "sldp_example_two_$(N).sof.json.gz");
            test_scenarios = 30,
        )
        return
    end
    SDDP.train(model; log_frequency = 10)
    bound = SDDP.calculate_bound(model)

    if N == 2
        Test.@test bound <= -57.0
    elseif N == 3
        Test.@test bound <= -59.33
    elseif N == 6
        Test.@test bound <= -61.22
    end
    return
end

sldp_example_two(N = 2)
sldp_example_two(N = 3)
sldp_example_two(N = 6)
-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
problem
  nodes           : 2
  state variables : 2
  scenarios       : 4.00000e+00
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [7, 11]
  AffExpr in MOI.EqualTo{Float64}         : [2, 2]
  AffExpr in MOI.LessThan{Float64}        : [2, 2]
  VariableRef in MOI.GreaterThan{Float64} : [5, 7]
  VariableRef in MOI.Integer              : [2, 2]
  VariableRef in MOI.LessThan{Float64}    : [4, 7]
  VariableRef in MOI.ZeroOne              : [4, 4]
numerical stability report
  matrix range     [1e+00, 6e+00]
  objective range  [1e+00, 3e+01]
  bounds range     [1e+00, 1e+02]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10  -9.800000e+01 -5.809615e+01  3.080893e-02        78   1
        20  -4.000000e+01 -5.809615e+01  9.104109e-02       148   1
        30  -4.000000e+01 -5.809615e+01  1.301620e-01       226   1
        40  -4.700000e+01 -5.809615e+01  1.640501e-01       296   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 1.640501e-01
total solves   : 296
best bound     : -5.809615e+01
simulation ci  : -5.381250e+01 ± 7.419419e+00
numeric issues : 0
-------------------------------------------------------------------

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
problem
  nodes           : 2
  state variables : 2
  scenarios       : 9.00000e+00
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [7, 11]
  AffExpr in MOI.EqualTo{Float64}         : [2, 2]
  AffExpr in MOI.LessThan{Float64}        : [2, 2]
  VariableRef in MOI.GreaterThan{Float64} : [5, 7]
  VariableRef in MOI.Integer              : [2, 2]
  VariableRef in MOI.LessThan{Float64}    : [4, 7]
  VariableRef in MOI.ZeroOne              : [4, 4]
numerical stability report
  matrix range     [1e+00, 6e+00]
  objective range  [1e+00, 3e+01]
  bounds range     [1e+00, 1e+02]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10  -4.000000e+01 -6.196125e+01  3.954792e-02       138   1
        20  -8.200000e+01 -6.196125e+01  7.803988e-02       258   1
        30  -4.700000e+01 -6.196125e+01  1.308970e-01       396   1
        40  -4.000000e+01 -6.196125e+01  1.681318e-01       516   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 1.681318e-01
total solves   : 516
best bound     : -6.196125e+01
simulation ci  : -5.381250e+01 ± 4.871988e+00
numeric issues : 0
-------------------------------------------------------------------

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
problem
  nodes           : 2
  state variables : 2
  scenarios       : 3.60000e+01
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [7, 11]
  AffExpr in MOI.EqualTo{Float64}         : [2, 2]
  AffExpr in MOI.LessThan{Float64}        : [2, 2]
  VariableRef in MOI.GreaterThan{Float64} : [5, 7]
  VariableRef in MOI.Integer              : [2, 2]
  VariableRef in MOI.LessThan{Float64}    : [4, 7]
  VariableRef in MOI.ZeroOne              : [4, 4]
numerical stability report
  matrix range     [1e+00, 6e+00]
  objective range  [1e+00, 3e+01]
  bounds range     [1e+00, 1e+02]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10  -6.300000e+01 -6.546793e+01  7.486200e-02       462   1
        20  -5.400000e+01 -6.546793e+01  1.337669e-01       852   1
        30  -7.000000e+01 -6.546793e+01  2.476978e-01      1314   1
        40  -5.900000e+01 -6.546793e+01  3.076539e-01      1704   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 3.076539e-01
total solves   : 1704
best bound     : -6.546793e+01
simulation ci  : -6.148750e+01 ± 5.461314e+00
numeric issues : 0
-------------------------------------------------------------------