search

odow / SDDP.jl

Stochastic Dual Dynamic Programming in Julia
https://sddp.dev
Other
289 stars 60 forks source link

Numerical stability issues #646

Closed bolbolim closed 11 months ago

bolbolim commented 1 year ago

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

odow commented 1 year ago

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).

bolbolim commented 1 year ago

unscaled problem that gives numerical instability warning

import Pkg
Pkg.add("JSON")
Pkg.add("Plots")
Pkg.add("DataFrames")
Pkg.add("XLSX")
Pkg.add("GLPK")
import XLSX
Pkg.add("CSV")
using SDDP, GLPK, Test
using Pkg
using DataFrames,XLSX,CSV
using JSON
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 = [6000, 7000, 11000, 8000, 7000, 11000]
inflow2 = [18000, 17000, 28000, 16000, 22000, 28000]
inflow3 = [18000, 18000, 22000, 19000, 18000, 24000]
inflow4 = [580, 0, 0, 0, 0, 600]
inflow5 = [80, 1550, 930, 0, 1230, 6540]
inflow6 = [3190, 4100, 1390, 3600, 2640, 2250]
load = [1871.12, 1601.09, 1678.08, 1624.6, 1645.42, 1697.42]
probability = [1.0, 0.68, 0.32, 0.442, 0.238, 0.32]
df = DataFrame(
    Node = node,
    Parent = parent,
    Prob = prob,
    Stage = stage,
    Inflow1 = inflow1,
    Inflow2 = inflow2,
    Inflow3 = inflow3,
    Inflow4 = inflow4,
    Inflow5 = inflow5,
    Inflow6 = inflow6,
    Load = load,
    probability = probability
)
n_stages=maximum(df.Stage)
graph = SDDP.Graph(:0)
for i = 1: length(df.Node)
    a= df.Parent[i]
    b= df.Node[i]
    c= Float64(df.Prob[i])
    exp  = :(SDDP.add_node(graph, $(QuoteNode(i))))
    exp2 = :(SDDP.add_edge(graph, $(QuoteNode(a)) => $(QuoteNode(b)), $c))
    eval(exp)
    eval(exp2)
end
graph
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
rows = ["initial_storage", "min_daily_ene_gen"]
fort_peck = [1.5208e7, 57.9]
garrison = [1.8195e7, 183.8]
oahe = [1.90485e7, 108.3]
big_bend = [1.718e6, 7.5]
fort_randall = [3.6975e6, 141.2]
gavins_point = [365500.0, 41.2]

input = DataFrame(Row = rows, FortPeck = fort_peck, Garrison = garrison, Oahe = oahe,
               BigBend = big_bend, FortRandall = fort_randall, GavinsPoint = gavins_point)
   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,.98,.85,.99]   
   Cascade_system = [
   Reservoir(4088000,18463000,input[1,2]*initial_storage_factor[1],[8800,7200,8800,7200,7200], 
                                       [43500,18250,43500,40000,40000],230000,6000, 8000, 5.84, 5, 12000,9000,input[2,2],109),
   Reservoir(4794000,21956000,input[1,3]*initial_storage_factor[2],[8200,8200,8200,8200,8200], 
                                       [121600,121600,121600,109250,109250],660000,9000,21000, 12.3, 5,12000, 9000,input[2,3],145),
   Reservoir(5315000,21876000,input[1,4]*initial_storage_factor[3],[7800,7800,7800,7800,7800,7800,7800], 
                                       [112290,112290,112290,112290,112290,112290,112290],80000,0,24000, 12.45, 7 , NaN, NaN,input[2,4],141),
   Reservoir(1631000,1749000,input[1,5]*initial_storage_factor[4],[12875,12875,12875,12875,12875,12875,12875,12875],
                                     [67275,67275,67275,67275,67275,67275,67275,67275],270000,0,27000, 4.86, 8, NaN, NaN,input[2,5],33),
   Reservoir(1469000,4307000,input[1,6]*initial_storage_factor[5],[5562,5562,5562,5562,5562,5562,5562,5562],
                                     [40000,40000,40000,40000,40000,40000,40000,40000],508000,9000,24000, 8.50, 8, 17000, 12000,input[2,6],43),
   Reservoir(295000 ,374000,input[1,7]*initial_storage_factor[6], [12000,12000,12000], 
                                   [44100,44100,44100],345000,9000,25000, 3.45, 3, 25000, 10000,input[2,7],13)]
 model = SDDP.PolicyGraph(
       graph,
       sense = :Min,
       bellman_function = SDDP.BellmanFunction(lower_bound = -1e+10, cut_type = SDDP.SINGLE_CUT),
       optimizer = GLPK.Optimizer,

       )do subproblem, node

           @variable(subproblem, Cascade_system[r].min_s <= storage[r=1:6] <= Cascade_system[r].max_s, SDDP.State, 
                       initial_value = Cascade_system[r].storage_initial)

           # state variable release, release at each stage and has lower bound of minimum release
           @variable(subproblem, release[r=1:6] >= Cascade_system[r].release_min, SDDP.State, initial_value = Cascade_system[r].release_initial)

           @variables(subproblem, begin

           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]

           inflow[r=1:6] == df[!,r+4][node]

           totaldemand == df[!,11][node]

           0 <= spill[r=1:6] <= Cascade_system[r].spill_cap

           total_gen[r=1:6] >=  Cascade_system[r].min_gen*1000  # (mw to kw)
           total_pflow[r=1:6] >=0
           end)

           @constraint(subproblem, demand_c, sum(total_gen[r]/1000  for r=1:6) == totaldemand)                                    #totalgen/1000 kw to mw 

           @constraint(subproblem, flow_gen[r=1:6,u=1:Cascade_system[r].n_unit], gen_output[r,u]  ==  power_flow[r,u]* Cascade_system[r].k    
           )

           @constraint(subproblem,gen[r=1:6], total_gen[r] == sum(gen_output[r,u] for u=1:Cascade_system[r].n_unit))

           @constraint(subproblem,pflow[r=1:6],  total_pflow[r] == sum(power_flow[r,u] for u=1:Cascade_system[r].n_unit))

            @constraint(subproblem,release_c1[r=1:6],release[r].out == total_pflow[r] + spill[r])

           @constraint(subproblem,release_c2[r=[1,5,6]],-Cascade_system[r].max_decrease <= release[r].out - release[r].in )                                          
           @constraint(subproblem,release_c3[r=[1,2]],  Cascade_system[r].max_increase >= release[r].out - release[r].in ) 

           @constraints(subproblem,begin
                   storage_equa[r = 1:6],
                   storage[r].out == storage[r].in + (inflow[r] + sum(release[j].out*A[r,j] for j=1:r))*1.9835
           end)

          @stageobjective(subproblem, sum(spill[r] for r=1:6))
   end

SDDP.train(model, stopping_rules = [SDDP.Statistical(; num_replications=5, iteration_period = 1, z_score = 1.96,
            verbose = true)],run_numerical_stability_report = true, add_to_existing_cuts=false)  #true if we want to add existing_cuts from last training 
SDDP.numerical_stability_report (

    model,
    by_node = true,
    print = true,
    warn = true,
)
bolbolim commented 1 year ago

scaled model that crashes

import Pkg
Pkg.add("JSON")
Pkg.add("Plots")
Pkg.add("DataFrames")
Pkg.add("XLSX")
Pkg.add("GLPK")
import XLSX
Pkg.add("CSV")
##Julia model
using SDDP, GLPK, Test
using Pkg
using DataFrames,XLSX,CSV
using JSON
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]
df = DataFrame(Node = node, Parent = parent, Prob = prob, Stage = stage,
               Inflow1 = inflow1, Inflow2 = inflow2, Inflow3 = inflow3,
               Inflow4 = inflow4, Inflow5 = inflow5, Inflow6 = inflow6,
               Load = load, probability = probability)
n_stages=maximum(df.Stage)
graph = SDDP.Graph(:0)
for i = 1: length(df.Node)
    a= df.Parent[i]
    b= df.Node[i]
    c= Float64(df.Prob[i])
    exp  = :(SDDP.add_node(graph, $(QuoteNode(i))))
    exp2 = :(SDDP.add_edge(graph, $(QuoteNode(a)) => $(QuoteNode(b)), $c))
    eval(exp)
    eval(exp2)
end
graph

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
rows = ["initial_storage", "min_daily_ene_gen"]
fort_peck = [1.5208e7, 57.9]
garrison = [1.8195e7, 183.8]
oahe = [1.90485e7, 108.3]
big_bend = [1.718e6, 7.5]
fort_randall = [3.6975e6, 141.2]
gavins_point = [365500.0, 41.2]
input = DataFrame(Row = rows, FortPeck = fort_peck, Garrison = garrison, Oahe = oahe,
               BigBend = big_bend, FortRandall = fort_randall, GavinsPoint = gavins_point)
input[2,7]
   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,.95,.85,.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,

       bellman_function = SDDP.BellmanFunction(lower_bound = -1e+10, cut_type = SDDP.SINGLE_CUT),
       optimizer = GLPK.Optimizer,
       )do subproblem, node

           @variable(subproblem, Cascade_system[r].min_s <= storage[r=1:6] <= Cascade_system[r].max_s, SDDP.State, 
                       initial_value = Cascade_system[r].storage_initial)

           @variable(subproblem, release[r=1:6] >= Cascade_system[r].release_min, SDDP.State, initial_value = Cascade_system[r].release_initial)

           @variables(subproblem, begin

           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]

           inflow[r=1:6] == df[!,r+4][node]

           totaldemand == df[!,11][node]

           0 <= spill[r=1:6] <= Cascade_system[r].spill_cap

           total_gen[r=1:6] >=  Cascade_system[r].min_gen*1000
           total_pflow[r=1:6] >=0
           end)

           @constraint(subproblem, demand_c, sum(total_gen[r]  for r=1:6) == totaldemand)                                    #totalgen/1000 kw to mw 

           @constraint(subproblem, flow_gen[r=1:6,u=1:Cascade_system[r].n_unit], gen_output[r,u]  ==  power_flow[r,u]* Cascade_system[r].k    
           )

           @constraint(subproblem,gen[r=1:6], total_gen[r] == sum(gen_output[r,u] for u=1:Cascade_system[r].n_unit))

           @constraint(subproblem,pflow[r=1:6],  total_pflow[r] == sum(power_flow[r,u] for u=1:Cascade_system[r].n_unit))

            @constraint(subproblem,release_c1[r=1:6],release[r].out == total_pflow[r] + spill[r])
           @constraint(subproblem,release_c2[r=[1,5,6]],-Cascade_system[r].max_decrease <= release[r].out - release[r].in )                                          
           @constraint(subproblem,release_c3[r=[1,2]],  Cascade_system[r].max_increase >= release[r].out - release[r].in ) 
            @constraints(subproblem,begin
                    storage_equa[r = 1:6],
                    storage[r].out == storage[r].in + (inflow[r] + sum(release[j].out*A[r,j] for j=1:r))*1.9835
            end)

          @stageobjective(subproblem, sum(spill[r] for r=1:6))

   end
SDDP.train(model, stopping_rules = [SDDP.Statistical(; num_replications=5, iteration_period = 1, z_score = 1.96,
            verbose = true)],run_numerical_stability_report = true, add_to_existing_cuts=false)  #true if we want to add existing_cuts from last training 
SDDP.numerical_stability_report(

    model,
    by_node = true,
    print = true,

    warn = true,
)
odow commented 12 months ago

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:

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)
bolbolim commented 12 months ago

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: @.***>

bolbolim commented 12 months ago

Dear Oscar Hi, I've corrected the inconsistencies in units of my model, used lower bound=0 and also used Gurobi optimizer. However, I still faced the same issues (unfeasibility in scaled model and numerical Warning in un_scaled model). I also tried different optimizer attributes with my model, like ScaleFlag etc On my unscaled problem (see code below) (see https://www.gurobi.com/documentation/9.5/refman/scaleflag.html for more information on scaleflag)

import MathOptInterface as MOI,
  model = SDDP.PolicyGraph(
        graph,
        sense = :Min,
        lower_bound = 0,
       optimizer= MOI.OptimizerWithAttributes(Gurobi.Optimizer,
        "ScaleFlag"=>2, #scale matrix
        "NumericFocus"=>3,
    )

but they didn't solve my problem and also after changing their values i didn't see any change in my matrix range

  Non-zero Matrix range     [1e-03, 1e+01]
  Non-zero Objective range  [1e+00, 1e+00]
  Non-zero Bounds range     [6e+03, 2e+07]
  Non-zero RHS range        [2e+02, 6e+04]

do you have any idea on how to proceed or why change in attributes won't effect matrix range?

odow commented 12 months ago

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.

bolbolim commented 11 months ago

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) but when I want to run it for the big tree I get numerical stability warning for some nodes(9,14,113,174,225,241,) any idea on what causes this or how to solve it? You can see a reproducible code below.

using SDDP, GLPK, Test
using Pkg
using DataFrames,XLSX,CSV
using JSON
using Gurobi
using JuMP

##defining factors
initial_storage_factor=[1,1,.99,.98,.85,.99]/10^6 
min_gen_factor=1
constant=1000

##define tree columns
node=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316]
parent=[0,1,1,2,2,3,4,5,5,6,6,7,8,9,10,10,11,15,12,12,13,14,15,16,17,20,18,19,20,21,22,23,24,25,26,32,27,28,29,32,30,31,32,33,34,35,36,37,38,39,40,41,42,43,43,44,45,46,47,48,50,49,50,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,69,70,71,72,73,74,77,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,93,96,97,98,99,100,100,101,102,103,104,105,106,107,108,109,110,117,111,112,113,114,115,116,117,118,119,120,136,121,122,123,136,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,148,153,154,155,156,157,158,159,160,161,162,169,163,164,169,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,203,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,228,209,210,211,212,213,214,215,216,217,218,219,220,221,218,222,223,224,225,226,227,228,229,230,231,232,233,244,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,283,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,280,283,284,285,286]
prob=[1,0.68,0.32,0.647058824,0.352941176,1,1,0.583333333,0.416666667,0.9375,0.0625,1,1,1,0.933333333,0.066666667,1,0.071428571,0.045454545,0.954545455,1,1,0.928571429,1,1,0.047619048,1,1,0.952380952,1,1,1,1,1,1,0.076923077,1,1,1,0.076923077,1,1,0.846153846,1,1,1,1,1,1,1,1,1,1,0.090909091,0.909090909,1,1,1,1,1,0.05,1,0.9,0.05,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.1,0.9,1,1,1,1,1,0.055555556,1,1,0.944444444,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.941176471,1,1,0.058823529,1,1,1,1,0.888888889,0.111111111,1,1,1,1,1,1,1,1,1,1,0.125,1,1,1,1,1,1,0.875,1,1,1,0.142857143,1,1,1,0.142857143,1,1,1,1,1,1,1,1,1,1,1,1,0.714285714,1,1,1,1,1,1,1,1,1,1,1,0.9375,1,1,1,1,0.0625,1,1,1,1,1,1,1,1,1,1,0.066666667,1,1,0.066666667,1,1,1,1,0.866666667,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.8,1,1,1,1,1,0.25,1,1,1,1,1,1,1,1,1,0.923076923,1,1,1,0.076923077,1,1,1,1,1,1,0.75,1,1,1,1,1,0.083333333,1,1,1,1,1,1,1,1,1,1,0.916666667,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.333333333,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.8,1,1,0.2,0.666666667,1,1,1,]
stage=[1,2,2,3,3,3,4,4,4,4,4,5,5,5,5,5,5,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10,10,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,]
inflow1=[6,7,11,8,7,11,8,7,5,11,14,9,6,5,12,14,14,16,7,9,7,5,12,14,13,12,18,8,9,7,6,10,14,13,13,19,25,8,9,9,6,8,10,13,13,11,17,25,8,9,11,6,9,8,10,13,13,11,15,25,10,8,8,9,13,7,9,7,10,15,14,12,20,33,9,8,9,8,11,6,6,8,13,10,25,15,13,12,25,9,9,7,9,6,10,6,5,10,11,10,36,16,15,9,20,9,10,7,9,7,11,10,6,5,9,13,10,9,34,16,17,13,18,9,12,7,8,9,22,11,10,6,5,11,13,10,9,31,17,15.5,18,13,16,18,9,12,7,8,10,21,10,10,6,4,13,16,10,9,32,18,15,24,13,25,18,8,13,7,7,8,24,10,8,4,7,4,10,14,11,9,40,18,14,27,7,29,25,8,16,8,12,7,7,10,20,10,8,4,8,4,9,13,12,10,45,16,14,15.6,19,6,9,25,6,15,8,11,8,8,10,18,10,9,4,9,5,9,13,13,11,39,15,15,14.7,9,18,6,12,25,6,14,8,10,8,8,10,18,11,9,9,5,8,5,8,14,13,11,36,14,16.5,14.3,24,8,17,6,15,23,6,12,10,10,7,9,10,15,15,12,9,5,8,6,11,12,13,11,34,14,12.3,17,14,8,8,22,7,20,22,8,11,12,11,7,9,10,17,17,11,9,5,7,5,9,13,10,12,11,30,14,]
inflow2=[18,17,28,16,22,28,13,18,13,28,55,12,15,11,31,46,46,40,25,12,14,9,34,39,43,20,40,26,13,15,10,35,36,43,30,40,35,24,16,28,15,9,35,35,38,25,38,40,24,17,32,17,12,32,34,33,36,25,38,45,15,20,18,13,33,18,12,31,35,33,34,21,38,80,17,18,20,12,38,16,11,31,67,34,36,34,21,38,75,20,17,18,19,11,29,17,10,32,70,30,38,30,20,35,70,20,17,18,20,9,32,24,19,8,30,57,27,25,41,28,35,32,70,20,17,18,21,10,34,36,25,16,10,27,66,26,26,44,28,37,35,31,70,59,15,17,18,15,27,32,28,25,13,8,29,49,27,28,48,28,38,39.5,30,70,45,14,17,18,12,20,35,22,14,14,14,8,32,55,27,28,52,28,38.5,46.5,17,45,90,10,43,14,21,16,14,15,35,25,12,12,13,10,31,53,26,28,63,27,38,30.5,46.5,15,45,160,10,46,14,21,14,14,12,35,29,16,15,14,7,30,52,28,28,68,26,35.5,31,30,48,14,35,160,10,50,14,21,13,13,12,31,28,16,18,14,17,7,29,46,27,28,56,25,33,30.5,29,32,59,14,30,100,10,44,15,20,12,13,15,27,29,16,18,15,16,9,32,44,27,28,53,22,28,32,30.5,29,25,60,14,30,85,10,32,15,20,11,13,17,25,32,14,19,12,13,9,32,49,16,27,28,53,21,]
inflow3=[18,18,22,19,18,24,18,17,13,25,58,14,17,13,25,31,56,32,59,17,19,17,26,31,55,24,32,36,16,18,19,31,32,51,24,35,36,35,15,24,18,18,32,31,49,24,50,48,34,15,22,19,20,42,41,30,46,24,30,40,23,35,15,26,18,18,19,55,40,32,43,24,40,43,27,31,14,19,21,19,21,61,48,40,35,41,30,15,40,23,30,26,13,18,25,18,22,91,55,35,37,38,38,55,45,33,28,21,14,17,30,38,18,19,121,48,35,30,36,35,44,30,51,30,25,16,16,30,37,34,38,19,19,139,41,35,30,34,32,43.5,40,35,70,59,25,22,12,14,40,41,43,35,19,17,147,38,35,30,34,35,44,40,15,45,50,24,19,10,14,38,37,51,30,16,19,17,160,45,34,31,36,50,44.5,40,20,40,40,15,45,20,37,10,13,36,40,67,29,14,18,17,153,37,33,31,37,49,47.5,42.5,40,16,30,40,15,40,20,58,10,12,22,43,67,33,22,19,19,105,38,31,32,37,46,63,42.5,20,44,16,35,80,15,34,20,68,10,14,27,39,78,18,33,22,20,19,75,35,31,34,37,40,69,42.5,38,28,46,16,30,83,15,21,17,76,10,15,31,40,92,18,38,22,19,21,54,28,30,35,36,34,37,66.5,42.5,35,24,44,16,20,70,15,19,16,32,10,13,22,36,74,18,38,19,19,18,49,34,30,28,34,35,32,]
inflow4=[0.58,0,0,0,0,0.6,0,0,0.04,0.6,0.96,0,0,0.02,0,0,5.14,2.6,11.08,0,0,0,0,0,0,3.64,2.16,5.66,0,0.622,0.06,0.08,0,3.72,6.74,0,0,4.96,0,0.32,0,0,0,0,0.54,10.2,1,0,6.2,0,0.38,0,0.52,0,0.1,0,1.26,1.82,3.14,0.74,0.5,0.48,1.64,0,0.4,0,0.78,0,0,0,4.78,3.98,1.56,7.74,0,0.62,0,1.32,0,0,0,0,0,0,0,3.16,0.88,0,5.14,3.88,0,1.76,0.62,0,0,0,0,0.88,0,1.08,0,0,0,0.28,0,0,5.12,1.32,0,1,0,3.6,0,0.68,1.94,0,0,0,3,0,3.44,0,0,2.8,0,0.8,2.4,4.14,0.22,0.48,0.82,0,0,0,0,0.76,0,0,0,0.58,0,0,0.28,0,1.1,0,1.36,1.24,13.82,0,5.94,2.92,0,0,0,0,0,0,0,1.96,0,3.64,0.42,0,4.04,1.68,0,1.6,0,9.74,0,1.74,0,0.02,2.56,0,5.7,0,0,0,0,2.54,0,3.08,2.06,7.68,0.34,0,0,1.7,2.5,0,0,12.24,2.8,2.42,0,0,0,0,1.72,0,0,0,0,9.06,0,0.38,3.18,3.96,5,0,0,0.52,0.74,17,0,0,6.42,0,4.46,1.96,0,0,0.04,0,1.2,0,0,0,5.56,3.56,1.28,0,1.66,2.46,2.26,6.4,0,0,0,15.78,1.66,0,4.3,0.28,1.42,1.7,0.72,0,0,0,0.36,0,0,0,0,3.84,9.04,0,0,1.1,0.68,0.22,3.48,0.56,0,0,0.76,11.8,0.5,0,11.18,0,0,0,0.64,1.38,0,1.46,0,1.06,0,0,0,1.2,2.64,11.64,0.64,0,1.88,0,2.3,2.16,0,0,0.24,0.98,6.1,1.6,0.94,11.38,0.58,5.46,0,2.06,1.5,0.14,1.24,0,0.18,0,0,0,0,0]
inflow5=[0.08,1.55,0.93,0,1.23,6.54,0,1.2,0.84,5.54,10.97,1.809,0,6.82,5.6,5.56,8.08,9.3,17.52,0,6.4177,1.76,0,4.81,5.76,3.14,0,10.22,1.03,0.7922,0.33,3.48,0,7.07,13.31,0.83,0,4.11,1.04,3.768,1.9111,0,0,1.22,12.75,13.51,0.34,0,7.65,4.21,0,3.15,7.65,5.61,4.95,2.94,9.88,4.67,3.34,3.4,0.38,2.64,0,5.66,13.8,9.79,2.08,5.39,1.72,9.77,17.78,3.25,4.4,3.89,2.08,3.09,0.53,5.84,10.41,0,3.39,0.39,1.86,5.81,3.42,12.78,0,0,9.41,14.94,0.85,4.07,2.08,0,0.16,0,2.61,5.25,6.91,15.44,7.1,6.77,3.82,4.66,0,16.45,7.91,8.23,0,0.59,0,13.21,0.03,0.88,12.81,4.03,0,5.25,10.2,5.03,0,0,0,0,0,1.73,1.57,2.03,0,14.48,14.43,8.07,0.96,5.47,9.95,0,8.75,10.91,5.71,3.09,9.14,0,0,10.36,3.08,2.26,0,0.89,9.79,0,7.49,4.08,5.25,1.39,6.87,11.16,4.12,7.34,0,9.81,3.06,4.96,0.65,0,2.51,1.78,0,1.29,1.57,12.51,0,4.42,0,2.07,1.95,0,17.23,2.36,1.56,2.1,0,16.2,0,9.41,0,4.59,0,4.62,0,3.57,1.56,0,0,16.33,7.68,4.44,0,0.96,0.47,0,10.46,2.99,0,1.86,8.5,20.73,4.42,10.49,10.56,10.44,22.27,1.2,0,0,2.69,14.45,1.84,0,14.77,0.32,8.66,0,2.24,4.31,6.24,1.6,10.71,0.47,4.79,8.65,17.73,1.77,1.81,3.13,0,14.29,0.41,16.77,4.14,15.75,0.05,15.03,12.26,0.86,9.54,1.58,7.23,2.5,0.55,5.72,4.89,5.14,0,6.18,8,4.62,10.61,13.71,6.6,2.93,6.45,8.91,0,3.44,7.86,4.95,0,0,3.44,10,3.22,0,12.28,0,5.75,0,5.6,6.12,0,3.44,4.77,2.37,5.52,4.98,10.26,13.31,0,13.44,1.5,5.23,4.21,0,0,0,0,10.79,0,0,0.87,2.4,0.6,19.83,0,7.89,1.25,4.39,3.42,2.6,9.37,4.52,9.03,9.4,0,1.66,7.63,15.31]
inflow6=[3.19,4.1,1.39,3.6,2.64,2.25,3.44,1.94,2.59,2.28,10.45,3.398,0.94,2.34,1.81,3.7,8.78,4.95,9.46,3.532,1.63,1.83,1.42,3.24,8.14,2.71,5.79,8.83,3.06,3.825,1.49,1.78,3.08,8.13,7.06,5.31,2.12,8.02,3.15,2.4289,2.641,1.35,0.85,3.17,7.81,8.94,3.53,3.51,9.52,3.9,4.4726,1.644,1.79,2.52,1.92,3.73,8.04,5.5,1.15,3.48,4.25,8.93,5.26,5.28,5.5,1.63,3.97,3.99,0.91,2.51,7.05,4.21,4.54,6.4,6.05,8.59,3.9,4.84,6.4,2.61,4.42,4.84,3.25,2.17,2.21,7.67,1.72,2.85,7.42,13.78,3.37,9.75,6.75,4.92,6.8,1.48,3.35,4.63,2.2,2.25,2.11,7.19,2.9,1.59,3.23,14.94,5.33,7.18,4.57,2.82,4.89,9.13,1.64,2.04,3.65,1.18,2.73,2.47,4.18,5.26,3.38,2.69,2.94,13.72,9.02,6.04,4.88,2.57,5.41,4.32,12.55,2.42,1.75,3.92,1.69,1.17,2.46,4.26,4.9,1.64,4.57,2.81,2.87,4.5,12.46,5.06,6.58,4.78,6.2,5.32,4.64,11.18,2.42,3.59,4.57,5.66,1.18,2.66,3.65,8.84,2.31,4.61,2.93,2.33,2.81,12.72,4.2,7.44,3.21,5.68,3.19,4.08,5.59,5.97,3.36,3.67,5.44,5.59,2.25,3.07,3.45,8.75,2.6,4.53,2.52,5.1,4.37,1.99,1.22,9.96,4.76,7.88,3.16,6.98,4.97,5.6,4.14,4.98,1.75,1.99,4.61,2.94,0.53,2.8,3.82,13.67,2.11,1.4,3.88,0,8.58,4.74,1.37,2.56,7.22,12.99,6.24,3.09,12.74,5.36,7.76,5.41,4,2.81,2.98,5.53,7.05,1.09,2.4,6.81,13.67,2.46,1.52,3.05,1.58,6.84,7.28,5.84,2.12,4.77,5.88,7.56,7.62,3.46,13.23,4.8,7.24,6.82,5.16,5,2.96,2.4,3.71,7.26,1.1,3.32,4.12,10.83,3.13,2.05,7.09,2.71,0.8,9.44,4.6,2.96,3.46,2.85,4.14,6.27,3.46,2.65,10.92,4.37,8.1,5.81,8.26,4.5,2.31,3.21,3.75,7.58,1.52,2.26,5.3,8.94,3.62,4.84,1.97,7.17,4.27,2.4,6.65,2.59,1.52,3.45,2.86,7.45,6.66,3.44,3.73,10.5,5.71,7.56,4.39,6.55,4.65,0.41,4.52,4.87,6.11,6.42,1.23,2.65,5.52,9]
load=[1.871122371,1.601088783,1.678076703,1.624601397,1.645421716,1.697422397,1.852998523,1.868662699,1.919092501,1.878320237,0.89931207,1.844848953,1.858355251,1.906552365,1.868729538,1.32824113,0.757056989,1.875976906,1.844568274,1.874171902,1.886767928,1.915066609,1.886739518,1.060228667,0.794517345,1.929511846,1.925816467,1.926541632,1.928918559,1.946076881,1.957033096,1.932055815,1.509715828,0.883226349,1.71917044,1.909475051,1.872593045,1.821596178,1.872305846,1.850795688,1.87652838,1.912242467,1.887261379,1.514076862,1.261459777,1.533107221,1.834635569,1.818650847,1.724980504,1.805650373,1.787321504,1.815190714,1.880662287,1.749922885,1.838739597,1.415536045,1.145129653,1.527612618,1.816391685,1.803405354,1.805299767,1.712941407,1.790771431,1.664844634,1.772881166,1.799383648,1.862216632,1.737031432,1.822598687,1.319634119,1.194596771,1.17706392,1.657688987,1.680160548,1.674518312,1.516887958,1.639424806,1.438095365,1.641902628,1.665741247,1.784273547,1.578748776,1.436263409,1.707152138,1.005969978,1.045445688,1.285363416,1.69377711,1.714994655,1.774823149,1.714424154,1.574689066,1.689844633,1.507131169,1.668791916,1.698784075,1.792804375,1.61680171,1.503781739,1.740309437,0.882710448,1.00702776,1.096494901,1.613250809,1.65728839,1.733033218,1.654236654,1.474564376,1.620245736,1.38941788,1.600320582,1.669696864,1.632899183,1.75555812,1.533614676,1.386168758,1.68800661,1.466736443,1.650939107,1.320606536,1.542771258,1.842622656,1.81995359,1.863863576,1.820317906,1.729357247,1.803243775,1.679645332,1.874451361,1.796197364,1.84035781,1.820022155,1.887807671,1.758337491,1.676853844,1.838980336,1.540213375,2,1.641426485,1.796703081,1.585492242,1.860317135,1.832002726,1.849692773,1.872239177,1.83278297,1.751381463,1.817687732,1.705754888,1.882208066,1.811361278,1.852652146,1.834226255,1.895645893,1.776838307,1.703006957,1.849692773,1.305842116,1.616358889,1.929147112,1.599751382,1.141622788,1.654851095,1.68640466,1.716537868,1.761645377,1.682707236,1.509078594,1.648121304,1.424608568,1.777204042,1.635890208,1.70357407,1.752725799,1.666710273,1.789589596,1.569375344,1.421664504,1.716537868,1.314014879,1.2350449,1.626962697,1.867309751,1.756966772,1.909845873,1.931167651,1.885828988,1.772830032,1.899167457,1.911259946,1.890098091,1.843002053,1.885213436,1.813541873,1.919818131,1.870182572,1.903257538,1.916434208,1.893375036,1.92631671,1.849796659,1.810198115,1.899167457,1.408094451,1.501197893,1.412808079,1.866418112,1.830035285,1.757657478,1.908353084,1.929633564,1.885258232,1.772670937,1.898598684,1.910358452,1.889778858,1.842779233,1.885336345,1.813651577,1.918988902,1.868733825,1.901946025,1.914760143,1.892335437,1.924370731,1.848729636,1.810220635,1.898598684,1.170290037,1.066560517,1.667065655,1.587163419,1.346999897,1.519453565,1.139888099,1.711322595,1.637037668,1.675184156,1.216075501,1.70501866,1.748045407,1.672748599,1.500673271,1.640668251,1.418685845,1.763796538,1.620221285,1.696723579,1.687880389,1.734764801,1.65271708,1.769928109,1.556244721,1.415347627,1.70501866,1.565744911,1.554818201,1.22311608,1.502667422,1.209107273,1.632581999,1.424830725,0.968261571,1.647887888,1.547400187,1.605425793,1.063554363,1.638438935,1.690363318,1.599495648,1.398291706,1.558286814,1.303224932,1.706877185,1.546436517,1.633464833,1.619597995,1.676177781,1.577163155,1.718612621,1.470682493,1.300648847,1.638438935,1.204692102,1.649306039,1.397089798,1.624268651,1.422424313,1.05257295,1.584573239,1.327084683,0.757056989,1.600762663,1.4678421,1.55204139,0.875209335,1.592239859,1.656297443,1.54419667,1.29365906,1.492917452,1.177063192,1.676229025,1.477106174,1.586994011,1.565862253,1.635663072,1.513511639,1.688013686,1.383907493,1.174141983,1.648179877,1.592239859,1.034229479,1.728439325,1.825072421]
pprobability=[1,0.68,0.32,0.44,0.24,0.32,0.44,0.14,0.1,0.3,0.02,0.44,0.14,0.1,0.28,0.02,0.02,0.02,0.02,0.42,0.14,0.1,0.26,0.02,0.02,0.02,0.02,0.02,0.4,0.14,0.1,0.26,0.02,0.02,0.02,0.02,0.02,0.02,0.4,0.02,0.14,0.1,0.22,0.02,0.02,0.02,0.02,0.02,0.02,0.4,0.02,0.14,0.1,0.02,0.2,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.36,0.02,0.02,0.14,0.1,0.02,0.2,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.36,0.02,0.02,0.14,0.1,0.02,0.02,0.18,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.34,0.02,0.02,0.14,0.1,0.02,0.02,0.18,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.32,0.02,0.02,0.02,0.14,0.1,0.02,0.02,0.16,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.32,0.02,0.02,0.02,0.02,0.14,0.1,0.02,0.02,0.14,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.32,0.02,0.02,0.02,0.02,0.14,0.1,0.02,0.02,0.1,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.3,0.02,0.02,0.02,0.02,0.02,0.14,0.1,0.02,0.02,0.1,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.26,0.02,0.02,0.02,0.02,0.02,0.14,0.1,0.02,0.02,0.1,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.26,0.02,0.02,0.02,0.02,0.02,0.14,0.1,0.02,0.02,0.08,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.24,0.02,0.02,0.02,0.02,0.02,0.02,0.14,0.1,0.02,0.02,0.06,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.22,0.02,0.02,0.02,0.02,0.02,0.02,0.14,0.1,0.02,0.02,0.06,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.22,0.02,0.02,0.02,0.02,0.02,0.02,0.14,0.08,0.02,0.02,0.02,0.04,0.02,0.02,0.02]

##construct tree
df = DataFrame(Node = node, Parent = parent, Prob = prob, Stage = stage,
               Inflow1 = inflow1, Inflow2 = inflow2, Inflow3 = inflow3,
               Inflow4 = inflow4, Inflow5 = inflow5, Inflow6 = inflow6,
               Load = load, probability = pprobability)

n_stages=maximum(df.Stage)

##Convert tree to sddp graph
graph = SDDP.Graph(:0)
for i = 1: length(df.Node)
    a= df.Parent[i]
    b= df.Node[i]
    c= Float64(df.Prob[i])
    exp  = :(SDDP.add_node(graph, $(QuoteNode(i))))
    exp2 = :(SDDP.add_edge(graph, $(QuoteNode(a)) => $(QuoteNode(b)), $c))
    eval(exp)
    eval(exp2)
end
graph

## Define Reservoir data structure
struct Reservoir
    min_s::Float64              ## minimum storage volume(af)
    max_s::Float64              ## maximum storage volume(af)
    storage_initial::Float64    ## initial storage volume(af)
                                ##!!!IMPORTANT: change each month
    arc_cap::Vector{Float64}    ## arc capacity (cfs)
    max_pwr::Vector{Float64}    ## max power generation (kw)
    spill_cap::Float64          ## spill capacity (cfs)
    release_min::Float64        ## min release (cfs) 
    release_initial::Float64    ## initial release (cfs) (from 50 year historical records 1970-2019)
    k::Float64                  ## k coefficient  (kw/cfs)     avg (energy(MWH)*1000/24) / avg(release)        
    n_unit::Int64               ## number of generators 
    max_decrease::Float64       ## maximum daily descrease rate of release change (cfs)
    max_increase::Float64       ## maximum daily inscrease rate of release change (cfs)
    min_gen::Float64            ## daily mean generation follow USACE forecast information (MW) (added 8/11/2022)
                                ##!!!IMPORTANT: change each month
    evap_avg_monthly::Float64   ## monthly avg evaporation (1000 af)
end

## Define System connectivity matrix and Engineering data

    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]
    ## Decomposed matrices for today and yesterday release
    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]    

##define Cascade system
Cascade_system = [
 Reservoir(4.088, 18.463, 15.208, [8.8, 7.2, 8.8, 7.2, 7.2], [43.5, 18.25, 43.5, 40.0, 40.0], 230.0, 6.0, 6.0, 5.84, 5, 12.0, 9.0, 57.9, 109.0)
 Reservoir(4.794, 21.956, 18.195, [8.2, 8.2, 8.2, 8.2, 8.2], [121.6, 121.6, 121.6, 109.25, 109.25], 660.0, 9.0, 9.0, 12.3, 5, 12.0, 9.0, 183.8, 145.0)
 Reservoir(5.315, 21.876, 18.858015, [7.8, 7.8, 7.8, 7.8, 7.8, 7.8, 7.8], [112.29, 112.29, 112.29, 112.29, 112.29, 112.29, 112.29], 80.0, 0.0, 0.0, 12.45, 7, NaN, NaN, 108.3, 141.0)
 Reservoir(1.631, 1.749, 1.6836399999999998, [12.875, 12.875, 12.875, 12.875, 12.875, 12.875, 12.875, 12.875], [67.275, 67.275, 67.275, 67.275, 67.275, 67.275, 67.275, 67.275], 270.0, 0.0, 0.0, 4.86, 8, NaN, NaN, 7.5, 33.0)
 Reservoir(1.469, 4.307, 3.142875, [5.562, 5.562, 5.562, 5.562, 5.562, 5.562, 5.562, 5.562], [40.0, 40.0, 40.0, 40.0, 40.0, 40.0, 40.0, 40.0], 508.0, 9.0, 9.0, 8.5, 8, 17.0, 12.0, 141.2, 43.0)
 Reservoir(0.295, 0.374, 0.361845, [12.0, 12.0, 12.0], [44.1, 44.1, 44.1], 345.0, 9.0, 9.0, 3.45, 3, 25.0, 10.0, 41.2, 13.0)
    ]

##define the model
  model = SDDP.PolicyGraph(
        graph,
        sense = :Min,
        lower_bound = 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]
        inflow[r=1:6] == df[!,r+4][node]
        totaldemand == df[!,11][node]
        0 <= spill[r=1:6] <= Cascade_system[r].spill_cap
        total_gen[r=1:6] >= Cascade_system[r].min_gen/constant
        total_pflow[r=1:6] >=0
        end)

    @constraints(subproblem,begin
        sum(total_gen[r]  for r=1:6) == totaldemand#1000 or 1 for scaled
        [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]*1000 == 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/constant   
        end)
    @stageobjective(subproblem, sum(spill[r] for r=1:6))
end 

##train the model
SDDP.train(model; iteration_limit = 5)

##numercila stability report
SDDP.numerical_stability_report(

    model,
    by_node = true,
    print = true,

    warn = true,
)
odow commented 11 months ago

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.

  rhs range        [4e-05, 2e+01]
WARNING: numerical stability issues detected
  - rhs range contains small coefficients
Very large or small absolute values of coefficients
can cause numerical stability issues. Consider
reformulating the model.

Things are probably still fine.

p.s. I still don't know why you need eval to construct the graph. Take a look at how I did it above.

bolbolim commented 11 months ago

Thanks, i'll take a look at it

odow commented 11 months ago

Any updates?

bolbolim commented 11 months ago

not really, we still get warnings for some nodes. And haven't been able to solve it yet.

odow commented 11 months ago

It doesn't matter if you get some warnings. They're warnings, not errors. It's a hint, not a requirement.

bolbolim commented 11 months ago

Yes. we probably going to ignore them. Thank you again!! Your comments were essential to improvement of our model!

odow commented 11 months ago

Closing as this seems resolved. Please re-open if you have further questions :smile: