Indexing problem SDDP package

Specific Domains Optimization (Mathematical)
1

I am trying to define a hydropower production problem in Julia using the SDDP package.

When running the code I get this error: “ERROR: KeyError: key Symbol(“w2[3]”) not found”. This error seems like it is being caused by the water balance constraint for reservoir 2. I am a little bit unsure of how the indexing in each subproblem is working. I thought each of the variables in one subproblem had to be initialized by the number of scenarios in that time step, but it seems like I am mistaken. I have also been experimenting with using the time and scenario variables extracted from the index variable without any luck.

I would appreciate any guidance in this matter!

The code I have written is the following:


# Read lattice datafiles
P = Any[Any[47.16], Any[50.783602, 50.783602, 61.125716, 61.125716], Any[62.898253, 62.898253, 50.222343, 50.222343], Any[64.728229, 64.728229, 49.976048, 49.976048], Any[49.512781, 49.512781, 65.75061, 65.75061]]                 # Prices
I = Any[Any[1.1788416e7], Any[3.3384999e7, 7.997725e6, 3.3384999e7, 7.997725e6], Any[4.9624374e7, 9.144343e6, 4.9624374e7, 9.144343e6], Any[5.9080642e7, 8.643971e6, 5.9080642e7, 8.643971e6], Any[3.6434978e7, 6.237361e6, 3.6434978e7, 6.237361e6]]                # Inflow
probs = [[0.031660448 0.559019552 0.021939552 0.38738044800000004], [0.07481887811568554 0.32106383333552324 0.11417365919774729 0.48994362935104396; 0.009880335634486002 0.38600237581672275 0.015077398938633187 0.589039889610158; 0.07439745807024682 0.3192554296459654 0.114595079243186 0.4917520330406018; 0.009824684285563117 0.3838282034306491 0.015133050287556075 0.5912140619962317], [0.0447310109816123 0.3413742628366524 0.07112084087023957 0.5427738853114957; 0.004371530862170055 0.38173374295609463 0.006950590741969661 0.6069441354397657; 0.04459979394961132 0.3403728520348853 0.07125205790224054 0.5437752961132628; 0.0043587071121039 0.38061393887239275 0.006963414492035815 0.6080639395234675], [0.17726648108050216 0.4555237218150404 0.10286819905417799 0.2643415980502795; 0.020548097843605957 0.6122421050519367 0.011924114509842761 0.35528568259461474; 0.17799569233068535 0.4573975843786121 0.1021389878039948 0.2624677354867078; 0.020632625409257907 0.6147606513000395 0.01183958694419081 0.3527671363465118]]       # Probabilities

# Define parameters
EC = 0.61748                    # energy coefficient
alpha = 0.02                    # discount rate
sigma = 604800                  # number of seconds the turbine is running per time period
phi = 0.5                       # percentage of inflow flowing into reservoir 2
W_V_bar = 44500000              # maximum capacity of reservoir 2
W_S_bar = 22500000              # maximum capacity of reservoir 1
X_bar = 10281600                # maximum capacity of turbine


W_S_underline = [15050000, 15050000, 0, 0, 0]

num_time_pers = 5           # The number of time periods
num_scen = 4                # The number of scenarios


# Create an array indicating the number of scenarios for each time period
scens = [1; fill(num_scen, num_time_pers - 1)]

graph = SDDP.MarkovianGraph(probs)


model = SDDP.PolicyGraph(
        graph;
        sense = :Max,
        lower_bound = 0.0,
        upper_bound = 1e100,
        optimizer = HiGHS.Optimizer,
    ) do subproblem, index
    # Unpack the time stage and Markov index (scenarios)
    time, scenario = index

    # State variables
    @variable(subproblem, W_S_underline[time] <= w1[i = 1:scens[time]] <= W_S_bar, SDDP.State, initial_value = 20000000)         # Water balance constraint reservoir 1  
    @variable(subproblem, 0 <= w2[i = 1:scens[time]] <= W_V_bar, SDDP.State, initial_value = 10000000)                           # Water balance constraint reservoir 2

    # Control variables
    @variable(subproblem, x[i = 1:scens[time]] >= 0)     # Variables for production volumes
    @variable(subproblem, s[i = 1:scens[time]] >= 0)     # Variables for spillage
    @variable(subproblem, fI[i = 1:scens[time]] >= 0)    # Variables for water flowing from reservoir 2 to reservoir 1
    @variable(subproblem, fO[i = 1:scens[time]] >= 0)    # Variables for water flowing from reservoir 1 to reservoir 2

    # Constraints
    @constraint(subproblem, [i = 1:scens[time]], w1[i].out == w1[i].in + (1-phi) * I[time][i] - fO[i] + fI[i])                   # Balance constraint for reservoir 1
                            
    @constraint(subproblem, [i = 1:scens[time]], w2[i].out == w2[i].in - x[i] - s[i] + phi * I[time][i] + fO[i] - fI[i])         # Balance constraint for reservoir 2

    @constraint(subproblem, [i = 1:scens[time]], x[i] <= X_bar)                                                            # Maximum capacity constraint turbine

    @stageobjective(subproblem, sum(EC * P[time][scen] * x[scen] * exp(-alpha*time) for scen in 1:scens[time]))                   # Objective for subproblem

end

SDDP.train(model) ```
2

Hi @heiwie,

You should start by taking a looking at the introductory tutorials for SDDP.jl:

There are also more complicated examples like:

You need to build the subproblem at each node. It doesn’t have variables within a subproblem indexed by scenario, etc.

You also have a major issue: the data you are using is too large. Solvers cannot optimize with numbers like 1e100.

See Words of warning · SDDP.jl

Here’s how I would write your model, after rescaling the inflow and reservoir data by 1e-6:

using SDDP, HiGHS

P = [
    [47.16],
    [50.783602, 50.783602, 61.125716, 61.125716],
    [62.898253, 62.898253, 50.222343, 50.222343],
    [64.728229, 64.728229, 49.976048, 49.976048],
    [49.512781, 49.512781, 65.75061, 65.75061],
]
I = [
    [1.1788416e1],
    [3.3384999e1, 7.997725, 3.3384999e1, 7.997725],
    [4.9624374e1, 9.144343, 4.9624374e1, 9.144343],
    [5.9080642e1, 8.643971, 5.9080642e1, 8.643971],
    [3.6434978e1, 6.237361, 3.6434978e1, 6.237361],
]
probs = [
    [1.0;;],
    [0.031660448 0.559019552 0.021939552 0.38738044800000004],
    [0.07481887811568554 0.32106383333552324 0.11417365919774729 0.48994362935104396; 0.009880335634486002 0.38600237581672275 0.015077398938633187 0.589039889610158; 0.07439745807024682 0.3192554296459654 0.114595079243186 0.4917520330406018; 0.009824684285563117 0.3838282034306491 0.015133050287556075 0.5912140619962317],
    [0.0447310109816123 0.3413742628366524 0.07112084087023957 0.5427738853114957; 0.004371530862170055 0.38173374295609463 0.006950590741969661 0.6069441354397657; 0.04459979394961132 0.3403728520348853 0.07125205790224054 0.5437752961132628; 0.0043587071121039 0.38061393887239275 0.006963414492035815 0.6080639395234675],
    [0.17726648108050216 0.4555237218150404 0.10286819905417799 0.2643415980502795; 0.020548097843605957 0.6122421050519367 0.011924114509842761 0.35528568259461474; 0.17799569233068535 0.4573975843786121 0.1021389878039948 0.2624677354867078; 0.020632625409257907 0.6147606513000395 0.01183958694419081 0.3527671363465118],
]
EC = 0.61748
alpha = 0.02
phi = 0.5
W_V_bar = 44.5
W_S_bar = 22.5
X_bar = 10.2816
W_S_underline = [15.05, 15.05, 0, 0, 0]

model = SDDP.MarkovianPolicyGraph(;
    transition_matrices = probs,
    sense = :Max,
    upper_bound = 1e6,
    optimizer = HiGHS.Optimizer,
) do subproblem, node
    t, i = node
    @variables(subproblem, begin
        W_S_underline[t] <= w1 <= W_S_bar, SDDP.State, (initial_value = 20)
        0 <= w2 <= W_V_bar, SDDP.State, (initial_value = 10)
        0 <= x <= X_bar
        s >= 0
        fI >= 0
        fO >= 0
    end)
    @constraints(subproblem, begin
        w1.out == w1.in + (1 - phi) * I[t][i] - fO + fI
        w2.out == w2.in - x - s + phi * I[t][i] + fO - fI
    end)
    @stageobjective(subproblem, EC * P[t][i] * x * exp(-alpha * t))
end

SDDP.train(model; log_every_iteration = true)
1 Like
3

Thank you, it seems to work fine now! I had to make a slight modification to the probs matrix, as it was being created as a vector instead of a matrix due to the first row and therefore crashed the code.

1 Like
4

What version of Julia are you using?

I think the [1.0 ;;] syntax was added in Julia 1.7 or 1.8. if you’re using 1.6 you’ll need reshape([1.0], 1, 1).

1 Like
5

Yes exactly, using 1.6 but managed to figure out the reshape!

1 Like