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