Hydro valleys
This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.
This problem is a version of the hydro-thermal scheduling problem. The goal is to operate two hydro-dams in a valley chain over time in the face of inflow and price uncertainty.
Turbine response curves are modelled by piecewise linear functions which map the flow rate into a power. These can be controlled by specifying the breakpoints in the piecewise linear function as the knots in the Turbine struct.
The model can be created using the hydro_valley_model function. It has a few keyword arguments to allow automated testing of the library. hasstagewiseinflows determines if the RHS noise constraint should be added. hasmarkovprice determines if the price uncertainty (modelled by a Markov chain) should be added.
In the third stage, the Markov chain has some unreachable states to test some code-paths in the library.
We can also set the sense to :Min or :Max (the objective and bound are flipped appropriately).
using SDDP, HiGHS, Test, Randomstruct Turbine flowknots::Vector{Float64} powerknots::Vector{Float64}endstruct Reservoir min::Float64 max::Float64 initial::Float64 turbine::Turbine spill_cost::Float64 inflows::Vector{Float64}endfunction hydro_valley_model(; hasstagewiseinflows::Bool = true, hasmarkovprice::Bool = true, sense::Symbol = :Max,) valley_chain = [ Reservoir( 0, 200, 200, Turbine([50, 60, 70], [55, 65, 70]), 1000, [0, 20, 50], ), Reservoir( 0, 200, 200, Turbine([50, 60, 70], [55, 65, 70]), 1000, [0, 0, 20], ), ] turbine(i) = valley_chain[i].turbine # Prices[t, Markov state] prices = [ 1 2 0 2 1 0 3 4 0 ] # Transition matrix if hasmarkovprice transition = Array{Float64,2}[[1.0]', [0.6 0.4], [0.6 0.4 0.0; 0.3 0.7 0.0]] else transition = [ones(Float64, (1, 1)) for t in 1:3] end flipobj = (sense == :Max) ? 1.0 : -1.0 lower = (sense == :Max) ? -Inf : -1e6 upper = (sense == :Max) ? 1e6 : Inf N = length(valley_chain) # Initialise SDDP Model return m = SDDP.MarkovianPolicyGraph(; sense = sense, lower_bound = lower, upper_bound = upper, transition_matrices = transition, optimizer = HiGHS.Optimizer, ) do subproblem, node t, markov_state = node # ------------------------------------------------------------------ # SDDP State Variables # Level of upper reservoir @variable( subproblem, valley_chain[r].min <= reservoir[r = 1:N] <= valley_chain[r].max, SDDP.State, initial_value = valley_chain[r].initial ) # ------------------------------------------------------------------ # Additional variables @variables( subproblem, begin outflow[r = 1:N] >= 0 spill[r = 1:N] >= 0 inflow[r = 1:N] >= 0 generation_quantity >= 0 # Total quantity of water # Proportion of levels to dispatch on 0 <= dispatch[r = 1:N, level = 1:length(turbine(r).flowknots)] <= 1 rainfall[i = 1:N] end ) # ------------------------------------------------------------------ # Constraints @constraints( subproblem, begin # flow from upper reservoir reservoir[1].out == reservoir[1].in + inflow[1] - outflow[1] - spill[1] # other flows flow[i = 2:N], reservoir[i].out == reservoir[i].in + inflow[i] - outflow[i] - spill[i] + outflow[i-1] + spill[i-1] # Total quantity generated generation_quantity == sum( turbine(r).powerknots[level] * dispatch[r, level] for r in 1:N for level in 1:length(turbine(r).powerknots) ) # ------------------------------------------------------------------ # Flow out turbineflow[r = 1:N], outflow[r] == sum( turbine(r).flowknots[level] * dispatch[r, level] for level in 1:length(turbine(r).flowknots) ) # Dispatch combination of levels dispatched[r = 1:N], sum( dispatch[r, level] for level in 1:length(turbine(r).flowknots) ) <= 1 end ) # rainfall noises if hasstagewiseinflows && t > 1 # in future stages random inflows @constraint(subproblem, inflow_noise[i = 1:N], inflow[i] <= rainfall[i]) SDDP.parameterize( subproblem, [ (valley_chain[1].inflows[i], valley_chain[2].inflows[i]) for i in 1:length(transition) ], ) do ω for i in 1:N fix(rainfall[i], ω[i]) end end else # in the first stage deterministic inflow @constraint( subproblem, initial_inflow_noise[i = 1:N], inflow[i] <= valley_chain[i].inflows[1] ) end # ------------------------------------------------------------------ # Objective Function if hasmarkovprice @stageobjective( subproblem, flipobj * ( prices[t, markov_state] * generation_quantity - sum(valley_chain[i].spill_cost * spill[i] for i in 1:N) ) ) else @stageobjective( subproblem, flipobj * ( prices[t, 1] * generation_quantity - sum(valley_chain[i].spill_cost * spill[i] for i in 1:N) ) ) end endendfunction test_hydro_valley_model() # For repeatability Random.seed!(11111) # deterministic deterministic_model = hydro_valley_model(; hasmarkovprice = false, hasstagewiseinflows = false, ) SDDP.train( deterministic_model; iteration_limit = 10, cut_deletion_minimum = 1, print_level = 0, ) @test SDDP.calculate_bound(deterministic_model) ≈ 835.0 atol = 1e-3 # stagewise inflows stagewise_model = hydro_valley_model(; hasmarkovprice = false) SDDP.train(stagewise_model; iteration_limit = 20, print_level = 0) @test SDDP.calculate_bound(stagewise_model) ≈ 838.33 atol = 1e-2 # Markov prices markov_model = hydro_valley_model(; hasstagewiseinflows = false) SDDP.train(markov_model; iteration_limit = 10, print_level = 0) @test SDDP.calculate_bound(markov_model) ≈ 851.8 atol = 1e-2 # stagewise inflows and Markov prices markov_stagewise_model = hydro_valley_model(; hasstagewiseinflows = true, hasmarkovprice = true) SDDP.train(markov_stagewise_model; iteration_limit = 10, print_level = 0) @test SDDP.calculate_bound(markov_stagewise_model) ≈ 855.0 atol = 1.0 # risk averse stagewise inflows and Markov prices riskaverse_model = hydro_valley_model() SDDP.train( riskaverse_model; risk_measure = SDDP.EAVaR(; lambda = 0.5, beta = 0.66), iteration_limit = 10, print_level = 0, ) @test SDDP.calculate_bound(riskaverse_model) ≈ 828.157 atol = 1.0 # stagewise inflows and Markov prices worst_case_model = hydro_valley_model(; sense = :Min) SDDP.train( worst_case_model; risk_measure = SDDP.EAVaR(; lambda = 0.5, beta = 0.0), iteration_limit = 10, print_level = 0, ) @test SDDP.calculate_bound(worst_case_model) ≈ -780.867 atol = 1.0 # stagewise inflows and Markov prices cutselection_model = hydro_valley_model() SDDP.train( cutselection_model; iteration_limit = 10, print_level = 0, cut_deletion_minimum = 2, ) @test SDDP.calculate_bound(cutselection_model) ≈ 855.0 atol = 1.0 # Distributionally robust Optimization dro_model = hydro_valley_model(; hasmarkovprice = false) SDDP.train( dro_model; risk_measure = SDDP.ModifiedChiSquared(sqrt(2 / 3) - 1e-6), iteration_limit = 10, print_level = 0, ) @test SDDP.calculate_bound(dro_model) ≈ 835.0 atol = 1.0 dro_model = hydro_valley_model(; hasmarkovprice = false) SDDP.train( dro_model; risk_measure = SDDP.ModifiedChiSquared(1 / 6), iteration_limit = 20, print_level = 0, ) @test SDDP.calculate_bound(dro_model) ≈ 836.695 atol = 1.0 # (Note) radius ≈ sqrt(2/3), will set all noise probabilities to zero except the worst case noise # (Why?): # The distance from the uniform distribution (the assumed "true" distribution) # to a corner of a unit simplex is sqrt(S-1)/sqrt(S) if we have S scenarios. The corner # of a unit simplex is just a unit vector, i.e.: [0 ... 0 1 0 ... 0]. With this probability # vector, only one noise has a non-zero probablity. # In the worst case rhsnoise (0 inflows) the profit is: # Reservoir1: 70 * $3 + 70 * $2 + 65 * $1 + # Reservoir2: 70 * $3 + 70 * $2 + 70 * $1 ### = $835endtest_hydro_valley_model()Test Passed