pythonpyomo

How to add param values along with variable values in the objective function using PYOMO?


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:

  1. Actual Revenue: Min( Facility Generation(Solar + Wind), Demand Profile) * PPA
  2. Excess Revenue: (Max(Facility Generation, Demand Profile) - Demand Profile) * 50% of PPA
  3. Shortfall Penalty (If Facility output is less than 90%(DFR) of Demand profile then this will apply): (Min(Facility output, DFR * Demand Profile) - Demand Profile * Shortfall Penalty

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

Solution

  • 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)