I have written a program to maximize the revenue, where there is a facility that contains Solar and Wind power plants. They serve a demand profile. To serve that demand profile they get some revenue based on PPA(Power Purchasing Agreement). There are three components which affect the revenue:
Total Revenue = Actual Revenue + Excess Revenue - Shortfall Penalty The program that I wrote for this:
import numpy as np
import pandas as pd
from pyomo.environ import *
total_instances = 8760
np.random.seed(total_instances)
load_profile = np.random.randint(80, 130, total_instances)
solar_profile = np.random.uniform(0, 0.9, total_instances)
wind_profile = np.random.uniform(0, 0.9, total_instances)
month_idx = pd.date_range('1/1/2023', periods=total_instances, freq='60min').month.tolist()
PPA = 10
shortfall_penalty_rate = 15
excess_rate = 5
DFR = 0.9
solar_capacity = 120
wind_capacity = 130
# Not in Use right now
load_profile_month = list(zip(month_idx, load_profile))
monthly_sum = [0] * 12
for i in range(len(load_profile_month)):
month_idx, load_value = load_profile_month[i]
monthly_sum[month_idx - 1] += load_value
model = ConcreteModel()
model.m_index = Set(initialize=list(range(len(load_profile))))
# variable
model.facility_output = Var(model.m_index, domain=NonNegativeReals)
# gen variable
model.solar_use = Var(model.m_index, domain=NonNegativeReals)
model.wind_use = Var(model.m_index, domain=NonNegativeReals)
# Load profile
model.load_profile = Param(model.m_index, initialize=load_profile)
model.solar_profile = Param(model.m_index, initialize=solar_profile)
model.wind_profile = Param(model.m_index, initialize=wind_profile)
model.lost_load = Var(model.m_index, domain=NonNegativeReals)
# Objective function
def revenue(model):
actual_revenue = sum(
min((model.solar_use[m] + model.wind_use[m]), model.load_profile[m]) * PPA
for m in model.m_index
)
excess_revenue = sum(
(max(model.solar_use[m] + model.wind_use[m], model.load_profile[m]) - model.load_profile[m]) * excess_rate
for m in model.m_index
)
shortfall_penalty = sum(
(min(model.solar_use[m] + model.wind_use[m], DFR * model.load_profile[m]) - model.load_profile[m] * DFR) * shortfall_penalty_rate
for m in model.m_index
)
total_revenue = actual_revenue + excess_revenue + shortfall_penalty
return total_revenue
model.obj = Objective(rule=revenue, sense=maximize)
def energy_balance(model, m):
return model.grid[m] == model.solar_use[m] + model.wind_use[m] + model.lost_load[m]
model.energy_balance = Constraint(model.m_index, rule=energy_balance)
def grid_limit(model, m):
return model.grid[m] >= model.load_profile[m]
model.grid_limit = Constraint(model.m_index, rule=grid_limit)
def max_solar(model, m):
eq = model.solar_use[m] <= solar_capacity * model.solar_profile[m]
return eq
model.max_solar = Constraint(model.m_index, rule=max_solar)
def max_wind(model, m):
eq = model.wind_use[m] <= wind_capacity * model.wind_profile[m]
return eq
model.max_wind = Constraint(model.m_index, rule=max_wind)
Solver = SolverFactory('gurobi')
Solver.options['LogFile'] = "gurobiLog"
# Solver.options['MIPGap'] = 0.0
print('\nConnecting to Gurobi Server...')
results = Solver.solve(model)
if (results.solver.status == SolverStatus.ok):
if (results.solver.termination_condition == TerminationCondition.optimal):
print("\n\n***Optimal solution found***")
print('obj returned:', round(value(model.obj), 2))
else:
print("\n\n***No optimal solution found***")
if (results.solver.termination_condition == TerminationCondition.infeasible):
print("Infeasible solution")
exit()
else:
print("\n\n***Solver terminated abnormally***")
exit()
grid_use = []
solar = []
wind = []
# e_in = []
# e_out = []
# soc = []
lost_load = []
load = []
for i in range(len(load_profile)):
grid_use.append(value(model.grid[i]))
solar.append(value(model.solar_use[i]))
wind.append(value(model.wind_use[i]))
lost_load.append(value(model.lost_load[i]))
load.append(value(model.load_profile[i]))
# e_in.append(value(model.e_in[i]))
# e_out.append(value(model.e_out[i]))
# soc.append(value(model.soc[i]))
pd.DataFrame({
'Grid': grid_use,
'Solar': solar,
'Wind': wind,
'Shortfall': lost_load,
'Load Profile': load,
}).to_excel('testing4.xlsx')
The error that I am getting is I think because I am using Load profile which is Param with the Variables Solar and wind:
ValueError: Error retrieving immutable Param value (load_profile[0]):
The Param value is undefined and no default value is specified.
ERROR: Rule failed when generating expression for objective obj: ValueError:
Error retrieving immutable Param value (load_profile[0]):
The Param value is undefined and no default value is specified.
ERROR: Constructing component 'obj' from data=None failed: ValueError: Error
retrieving immutable Param value (load_profile[0]):
The Param value is undefined and no default value is specified.
Process finished with exit code 1
All of the data items that you are using to initialize your parameters should be dictionaries (or other data types like panda series) with key:value pairs. Your parameter is indexed so the data that you are feeding it needs to be keyed by the same values
favorability = {'dog': 10, 'cat': 2, 'pig': 4}
random_numbers = { 1: 22.5, 2: 52.9, 3: 4.1}
in the model, something like:
m.animals = m.Set(initialize=favorability.keys())
m.case = m.Set(initialize=random_numbers.keys())
m.favorability = m.Param(m.animals, initialize=favorability)
m.cost = m.Param(m.case, initialize=random_numbers)