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 SDDP
using Test

import HiGHS
import Random

function 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::JuMP.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 -> JuMP.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))
    return
end

test_inner_hydro_1d()

Building and solving primal outer model for lower bounds
-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-25
-------------------------------------------------------------------
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  9.738922e-03        35   1
        10   7.500000e+02  2.935390e+03  4.066491e-02       350   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 4.066491e-02
total solves   : 350
best bound     :  2.935390e+03
simulation ci  :  1.544902e+03 ± 5.533339e+02
numeric issues : 0
-------------------------------------------------------------------

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

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