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 pyoload_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 scenariovector = [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 timestepsm.N = horizon_lengthm.S = len(scenario_length)# Define the horizon set starting at hour 1 until horizon length +1m.HORIZON = pyo.Set(initialize=range(1, m.N + 1))# Define scenario setm.SCENARIO = pyo.Set(initialize=range(1, m.S + 1))## Parameters# Round trip efficiencym.teta = pyo.Param(initialize=0.95)# Energy [MWh] in battery at t=0m.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 powerm.d_max = pyo.Param(initialize=5)# Maximum charge powerm.c_max = pyo.Param(initialize=5)# Maximum export powerm.im_max = pyo.Param(initialize=10)# Maximum import powerm.ex_max = pyo.Param(initialize=100)## CREATE DICTS FOR DATA: Price, Load & Calc# Create empty dictionaryprice_data_dict = {}# Loop over Price data elements of numpy arrayfor i in range(0, len(price_data)):# Add element to data_dictprice_data_dict[i + 1] = price_data[i]# Create empty dictionaryload_data_dict = {}# Loop over Load data elements of numpy arrayfor i in range(0, len(load_calc)):# Add element to data_dictload_data_dict[i + 1] = load_calc[i]# Create empty dictionarypark_data_dict = {}# Loop over Wind park data elements of numpy arrayfor i in range(0, len(park_calc)):# Add element to data_dictpark_data_dict[i + 1] = park_calc[i]# Create empty dictionaryprob_dict = {}# Loop over probability data elements of numpy arrayfor i in range(0, len(vector)):# Add element to prob_dictprob_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 datam.price = pyo.Param(m.HORIZON, initialize=price_data_dict, domain=pyo.Reals, mutable=True)# Load datam.Load = pyo.Param(m.HORIZON, initialize=load_data_dict, domain=pyo.Reals, mutable=True)# Wind park datam.wind = pyo.Param(m.SCENARIO, m.HORIZON, initialize=wind_1, mutable=True) #park_data_dict# Scenario probabilitym.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 powerm.Bat = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.NonNegativeReals) # initialize=0.0,# Binary variables charging and gridm.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 timestepif t == 1:return m.E[t] == m.E0 + m.c[t] * m.teta - m.d[t] / m.teta# Subsequent timestepselse:return m.E[t] == m.E[t - 1] + m.c[t] * m.teta - m.d[t] / m.tetadef 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 <= 0def ExMax(m, t):return m.G_im[t] + m.v[t] * m.im_max <= m.im_maxdef ChargeMax(m, t):return m.d[t] - m.u[t] * m.d_max <= 0def DischargeMax(m, t):return m.c[t] + m.u[t] * m.c_max <= m.c_maxm.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 SOLUTIONreturn mdef extract_solution(m):'''Function to extract the solution from the solved pyomo modelInputs:m: Solved pyomo modelOutputs:c_control: numpy array of charging power values, MWd_control: numpy array of discharging power values, MWE_control: numpy array of energy stored values, MWht: 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 = 0for i in m.HORIZON:t[count] = pyo.value(i)# Use value( ) function to extract the solution for each variable and append to the results listsc_control[count] = pyo.value(m.c[i])# Adding negative sign to discharge for plottingd_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, tdef plot_solution(c_control, d_control, Load_control, Grid_control, G_im_control, G_ex_control, t=None):plt.figure()# Make a plotplt.plot(t, Grid_control, 'm.-', label='Cold ironing power') # Energy level Was P_ciplt.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_windplt.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 plotplt.step(t, c_control, 'b.-') # Energy levelplt.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 plotax2 = ax.twinx()ax2.step(t, three_days, 'r.-') # 96 was three_daysax2.set_ylabel("Day-ahead price [EUR/MWh]", color="red")# Make sure axis label fitplt.tight_layout()plt.show()fig, ax = plt.subplots()# Make a plotplt.step(t, c_control, 'b.-') # Energy levelplt.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 plotax2 = ax.twinx()ax2.step(t, park_data, 'k.-') # 96 was three_daysax2.set_ylabel("Wind power production [MW]", color="black")# Make sure axis label fitplt.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)