Hydro-thermal with inner approximation

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

This is a simple version of the hydro-thermal scheduling problem. The goal is to operate one hydro-dam and two thermal plants over time in the face of inflow uncertainty.

The purpose of this tutorial is to provide a demonstration of the experimental SDDP.Inner submodule.

using SDDPusing Testimport HiGHSimport Randomfunction test_inner_hydro_1d()    Random.seed!(2)    stages = 4    Ω = max.(40 .+ 20.0 * randn(10), 0.0)    risk_measure = SDDP.EAVaR(; lambda = 0.5, beta = 0.2)    function build_subproblem(sp::Model, t::Int)        @variable(sp, 0 <= x_vol <= 100, SDDP.State, initial_value = 83.222)        @variable(sp, 0 <= u_gh <= 60)        @variable(sp, 0 <= u_spill <= 200)        @variable(sp, 0 <= u_gt[1:2] <= 15)        @variable(sp, 0 <= u_def <= 75)        @variable(sp, w_inflow == 0.0)        @constraint(sp, sum(u_gt) + u_gh + u_def == 75)        @constraint(sp, x_vol.in + w_inflow - u_gh - u_spill == x_vol.out)        if t > 1            SDDP.parameterize(w -> fix(w_inflow, w), sp, Ω)        end        @stageobjective(sp, u_spill + 5 * u_gt[1] + 10 * u_gt[2] + 50 * u_def)    end    println("Building and solving primal outer model for lower bounds")    model = SDDP.LinearPolicyGraph(        build_subproblem;        stages,        sense = :Min,        optimizer = HiGHS.Optimizer,        lower_bound = 0.0,    )    SDDP.train(model; iteration_limit = 10, risk_measure)    lower_bound = SDDP.calculate_bound(model)    simulations =        SDDP.simulate(model, 500, [:x_vol, :u_gh, :u_spill, :u_gt, :u_def])    objs = [sum(data[:stage_objective] for data in sim) for sim in simulations]    μ, ci = round.(SDDP.confidence_interval(objs); digits = 2)    println("Building and solving inner model for upper bounds:")    inner_model, upper_bound = SDDP.Inner.inner_dp(        build_subproblem,        model;        stages,        sense = :Min,        optimizer = HiGHS.Optimizer,        lower_bound = 0.0,        risk_measure,        bellman_function = SDDP.Inner.InnerBellmanFunction(            t -> 50.0 * (stages - t);            upper_bound = t -> 75.0 * 50.0 * (stages - t),            vertex_type = SDDP.SINGLE_CUT,        ),    )    @test lower_bound <= upper_bound    println()    println("Bounds:")    println("  Risk-neutral confidence interval: ", μ, " ± ", ci)    println("  Risk-adjusted lower bound: ", round(lower_bound; digits = 2))    println("  Risk-adjusted upper bound: ", round(upper_bound; digits = 2))    returnendtest_inner_hydro_1d()

Building and solving primal outer model for lower bounds
-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-26
-------------------------------------------------------------------
problem
  nodes           : 4
  state variables : 1
  scenarios       : 1.00000e+03
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : A convex combination of 0.5 * SDDP.Expectation() + 0.5 * SDDP.AVaR(0.2)
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [9, 9]
  AffExpr in MOI.EqualTo{Float64}         : [2, 2]
  VariableRef in MOI.EqualTo{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  [1e+00, 5e+01]
  bounds range     [2e+01, 2e+02]
  rhs range        [8e+01, 8e+01]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
         1   1.948878e+03  2.847167e+03  2.048397e-02        35   1
        10   7.500000e+02  2.935390e+03  1.234879e-01       350   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 1.234879e-01
total solves   : 350
best bound     :  2.935390e+03
numeric issues : 0
-------------------------------------------------------------------

Building and solving inner model for upper bounds:
Node: 3 - elapsed time: 0.05 plus 0.02 for vertex selection.
Node: 2 - elapsed time: 0.01 plus 0.0 for vertex selection.
Node: 1 - elapsed time: 0.0 plus 0.0 for vertex selection.
First-stage upper bound: 2969.680973503913
Total time for upper bound: 0.089246535

Bounds:
  Risk-neutral confidence interval: 1411.99 ± 82.02
  Risk-adjusted lower bound: 2935.39
  Risk-adjusted upper bound: 2969.68