Historical data - SDDP.jl

Specific Domains Optimization (Mathematical)
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