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, Random
struct Turbine
flowknots::Vector{Float64}
powerknots::Vector{Float64}
end
struct Reservoir
min::Float64
max::Float64
initial::Float64
turbine::Turbine
spill_cost::Float64
inflows::Vector{Float64}
end
function 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
JuMP.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
end
end
function 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
### = $835
end
test_hydro_valley_model()Test Passed