# Historical data - SDDP.jl

**URL:** <https://discourse.julialang.org/t/historical-data-sddp-jl/117207>\
**Category:** Optimization (Mathematical)\
**Tags:** sddp\
**Created:** [July 19, 2024, 1:03am UTC](https://discourse.julialang.org/t/historical-data-sddp-jl/117207 "2024-07-19T01:03:42Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sorooshsa](https://avatars.discourse-cdn.com/v4/letter/s/977dab/32.png) [@Sorooshsa](https://discourse.julialang.org/u/Sorooshsa)\
**Post date:** [July 19, 2024, 1:03am UTC](https://discourse.julialang.org/t/historical-data-sddp-jl/117207/1 "2024-07-19T01:03:42Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [July 19, 2024, 1:15am UTC](https://discourse.julialang.org/t/historical-data-sddp-jl/117207/2 "2024-07-19T01:15:09Z")

</div>

Hard to say what is happening without a reproducible example 😄

How close is “close together”?

---

<div class="post-metadata">

**Author:** ![Sorooshsa](https://avatars.discourse-cdn.com/v4/letter/s/977dab/32.png) [@Sorooshsa](https://discourse.julialang.org/u/Sorooshsa)\
**Post date:** [July 19, 2024, 1:40am UTC](https://discourse.julialang.org/t/historical-data-sddp-jl/117207/3 "2024-07-19T01:40:18Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [July 19, 2024, 2:47am UTC](https://discourse.julialang.org/t/historical-data-sddp-jl/117207/4 "2024-07-19T02:47:52Z")

</div>

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:

```julia
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

```julia
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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Sorooshsa](https://avatars.discourse-cdn.com/v4/letter/s/977dab/32.png) [@Sorooshsa](https://discourse.julialang.org/u/Sorooshsa)\
**Post date:** [July 19, 2024, 4:00am UTC](https://discourse.julialang.org/t/historical-data-sddp-jl/117207/5 "2024-07-19T04:00:32Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Sorooshsa](https://avatars.discourse-cdn.com/v4/letter/s/977dab/32.png) [@Sorooshsa](https://discourse.julialang.org/u/Sorooshsa)\
**Post date:** [August 6, 2024, 10:53pm UTC](https://discourse.julialang.org/t/historical-data-sddp-jl/117207/6 "2024-08-06T22:53:07Z")

</div>

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.

```julia
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])

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 6, 2024, 11:42pm UTC](https://discourse.julialang.org/t/historical-data-sddp-jl/117207/7 "2024-08-06T23:42:24Z")

</div>

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](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:

```julia
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:

```julia
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)]

```