-
Notifications
You must be signed in to change notification settings - Fork 61
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Numerical stability issues #646
Comments
Can you provide a reproducible example (or summary) in text format as a comment here? I don't want to download the random link, and it's also good practice for future readers (the link might break). |
##unscaled problem that gives numerical instability warning
|
scaled model that crashes
|
The second problem is infeasible because there is a mis-match between the units of load and Here are a few other suggestions:
This code doesn't run because of the unit issue, but it should point you in the right direction: using SDDP
using GLPK
using DataFrames
struct Reservoir
min_s::Float64
max_s::Float64
storage_initial::Float64
arc_cap::Vector{Float64}
max_pwr::Vector{Float64}
spill_cap::Float64
release_min::Float64
release_initial::Float64
k::Float64
n_unit::Int64
max_decrease::Float64
max_increase::Float64
min_gen::Float64
evap_avg_monthly::Float64
end
df = DataFrame(
Node = [1, 2, 3, 4, 5, 6],
Parent = [0, 1, 1, 2, 2, 3],
Prob = [1.0, 0.68, 0.32, 0.65, 0.35, 1.0],
Stage = [1, 2, 2, 3, 3, 3],
Inflow1 = [6, 7, 11, 8, 7, 11],
Inflow2 = [18, 17, 28, 16, 22, 28],
Inflow3 = [18, 18, 22, 19, 18, 24],
Inflow4 = [0.58, 0.0, 0.0, 0.0, 0.0, 0.6],
Inflow5 = [0.08, 1.55, 0.93, 0.0, 1.23, 6.54],
Inflow6 = [3.19, 4.1, 1.39, 3.6, 2.64, 2.25],
Load = [1.87112, 1.60109, 1.67808, 1.6246, 1.64542, 1.69742],
probability = [1.0, 0.68, 0.32, 0.442, 0.238, 0.32],
)
graph = SDDP.Graph(0)
for i in 1:length(df.Node)
SDDP.add_node(graph, i)
SDDP.add_edge(graph, df.Parent[i] => df.Node[i], df.Prob[i])
end
input = DataFrame(
Row = ["initial_storage", "min_daily_ene_gen"],
FortPeck = [1.5208e7, 57.9],
Garrison = [1.8195e7, 183.8],
Oahe = [1.90485e7, 108.3],
BigBend = [1.718e6, 7.5],
FortRandall = [3.6975e6, 141.2],
GavinsPoint = [365500.0, 41.2],
)
A=[
-1 0 0 0 0 0
1 -1 0 0 0 0
0 1 -1 0 0 0
0 0 1 -1 0 0
0 0 0 1 -1 0
0 0 0 0 1 -1
]
A1 = [
-1 0 0 0 0 0
1 -1 0 0 0 0
0 1 -1 0 0 0
0 0 0.8 -1 0 0
0 0 0 0.7 -1 0
0 0 0 0 1 -1
]
A2 = [
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0.2 0 0 0
0 0 0 0.3 0 0
0 0 0 0 0 0
]
initial_storage_factor = [1, 1, 1, 0.95, 0.85, 0.98]
Cascade_system = [
Reservoir(
4088000/10^6,
18463000/10^6,
input[1,2]*initial_storage_factor[1]/10^6,
[8800/10^3,7200/10^3,8800/10^3,7200/10^3,7200/10^3],
[43500/10^3,18250/10^3,43500/10^3,40000/10^3,40000/10^3],
230000/10^3,
6000/10^3,
6300/10^3,
5.84,
5,
12000/10^3,
9000/10^3,
input[2,2],
109,
),
Reservoir(
4794000/10^6,
21956000/10^6,
input[1,3]*initial_storage_factor[2]/10^6,
[8200/10^3,8200/10^3,8200/10^3,8200/10^3,8200/10^3],
[121600/10^3,121600/10^3,121600/10^3,109250/10^3,109250/10^3],
660000/10^3,
9000/10^3,
16100/10^3,
12.3,
5,
12000/10^3,
9000/10^3,
input[2,3],
145,
),
Reservoir(
5315000/10^6,
21876000/10^6,
input[1,4]*initial_storage_factor[3]/10^6,
[7800/10^3,7800/10^3,7800/10^3,7800/10^3,7800/10^3,7800/10^3,7800/10^3],
[112290/10^3,112290/10^3,112290/10^3,112290/10^3,112290/10^3,112290/10^3,112290/10^3],
80000/10^3,0,19550/10^3, 12.45, 7 , NaN, NaN,input[2,4],141
),
Reservoir(
1631000/10^6,
1749000/10^6,
input[1,5]*initial_storage_factor[4]/10^6,
[12875/10^3,12875/10^3,12875/10^3,12875/10^3,12875/10^3,12875/10^3,12875/10^3,12875/10^3],
[67275/10^3,67275/10^3,67275/10^3,67275/10^3,67275/10^3,67275/10^3,67275/10^3,67275/10^3],
270000/10^3,
0,
18950/10^3,
4.86,
8,
NaN,
NaN,
input[2,5],
33,
),
Reservoir(
1469000/10^6,
4307000/10^6,
input[1,6]*initial_storage_factor[5]/10^6,
[5562/10^3,5562/10^3,5562/10^3,5562/10^3,5562/10^3,5562/10^3,5562/10^3,5562/10^3],
[40000/10^3,40000/10^3,40000/10^3,40000/10^3,40000/10^3,40000/10^3,40000/10^3,40000/10^3],
508000/10^3,
9000/10^3,
21700/10^3,
8.50,
8,
17000/10^3,
12000/10^3,
input[2,6],
43,
),
Reservoir(
295000/10^6,
374000/10^6,
input[1,7]*initial_storage_factor[6]/10^6,
[12000/10^3,12000/10^3,12000/10^3],
[44100/10^3,44100/10^3,44100/10^3],
345000/10^3,
9000/10^3,
24500/10^3,
3.45,
3,
25000/10^3,
10000/10^3,
input[2,7],
13,
)
]
model = SDDP.PolicyGraph(
graph;
sense = :Min,
lower_bound = 0.0,
optimizer = Gurobi.Optimizer,
) do subproblem, node
@variables(subproblem, begin
Cascade_system[r].min_s <= storage[r = 1:6] <= Cascade_system[r].max_s, SDDP.State, (initial_value = Cascade_system[r].storage_initial)
release[r = 1:6] >= Cascade_system[r].release_min, SDDP.State, (initial_value = Cascade_system[r].release_initial)
0 <= gen_output[r = 1:6, u = 1:Cascade_system[r].n_unit] <= Cascade_system[r].max_pwr[u]
0 <= power_flow[r = 1:6, u = 1:Cascade_system[r].n_unit] <= Cascade_system[r].arc_cap[u]
0 <= spill[r = 1:6] <= Cascade_system[r].spill_cap
total_gen[r = 1:6] >= 1000 * Cascade_system[r].min_gen
total_pflow[1:6] >= 0
end)
@constraints(subproblem, begin
# sum(total_gen) == df[node, :Load]
[r=1:6,u=1:Cascade_system[r].n_unit], gen_output[r,u] == power_flow[r,u]* Cascade_system[r].k
[r=1:6], total_gen[r] == sum(gen_output[r,u] for u=1:Cascade_system[r].n_unit)
[r=1:6], total_pflow[r] == sum(power_flow[r,u] for u=1:Cascade_system[r].n_unit)
[r = 1:6], release[r].out == total_pflow[r] + spill[r]
[r = [1, 5, 6]], -Cascade_system[r].max_decrease <= release[r].out - release[r].in
[r = [1, 2]], Cascade_system[r].max_increase >= release[r].out - release[r].in
[r = 1:6], storage[r].out == storage[r].in + (df[node, Symbol("Inflow$r")] + sum(release[j].out*A[r,j] for j=1:r))*1.9835
end)
@stageobjective(subproblem, sum(spill))
end
SDDP.train(model; iteration_limit = 5) |
Hi Oscar,
Thank you for the feedback. I'll implement your suggestions.
*Kind Regards,*
*Mohammad Bolboli Zadeh**Ph.D. *
*studentCivil and Environmental Engineering DepartmentCenter for
Hydrometeorology and Remote Sensing (CHRS)University of California, Irvine*
…On Tue, Aug 1, 2023 at 7:07 PM Oscar Dowson ***@***.***> wrote:
The second problem is infeasible because there is a mis-match between the
units of load and total_gen. I'd double check your data inputs.
Here are a few other suggestions:
- Don't use GLPK. It is a very poor solver for SDDP. Use Gurobi, or at
minimum, HiGHS.
- Don't use large coefficients. The lower bound of 1e10 on the value
function is too large. Since the spill is non-negative, use a lower bound
of 0.
- Format your code for better readability
- I don't understand why you needed to use eval and QuoteNode.
This code doesn't run because of the unit issue, but it should point you
in the right direction:
using SDDPusing GLPKusing DataFrames
struct Reservoir
min_s::Float64
max_s::Float64
storage_initial::Float64
arc_cap::Vector{Float64}
max_pwr::Vector{Float64}
spill_cap::Float64
release_min::Float64
release_initial::Float64
k::Float64
n_unit::Int64
max_decrease::Float64
max_increase::Float64
min_gen::Float64
evap_avg_monthly::Float64end
df = DataFrame(
Node = [1, 2, 3, 4, 5, 6],
Parent = [0, 1, 1, 2, 2, 3],
Prob = [1.0, 0.68, 0.32, 0.65, 0.35, 1.0],
Stage = [1, 2, 2, 3, 3, 3],
Inflow1 = [6, 7, 11, 8, 7, 11],
Inflow2 = [18, 17, 28, 16, 22, 28],
Inflow3 = [18, 18, 22, 19, 18, 24],
Inflow4 = [0.58, 0.0, 0.0, 0.0, 0.0, 0.6],
Inflow5 = [0.08, 1.55, 0.93, 0.0, 1.23, 6.54],
Inflow6 = [3.19, 4.1, 1.39, 3.6, 2.64, 2.25],
Load = [1.87112, 1.60109, 1.67808, 1.6246, 1.64542, 1.69742],
probability = [1.0, 0.68, 0.32, 0.442, 0.238, 0.32],
)
graph = SDDP.Graph(0)for i in 1:length(df.Node)
SDDP.add_node(graph, i)
SDDP.add_edge(graph, df.Parent[i] => df.Node[i], df.Prob[i])end
input = DataFrame(
Row = ["initial_storage", "min_daily_ene_gen"],
FortPeck = [1.5208e7, 57.9],
Garrison = [1.8195e7, 183.8],
Oahe = [1.90485e7, 108.3],
BigBend = [1.718e6, 7.5],
FortRandall = [3.6975e6, 141.2],
GavinsPoint = [365500.0, 41.2],
)
A=[
-1 0 0 0 0 0
1 -1 0 0 0 0
0 1 -1 0 0 0
0 0 1 -1 0 0
0 0 0 1 -1 0
0 0 0 0 1 -1
]
A1 = [
-1 0 0 0 0 0
1 -1 0 0 0 0
0 1 -1 0 0 0
0 0 0.8 -1 0 0
0 0 0 0.7 -1 0
0 0 0 0 1 -1
]
A2 = [
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0.2 0 0 0
0 0 0 0.3 0 0
0 0 0 0 0 0
]
initial_storage_factor = [1, 1, 1, 0.95, 0.85, 0.98]
Cascade_system = [
Reservoir(
4088000/10^6,
18463000/10^6,
input[1,2]*initial_storage_factor[1]/10^6,
[8800/10^3,7200/10^3,8800/10^3,7200/10^3,7200/10^3],
[43500/10^3,18250/10^3,43500/10^3,40000/10^3,40000/10^3],
230000/10^3,
6000/10^3,
6300/10^3,
5.84,
5,
12000/10^3,
9000/10^3,
input[2,2],
109,
),
Reservoir(
4794000/10^6,
21956000/10^6,
input[1,3]*initial_storage_factor[2]/10^6,
[8200/10^3,8200/10^3,8200/10^3,8200/10^3,8200/10^3],
[121600/10^3,121600/10^3,121600/10^3,109250/10^3,109250/10^3],
660000/10^3,
9000/10^3,
16100/10^3,
12.3,
5,
12000/10^3,
9000/10^3,
input[2,3],
145,
),
Reservoir(
5315000/10^6,
21876000/10^6,
input[1,4]*initial_storage_factor[3]/10^6,
[7800/10^3,7800/10^3,7800/10^3,7800/10^3,7800/10^3,7800/10^3,7800/10^3],
[112290/10^3,112290/10^3,112290/10^3,112290/10^3,112290/10^3,112290/10^3,112290/10^3],
80000/10^3,0,19550/10^3, 12.45, 7 , NaN, NaN,input[2,4],141
),
Reservoir(
1631000/10^6,
1749000/10^6,
input[1,5]*initial_storage_factor[4]/10^6,
[12875/10^3,12875/10^3,12875/10^3,12875/10^3,12875/10^3,12875/10^3,12875/10^3,12875/10^3],
[67275/10^3,67275/10^3,67275/10^3,67275/10^3,67275/10^3,67275/10^3,67275/10^3,67275/10^3],
270000/10^3,
0,
18950/10^3,
4.86,
8,
NaN,
NaN,
input[2,5],
33,
),
Reservoir(
1469000/10^6,
4307000/10^6,
input[1,6]*initial_storage_factor[5]/10^6,
[5562/10^3,5562/10^3,5562/10^3,5562/10^3,5562/10^3,5562/10^3,5562/10^3,5562/10^3],
[40000/10^3,40000/10^3,40000/10^3,40000/10^3,40000/10^3,40000/10^3,40000/10^3,40000/10^3],
508000/10^3,
9000/10^3,
21700/10^3,
8.50,
8,
17000/10^3,
12000/10^3,
input[2,6],
43,
),
Reservoir(
295000/10^6,
374000/10^6,
input[1,7]*initial_storage_factor[6]/10^6,
[12000/10^3,12000/10^3,12000/10^3],
[44100/10^3,44100/10^3,44100/10^3],
345000/10^3,
9000/10^3,
24500/10^3,
3.45,
3,
25000/10^3,
10000/10^3,
input[2,7],
13,
)
]
model = SDDP.PolicyGraph(
graph;
sense = :Min,
lower_bound = 0.0,
optimizer = Gurobi.Optimizer,
) do subproblem, node
@variables(subproblem, begin
Cascade_system[r].min_s <= storage[r = 1:6] <= Cascade_system[r].max_s, SDDP.State, (initial_value = Cascade_system[r].storage_initial)
release[r = 1:6] >= Cascade_system[r].release_min, SDDP.State, (initial_value = Cascade_system[r].release_initial)
0 <= gen_output[r = 1:6, u = 1:Cascade_system[r].n_unit] <= Cascade_system[r].max_pwr[u]
0 <= power_flow[r = 1:6, u = 1:Cascade_system[r].n_unit] <= Cascade_system[r].arc_cap[u]
0 <= spill[r = 1:6] <= Cascade_system[r].spill_cap
total_gen[r = 1:6] >= 1000 * Cascade_system[r].min_gen
total_pflow[1:6] >= 0
end)
@Constraints(subproblem, begin
# sum(total_gen) == df[node, :Load]
[r=1:6,u=1:Cascade_system[r].n_unit], gen_output[r,u] == power_flow[r,u]* Cascade_system[r].k
[r=1:6], total_gen[r] == sum(gen_output[r,u] for u=1:Cascade_system[r].n_unit)
[r=1:6], total_pflow[r] == sum(power_flow[r,u] for u=1:Cascade_system[r].n_unit)
[r = 1:6], release[r].out == total_pflow[r] + spill[r]
[r = [1, 5, 6]], -Cascade_system[r].max_decrease <= release[r].out - release[r].in
[r = [1, 2]], Cascade_system[r].max_increase >= release[r].out - release[r].in
[r = 1:6], storage[r].out == storage[r].in + (df[node, Symbol("Inflow$r")] + sum(release[j].out*A[r,j] for j=1:r))*1.9835
end)
@stageobjective(subproblem, sum(spill))end
SDDP.train(model; iteration_limit = 5)
—
Reply to this email directly, view it on GitHub
<https://urldefense.com/v3/__https://github.com/odow/SDDP.jl/issues/646*issuecomment-1661378421__;Iw!!CzAuKJ42GuquVTTmVmPViYEvSg!PIXGgpdoYNq0zRuFe6Aj366_NsFJQd-ucC4HlDtfk-wgZpZX_YEXseXtmiCscf9S-EMSHe0lf-NthuAoBgnzPGM6$>,
or unsubscribe
<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/BBR3HAIFV2O4VOGBLBRMBFTXTGY7RANCNFSM6AAAAAA22RPRZQ__;!!CzAuKJ42GuquVTTmVmPViYEvSg!PIXGgpdoYNq0zRuFe6Aj366_NsFJQd-ucC4HlDtfk-wgZpZX_YEXseXtmiCscf9S-EMSHe0lf-NthuAoBqlPQcUD$>
.
You are receiving this because you authored the thread.Message ID:
***@***.***>
|
Dear Oscar
but they didn't solve my problem and also after changing their values i didn't see
do you have any idea on how to proceed or why change in attributes won't effect matrix range? |
Please provide a reproducible example of your new code. If you have infeasibilities in your scaled model, it's likely because you still have a modeling error somewhere. Try commenting out the constraints and re-add them one-by-one. Which one causes the problem? Does it happen on the first iteration, or after some time? Your model must be feasible for any set of releases chosen in the previous time-step. |
Thanks for your comments. i've worked on the model and solved its inconsistencies. The model now runs with no problem on the small tree (the previous code)
|
The warning is just a warning. In general, you want coefficients that are between 1e-04 and 1e+04. This isn't a requirement though.
Things are probably still fine. p.s. I still don't know why you need |
Thanks, i'll take a look at it |
Any updates? |
not really, we still get warnings for some nodes. And haven't been able to solve it yet. |
It doesn't matter if you get some warnings. They're warnings, not errors. It's a hint, not a requirement. |
Yes. we probably going to ignore them. |
Closing as this seems resolved. Please re-open if you have further questions 😄 |
I am working on a Multistage Stochastic Optimization problem, for dynamic cascade reservoir optimization. Same as post #645. But I'm having numerical stability issues. Sometimes it's just a warning, and sometimes it leads to an error. I've tried to scale the problem, but I'm still having issues.
You can find the codes and necessary files in this folder. The summary.xlsx contains the summary of runs.
!C-1_not_scaled is code for the problem without scaling.
!C-1_scaled is code for the problem with scaling.
https://www.dropbox.com/scl/fo/jjbjnwmruo8upacamuhh2/h?rlkey=kq06gwbxvslk19sp2k9ifiwkv&dl=0
The text was updated successfully, but these errors were encountered: