User question: jarandh

[jarandh]

I'm trying to add multi-dimensional noise to my model. The noise consists of an inflow and a demandscenario. Somehow, the SDDP.parameterize doesn't assign the noise to my variables through JuMP.fix because my \omega is of type Nothing.

Here is my code, it is based on the multi-dimensional noise term how-to-guide:

        Ω = [(inf = v, demscen = c) for v in [0, 10000, 20000] for c in [1, 2, 3, 4, 5]]
        P = [v * c for v in [1/3, 1/3, 1/3] for c in [1/5, 1/5, 1/5, 1/5, 1/5]]

        SDDP.parameterize(subproblem, Ω, P) do ω
            JuMP.fix(inflow, ω.inf)
            for i in 1:numloadblocks
                index = Int(1)
                JuMP.fix(demand[i], mean(demandvec[ω.demscen][index:Int(round((index + rounded_counts[i]-1)))]); force = true)
                index += Int(rounded_counts[i])
            end
        end

I end up with the Error: type Nothing has no field inf. When trying to evaluate or simulate the policy. I am able to train the model, but the lower bound is just 0, which suggests that demand is never assigned from the \Omega to the variable demand.

It feels like I've tried everything. Do anyone have any suggestions for what I can do, or is it something with SDDP.parameterize?

[odow]

Can you provide a reproducible example of the full that I can copy and paste?

It'd also be helpful to have the full text of the error.

The code as written looks correct, so something else must be happening.

[jarandh]

Yes, I will write a reproducible example here right away. The full error looks like this:

ERROR: type Nothing has no field inf
Stacktrace:
 [1] getproperty(x::Nothing, f::Symbol)
   @ Base .\Base.jl:38
 [2] (::var"#40#52"{Vector{VariableRef}, VariableRef})(ω::Nothing)
   @ Main c:\Users\jho\OneDrive - Norges vassdrags- og energidirektorat\Dokumenter\jaju\nowweareinvestingV6.jl:166
 [3] parameterize(node::SDDP.Node{Int64}, noise::Nothing)
   @ SDDP C:\Users\jho\.julia\packages\SDDP\2Jypx\src\algorithm.jl:255
 [4] solve_subproblem(model::SDDP.PolicyGraph{Int64}, node::SDDP.Node{Int64}, state::Dict{Symbol, Float64}, noise::Nothing, scenario_path::Vector{Tuple{Int64, Any}}; duality_handler::Nothing)
   @ SDDP C:\Users\jho\.julia\packages\SDDP\2Jypx\src\algorithm.jl:375
 [5] evaluate(rule::SDDP.DecisionRule{Int64}; incoming_state::Dict{Symbol, Float64}, noise::Nothing, controls_to_record::Vector{Symbol})
   @ SDDP C:\Users\jho\.julia\packages\SDDP\2Jypx\src\algorithm.jl:1327
 [6] top-level scope
   @ c:\Users\jho\OneDrive - Norges vassdrags- og energidirektorat\Dokumenter\jaju\nowweareinvestingV6.jl:187

[odow]

Ah. You're calling SDDP.evaluate? You need to pass a noise keyword argument.

SDDP.evaluate(
    rule;
    noise = (inf = 10000, demscen = 3),
    incoming_state = ... TODO ...
)

[jarandh]

Oh, sorry about that. That's not really the core of my problem though. I have a model where I don't need noise in my first node, but in the first node, the controls are investments in plant capacities, which are state variables in the model. I will show you in the example, but the problem is that lower bound is 0 throughout the training, which I think means that adding demand as noise in the code above does not work. The model worked fine when demand was just constant defined in the model constraint. I'm trying to refine a understandable example to paste here.

[jarandh]

Below is an example that should be ok to copy and paste:

using SDDP, HiGHS, Plots, Statistics

# Model setup
numstages = 52 #weeks
β = 0.95 # discount factor
i = 1 - β # interest rate
β_w = β^(1/52) # weekly discount factor
N = 20 # lifetime of investments

α = (i*(1+i)^N)/((1+i)^N-1)  # annuity factor

# Costs 
shortagecost = 100000 #$/MWh

c_peaker = 50000 #$/MW
c_peaker_infinite = c_peaker*α/(1-β)

c_baseload = 250000 #$/MW
c_baseload_infinite = c_baseload*α/(1-β)

fuelcost_peaker = 200 #$/MWh
fuelcost_baseload = 100 #$/MWh

# Demand scenarios
numdemandscenarios = 5
peakhours = 33
offpeakhours = 168- peakhours 
demandvec = Vector{Vector{Float64}}(undef, numdemandscenarios)

for i = 1:numdemandscenarios
    dp = rand(200:300, peakhours)
    db = rand(100:200, offpeakhours)
    d = sort(vcat(dp,db), rev=true)
    demandvec[i] = d
end

# Demand to loadblocks
numloadblocks = 5
rounded_counts = [19, 25, 41, 36, 47]

# Defining the graph. The first node is only visited once to make the investment in thermal capacity. Weekly discounting.

graph = SDDP.Graph(0)

for i = 1:numstages+1
    SDDP.add_node(graph, i)
end

SDDP.add_edge(graph, 0 => 1, 1.0)
SDDP.add_edge(graph, 1 => 2, 1.0)

for i = 2:(numstages)
    SDDP.add_edge(graph, i => i+1, β_w^(i-1))
end

SDDP.add_edge(graph, numstages+1 => 2, β)

# Building SDDP-model
model = SDDP.PolicyGraph(
    graph,
    sense = :Min,
    lower_bound = 0.0,
    optimizer = HiGHS.Optimizer,
) do subproblem, node

# State variables, the thermal capacity is following the problem in every stage
@variable(subproblem, 0 <= volume <= 30000, SDDP.State, initial_value = 20000) #Reservoir level in MWh   
@variable(subproblem, 0 <= peaker_capacity, SDDP.State, initial_value = 0) #MW
@variable(subproblem, 0 <= baseload_capacity, SDDP.State, initial_value = 0) #MW

# Special treatment for the investment stage: Just investment, no demand to meet and no noise to observe
    if node == 1
        # Control in the first node
        @variables(subproblem, begin
        peaker_investment >= 0
        baseload_investment >= 0
        end)

        # Constraint in the first node
        @constraints(subproblem, begin
        added_peaker_investment, peaker_capacity.out == peaker_capacity.in + peaker_investment
        added_baseload_investment, baseload_capacity.out == baseload_capacity.in + baseload_investment
        firststage_water_balance, volume.out == volume.in
        end)

        # Stage objective in the first node
        @stageobjective(subproblem, c_peaker_infinite*peaker_investment + c_baseload_infinite*baseload_investment)

# The rest of the nodes are implemented with a very simple power market model    
    else

        # Controls in in the rest of the nodes. First block is peak hours, and second block is off-peak hours
        # Arrays with generation in MWh, range from higehst load center to lowest load center
        @variables(subproblem, begin
        peaker_generation[1:numloadblocks] >= 0
        baseload_generation[1:numloadblocks] >= 0
        hydro_generation[1:numloadblocks] >= 0
        hydro_spill[1:numloadblocks] >= 0
        inflow
        shortage[1:numloadblocks] >= 0
        demand[1:numloadblocks] >=0
        end
        )

        # Constraints in the rest of the nodes 
        @constraints(subproblem, begin
        water_balance, volume.out == volume.in + inflow - sum(hydro_generation) - sum(hydro_spill)
        demand_constraint, peaker_generation + baseload_generation + hydro_generation + shortage .== demand
        peaker_cap_constraint, peaker_generation .<= peaker_capacity.in.*rounded_counts
        baseload_cap_constraint, baseload_generation .<= baseload_capacity.in.*rounded_counts
        no_more_peaker, peaker_capacity.out == peaker_capacity.in
        no_more_baseload, baseload_capacity.out == baseload_capacity.in
        end)

        # Noise
        Ω = [(inf = v, demscen = c) for v in [0, 10000, 20000] for c in [1, 2, 3, 4, 5]]
        P = [v * c for v in [1/3, 1/3, 1/3] for c in [1/5, 1/5, 1/5, 1/5, 1/5]]

        SDDP.parameterize(subproblem, Ω, P) do ω
            JuMP.fix(inflow, ω.inf)
            for i in 1:numloadblocks
                index = Int(1)
                JuMP.fix(demand[i], mean(demandvec[ω.demscen][index:Int(round((index + rounded_counts[i]-1)))]); force = true)
                index += Int(rounded_counts[i])
            end
        end

        # Stage objective in the rest of the nodes
        @stageobjective(
            subproblem,
            fuelcost_peaker*sum(peaker_generation) + fuelcost_baseload*sum(baseload_generation) + shortagecost*sum(shortage)) 
    end
end

SDDP.train(model; iteration_limit = 20)

[odow]

Overall, nice model. You're heading in the right direction and nearly there. But there's quite a few different problems that are combining to confuse you :smile:

Bug in parameterize

I'm guessing you don't want to restore index = 1 in every iteration of the for loop.

        SDDP.parameterize(subproblem, Ω, P) do ω
            JuMP.fix(inflow, ω.inf)
            for i in 1:numloadblocks
                index = Int(1)
                JuMP.fix(demand[i], mean(demandvec[ω.demscen][index:Int(round((index + rounded_counts[i]-1)))]); force = true)
                index += Int(rounded_counts[i])
            end
        end

should be

        SDDP.parameterize(subproblem, Ω, P) do ω
            JuMP.fix(inflow, ω.inf)
            index = 1
            for i in 1:numloadblocks
                JuMP.fix(demand[i], mean(demandvec[ω.demscen][index:(index+rounded_counts[i]-1)]); force = true)
                index += rounded_counts[i]
            end
        end

Discount factors

It's also cool to see that you're constructing a non-standard graph. This is where SDDP.jl shines compared to alternatives. But your graph also needs some tweaking.

If I set numstages = 4 to make it a bit easier to understand, you currently have:

julia> numstages = 4
4

julia> graph = SDDP.Graph(0)
Root
 0
Nodes
 {}
Arcs
 {}

julia> for i = 1:numstages+1
           SDDP.add_node(graph, i)
       end

julia> SDDP.add_edge(graph, 0 => 1, 1.0)

julia> SDDP.add_edge(graph, 1 => 2, 1.0)

julia> for i = 2:(numstages)
           SDDP.add_edge(graph, i => i+1, β_w^(i-1))
       end

julia> SDDP.add_edge(graph, numstages+1 => 2, β)

julia> graph
Root
 0
Nodes
 1
 2
 3
 4
 5
Arcs
 0 => 1 w.p. 1.0
 1 => 2 w.p. 1.0
 2 => 3 w.p. 0.9990140768344814
 3 => 4 w.p. 0.9980291257134511
 4 => 5 w.p. 0.997045145678548
 5 => 2 w.p. 0.95

This has a weekly discount factor from week 2-3, 3-4, and 4-5 that increases over time! Then you have a bit annual discount from weeks 5-2. You either need a constant weekly discount, or you just need a single annual discount in in the 5-2 arc. Currently you're mixing things together.

I'd do instead:

julia> numstages = 4
4

julia> graph = SDDP.LinearGraph(numstages+1);

julia> SDDP.add_edge(graph, numstages+1 => 2, β)

julia> graph
Root
 0
Nodes
 1
 2
 3
 4
 5
Arcs
 0 => 1 w.p. 1.0
 1 => 2 w.p. 1.0
 2 => 3 w.p. 1.0
 3 => 4 w.p. 1.0
 4 => 5 w.p. 1.0
 5 => 2 w.p. 0.95

Too large costs

You'll also likely run into some numerical issues because your costs are too large. SDDP.jl depends on solvers like HiGHS, and they don't like it when coefficients are much larger than 1e6 or so. So dealing in dollars/MW can lead, in the first stage, to very high costs. Cut things down to thousands of dollars per MW.

shortagecost = 100 #'000$/MWh

c_peaker = 50 #'000$/MW
c_peaker_infinite = c_peaker*α/(1-β)

c_baseload = 250 #'000$/MW
c_baseload_infinite = c_baseload*α/(1-β)

fuelcost_peaker = 0.200 #'000$/MWh
fuelcost_baseload = 0.100 #'000$/MWh

Related, there's nothing to stop SDDP.jl from exploring very large investments in the first-stage. Although not optimal, this can also lead to numerical issues. You'll likely get better performance by adding something like:

        @variables(subproblem, begin
            0 <= peaker_investment <= 1000
            0 <= baseload_investment <= 1000
        end)

Too high inflows

Finally, when I run the model, it seems that it can run everything via hydro without buying any peaker or baseload capacity.

Your inflows currently have a mean of 10_000 units, and your reservoir capacity is only 30,000, so you can fill 1/3 of the reservoir every week!

[jarandh]

Thanks so much Oscar! Very much appreciated. Yes, it's been really fun working with SDDP.jl so far!

From your last sentence, I realize what is wrong with this implementation compared to the previous models. The demand that is to be met in the constraint isn't multiplied with the number of hours. I must have been confused by the error combined with the lack of the noise-argument in SDDP.evaluate. But I'll fix all your suggestions, and hopefully it will all make sense.

Question about weekly discounting: don't you need to also discount the cyclic edge from the last stage?

[odow]

don't you need to also discount the cyclic edge from the last stage?

On a technical level, you need the product of the probabilities around every loop to be strictly less than one. It doesn't really matter if this happen on a single edge, or over lots of edges.

In this example, the edge from 5 to 2 is discounted:

julia> numstages = 4
4

julia> graph = SDDP.LinearGraph(numstages+1);

julia> SDDP.add_edge(graph, numstages+1 => 2, β)

julia> graph
Root
 0
Nodes
 1
 2
 3
 4
 5
Arcs
 0 => 1 w.p. 1.0
 1 => 2 w.p. 1.0
 2 => 3 w.p. 1.0
 3 => 4 w.p. 1.0
 4 => 5 w.p. 1.0
 5 => 2 w.p. 0.95

and the product around the 2-3-4-5-2 loop is 0.95.

[jarandh]

Oh, I see my mistake now. Easy to see that I'm not an economist. The weekly discount factors from a node to another node should off course be the same, and not refer to the first node. Anyways, thank you so much Oscar for helping out. I'm staying at EPOC with Andy for a year, so you'll probably hear from me again.

[odow]

Here's how I'd write your model:

using SDDP
import Gurobi
import Statistics
const GRB_ENV = Gurobi.Env()

β = 0.95

function create_demand_scenario(load_blocks)
    B = length(load_blocks)
    peak_hours = 33
    off_peak_hours = 168 - 33
    demand = [rand(200:300, peak_hours); rand(100:200, off_peak_hours)]
    sort!(demand; rev = true)
    block_demand = zeros(B)
    index = 1
    for (i, block) in enumerate(load_blocks)
        block_demand[i] = block * Statistics.mean(demand[index:(index+block-1)])
        index += block
    end
    return block_demand
end

graph = SDDP.LinearGraph(53)
SDDP.add_edge(graph, 53 => 2, β)

model = SDDP.PolicyGraph(
    graph,
    sense = :Min,
    lower_bound = 0.0,
    optimizer = () -> Gurobi.Optimizer(GRB_ENV),
) do sp, node
    α = (2 - β)^20 /((2 - β)^20 - 1)
    c_shortage = 100.0
    c_cap_peaker = 50 * α
    c_cap_baseload = 250 * α
    c_gen_peaker = 0.2
    c_gen_baseload = 0.1
    load_blocks = [19, 25, 41, 36, 47]
    B = length(load_blocks)
    @variables(sp, begin
        0 <= x_hydro <= 30_000, SDDP.State, (initial_value = 20_000)
        0 <= x_peaker <= 1_000, SDDP.State, (initial_value = 0)
        0 <= x_baseload <= 1_000, SDDP.State, (initial_value = 0)
        u_cap_peaker >= 0
        u_cap_baseload >= 0
        u_gen_peaker[1:B] >= 0
        u_gen_baseload[1:B] >= 0
        u_gen_hydro[1:B] >= 0
        u_shortage[1:B] >= 0
        w_inflow
        w_demand[1:B]
    end)
    @constraints(sp, begin
        x_peaker.out == x_peaker.in + u_cap_peaker
        x_baseload.out == x_baseload.in + u_cap_baseload
        x_hydro.out <= x_hydro.in + w_inflow - sum(u_gen_hydro)
        [i=1:B], u_gen_peaker[i] + u_gen_baseload[i] + u_gen_hydro[i] + u_shortage[i] == w_demand[i]
        [i=1:B], u_gen_peaker[i] <= load_blocks[i] * x_peaker.in
        [i=1:B], u_gen_baseload[i] <= load_blocks[i] * x_baseload.in
    end)
    @stageobjective(
        sp,
        c_cap_peaker * u_cap_peaker +
        c_cap_baseload * u_cap_baseload +
        c_gen_peaker * sum(u_gen_peaker) +
        c_gen_baseload * sum(u_gen_baseload) +
        c_shortage * sum(u_shortage),
    )
    if node == 1
        # Node 1 is special. All operational variables should go to zero, and we
        # should install baseload and peaker capacity.
        fix(w_inflow, 0)
        fix.(w_demand, 0)
    else
        # Cannot install capacity in operation nodes
        fix(u_cap_peaker, 0; force = true)
        fix(u_cap_baseload, 0; force = true)
        # Parameterize randomness
        Ω_inflow = [0, 10_000, 20_000]
        Ω_demand = [create_demand_scenario(load_blocks) for _ in 1:5]
        Ω = [(inflow = i, demand = d) for i in Ω_inflow for d in Ω_demand]
        SDDP.parameterize(sp, Ω) do ω
            JuMP.fix(w_inflow, ω.inflow)
            JuMP.fix.(w_demand, ω.demand)
            return
        end
    end
    return
end

SDDP.train(model; iteration_limit = 100)

Creating a single subproblem is a little weird, but it will simplify SDDP.simulate later on. I also try to put all data inside the do ... end block, which simplifies things if you want to use parallelism, and I find forces you to write cleaner models and data.

I'm staying at EPOC with Andy for a year

Let me know if you ever want a coffee/beer to discuss this. o.dowson@gmail.com.