Infinite horizon hydro-thermal

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

using SDDP, HiGHS, Test, Statistics

function infinite_hydro_thermal(; cut_type)
    Ω = [
        (inflow = 0.0, demand = 7.5),
        (inflow = 5.0, demand = 5),
        (inflow = 10.0, demand = 2.5),
    ]
    graph = SDDP.Graph(
        :root_node,
        [:week],
        [(:root_node => :week, 1.0), (:week => :week, 0.9)],
    )
    model = SDDP.PolicyGraph(
        graph;
        lower_bound = 0,
        optimizer = HiGHS.Optimizer,
    ) do subproblem, node
        @variable(
            subproblem,
            5.0 <= reservoir <= 15.0,
            SDDP.State,
            initial_value = 10.0
        )
        @variables(subproblem, begin
            thermal_generation >= 0
            hydro_generation >= 0
            spill >= 0
            inflow
            demand
        end)
        @constraints(
            subproblem,
            begin
                reservoir.out == reservoir.in - hydro_generation - spill + inflow
                hydro_generation + thermal_generation == demand
            end
        )
        @stageobjective(subproblem, 10 * spill + thermal_generation)
        SDDP.parameterize(subproblem, Ω) do ω
            JuMP.fix(inflow, ω.inflow)
            return JuMP.fix(demand, ω.demand)
        end
    end
    SDDP.train(
        model;
        cut_type = cut_type,
        log_frequency = 100,
        sampling_scheme = SDDP.InSampleMonteCarlo(; terminate_on_cycle = true),
        parallel_scheme = SDDP.Serial(),
        cycle_discretization_delta = 0.1,
    )
    @test SDDP.calculate_bound(model) ≈ 119.167 atol = 0.1

    results = SDDP.simulate(model, 500)
    objectives =
        [sum(s[:stage_objective] for s in simulation) for simulation in results]
    sample_mean = round(Statistics.mean(objectives); digits = 2)
    sample_ci = round(1.96 * Statistics.std(objectives) / sqrt(500); digits = 2)
    println("Confidence_interval = $(sample_mean) ± $(sample_ci)")
    @test sample_mean - sample_ci <= 119.167 <= sample_mean + sample_ci
    return
end

infinite_hydro_thermal(; cut_type = SDDP.SINGLE_CUT)
infinite_hydro_thermal(; cut_type = SDDP.MULTI_CUT)

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
  nodes           : 1
  state variables : 1
  scenarios       : Inf
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [8, 8]
  AffExpr in MOI.EqualTo{Float64}         : [2, 2]
  VariableRef in MOI.GreaterThan{Float64} : [5, 5]
  VariableRef in MOI.LessThan{Float64}    : [1, 1]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [1e+00, 1e+01]
  bounds range     [5e+00, 2e+01]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
       100   2.500000e+01  1.188965e+02  8.096220e-01      1946   1
       200   2.500000e+01  1.191634e+02  1.015996e+00      3920   1
       300   0.000000e+00  1.191666e+02  1.228354e+00      5902   1
       330   2.500000e+01  1.191667e+02  1.270023e+00      6224   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 1.270023e+00
total solves   : 6224
best bound     :  1.191667e+02
simulation ci  :  2.158333e+01 ± 3.290252e+00
numeric issues : 0
-------------------------------------------------------------------

Confidence_interval = 120.33 ± 12.06
-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
  nodes           : 1
  state variables : 1
  scenarios       : Inf
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [8, 8]
  AffExpr in MOI.EqualTo{Float64}         : [2, 2]
  VariableRef in MOI.GreaterThan{Float64} : [5, 5]
  VariableRef in MOI.LessThan{Float64}    : [1, 1]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [1e+00, 1e+01]
  bounds range     [5e+00, 2e+01]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
       100   0.000000e+00  1.191285e+02  2.950759e-01      2874   1
       200   2.500000e+00  1.191666e+02  5.294189e-01      4855   1
       282   7.500000e+00  1.191667e+02  6.648269e-01      5733   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 6.648269e-01
total solves   : 5733
best bound     :  1.191667e+02
simulation ci  :  2.104610e+01 ± 3.492245e+00
numeric issues : 0
-------------------------------------------------------------------

Confidence_interval = 116.06 ± 13.65