Hello everyone,
I was experimenting with SDDP.jl, changing parameters to see how it behaved.
I defined a random variable but parametrized it equally (9 equal realizations), so the problem became deterministic and was solved with SDDP.jl.
Then I used JuMP to directly formulate the problem, which I then solved using Gurobi.
However, the optimal values for these two problems were surprisingly different (by far).
Moreover, the lower bound in SDDP.jl did not change after the first steps (which is reasonable since we do not have different trajectories in the policy graph), however, the gap between Bound (order of 10^7) and Simulation columns (order of 10^8) was significant.
Is there any explanation besides the possibility that I made a mistake while formulating them?
Is there any explanation besides the possibility that I made a mistake while formulating them?
I guess there could be a bug in SDDP.jl or Gurobi, but the far likelier explanation is that there is a bug in your code.
Can you share a reproducible example?
Yes. In the following, I produce small expansion + production problem, and solve once by Gurobi, and again by SDDP.jl. There is no uncertainty in the problem.
There are two hydropower plants, the first one has reservoir, second one is ROR.
Hydro 1 locates above hydro 2, so the water discharged (wd) and spilled (ww) of hydro 1 goes to hydro.
import Pkg
Pkg.add("JuMP")
Pkg.add("Gurobi")
Pkg.add("SDDP")
using Pkg
using JuMP
using Gurobi
using SDDP
Hydropowers = [1, 2]
DEMAND = 0.03.*[33541.59,33407.22, 31573.49,29763.47, 32625.51,34560.82, 33225.80, 31276.97, 27555.82, 28167.11, 24756.90]
INFLOW_A = [ 1.5*328.421,1.5*323.052,1.5*318.552,1.5*314.842,1.5*309.868,1.5*305.526,1.5*301.421,1.5*296.763,1.5*293.131,1.5*289.342,1.5*285.315]
# Inflow of the first hydro
INFLOW_B = [145.706, 143.390, 141.431,139.518,137.240,135.45,133.303,131.559,129.834,127.771,126.103]
# Inflow of the second hydro
INIT_CAP_HYDRO = [0.469*1000, 0.319*1000]
# Initial capacity of generation for hydro 1 and hydro 2
CONV_FAC_DAILY = (3600*24)/(100^3)
# m^3/s to hm^3
MAX_SPILL = [0,5000]
MAX_SPILL_OUT = [3340, 5000]
PROD_FAC2 = [0.415044, 0.407028]
# Productivity factor. m^3/s to MW
INIT_WATER = 49880.0
#Initial water in the reservoir
OMEG = Dict{Any, Any}()
for i in 1:11
OMEG[i] = [(demand = h) for h in [1*DEMAND[i], 1*DEMAND[i], 1*DEMAND[i]]]
end
P = [1/3, 1/3, 1/3]
graph = SDDP.LinearGraph(1 + 11)
function builder(model, node)
@variable(model, cap_hyd[i in Hydropowers], SDDP.State, lower_bound = INIT_CAP_HYDRO[i], initial_value = INIT_CAP_HYDRO[i])
@variable(model, ws, SDDP.State, lower_bound = 0, initial_value = INIT_WATER)
# Stored water in the reservoir at the beginning of day j, Dimension [hm^3]
if node == 1
#____________________________Expanion node_______________________________
@constraint(model, [i = Hydropowers], cap_hyd[i].out <= 1.2*INIT_CAP_HYDRO[i])
@constraint(model, ws.out == ws.in)
@stageobjective(model, sum((cap_hyd[i].out - INIT_CAP_HYDRO[i]) for i in Hydropowers)*1000)
else
#______________________Production node: hydro____________________________
@variable(model, gen_hyd[i in Hydropowers], lower_bound = 0)
# Hourly electricity generation by hydro i, Dimension [MW]
@variable(model, wd[i in Hydropowers], lower_bound = 0)
# Water discharge from hydropower plant i, Dimension [m^3/s]
@variable(model, wo[i in Hydropowers], lower_bound = 0)
# Out of system water from hydropower plant i, Dimension [m^3/s]
@variable(model, ww[i in Hydropowers], lower_bound = 0)
# Water spilled from hydropower plant i, Dimension [m^3/s]
@variable(model, pe, lower_bound = 0)
# Amount of electricity purchased. Dimesnion [MW]
@variable(model, se, lower_bound = 0)
# Amount of electricity sold. Dimension [MW]
@variable(model, penalty, lower_bound = 0)
@variable(model, demand)
SDDP.parameterize(model, OMEG[node-1], P) do ω
JuMP.fix(demand, ω)
end
@constraint(model, [i = Hydropowers], cap_hyd[i].out == cap_hyd[i].in)
@constraint(model, [i = Hydropowers], gen_hyd[i] <= cap_hyd[i].in)
@constraint(model, ws.out == ws.in + CONV_FAC_DAILY*(INFLOW_A[node-1] - ww[1] - wo[1] - wd[1]))
if node == 11
@constraint(model, 0.8*INIT_WATER[1] - ws.out == penalty)
end
@constraint(model, 0 == CONV_FAC_DAILY*(INFLOW_B[node-1] - ww[2] - wo[2] -wd[2] + ww[1] + wd[1]))
@constraint(model, [i = Hydropowers], ww[i] <= MAX_SPILL[i])
@constraint(model, [i = Hydropowers], wo[i] <= MAX_SPILL_OUT[i])
@constraint(model, [i = Hydropowers], gen_hyd[i] <= PROD_FAC2[i]*wd[i])
@constraint(model, pe <= 100000)
@constraint(model, se <= 2000)
@constraint(model, sum(gen_hyd[i] for i in Hydropowers) + pe - se == demand)
@stageobjective(model, sum(gen_hyd[i]*0.001*24 for i in Hydropowers) + 30*24*pe - 5*24*se + penalty*200000)
end
return model
end
main_model = SDDP.PolicyGraph(
builder,
graph;
sense = :Min,
lower_bound = 0,
optimizer = Gurobi.Optimizer,
)
SDDP.numerical_stability_report(main_model)
SDDP.train(main_model; iteration_limit = 100, cut_type = SDDP.SINGLE_CUT)
simulations = SDDP.simulate(main_model, 100, [:gen_wind, :cap_hyd, :wo,:ww,:ws,:wd,:se ,:pe,:demand, :gen_hyd], skip_undefined_variables=true)
for i in 2:12
println(i," " ,simulations[1][i][:se], " ",simulations[1][i][:pe], " ", simulations[1][i][:gen_hyd][1], " ", simulations[1][i][:gen_hyd][2])
end
Now, we formulate the problem directly using JuMP and solving by Gurobi.
T = 1:11
model = Model(Gurobi.Optimizer)
@variable(model, cap_hyd[i in Hydropowers , t in T], lower_bound = INIT_CAP_HYDRO[i])
@variable(model, ws[t in 1:12], lower_bound = 0)
# Level of the water at the beginning of time t
@variable(model, gen_hyd[i in Hydropowers, t in 1:11], lower_bound = 0)
@variable(model, wd[i in Hydropowers, t in 1:11], lower_bound = 0)
@variable(model, wo[i in Hydropowers, t in 1:11], lower_bound = 0)
@variable(model, ww[i in Hydropowers, t in 1:11], lower_bound = 0)
@variable(model, pe[t in 1:11], lower_bound = 0)
@variable(model, se[t in 1:11], lower_bound = 0)
@variable(model, penalty, lower_bound = 0)
@constraint(model, [i = Hydropowers], cap_hyd[i,1] <= 1.2*INIT_CAP_HYDRO[i])
@constraint(model, ws[1] == INIT_WATER)
for t in 1:11
if t in 1:10
@constraint(model, [i = Hydropowers], cap_hyd[i,t+1] == cap_hyd[i,t])
end
@constraint(model, [i = Hydropowers], gen_hyd[i,t] <= cap_hyd[i,t])
@constraint(model, ws[t+1] == ws[t] + CONV_FAC_DAILY*(INFLOW_A[t] - ww[1,t] - wo[1,t] - wd[1,t]))
@constraint(model, 0.8*INIT_WATER[1] - ws[12] == penalty)
@constraint(model, 0 == CONV_FAC_DAILY*(INFLOW_B[t] - ww[2,t] - wo[2,t] - wd[2,t] + ww[1,t] + wd[1,t]))
@constraint(model, [i = Hydropowers], ww[i,t] <= MAX_SPILL[i])
@constraint(model, [i = Hydropowers], wo[i,t] <= MAX_SPILL_OUT[i])
@constraint(model, [i = Hydropowers], gen_hyd[i,t] <= PROD_FAC2[i]*wd[i,t])
@constraint(model, pe[t] <= 100000)
@constraint(model, se[t] <= 2000)
@constraint(model, sum(gen_hyd[i,t] for i in Hydropowers) + pe[t] - se[t] == DEMAND[t])
end
@objective(model, Min, sum((cap_hyd[i,1] - INIT_CAP_HYDRO[i]) for i in Hydropowers)*1000 +
sum(gen_hyd[i,t]*0.001*24 for i in Hydropowers for t in 1:11) + sum(30*24*pe[t]- 5*24*se[t] for t in 1:11) + penalty*200000)
optimize!(model)
for t in 1:11
println(value(se[t]), " ", value(pe[t]) ," ", value(gen_hyd[1,t])," " ,value(gen_hyd[2,t]))
end
As you can see:
- The optimal value and solution obtained by SDDP.jl and JuMP is different (see the last four rows of the results I printed in both problems).
- The simulation column is not equal to the bound column while using SDDP.jl
One thing that might help us figure out what’s causing this error is that if we change the coefficient in the DEMAND parameter from 0.03 to 0.04 in both problems, we’ll have more demand in each node and won’t be able to sell any electricity. SDDP.jl and JuMP will have the same optimal solution as a result of this change. In addition, the simulation and bound columns in SDDP.jl will be the same.
It is very likely that I had a bug in my code or formulation that I was unable to identify.
1 Like
Im not at my computer right now, but Ill take a look when I get a chance.
So this was actually slightly tricky to track down. I ended up re-writing your code to identify the typo between the two models…and there wasn’t one. (Although you’ll see a few minor changes to the code, like using variable bounds instead of constraints. These are mathematically the same problem, but the variable bound version is “nicer” for the solvers.)
The problem turns out to be the lower bound that you set for SDDP. At first glance, a lower bound of 0 seemed reasonable, because the true cost is much higher than that, and most terms are non-negative.
The exception is the se term, which has a coefficient of -120. The lower bound needs to be a valid lower bound, not just for the full solution, but for the cost to go of every stage. Setting different bound fixes this.
I need to come up with a strategy for detecting when the bound choice is wrong. At the moment, the only real approach is to try different values and see if things change.
using JuMP
using Gurobi
using SDDP
Hydropowers = [1, 2]
DEMAND = 0.03 .* [33541.59,33407.22, 31573.49,29763.47, 32625.51,34560.82, 33225.80, 31276.97, 27555.82, 28167.11, 24756.90]
INFLOW_A = [1.5*328.421,1.5*323.052,1.5*318.552,1.5*314.842,1.5*309.868,1.5*305.526,1.5*301.421,1.5*296.763,1.5*293.131,1.5*289.342,1.5*285.315]
INFLOW_B = [145.706, 143.390, 141.431,139.518,137.240,135.45,133.303,131.559,129.834,127.771,126.103]
INIT_CAP_HYDRO = [0.469*1000, 0.319*1000]
CONV_FAC_DAILY = (3600*24) / (100^3)
MAX_SPILL = [0, 5000]
MAX_SPILL_OUT = [3340, 5000]
PROD_FAC2 = [0.415044, 0.407028]
INIT_WATER = 49880.0
T = 1:11
_string(s::Int) = lpad(s, 2)
_string(s) = lpad(round(s; digits = 2), 6)
begin
main_model = SDDP.LinearPolicyGraph(;
stages = 12,
sense = :Min,
lower_bound = -100_000,
optimizer = Gurobi.Optimizer,
) do model, node
t = node - 1
@variable(
model,
INIT_CAP_HYDRO[i] <= cap_hyd[i in Hydropowers] <= 1.2 * INIT_CAP_HYDRO[i],
SDDP.State,
initial_value = INIT_CAP_HYDRO[i],
)
@variable(model, ws >= 0, SDDP.State, initial_value = INIT_WATER)
if node == 1
@constraint(model, ws.out == ws.in)
@stageobjective(
model,
1_000 * sum(cap_hyd[i].out - INIT_CAP_HYDRO[i] for i in Hydropowers),
)
else
@variables(model, begin
gen_hyd[i in Hydropowers] >= 0
wd[i in Hydropowers] >= 0
0 <= wo[i in Hydropowers] <= MAX_SPILL_OUT[i]
0 <= ww[i in Hydropowers] <= MAX_SPILL[i]
0 <= pe <= 100_000
0 <= se <= 2_000
penalty >= 0
demand == DEMAND[t]
end)
@constraints(model, begin
[i in Hydropowers], cap_hyd[i].out == cap_hyd[i].in
[i in Hydropowers], gen_hyd[i] <= cap_hyd[i].in
[i in Hydropowers], gen_hyd[i] <= PROD_FAC2[i] * wd[i]
ws.out == ws.in + CONV_FAC_DAILY*(INFLOW_A[t] - ww[1] - wo[1] - wd[1])
0 == CONV_FAC_DAILY*(INFLOW_B[t] - ww[2] - wo[2] - wd[2] + ww[1] + wd[1])
sum(gen_hyd) + pe - se == demand
end)
if t == 11
@constraint(model, 0.8 * INIT_WATER - ws.out <= penalty)
end
@stageobjective(
model,
24 * (0.001sum(gen_hyd) + 30pe - 5se) + 200_000 * penalty
)
end
end
SDDP.train(main_model)
simulations = SDDP.simulate(
main_model,
1,
[:gen_wind, :cap_hyd, :wo,:ww,:ws,:wd,:se ,:pe,:demand, :gen_hyd];
skip_undefined_variables = true,
);
data = simulations[1];
for t in T
println(
_string(t),
" ",
_string(simulations[1][1+t][:cap_hyd][1].out),
" ",
_string(simulations[1][1+t][:cap_hyd][2].out),
" ",
_string(simulations[1][1+t][:ws].in),
" ",
_string(simulations[1][1+t][:se]),
" ",
_string(simulations[1][1+t][:pe]),
" ",
_string(simulations[1][1+t][:gen_hyd][1]),
" ",
_string(simulations[1][1+t][:gen_hyd][2]),
)
end
end
begin
model = Model(Gurobi.Optimizer)
@variables(model, begin
INIT_CAP_HYDRO[i] <= cap_hyd[i in Hydropowers, T] <= 1.2 * INIT_CAP_HYDRO[i]
ws[1:12] >= 0
gen_hyd[Hydropowers, T] >= 0
wd[Hydropowers, T] >= 0
0 <= wo[i in Hydropowers, T] <= MAX_SPILL_OUT[i]
0 <= ww[i in Hydropowers, T] <= MAX_SPILL[i]
0 <= pe[T] <= 100_000
0 <= se[T] <= 2_000
penalty >= 0
end)
fix(ws[1], INIT_WATER; force = true)
@constraints(model, begin
[i in Hydropowers, t in 2:11], cap_hyd[i,t] == cap_hyd[i,t-1]
[i in Hydropowers, t in T], gen_hyd[i,t] <= cap_hyd[i,t]
[i in Hydropowers, t in T], gen_hyd[i,t] <= PROD_FAC2[i] * wd[i,t]
[t in T], ws[t+1] == ws[t] + CONV_FAC_DAILY*(INFLOW_A[t] - ww[1,t] - wo[1,t] - wd[1,t])
[t in T], 0 == CONV_FAC_DAILY*(INFLOW_B[t] - ww[2,t] - wo[2,t] - wd[2,t] + ww[1,t] + wd[1,t])
[t in T], sum(gen_hyd[:, t]) + pe[t] - se[t] == DEMAND[t]
0.8 * INIT_WATER - ws[12] <= penalty
end)
@objective(
model,
Min,
1_000 * sum(cap_hyd[i,1] - INIT_CAP_HYDRO[i] for i in Hydropowers) +
24 * (0.001sum(gen_hyd) + 30sum(pe) - 5sum(se)) +
200_000 * penalty,
)
optimize!(model)
for t in 1:11
println(
_string(t),
" ",
_string(value(cap_hyd[1, t])),
" ",
_string(value(cap_hyd[2, t])),
" ",
_string(value(ws[t])),
" ",
_string(value(se[t])),
" ",
_string(value(pe[t])),
" ",
_string(value(gen_hyd[1,t])),
" ",
_string(value(gen_hyd[2,t])),
)
end
end
Great.
Thank you very much for reproducing the problem and solving the issue.
1 Like