Historical data - SDDP.jl

Specific Domains Optimization (Mathematical)
1

Hi everyone,

I trained a big model using SDDP.jl. Then, I used historical data to simulate the policy multiple times. However, when I examine the results of different simulations with the same historical data, I notice that the optimal solutions and optimal values are different (though close together).

When I check the random noises, they are the same across different runs of simulations, as expected. Shouldn’t the results be the same for the same historical data in different runs? Or am I missing something?

2

Hard to say what is happening without a reproducible example :smile:

How close is “close together”?

1 Like
3

For instance, the optimal values for different runs vary by less than 4%. I will provide a reproducible example shortly. So, the short answer is that they should be the same, and I might be doing something wrong?

4

Are you just simulating the same policy multiple times and getting different results? This would be unexpected and we should investigate.

Or are you training multiple policies (with different random seeds) and then evaluating their historical scenarios? This could be expected, particularly if the historical scenarios are not drawn from the in-sample scenarios used during training.

In pseudo-code:

model = create_model()
SDDP.train(model)
results_1 = SDDP.simulate(model; sampling_scheme = historical)
results_2 = SDDP.simulate(model; sampling_scheme = historical)
compare(results_1, results_2)  # 4%

or

model = create_model()
SDDP.train(model)
results_1 = SDDP.simulate(model; sampling_scheme = historical)

model = create_model()
SDDP.train(model)
results_2 = SDDP.simulate(model; sampling_scheme = historical)

compare(results_1, results_2)  # 4%

If it is the latter, you can make the process deterministic by setting the random seed:

import Random
model = create_model()
Random.seed!(1234)
SDDP.train(model)
results_1 = SDDP.simulate(model; sampling_scheme = historical)

model = create_model()
Random.seed!(1234)
SDDP.train(model)
results_2 = SDDP.simulate(model; sampling_scheme = historical)

compare(results_1, results_2)  # exactly the same
1 Like
5

I exactly did the first one.
I will get back with a reproducible example.
Thank you.

6

Hi,
I have addressed the issue again with a reproducible code.
I should clarify that when simulating the same policy multiple times with the same historical data (Hist_1), the results will be same. However, once we simulate the model with another set of historical data (Hist_2) and then return to Hist_1 for simulation, the results will differ.

To explain, consider the code below. I first trained a model and then simulated it with Hist_1. Next, I simulated it with Hist_2, and finally, I simulated it again with Hist_1. I expected that the results of the first and third simulations would be the same, but they are not.

As we can see, both the optimal solutions and values are different, when we simulate with Hist_1 the first and second time.

import Pkg
using Pkg
Pkg.add("SDDP")
Pkg.add("Gurobi")
using SDDP
using Gurobi
using Random

Hydropower_plant = ["H1", "H2", "H3"]
Demand = Dict(vcat(collect(1:20) .=> 4000, collect(21:30) .=> 5000, collect(31:150) .=> 3000)...)
Product_factor = Dict("H1" => 0.5, "H2" => 0.5, "H3" => 0.5)
Probability = [1/4, 1/4, 1/4, 1/4]
Omega = [1000.0, 750.0, 500.0, 250]

graph = SDDP.LinearGraph(150)
function builder(model, node)
    @variables(model, begin
        0 <= gen_hyd[i in Hydropower_plant] 
        water_inf[i in Hydropower_plant]
        0 <= shedding <= 10000 # Big M for feasibility
        0 <= water_discharge[i in Hydropower_plant]
        0 <= purchase <= 1000 
    end)
    @variable(model, 1000 <= water_level[i in Hydropower_plant] <= 4000, SDDP.State, initial_value = 1000)

    
    @constraints(model, begin
        [i in Hydropower_plant], gen_hyd[i] <= 4000
        [i in Hydropower_plant], gen_hyd[i] <= Product_factor[i]*water_discharge[i]
    end)

    @constraint(model, sum(gen_hyd[i] for i in Hydropower_plant) + shedding + purchase  == Demand[node]) 

    @constraints(model, begin
        water_level["H1"].out ==  water_level["H1"].in + water_inf["H1"] - water_discharge["H1"] 
        water_level["H2"].out ==  water_level["H2"].in + water_inf["H2"] + water_discharge["H1"] - water_discharge["H2"]
        water_level["H3"].out ==  water_level["H3"].in + water_inf["H3"] + water_discharge["H2"] - water_discharge["H3"]
    end)

    SDDP.parameterize(model, Omega, Probability) do ω
        JuMP.fix(water_inf["H1"], ω)
        JuMP.fix(water_inf["H2"], ω)
        JuMP.fix(water_inf["H3"], ω)
    end

    @stageobjective(model, sum(gen_hyd[i]*20 for i in Hydropower_plant) + 25*purchase + 50*shedding)
end

main_model = SDDP.PolicyGraph(
builder,
graph;
sense = :Min,
lower_bound = 0,
optimizer = Gurobi.Optimizer,
)
SDDP.train(main_model; stopping_rules = [ SDDP.TimeLimit(30)], cut_type = SDDP.MULTI_CUT)
#__________________

Hist_1 = Vector{Tuple{Int, Int}}()
j = 1
for t in 1:150
    Random.seed!(t)
    push!(Hist_1, (j, rand(1:2000)))
    j += 1
end

Sim_1 = SDDP.simulate(main_model, 1, [:gen_hyd, :shedding, :purchase], sampling_scheme = SDDP.Historical(Hist_1))
Sol_1 = [sum(Sim_1[1][i][:gen_hyd][j] for j in Hydropower_plant) for i in 1:150]
Obj_1 = sum(Sim_1[1][i][:stage_objective] for i in 1:150)


Hist_2 = Vector{Tuple{Int, Int}}()
j = 1
for t in 1:150
    Random.seed!(t + 10)
    push!(Hist_2, (j, rand(1:2000)))
    j += 1
end
Sim_2 = SDDP.simulate(main_model, 1, [:gen_hyd, :shedding, :purchase], sampling_scheme = SDDP.Historical(Hist_2))


Sim_3 = SDDP.simulate(main_model, 1, [:gen_hyd, :shedding, :purchase], sampling_scheme = SDDP.Historical(Hist_1))
Sol_3 = [sum(Sim_3[1][i][:gen_hyd][j] for j in Hydropower_plant) for i in 1:150]
Obj_3 = sum(Sim_3[1][i][:stage_objective] for i in 1:150)

#___________Solutions comparison
println(Obj_3 - Obj_1)
println(Sol_3 - Sol_1)
println([Sim_3[1][i][:purchase] for i in 1:150] - [Sim_1[1][i][:purchase] for i in 1:150])
7

Aha! Dual degeneracy strikes again. @joaquimg will be interested to see this.

Take a read of section 2.7 from https://researchspace.auckland.ac.nz/bitstream/handle/2292/37700/whole.pdf?sequence=2&isAllowed=y

The issue is that the Sim_2 solve results in the Sim_3 models being warm-started from a different point than the Sim_1 models, and Gurobi finds a different solution that has the same objective value. From there we depart, and the two out-of-sample simulations have different costs.

This is pretty sensitive to the time limit, etc. It depends upon one of the stages having multiple optimal solutions in final policy. I was able to reproduce this a few times, and then not able with different seeds.

If you find a case that differs, you’ll probably find that Sim_1[1][1] == Sim_3[1][1] but then some other element is not equal.

You should be able to fix this by calling:

for node in values(model.nodes)
    MOI.Utilities.reset_optimizer(node.subproblem)
end

before each simulation, which resets the model into a clean state and doesn’t warm-start from previous solutions.

As a small aside, here is how I would write your model:

using SDDP
import Gurobi
import Random
c_demand = Dict((1:20 .=> 4000)..., (21:30 .=> 5000)..., (31:150 .=> 3000)...)
c_factor = Dict("H1" => 0.5, "H2" => 0.5, "H3" => 0.5)
s_hydros = collect(keys(c_factor))
c_upstream = Dict("H1" => [], "H2" => ["H1"], "H3" => ["H2"])
function builder(model, node)
    @variables(model, begin
        1_000 <= x_water_level[s_hydros] <= 4_000, SDDP.State, (initial_value = 1_000)
        0 <= u_gen_hyd[s_hydros] <= 4_000
        0 <= u_shedding
        0 <= u_water_discharge[s_hydros]
        0 <= u_purchase <= 1_000 
        w_water_inf[s_hydros]
    end)
    @constraints(model, begin
        [i in s_hydros], u_gen_hyd[i] <= c_factor[i] * u_water_discharge[i]
        [i in s_hydros], x_water_level[i].out == x_water_level[i].in + w_water_inf[i] - u_water_discharge[i] + sum(u_water_discharge[j] for j in c_upstream[i])
        sum(u_gen_hyd[i] for i in s_hydros) + u_shedding + u_purchase == c_demand[node]
    end)
    Ω = [1000.0, 750.0, 500.0, 250]
    SDDP.parameterize(model, Ω) do ω
        JuMP.fix.(w_water_inf, ω)
        return
    end
    @stageobjective(
        model,
        20 * sum(u_gen_hyd) + 25 * u_purchase + 50 * u_shedding,
    )
    return
end
model = SDDP.LinearPolicyGraph(
    builder;
    stages = 150,
    sense = :Min,
    lower_bound = 0.0,
    optimizer = Gurobi.Optimizer,
)
SDDP.train(model; time_limit = 30, cut_type = SDDP.MULTI_CUT)
Random.seed!(1234)
hist_1 = [(t, rand(1:2_000)) for t in 1:150]
sim_1 = SDDP.simulate(
    model, 
    1, 
    [:u_gen_hyd, :u_shedding, :u_purchase, :u_water_discharge, :w_water_inf]; 
    sampling_scheme = SDDP.Historical(hist_1),
)
obj_1 = sum(data[:stage_objective] for data in first(sim_1))
sol_1 = [sum(data[:u_gen_hyd]) for data in first(sim_1)]
1 Like