How to properly index scenarios in a two stage stochastic programming in Pyomo?

Question

Main question: How to define the wind scenario parameter so that it can provide different wind scenarios for different probabilities in scenario set 'S'.

When running the stochastic code with scenario set S = [1,2,3,...10] and the corresponding probabilities Prob = [0.1, 0.1, ..., 0.1]. With t = [1, 2, ..., 48] with corresponding windspeeds m.wind.

I tried the following options to create a (s*t) matrix for the objective function below.

Objective function

Option 1: m.wind = dict{:48}, this is the dict with the wind speeds at 48 timesteps used in the deterministic model

Option 2: m. wind = {ndarray:(10,48)} = array([[...],..,[...]]), In this I constructed an array with ten scenarios consisting of wind speeds at 48 timesteps for 10 scenarios.

Option 3: m.windd = {list:10} = array[(...),...,(...)], In this option I put the different scenarios between parentheses after having read that Pyomo sometimes doesn't recognize square brackets.

Option 4: The last method for constructing a scenario set for wind was by creating a (x*t) dictionary.

All options resulted in the error: ERROR: index '0' not valid for indeed component 'wind'

Do you know how this error can be resolved and how the wind should be properly indexed?

`

def build_model(price_data, horizon_length, scenario_length, load_calc, park_calc):

m = pyo.ConcreteModel()

### BEGIN SOLUTION

# test vector
vector = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
vector = vector.reshape(10,1)

## Sets

# Save the number of timesteps
m.N = horizon_length
m.S = len(scenario_length)

# Define the horizon set starting at hour 1 until horizon length +1
m.HORIZON = pyo.Set(initialize=range(1, m.N + 1))

# Define scenario set
m.SCENARIO = pyo.Set(initialize=range(1, m.S + 1))

## Parameters

# Round trip efficiency
m.teta = pyo.Param(initialize=0.95)

# Energy [MWh] in battery at t=0
m.E0 = pyo.Param(initialize=2.0, mutable=True)

# Guarantee of origin for local wind [€/MWh]
m.goNL = pyo.Param(initialize=5)

# Guarantee of origin for grid power [€/MWh]
m.goBE = pyo.Param(initialize=150)

# Maximum discharge power
m.d_max = pyo.Param(initialize=5)

# Maximum charge power
m.c_max = pyo.Param(initialize=5)

# Maximum export power
m.im_max = pyo.Param(initialize=10)

# Maximum import power
m.ex_max = pyo.Param(initialize=100)

## CREATE DICTS FOR DATA: Price, Load & Calc
# Create empty dictionary
price_data_dict = {}

# Loop over Price data elements of numpy array
for i in range(0, len(price_data)):
    # Add element to data_dict
    price_data_dict[i + 1] = price_data[i]

# Create empty dictionary
load_data_dict = {}

# Loop over Load data elements of numpy array
for i in range(0, len(load_calc)):
    # Add element to data_dict
    load_data_dict[i + 1] = load_calc[i]

# Create empty dictionary
park_data_dict = {}

# Loop over Wind park data elements of numpy array
for i in range(0, len(park_calc)):
    # Add element to data_dict
    park_data_dict[i + 1] = park_calc[i]

# Create empty dictionary
prob_dict = {}

# Loop over probability data elements of numpy array
for i in range(0, len(vector)):
    # Add element to prob_dict
    prob_dict[i + 1] = vector[i]

# Repeat the wind data to a matrix for 10 similar scenario's
wind_matrix = np.tile(park_calc, (10, 1))
# wind_matrix = np.tile(park_calc, (10, 1)) * vector

park_data_dict_2 = {1: park_data_dict, 2: park_data_dict, 3: park_data_dict, 4: park_data_dict, 5: park_data_dict,
                    6: park_data_dict, 7: park_data_dict, 8: park_data_dict, 9: park_data_dict, 10: park_data_dict}

# Price data
m.price = pyo.Param(m.HORIZON, initialize=price_data_dict, domain=pyo.Reals, mutable=True)

# Load data
m.Load = pyo.Param(m.HORIZON, initialize=load_data_dict, domain=pyo.Reals, mutable=True)

# Wind park data
m.wind = pyo.Param(m.SCENARIO, m.HORIZON, initialize=park_data_dict_2, mutable=True) #park_data_dict

# Scenario probability
m.prob = pyo.Param(m.SCENARIO, initialize=vector)  # Was Scen_prob

# # New description of wind in 10 different scenarios
# m.wind = pyo.Param(m.SCENARIO, m.HORIZON, initialize=wind_matrix_2) #  initialize=wind_matrix_2

## Variables

## Battery related variables
# Charging rate [MW]
m.c = pyo.Var(m.HORIZON, initialize=0.0, bounds=(0, 10), domain=pyo.NonNegativeReals)

# Discharging rate [MW]
m.d = pyo.Var(m.HORIZON, initialize=0.0, bounds=(0, 10), domain=pyo.NonNegativeReals)

# Battery power
m.Bat = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.NonNegativeReals)

# Binary variables charging and grid
m.u = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.Binary)
m.v = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.Binary)

# Energy (state-of-charge) [MWh]
m.E = pyo.Var(m.HORIZON, initialize=2.0, bounds=(0, 5), domain=pyo.NonNegativeReals)
m.G_im = pyo.Var(m.HORIZON, initialize=0, bounds=(0, 10), domain=pyo.NonNegativeReals)
m.G_ex = pyo.Var(m.HORIZON, initialize=0, bounds=(0, 100), domain=pyo.NonNegativeReals)
m.grid = pyo.Var(m.HORIZON, initialize=m.Load, bounds=(0, 10), domain=pyo.NonNegativeReals)

# Objective function

# def objfun(model):
#     return sum((m.price[t] + m.goNL) * m.wind[t] + (m.price[t] + m.goBE) * m.G_im[t] for t in m.HORIZON)

def objfun(model):
    return sum((m.price[t] + m.goBE) * m.G_im[t] + (m.price[t] + m.goNL) * sum(m.prob[s] * m.wind[s, t] for s in m.SCENARIO) for t in m.HORIZON)

m.OBJ = pyo.Objective(rule=objfun, sense=pyo.minimize)

def PowerBalance(m, t):
    return m.Load[t] + m.c[t] == m.grid[t] + m.d[t]

# Define Energy Balance constraints. [MWh] = [MW]*[1 hr]
# Note: assume 1-hour timestep in price data and control actions.
def EnergyBalance(m, t):
    # First timestep
    if t == 1:
        return m.E[t] == m.E0 + m.c[t] * m.teta - m.d[t] / m.teta

        # Subsequent timesteps
    else:
        return m.E[t] == m.E[t - 1] + m.c[t] * m.teta - m.d[t] / m.teta

# def ColdIroning(m, t):
#     return m.c[t] + m.d[t] + m.Load[t] <= m.CI

def GridBalance(m, t, s):
    return m.grid[t] == m.wind[t, s] + m.G_im[t] - m.G_ex[t]

def ImMax(m, t):
    return m.G_ex[t] - m.v[t] * m.ex_max <= 0

def ExMax(m, t):
    return m.G_im[t] + m.v[t] * m.im_max <= m.im_max

# def BatteryBalance(m, t):
#     return m.Bat[t] - m.d[t] + m.c[t] == 0
#
def ChargeMax(m, t):
    return m.d[t] - m.u[t] * m.d_max <= 0

def DischargeMax(m, t):
    return m.c[t] + m.u[t] * m.c_max <= m.c_max

m.EnergyBalance_Con = pyo.Constraint(m.HORIZON, rule=EnergyBalance)
m.PowerBalance_Con = pyo.Constraint(m.HORIZON, rule=PowerBalance)
# m.ColdIroning_Con = pyo.Constraint(m.HORIZON, rule=ColdIroning)
m.GridBalance_Con = pyo.Constraint(m.HORIZON, m.SCENARIO, rule=GridBalance)
# m.BatteryBalance_Con = pyo.Constraint(m.HORIZON, rule=BatteryBalance)
m.ChargeMax_Con = pyo.Constraint(m.HORIZON, rule=ChargeMax)
m.DischargeMax_Con = pyo.Constraint(m.HORIZON, rule=DischargeMax)
m.ImMax_Con = pyo.Constraint(m.HORIZON, rule=ImMax)
m.ExMax_Con = pyo.Constraint(m.HORIZON, rule=ExMax)
## END SOLUTION

return m`

Answer 1

Two quick things on your model:

You don't need to (and probably shouldn't) initialize variables, just let the solver do its work.

You have an egregious typo wind[t, s] (instead of wind[s, t]) in your GridBalance constraint, which would be a devil to find if |T| == |S|.

You didn't say what format the wind data came to you in. Perhaps you are just hand-jamming in some table data. So let's start with what pyomo wants. It wants a tuple-indexed dictionary to initialize the parameter. Meaning that the key values of the dictionary are tuples of (s, t) values. This is also known as a "flat" data structure where all of the keys are enumerated and the data value of interest is in 1 column (vice a matrix format or something).

So you want to initialize your parameter from something like this:

import pyomo.environ as pyo

# what we want: a "flat" dictionary

#         s, t : w
wind = { (1, 1): 12,
         (1, 2): 11,
         (1, 3): 10,
         (2, 1): 9,
         (2, 2): 13,
         (2, 3): 14}

m = pyo.ConcreteModel()

# SETS
m.S = pyo.Set(initialize=[1, 2])
m.T = pyo.Set(initialize=[1, 2, 3])

# PARAMS

m.wind = pyo.Param(m.S, m.T, initialize=wind)

m.pprint()

Output:

3 Set Declarations
    S : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}
    T : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {1, 2, 3}
    wind_index : Size=1, Index=None, Ordered=True
        Key  : Dimen : Domain : Size : Members
        None :     2 :    S*T :    6 : {(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3)}

1 Param Declarations
    wind : Size=6, Index=wind_index, Domain=Any, Default=None, Mutable=False
        Key    : Value
        (1, 1) :    12
        (1, 2) :    11
        (1, 3) :    10
        (2, 1) :     9
        (2, 2) :    13
        (2, 3) :    14

4 Declarations: S T wind_index wind

There are several techniques to make a tuple-indexed dictionary from data, dependent on the structure of the source data, obviously. Before that though, I see you converting lists to dicts with loops. Certainly doable, or you could use a shortcut with enumerate, which generates index:value pairs and just pass that to a dictionary constructor. Note that enumerate takes an optional start argument, if you like your data 1-indexed vice 0.

prices = [3.5, 4.2, 9.8]
price_dict = dict(enumerate(prices, start=1))

print(price_dict)

# {1: 3.5, 2: 4.2, 3: 9.8}

So if you are purely hand-jamming in the values, you could type out a flat dictionary as shown above, or if you have a list-of-lists (aka matrix) of wind data you can convert it several ways, depending on your comfort level with dictionary comprehensions and such. All 3 of these below generate the same flat dictionary usable in your model:

raw_wind_data = [[4, 5, 9],
                 [3, 0, 12]]

wind_1 = {}
for i in range(len(raw_wind_data)):
    for j in range(len(raw_wind_data[0])):
        wind_1[(i+1, j+1)] = raw_wind_data[i][j]

wind_2 = { (r+1, c+1) : raw_wind_data[r][c] 
                    for r in range(len(raw_wind_data)) 
                    for c in range(len(raw_wind_data[0]))}

wind_3 = {(r_idx, c_idx): w 
            for r_idx, row in enumerate(raw_wind_data, 1)
            for c_idx, w   in enumerate(row, 1)}

print(wind_1)
print(wind_2)
print(wind_3)

# {(1, 1): 4, (1, 2): 5, (1, 3): 9, (2, 1): 3, (2, 2): 0, (2, 3): 12}
# {(1, 1): 4, (1, 2): 5, (1, 3): 9, (2, 1): 3, (2, 2): 0, (2, 3): 12}
# {(1, 1): 4, (1, 2): 5, (1, 3): 9, (2, 1): 3, (2, 2): 0, (2, 3): 12}