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.
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`
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()
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}