how to model if else condition in the objective function in PYOMO?

46 Views Asked by At

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:

  1. Two generators gen1 and gen2 which will get revenue if they supply continuous energy according to PPA.
  2. If there is some shortfall then there is some penalty associated with it.
  3. if there is an excess supply then there is revenue but with 50% of PPA only.
  4. 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?

0

There are 0 best solutions below