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