I am trying to understand how can model if-else conditions in the objective function. I know how to do it in a constraint using binary variable but I am not sure how to do that in the objective function itself. For example, I have created a very simple energy supply model, here are its configurations:
- Two generators gen1 and gen2 which will get revenue if they supply continuous energy according to PPA.
- If there is some shortfall then there is some penalty associated with it.
- if there is an excess supply then there is revenue but with 50% of PPA only.
- All these calculations are happening monthly basis for every hour( For example if the model is running for 2 months then there will be a total of 60 Days * 24 Hours = 1440 Hours)
Here is my approach:
import numpy as np
from pyomo.environ import *
no_of_days = 60
hours = 24
total_instance = no_of_days * hours
np.random.seed(total_instance)
load_profile = np.array([100] * total_instance)
DFR = 0.9
interval_freq = 60
gen1_capacity = 200
gen2_capacity = 310
PPA = 10
shortfall_penalty_rate = 15
excess_revenue_rate = 0.5 * PPA
model = ConcreteModel()
model.month_index = Set(initialize=list(range(1, 3)))
model.m_index = Set(initialize=list(range(len(load_profile))))
# variable
model.grid = Var(model.month_index, model.m_index, domain=NonNegativeReals)
# gen variable
model.gen1_use = Var(model.month_index, model.m_index, domain=NonNegativeReals)
model.gen2_use = Var(model.month_index, model.m_index, domain=NonNegativeReals)
# Load profile
model.load_profile = Param(model.month_index, model.m_index, initialize=lambda model, month, m: load_profile[month] * interval_freq / 60.0)
model.lost_load = Var(model.month_index, model.m_index, domain=NonNegativeReals)
# Objective function
def revenue(model):
monthly_demand = sum(
model.load_profile[month, m] for month in model.month_index for m in model.m_index
)
monthly_energy = sum(
model.gen1_use[month, m] + model.gen2_use[month, m]
for month in model.month_index for m in model.m_index
)
# Revenue Calculations
monthly_revenue = min(monthly_energy, monthly_demand) * PPA
# Excess Revenue Calculations
if monthly_energy > monthly_demand:
excess_revenue = (monthly_energy - monthly_demand) * excess_revenue_rate
else:
excess_revenue = 0
# Shortfall Penalty Calculations
actual_DFR = monthly_energy / monthly_demand
if actual_DFR < DFR:
shortfall_penalty = (DFR * monthly_demand - monthly_energy) * -shortfall_penalty_rate
else:
shortfall_penalty = 0
total_cost = monthly_revenue + excess_revenue + shortfall_penalty
return total_cost
model.obj = Objective(rule=revenue, sense=maximize)
def energy_balance(model, month, m):
return model.grid[month, m] == model.gen1_use[month, m] + model.gen2_use[month, m] + model.lost_load[month, m]
model.energy_balance = Constraint(model.month_index, model.m_index, rule=energy_balance)
def grid_limit(model, month, m):
return model.grid[month, m] >= model.load_profile[m]
model.grid_limit = Constraint(model.month_index, model.m_index, rule=grid_limit)
def max_gen1(model, month, m):
eq = model.gen1_use[month, m] <= gen1_capacity
return eq
model.max_gen1 = Constraint(model.month_index, model.m_index, rule=max_gen1)
def max_gen2(model, month, m):
eq = model.gen2_use[month, m] <= gen2_capacity
return eq
model.max_gen2 = Constraint(model.month_index, model.m_index, rule=max_gen2)
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 = []
gen1 = []
gen2 = []
lost_load = []
load = []
for month in range(1, 13):
for i in range(len(load_profile)):
grid_use.append(value(model.grid[month, i]) * (60 / interval_freq))
gen1.append(value(model.gen1_use[month, i]) * (60 / interval_freq))
gen2.append(value(model.gen2_use[month, i]) * (60 / interval_freq))
lost_load.append(value(model.lost_load[month, i]) * (60 / interval_freq))
load.append(value(model.load_profile[month, i]) * (60 / interval_freq))
I am getting error in the objective function itself:
if monthly_energy > monthly_demand:
File "C:\Users\nvats\AppData\Local\Continuum\anaconda3\envs\microgrid\lib\site-packages\pyomo\core\expr\numvalue.py", line 730, in __gt__
return _generate_relational_expression(_lt, other, self)
File "C:\Users\nvats\AppData\Local\Continuum\anaconda3\envs\microgrid\lib\site-packages\pyomo\core\expr\logical_expr.py", line 380, in _generate_relational_expression
"%s < %s" % ( prevExpr.to_string(), lhs, rhs ))
TypeError: Cannot create a compound inequality with identical upper and lower
bounds using strict inequalities: constraint infeasible:
288000.0 < gen1_use[1,0] + gen2_use[1,0] + gen1_use[1,1] + gen2_use[1,1] + gen1_use[1,2] + gen2_use[1,2] + gen1_use[1,3] + gen2_use[1,3] + gen1_use[1,4] + gen2_use[1,4] + .... so on
Can someone please help and explain what's wrong with the objective function and how can i correct it?