When running a stochastic programming problem in Pyomo, the resulting solution works only when running 10 precisely the same scenarios but the results remain zero when running different scenarios.
I aim to run 10 different scenarios with corresponding probabilities in a stochastic fashion. To provide some context, the problem concerns a cost-minimization battery dispatch problem with uncertainty regarding wind power production.
As said before, when running the 10 similar scenarios with corresponding probabilities 0.1 this results in the following plot. The wind scenarios are described in a [s, t] tuple-indexed dictionary. The wind parameter is indexed similarly in the objective function and constraints.
Battery dispatch with 10 similar scenarios
However, when running 10 different scenarios all with probability 0.1 the resulting plot ends up empty. The resulting plot is provided below. It says that no primal feasible solution can be found, but I'm not sure why this is the case.
Empty battery dispatch plot 10 different scenarios
Hence my question, how do I retrieve indexed results?
## Import packages
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pyomo.environ as pyo
load_data = ([5.76186621, 5.37067822, 5.29351123, 5.52940967, 5.39076172,
5.43196387, 5.48100537, 5.81333301, 5.69059522, 5.48861719,
6.69299756, 5.19573291, 5.09057129, 5.21503565, 5.21353174,
5.39826172, 5.23854395, 5.54866699, 5.22932715, 5.57071094,
5.68876563, 5.39178857, 5.37514307, 5.63663672, 5.73555957,
5.81333301, 5.32496875, 5.8409707 , 5.1021416 , 5.44398438,
5.30204785, 5.63930811, 5.50765918, 5.17532715, 5.24306397,
5.25783154, 5.2048418 , 5.91519141, 5.90498242, 5.44157129,
5.84772656, 5.39375586, 5.50630371, 5.41415381, 5.34967041,
5.58918799, 4.95105859, 5.34705859]) # Electrical load of the machine [MW]
three_days = ([ 66.75, 50.74, 50. , 50.72, 51.3 , 48.6 , 64.46, 85.66,
100. , 67.8 , 50.28, 48.12, 45.57, 40.26, 39.09, 38.94,
40.71, 44.84, 70.72, 80. , 80.1 , 78.16, 55. , 47.92,
46.65, 48.7 , 47.85, 50.08, 50. , 46. , 58. , 72.97,
75.12, 51.47, 45.91, 45.4 , 42.9 , 42.71, 41.26, 41.84,
44. , 66.36, 54.09, 69.83, 88.32, 63.6 , 58.88, 47.56]) # Electricity price data [$/MWh]
park_data = ([10.548, 15.84 , 10.548, 3.636, 6.525, 0. , 0.315, 10.548,
6.525, 0.315, 0.315, 6.525, 6.525, 10.548, 10.548, 6.525,
6.525, 10.548, 10.548, 10.548, 6.525, 6.525, 6.525, 6.525,
6.525, 3.636, 6.525, 6.525, 6.525, 6.525, 3.636, 3.636,
0.315, 0.315, 0. , 0.315, 0.315, 1.656, 3.636, 1.656,
3.636, 1.656, 0.315, 0. , 0. , 0. , 0. , 0. ]) # Wind power production [MW]
park_data_2 = [x*0.9 for x in park_data]
Scen_prob = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] # Scenario probabilities
Scen_prob = np.array(Scen_prob) # List to array
Scen_prob = Scen_prob.reshape(10,1) # List to array
# 10 similar scenarios
raw_wind_data = [park_data,
park_data,
park_data,
park_data,
park_data,
park_data,
park_data,
park_data,
park_data,
park_data]
# 10 different scenarios
# raw_wind_data = [park_data,
# park_data,
# park_data,
# park_data,
# park_data,
# park_data_2,
# park_data_2,
# park_data_2,
# park_data_2,
# park_data_2]
#Transform wind data to dataframe
raw_wind_df = pd.DataFrame(raw_wind_data)
raw_wind_df = raw_wind_df.transpose()
raw_wind_df.columns = ['Scenario 1', 'Scenario 2', 'Scenario 3', 'Scenario 4', 'Scenario 5', 'Scenario 6', 'Scenario 7', 'Scenario 8', 'Scenario 9', 'Scenario 10']
ax = raw_wind_df.plot()
plt.show()
def build_model(price_data, horizon_length, scenario_length, load_calc, park_calc, raw_wind_data):
""" Builds HORIZON pyomo model
:param raw_wind_data:
:param scenario_length:
:param price_data: dictionary of price data
:param horizon_length: number of timesteps to evaluate
:param include_pbc: Boolean, include periodic boundary condition constraint
:param load_calc: dictionary of laod data
:param park_calc: dictionary of windpark power production
:return: m: pyomo model
"""
m = pyo.ConcreteModel()
### BEGIN SOLUTION
# Probability distribution per scenario
vector = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.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]
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]
# 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=wind_1, mutable=True) #park_data_dict
# Scenario probability
m.prob = pyo.Param(m.SCENARIO, initialize=prob_dict) # Was Scen_prob
## Variables
## Battery related variables
# Charging rate [MW]
m.c = pyo.Var(m.HORIZON, initialize=0.0, bounds=(0, 10), domain=pyo.NonNegativeReals) # initialize=0.0,
# Discharging rate [MW]
m.d = pyo.Var(m.HORIZON, initialize=0.0, bounds=(0, 10), domain=pyo.NonNegativeReals) # initialize=0.0,
# Battery power
m.Bat = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.NonNegativeReals) # initialize=0.0,
# Binary variables charging and grid
m.u = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.Binary) # initialize=0.0,
m.v = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.Binary) # initialize=0.0,
# Energy (state-of-charge) [MWh]
m.E = pyo.Var(m.HORIZON, initialize=2.0, bounds=(0, 5), domain=pyo.NonNegativeReals) # initialize=2.0,
m.G_im = pyo.Var(m.HORIZON, initialize=0.0, bounds=(0, 10), domain=pyo.NonNegativeReals) # initialize=0,
m.G_ex = pyo.Var(m.HORIZON, initialize=0.0, bounds=(0, 100), domain=pyo.NonNegativeReals) # initialize=0,
m.grid = pyo.Var(m.HORIZON, initialize=m.Load, bounds=(0, 10), domain=pyo.NonNegativeReals) # initialize=m.Load,
def objfun(m):
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]
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 GridBalance(m, s, t):
return m.grid[t] == m.wind[s, t] + m.G_im[t] - m.G_ex[t] # for s in m.SCENARIO for t in m.HORIZON)
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 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.GridBalance_Con = pyo.Constraint(m.SCENARIO, m.HORIZON, rule=GridBalance)
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
def extract_solution(m):
'''
Function to extract the solution from the solved pyomo model
Inputs:
m: Solved pyomo model
Outputs:
c_control: numpy array of charging power values, MW
d_control: numpy array of discharging power values, MW
E_control: numpy array of energy stored values, MWh
t: numpy array of time indices from the Pyomo model, hr
'''
c_control = np.empty(m.N)
d_control = np.empty(m.N)
E_control = np.empty(m.N)
price_control = np.empty(m.N)
Load_control = np.empty(m.N)
Grid_control = np.empty(m.N)
G_im_control = np.empty(m.N)
G_ex_control = np.empty(m.N)
t = np.empty(m.N)
### BEGIN SOLUTION
# Loop over elements of HORIZON set.
count = 0
for i in m.HORIZON:
t[count] = pyo.value(i)
# Use value( ) function to extract the solution for each variable and append to the results lists
c_control[count] = pyo.value(m.c[i])
# Adding negative sign to discharge for plotting
d_control[count] = -pyo.value(m.d[i])
E_control[count] = pyo.value(m.E[i])
price_control[count] = pyo.value(m.price[i])
Load_control[count] = pyo.value(m.Load[i])
Grid_control[count] = pyo.value(m.grid[i])
G_im_control[count] = pyo.value(m.G_im[i])
G_ex_control[count] = pyo.value(m.G_ex[i])
count = count + 1
### END SOLUTION
## Calculate revenue
# include minus sign in order to have proper values [Charge = Neg][Discharge = Pos]
all_control_sum = np.add(-c_control, -d_control)
# original: [Charge = Pos][Discharge = Neg]
all_control_norm = np.add(c_control, d_control)
my_data_list = list(price_control)
revenue_hour = [a * b for a, b in zip(all_control_sum, my_data_list)]
revenue_hour_array = np.array(revenue_hour)
CumRev = revenue_hour_array.cumsum()
return c_control, d_control, Load_control, Grid_control, G_im_control, G_ex_control, t
def plot_solution(c_control, d_control, Load_control, Grid_control, G_im_control, G_ex_control, t=None):
plt.figure()
# Make a plot
plt.plot(t, Grid_control, 'm.-', label='Cold ironing power') # Energy level Was P_ci
plt.plot(t, G_im_control, 'k.-', label='Grid import')
plt.plot(t, G_ex_control, 'b.-', label='Grid export')
plt.plot(t, Load_control, 'g.-', label='Load Power')
# plt.plot(t, P_grid, 'm.-', label='Grid Power')
plt.plot(t, park_data, 'y.-', label='Wind Power') # was P_wind
plt.plot(t, c_control, 'c.-', label='Charge')
plt.plot(t, d_control, 'c.-', label='Discharge')
plt.title('Check: P_ci is ideally lower than P_Load')
plt.xlabel('Time (hr)')
plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15),
ncol=3, fancybox=True, shadow=True)
plt.ylabel('Cold ironing and load energy [MWh]') # Cumulative revenue [€]
plt.xticks() # 72 was len(three_days) and 10 was int(step_size)
# plt.legend()
plt.grid(True)
plt.show()
fig, ax = plt.subplots()
# Make a plot
plt.step(t, c_control, 'b.-') # Energy level
plt.step(t, d_control, 'b.-')
plt.xlabel('Time (hr)')
plt.ylabel('Energy in battery [MWh]', color="blue") # Cumulative revenue [€]
plt.xticks(range(0, 50, 10)) # 72 was len(three_days) and 10 was int(step_size)
plt.grid(True)
plt.title('Charging based on day-ahead price')
# insert the second plot
ax2 = ax.twinx()
ax2.step(t, three_days, 'r.-') # 96 was three_days
ax2.set_ylabel("Day-ahead price [EUR/MWh]", color="red")
# Make sure axis label fit
plt.tight_layout()
plt.show()
fig, ax = plt.subplots()
# Make a plot
plt.step(t, c_control, 'b.-') # Energy level
plt.step(t, d_control, 'b.-')
plt.xlabel('Time (hr)')
plt.ylabel('Energy in battery [MWh]', color="blue") # Cumulative revenue [€]
plt.xticks(range(0, 50, 10)) # 72 was len(three_days) and 10 was int(step_size)
plt.grid(True)
plt.title('Charging and Wind speeds')
# insert the second plot
ax2 = ax.twinx()
ax2.step(t, park_data, 'k.-') # 96 was three_days
ax2.set_ylabel("Wind power production [MW]", color="black")
# Make sure axis label fit
plt.tight_layout()
plt.show()
return
## BEGIN SOLUTION
# Build model
m_og = build_model(three_days, len(three_days), Scen_prob, load_data, park_data, raw_wind_data)
# Specify the solver
solver = pyo.SolverFactory('glpk')
results = solver.solve(m_og, tee=True)
# Extract solution
c_control, d_control, Load_control, Grid_control, G_im_control, G_ex_control, t = extract_solution(m_og)
plot_solution(c_control, d_control, Load_control, Grid_control, G_im_control, G_ex_control, t)
