Markov chain leads to memory saturation

[AngeBlanchard]

Hi,

I'm struggling with my code. I had a model with 365 nodes that I used for a power dispatch problem on an electricity system. I wanted to add a Markov variable accounting for a supply disruption in gas (ie state 1 when imports are available, state 2 when they are not).

Since I know my laptop is not a super calculator (32GB RAM), I decided to reduce the problem size by representing the system with only 52 temporal nodes, each one of them with 2 possible states for the Markov variable, ending up with 104 nodes in total.

Still, my machine seems unable to simulate more than 10 scenarios, while it could simulate more than 500 with the 365 nodes version of my problem. Here Gurobi sends me "Error 10001"...

Is it normal ? Is the only solution to find a computer with more RAM ?
Thanks a lot for your help,
Ange

PS: here is the error message I get

ERROR: Gurobi Error 10001:
Stacktrace:
[1] _check_ret
@ C:\Users\Ange.blanchard.julia\packages\Gurobi\FliRK\src\MOI_wrapper\MOI_wrapper.jl:393 [inlined]
[2] Gurobi.Env(; output_flag::Int64, memory_limit::Nothing, started::Bool)
@ Gurobi C:\Users\Ange.blanchard.julia\packages\Gurobi\FliRK\src\MOI_wrapper\MOI_wrapper.jl:103
[3] Env
@ C:\Users\Ange.blanchard.julia\packages\Gurobi\FliRK\src\MOI_wrapper\MOI_wrapper.jl:95 [inlined]
[4] Gurobi.Optimizer(env::Nothing; enable_interrupts::Bool)
@ Gurobi C:\Users\Ange.blanchard.julia\packages\Gurobi\FliRK\src\MOI_wrapper\MOI_wrapper.jl:324
[5] Optimizer (repeats 2 times)
@ C:\Users\Ange.blanchard.julia\packages\Gurobi\FliRK\src\MOI_wrapper\MOI_wrapper.jl:318 [inlined]
[6] _instantiate_and_check(optimizer_constructor::Any)
@ MathOptInterface C:\Users\Ange.blanchard.julia\packages\MathOptInterface\cl3eR\src\instantiate.jl:94
[7] _instantiate_and_check(optimizer_constructor::MathOptInterface.OptimizerWithAttributes)
@ MathOptInterface C:\Users\Ange.blanchard.julia\packages\MathOptInterface\cl3eR\src\instantiate.jl:117
[8] instantiate(optimizer_constructor::Any; with_bridge_type::Type{Float64})
@ MathOptInterface C:\Users\Ange.blanchard.julia\packages\MathOptInterface\cl3eR\src\instantiate.jl:149
[9] set_optimizer(model::Model, optimizer_constructor::Any; add_bridges::Bool)
@ JuMP C:\Users\Ange.blanchard.julia\packages\JuMP\yYfHy\src\optimizer_interface.jl:405
[10] set_optimizer
@ C:\Users\Ange.blanchard.julia\packages\JuMP\yYfHy\src\optimizer_interface.jl:398 [inlined]
[11] SDDP.ValueFunction(node::SDDP.Node{Tuple{Int64, Int64}})
@ SDDP C:\Users\Ange.blanchard.julia\packages\SDDP\OoTfR\src\visualization\value_functions.jl:113
[12] #ValueFunction#248
@ C:\Users\Ange.blanchard.julia\packages\SDDP\OoTfR\src\visualization\value_functions.jl:105 [inlined]
[13] top-level scope
@ c:\Users\Ange.blanchard\Nextcloud\H2 SDDP\GITHUB_H2_SDDP\CODE\TRAINING_markov.jl:87

[odow]

I had a model with 365 nodes

Keep in mind that this will create 365 Gurobi models. If each of these has thousands of variables and constraints then yes... it can use a lot of RAM.

But I wouldn't say it was normal. I've solved bigger problems on a smaller machine. How big is each subproblem? (SDDP.jl prints out some summary statistics are the start of each run.)

It also looks like the error is when you are trying to visualize the value function after a training? One approach could be to use write_cuts_to_file after training, then close Julia to free RAM, restart, read_cuts_from_file, and then look at the value function. Otherwise, just use SDDP.simulate. I don't find the ValueFunction that helpful.

[AngeBlanchard]

Hi Oscar,
Thank you for your insights. Well, this is the summary I got when I implement the Markov Chain:

Problem
  Nodes           : 103
  State variables : 3
  Scenarios       : 5.00000e+51
  Existing cuts   : false
  Subproblem structure                      : (min, max)
    Variables                               : (8246, 8246)
    VariableRef in MOI.LessThan{Float64}    : (5715, 5716)
    AffExpr in MOI.LessThan{Float64}        : (336, 336)
    VariableRef in MOI.GreaterThan{Float64} : (7739, 7739)
    AffExpr in MOI.EqualTo{Float64}         : (2692, 3367)
Options
  Solver          : serial mode
  Risk measure    : SDDP.Expectation()
  Sampling scheme : SDDP.InSampleMonteCarlo

Actually I use ValueFunction for extracting the "water values" of my problem, and those are very valuable to me. I use this function inside SDDP.simulate. I think the problem is linked to that simulate process indeed, since the training phase is working well, but the simulate phase saturates the RAM after more than 10 runs...

I just ran my code without calling ValueFunction and it worked perfectly.
So I'm wondering now, how to retrieve the water values and filling levels without calling ValueFunction ? I tried to save the cuts as you told me, using write_cuts_to_file but the issue remains. The problem really is due to an over use of RAM inside the SDDP.simulate, and is not linked with the training phase at all, apparently.

Thanks a lot for your help
Ange

PS: here is the block of code that creates the error 10001, inside the SDDP.simulate

op_cost_H2 = Array{Float64}(undef, N_red, Nbsimus)

  for t in 1:N_red # looping through the stages
      for sim in 1:Nbsimus # looping through the simulations
          
          markov_2 = Int(sims[sim][t][:markov])
          V = SDDP.ValueFunction(model, node = (t, markov_2)) # node identification
  
          volume_PHS = Dict("volume_PHS" => sims[sim][t][:volume_PHS].out)
          volume_H2 = Dict("volume_H2" => sims[sim][t][:volume_H2].out)
          volume_BAT = Dict("volume_BAT" => sims[sim][t][:volume_BAT].out) # the 3 reservoirs levels
  
          levels = merge(merge(merge(volume_PHS), volume_H2), volume_BAT)
  
          c, p = SDDP.evaluate(V, levels)
  
          op_cost_H2[t, sim] = -p[Symbol("volume_H2")] # water values
      end
  end
  
  op_cost_H2_df = DataFrame(op_cost_H2, :auto)
  

[odow]

 for t in 1:N_red # looping through the stages
      for sim in 1:Nbsimus # looping through the simulations
          
          markov_2 = Int(sims[sim][t][:markov])
          V = SDDP.ValueFunction(model, node = (t, markov_2)) # node identification

This is creating a new Gurobi model for every t and every sim. It doesn't re-use the existing models that were trained, so you're creating a loooooooot of models!!!

One quick fix would be something like

value_functions = Dict(
    node => SDDP.ValueFunction(model; node = node) for node in keys(model.nodes)
)
op_cost_H2 = Array{Float64}(undef, N_red, Nbsimus)
for t in 1:N_red # looping through the stages
    for sim in 1:Nbsimus # looping through the simulations
        markov_2 = Int(sims[sim][t][:markov])
        V = value_functions[(t, markov_2)] # node identification
        levels = Dict(
            "volume_PHS" => sims[sim][t][:volume_PHS].out,
            "volume_H2" => sims[sim][t][:volume_H2].out,
            "volume_BAT" => sims[sim][t][:volume_BAT].out,
        )
        c, p = SDDP.evaluate(V, levels)
        op_cost_H2[t, sim] = -p[Symbol("volume_H2")] # water values
    end
end
op_cost_H2_df = DataFrame(op_cost_H2, :auto)

[odow]

See also https://odow.github.io/SDDP.jl/stable/tutorial/first_steps/#Custom-recorders.

[AngeBlanchard]

Thank you, I'll try it as soon as I can.
Something that puzzles me a bit is that my code seems to work just fine with CPLEX, while it leads to this 10001 error with Gurobi... Maybe the two solvers handle this differently.
But I hope it works with Gurobi anyway, since I don't have easy access to a laptop with CPLEX installed.
Thanks again,
Ange

[odow]

my code seems to work just fine with CPLEX

Even if it works, it'll be much slower than necessary

Maybe the two solvers handle this differently.

Yep. There are quite a few differences behind the scenes.

[odow]

This should point you in the right direction:

using SDDP
import DataFrames
import HiGHS

model = SDDP.LinearPolicyGraph(
    stages = 3,
    sense = :Min,
    lower_bound = 0.0,
    optimizer = HiGHS.Optimizer,
) do sp, t
    @variable(sp, 0 <= x <= 200, SDDP.State, initial_value = 200)
    @variable(sp, g_t >= 0)
    @variable(sp, g_h >= 0)
    @variable(sp, s >= 0)
    SDDP.parameterize(ω -> set_normalized_rhs(balance, ω), sp, [0, 50, 100])
    @constraint(sp, balance, x.out == x.in - g_h - s + 0.0)
    @constraint(sp, demand_constraint, g_h + g_t == 150)
    @stageobjective(sp, 50 * t * g_t)
end
SDDP.train(model; iteration_limit = 10)

simulations = SDDP.simulate(
    model,
    1,
    [:x];
    custom_recorders = Dict{Symbol,Function}(
        :x_price_out => sp -> reduced_cost(sp[:x].out),
    ),
)
results_simulate = DataFrames.DataFrame(vec([
    (
        simulation = i, 
        node = data[:node_index], 
        x = data[:x].out,
        x_price = data[:x_price_out],
    )
    for (i, sim) in enumerate(simulations) for data in sim
]))

V = Dict(node => SDDP.ValueFunction(model; node = node) for node in keys(model.nodes))
results = Any[]
for (i, sim) in enumerate(simulations), data in sim
    cost, price = SDDP.evaluate(V[data[:node_index]], Dict("x" => data[:x].out))
    push!(
        results, 
        (
            simulation = i, 
            node = data[:node_index], 
            x = data[:x].out,
            x_price = price[:x],
        ),
    )
end
results_value_function = DataFrames.DataFrame(results)

Which method to use will depend on what price you want to compute. The two methods differ by the cheapest marginal cost of thermal ($50) because the simulate method includes the first-stage objective, and the ValueFunction uses only the first-stage value function.

julia> results_simulate
3×4 DataFrame
 Row │ simulation  node   x        x_price  
     │ Int64       Int64  Float64  Float64  
─────┼──────────────────────────────────────
   1 │          1      1    200.0  -33.3333
   2 │          1      2    150.0    0.0
   3 │          1      3     -0.0    0.0

julia> results_value_function
3×4 DataFrame
 Row │ simulation  node   x        x_price  
     │ Int64       Int64  Float64  Float64  
─────┼──────────────────────────────────────
   1 │          1      1    200.0  -83.3333
   2 │          1      2    150.0    0.0
   3 │          1      3     -0.0    0.0

[odow]

Is there anything actionable to do here? If not, I will close.

[AngeBlanchard]

Yes, sorry, it's all good to me.
Thanks again !
Ange

[odow]

No problem. Let me know if you have other questions