""" Energy Optimization with Pyomo. Refactored from Navid Shirzadi code based on his research work. SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2023 Project Author Alireza Adli alireza.adli@concordia.ca """ import matplotlib.pyplot as plt import pyomo.environ as po from pyomo.opt import SolverFactory import pandas as pd import numpy as np from typing import Dict def pyomo_energy_optimization(user_input: Dict): # Constraints of the model def storage_1(model, t): if t == 1: return model.battery_energy[t] == model.battery_capacity else: return model.battery_energy[t] == model.battery_energy[t - 1] + \ model.charging_capacity[t] * model.battery_charging - \ model.discharging_capacity[t] / model.battery_discharging def storage_2(model, t): return model.battery_energy[t] >= model.battery_capacity * \ model.charge_state_min def storage_3(model, t): return model.battery_energy[t] <= model.battery_capacity def storage_4(model, t): return model.charging_capacity[t] * model.battery_charging + \ model.discharging_capacity[t] / model.battery_discharging <= \ model.charge_discharge_max * model.battery_capacity def storage_5(model, t): return model.charging_capacity[t] <= \ model.battery_capacity * model.gamma[t] def storage_6(model, t): return model.discharging_capacity[t] <= \ model.battery_capacity * model.teta[t] def storage_7(model, t): return model.gamma[t] + model.teta[t] == 1 def grid_1(model, t): return model.sold_grid_capacity[t] <= \ (model.wind_power_surplus[t] + model.pv_power_surplus[t]) * \ model.electricity_to_grid * model.eta[t] def grid_2(model, t): return model.purchased_grid_capacity[t] <= model.grid_purchase_max * \ model.lambdaa[t] def grid_3(model, t): return model.eta[t] + model.lambdaa[t] == 1 def wind_surplus(model, t): return model.wind_turbines_number * model.wind[t] == \ (model.wind_used_energy[t] + model.wind_power_surplus[t]) def pv_surplus(model, t): return model.pv_panels_number * model.pv[t] == model.pv_used_power[t] \ + model.pv_power_surplus[t] def balance(model, t): return model.wind_used_energy[t] + model.pv_used_power[t] + \ model.discharging_capacity[t] + model.purchased_grid_capacity[t] + \ model.loss[t] >= model.Load[t] + model.charging_capacity[t] + \ model.sold_grid_capacity[t] def loss_constraint(model, t): return sum(model.loss[t] for t in model.t) <= 0.01 * \ sum(model.Load[t] for t in model.t) def objective_rule(model): return model.pv_panels_number * \ (model.pv_cost + model.dc_dc_converter_cost) * \ (1 + model.pv_operating_cost / model.capital_recovery_factor) + \ (model.wind_turbines_number * model.wind_turbine_cost * (1 + model.wind_turbine_operating_cost / model.capital_recovery_factor)) + \ model.wind_turbines_number * model.ac_dc_rectifier_cost + \ model.battery_capacity * (model.battery_kwh_price + model.dc_ac_inverter_cost) + \ model.battery_operating_cost * \ (sum((model.discharging_capacity[t] + model.charging_capacity[t]) / model.capital_recovery_factor for t in model.t)) + \ ((model.battery_maintenance_cost * model.charge_discharge_max * model.battery_capacity) + (model.battery_replacement_cost * model.battery_capacity)) * \ model.capital_recovery_factor_battery + \ sum((model.grid_purchase_price * model.purchased_grid_capacity[t]) / model.capital_recovery_factor for t in model.t) - \ sum((model.grid_selling_price * model.sold_grid_capacity[t]) / model.capital_recovery_factor for t in model.t) + \ sum(model.loss[t] * model.loss_coefficient for t in model.t) + \ sum(model.purchased_grid_capacity[t] * model.nonrenewable_grid_portion * model.environmental_coefficient_penalty for t in model.t) default_values = dict(building_load=None, battery_charging=0.95, battery_discharging=0.95, charge_discharge_max=0.35, electricity_to_grid=1, pv_cost=1246, dc_dc_converter_cost=100, grid_purchase_max=200000, wind_turbine_operating_cost=0.002, wind_turbine_cost=2077, capital_recovery_factor=0.064, capital_recovery_factor_battery=1.7743, ac_dc_rectifier_cost=100, dc_ac_inverter_cost=100, battery_kwh_price=670, battery_operating_cost=0.00054, battery_maintenance_cost=13.2, battery_replacement_cost=670, grid_purchase_price=0.08, grid_selling_price=0.08, charge_state_min=0.2, loss_coefficient=0.8, pv_operating_cost=0.001, nonrenewable_grid_portion=0, environmental_coefficient_penalty=1, area=5020.0, wind_turbines_number=1200, battery_capacity=None, battery_energy=None, charging_capacity=None, discharging_capacity=None, wind_used_energy=None, wind_power_surplus=None, pv_used_power=None, pv_power_surplus=None, purchased_grid_capacity=2000, sold_grid_capacity=2000, loss=100) for each in user_input: if user_input[each] is None: user_input[each] = default_values[each] input_daily_sum = pd.read_excel('data/Input_Daily_Sum.xlsx') model = po.ConcreteModel() model.t = po.RangeSet(1, 365) if user_input['building_load']: model.Load = user_input['building_load'] else: model.Load = \ po.Param(model.t, initialize=dict(zip(input_daily_sum.time, input_daily_sum.Load))) model.wind = \ po.Param(model.t, initialize=dict(zip(input_daily_sum.time, input_daily_sum.Wind))) model.pv = \ po.Param(model.t, initialize=dict(zip(input_daily_sum.time, input_daily_sum.PV))) # Parameters of the model model.battery_charging = po.Param( initialize=user_input['battery_charging']) model.battery_discharging = po.Param( initialize=user_input['battery_discharging']) model.charge_discharge_max = po.Param( initialize=user_input['charge_discharge_max']) model.electricity_to_grid = po.Param( initialize=user_input['electricity_to_grid']) model.pv_cost = po.Param( initialize=user_input['pv_cost']) model.dc_dc_converter_cost = po.Param( initialize=user_input['dc_dc_converter_cost']) model.grid_purchase_max = po.Param( initialize=user_input['grid_purchase_max']) model.wind_turbine_cost = po.Param( initialize=user_input['wind_turbine_cost']) model.wind_turbine_operating_cost = po.Param( initialize=user_input['wind_turbine_operating_cost']) model.capital_recovery_factor = po.Param( initialize=user_input['capital_recovery_factor']) model.capital_recovery_factor_battery = po.Param( initialize=user_input['capital_recovery_factor_battery']) model.ac_dc_rectifier_cost = po.Param( initialize=user_input['ac_dc_rectifier_cost']) model.dc_ac_inverter_cost = po.Param( initialize=user_input['dc_ac_inverter_cost']) model.battery_kwh_price = po.Param( initialize=user_input['battery_kwh_price']) model.battery_operating_cost = po.Param( initialize=user_input['battery_operating_cost']) model.battery_maintenance_cost = po.Param( initialize=user_input['battery_maintenance_cost']) model.battery_replacement_cost = po.Param( initialize=user_input['battery_replacement_cost']) model.grid_purchase_price = po.Param( initialize=user_input['grid_purchase_price']) model.grid_selling_price = po.Param( initialize=user_input['grid_selling_price']) model.charge_state_min = po.Param( initialize=user_input['charge_state_min']) model.loss_coefficient = po.Param( initialize=user_input['loss_coefficient']) model.pv_operating_cost = po.Param( initialize=user_input['pv_operating_cost']) model.nonrenewable_grid_portion = po.Param( initialize=user_input['nonrenewable_grid_portion']) model.environmental_coefficient_penalty = po.Param( initialize=user_input['environmental_coefficient_penalty']) # Variables of the model model.pv_panels_number = po.Var( within=po.Integers, bounds=(0, user_input['area'] / 1.32)) model.wind_turbines_number = po.Var( within=po.Integers, bounds=(0, user_input['wind_turbines_number'])) model.battery_capacity = po.Var( bounds=(0, user_input['battery_capacity'])) model.battery_energy = po.Var( model.t, bounds=(0, user_input['battery_energy'])) model.charging_capacity = po.Var( model.t, bounds=(0, user_input['charging_capacity'])) model.discharging_capacity = po.Var( model.t, bounds=(0, user_input['discharging_capacity'])) model.wind_used_energy = po.Var( model.t, bounds=(0, user_input['wind_used_energy'])) model.wind_power_surplus = po.Var( model.t, bounds=(0, user_input['wind_power_surplus'])) model.pv_used_power = po.Var( model.t, bounds=(0, user_input['pv_used_power'])) model.pv_power_surplus = po.Var( model.t, bounds=(0, user_input['pv_power_surplus'])) model.purchased_grid_capacity = po.Var( model.t, bounds=(0, user_input['purchased_grid_capacity'])) model.sold_grid_capacity = po.Var( model.t, bounds=(0, user_input['sold_grid_capacity'])) model.loss = po.Var( model.t, domain=po.NonNegativeReals, bounds=(0, user_input['loss'])) model.gamma = po.Var( model.t, within=po.Binary) model.teta = po.Var( model.t, within=po.Binary) model.lambdaa = po.Var( model.t, within=po.Binary) model.eta = po.Var( model.t, within=po.Binary) # Assigning the constraints functions model.constraint_1 = po.Constraint(model.t, rule=storage_1) model.constraint_2 = po.Constraint(model.t, rule=storage_2) model.constraint_3 = po.Constraint(model.t, rule=storage_3) model.constraint_4 = po.Constraint(model.t, rule=storage_4) model.constraint_5 = po.Constraint(model.t, rule=storage_5) model.constraint_6 = po.Constraint(model.t, rule=storage_6) model.constraint_7 = po.Constraint(model.t, rule=storage_7) model.constraint_8 = po.Constraint(model.t, rule=grid_1) model.constraint_9 = po.Constraint(model.t, rule=grid_2) model.constraint_10 = po.Constraint(model.t, rule=grid_3) model.constraint_11 = po.Constraint(model.t, rule=wind_surplus) model.constraint_12 = po.Constraint(model.t, rule=pv_surplus) model.constraint_13 = po.Constraint(model.t, rule=balance) model.constraint_14 = po.Constraint(model.t, rule=loss_constraint) model.objective = po.Objective(rule=objective_rule, sense=po.minimize) model.write('data/Opt_New.lp') optmization = SolverFactory('scip') results = optmization.solve(model, tee=True) load_sum = input_daily_sum['Load'].sum() object_fu = po.value(model.objective) cost_of_energy = po.value( model.objective) * model.capital_recovery_factor / load_sum number_of_pv = po.value(model.pv_panels_number) number_of_wind_turbines = po.value(model.wind_turbines_number) battery_storage_capacity = po.value(model.battery_capacity) gs_values = pd.DataFrame(list(model.sold_grid_capacity[:].value), columns=['Gs']) gp_values = pd.DataFrame(list(model.purchased_grid_capacity[:].value), columns=['Gp']) pd_values = pd.DataFrame(list(model.discharging_capacity[:].value), columns=['Pd']) pc_values = pd.DataFrame(list(model.charging_capacity[:].value), columns=['Pc']) eb_values = pd.DataFrame(list(model.battery_energy[:].value), columns=['Eb']) e_wind_values = pd.DataFrame(list(model.wind_used_energy[:].value), columns=['EWind']) e_pv_values = pd.DataFrame(list(model.pv_used_power[:].value), columns=['EPV']) es_wind_values = pd.DataFrame(list(model.wind_power_surplus[:].value), columns=['ESWind']) es_pv_values = pd.DataFrame(list(model.pv_power_surplus[:].value), columns=['ESPV']) unmet_load = pd.DataFrame(list(model.loss[:].value), columns=['Loss']) indexed_result = pd.concat([gs_values, gp_values, pd_values, pc_values, eb_values, e_wind_values, e_pv_values, es_wind_values, es_pv_values, unmet_load], axis=1) total_capital_cost = po.value(model.pv_panels_number) * \ (model.pv_cost + model.dc_dc_converter_cost) + \ (po.value(model.wind_turbines_number) * model.wind_turbine_cost) + \ (po.value(model.battery_capacity) * (model.battery_kwh_price + model.dc_ac_inverter_cost)) renewable_penetration = (((np.array(e_pv_values.sum())) + (np.array(e_wind_values.sum())) - (gs_values.sum())) / load_sum).iloc[0] benefit = ((np.array(load_sum) + np.array(gs_values.sum()) - np.array(gp_values.sum())) * 0.07) - 0.02 * total_capital_cost payback_period = np.round(np.array(total_capital_cost)/(benefit * 1.7743)) indexed_result.to_csv('data/results_local_new.csv') results_local_new = pd.read_csv('data/results_local_new.csv') plt.plot(results_local_new['Loss']) return object_fu, cost_of_energy, number_of_pv, \ number_of_wind_turbines, battery_storage_capacity, \ total_capital_cost, renewable_penetration, payback_period[0]