How to run yearly optimization problem in monthly format?

Question

I have created a RTC (Round the clock) model using PYOMO with solar, wind, and battery as the main component and lost_load(Shortfall). It's a simple model where these three main sources of energy try to fulfill the demand of 100MW throughout the year for which they get revenue according to PPA ( which is 10 INR/kWh). And if they fail to do so there is a penalty (5 INR/kWh). This model is running on 15min intervals for a year, so there are a total of 35040 instances. I want to run this model monthly. It should be something like the model should try to fulfill monthly demand and if the model is not able to fulfill 90% of the total monthly demand then a shortfall penalty would get applied to it. My Approach :

import datetime
import pandas as pd
import numpy as np
from pyomo.environ import *


total_instances = 35040
np.random.seed(total_instances)
load_profile = np.array([100]*total_instances)
solar_profile = np.random.rand(total_instances)
wind_profile = np.random.rand(total_instances)

DFR = 0.9
interval_freq = 15
solar_capacity = 120
wind_capacity = 140
power = 50
energy = 200
round_trip_eff = 0.8
PPA = 10
shortfall_penalty = 15
excess_revenue_rate = 0.5 * PPA

model = ConcreteModel()
model.m_index = Set(initialize=list(range(len(load_profile))))

# variable
model.grid = Var(model.m_index, domain=NonNegativeReals)

# storage
# variable
model.storage_power = Param(initialize=power)
model.storage_energy = Param(initialize=energy)
model.chargeEff = Param(initialize=np.sqrt(round_trip_eff))
model.dchargeEff = Param(initialize=np.sqrt(round_trip_eff))
model.e_in = Var(model.m_index, domain=NonNegativeReals)
model.e_out = Var(model.m_index, domain=NonNegativeReals)
model.soc = Var(model.m_index, domain=NonNegativeReals)
model.ein_eout_lmt = Var(model.m_index, domain=Binary)

#  Solar variable
model.solar_capacity = Param(initialize=solar_capacity)
model.solar_use = Var(model.m_index, domain=NonNegativeReals)

#  Wind variables
model.wind_capacity = Param(initialize=wind_capacity)
model.wind_use = Var(model.m_index, domain=NonNegativeReals)

# Solar profile
model.solar_avl = Param(model.m_index, initialize=dict(zip(model.m_index, solar_profile * solar_capacity * interval_freq / 60.0)))

# Wind profile
model.wind_avl = Param(model.m_index, initialize=dict(zip(model.m_index, wind_profile * wind_capacity * interval_freq / 60.0)))

# Load profile
model.load_profile = Param(model.m_index, initialize=dict(zip(model.m_index, load_profile * interval_freq / 60.0)))

model.lost_load = Var(model.m_index, domain=NonNegativeReals)


# Objective function
def revenue(model):
    total_revenue = sum(
        (model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m]) * PPA
        for m in model.m_index
    )
    total_penalty = sum(
        model.lost_load[m] * -shortfall_penalty
        for m in model.m_index
    )

    total_cost = total_revenue + total_penalty

    return total_cost


model.obj = Objective(rule=revenue, sense=maximize)


def energy_balance1(model, m):
    return model.grid[m] == model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m] + model.lost_load[m]


model.energy_balance1 = Constraint(model.m_index, rule=energy_balance1)


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_gen(model, m):
    eq = model.solar_use[m] <= model.solar_avl[m]
    return eq


model.max_solar_gen = Constraint(model.m_index, rule=max_solar_gen)


def max_wind_gen(model, m):
    eq = model.wind_use[m] <= model.wind_avl[m]
    return eq


model.max_wind_gen = Constraint(model.m_index, rule=max_wind_gen)


# Charging and discharging controller
def ein_limit1(model, m):
    eq = model.e_in[m] + (1 - model.ein_eout_lmt[m]) * 1000000 >= 0

    return eq


model.ein_limit1 = Constraint(model.m_index, rule=ein_limit1)


def eout_limit1(model, m):
    eq = model.e_out[m] + (1 - model.ein_eout_lmt[m]) * -1000000 <= 0
    return eq


model.eout_limit1 = Constraint(model.m_index, rule=eout_limit1)


def ein_limit2(model, m):
    eq = model.e_out[m] + (model.ein_eout_lmt[m]) * 1000000 >= 0
    return eq


model.ein_limit2 = Constraint(model.m_index, rule=ein_limit2)


def eout_limit2(model, m):
    eq = model.e_in[m] + (model.ein_eout_lmt[m]) * -1000000 <= 0
    return eq


model.eout_limit2 = Constraint(model.m_index, rule=eout_limit2)


# max charging
def max_charge(model, m):
    return model.e_in[m] <= model.storage_power * interval_freq / 60


model.max_charge = Constraint(model.m_index, rule=max_charge)


# max discharging
def max_discharge(model, m):
    return model.e_out[m] <= model.storage_power * interval_freq / 60


model.max_max_discharge = Constraint(model.m_index, rule=max_discharge)


def soc_update(model, m):
    if m == 0:
        eq = model.soc[m] == \
             (
                     (model.storage_energy * 1) +
                     (model.e_in[m] * model.chargeEff) -
                     (model.e_out[m] / model.dchargeEff)
             )
    else:
        eq = model.soc[m] == \
             (
                     model.soc[m - 1] +
                     (model.e_in[m] * model.chargeEff) -
                     (model.e_out[m] / model.dchargeEff)
             )

    return eq


model.soc_update = Constraint(model.m_index, rule=soc_update)


def soc_min(model, m):
    eq = model.soc[m] >= model.storage_energy * 0.2

    return eq


model.soc_min = Constraint(model.m_index, rule=soc_min)


def soc_max(model, m):
    eq = model.soc[m] <= model.storage_energy * 1

    return eq


model.soc_max = Constraint(model.m_index, rule=soc_max)


def throughput_limit(model):
    energy_sum = sum(model.e_out[idx] for idx in model.m_index)
    return energy_sum <= model.storage_energy * 365


model.throughput_limit = Constraint(rule=throughput_limit)

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 = []
insufficient_dispatch = []
load = []
idx = []

for i in range(len(load_profile)):
    grid_use.append(value(model.grid[i]) * (60 / interval_freq))
    solar.append(value(model.solar_use[i]) * (60 / interval_freq))
    wind.append(value(model.wind_use[i]) * (60 / interval_freq))
    lost_load.append(value(model.lost_load[i]) * (60 / interval_freq))
    e_in.append(value(model.e_in[i]) * (60 / interval_freq))
    e_out.append(value(model.e_out[i]) * (60 / interval_freq))
    soc.append(value(model.soc[i]))
    load.append(value(model.load_profile[i]) * (60 / interval_freq))

time_stamp = pd.date_range('1/1/2022', periods=len(load_profile), freq='15min')
df_out = pd.DataFrame()

# Yearly Dispatch
df_out["Time Stamp"] = time_stamp
df_out["Month"] = time_stamp.month
df_out["Load"] = load
df_out["Solar"] = solar
df_out["Wind"] = wind
df_out["Charge"] = e_in
df_out["Discharge"] = e_out
df_out["SoC"] = np.array(soc)
df_out["Facility Output"] = df_out["Solar"] + df_out["Wind"] + df_out["Discharge"] - df_out["Charge"]
df_out["Shortfall"] = lost_load
df_out["DFR"] = df_out["Facility Output"] / df_out["Load"]

# Monthly Energy Outputs
summary_table = df_out.groupby('Month').sum() / 4
soc_monthly_avg = df_out.groupby('Month').mean()
selected_df = summary_table[['Load', 'Solar', 'Wind', 'Charge', 'Discharge', 'Facility Output']]
selected_df.loc[:, 'SoC'] = soc_monthly_avg['SoC']
selected_df.loc[:, 'Cycles'] = selected_df["Discharge"] / energy
selected_df.loc[:, 'DFR'] = selected_df["Facility Output"] / selected_df["Load"]
selected_df["Excess"] = selected_df.apply(lambda x: (x['Facility Output'] - x['Load']) * excess_revenue_rate if x['Facility Output'] > x['Load'] else 0, axis=1)
selected_df.loc[:, 'Shortfall'] = selected_df.apply(lambda x: DFR * x['Load'] - x['Facility Output'] if x['DFR'] < DFR else 0, axis=1)
selected_df.loc[:, 'Revenue'] = selected_df[['Facility Output', 'Load']].min(axis=1) * PPA
selected_df.loc[:, 'Shortfall Penalty'] = selected_df["Shortfall"] * -shortfall_penalty
selected_df.loc[:, 'Excess Revenue'] = selected_df["Excess"] * excess_revenue_rate
selected_df.loc[:, 'Total '] = selected_df[['Revenue', 'Shortfall Penalty', 'Excess Revenue']].sum(axis=1)
timestamp = datetime.datetime.now().strftime('%y-%m-%d_%H-%M')
with pd.ExcelWriter('sample output.xlsx') as writer:
    df_out.to_excel(writer, sheet_name='dispatch', index=False)
    selected_df.to_excel(writer, sheet_name='Summy')

this is the output from this code:

These are the post calculations. I have already done these calculations in the last few lines of the code after I got the yearly output. I want to do these calculations in the code itself. As I mentioned above, if the monthly facility output is less than 90% of the Total monthly demand then only apply the penalty to it. If the same monthly facility output is greater than the monthly demand profile then the facility should get the revenue based on excess_revenue_rate (50% of PPA).

NOTE: Facility output = Solar + Wind + Battery Discharge - Battery charge