From c2b9f7bec56ae64d080741f59e1483530ce12110 Mon Sep 17 00:00:00 2001 From: Sadam Hussain Date: Sat, 6 Apr 2024 17:08:29 -0400 Subject: [PATCH] Add HEMS --- HEMS | 919 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 919 insertions(+) create mode 100644 HEMS diff --git a/HEMS b/HEMS new file mode 100644 index 0000000..523b26d --- /dev/null +++ b/HEMS @@ -0,0 +1,919 @@ +# NOTE: for two homes + +# +import itertools + +import numpy as np +import pandas as pd +import pylab as pl +import pyomo.environ as py + +# NOTE: model + +from matplotlib import pyplot as plt, gridspec +from matplotlib.pyplot import legend +from pyomo.core import value, RangeSet + +# data = pd.ExcelFile('input_new_cluster.xlsx') #it will read the excel one time you do not need to read again and gain +# Parse the different tab +sheet_data = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='data') +sheet_EV = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='EV') +sheet_EWH = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='EWH') +sheet_wateruse = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='Water_Use') +sheet_watertemp = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='Water_Temp') +sheet_home = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='Home') +sheet_ESS = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='ESS') +sheet_PV = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='PV') +sheet_ambient = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='Ambient_Temp') +sheet_load = pd.read_excel("input_new_cluster_trying.xlsx", sheet_name='Base_Load') + +# +# +model = py.ConcreteModel() +# # i think it should be outside the loop +# +model.h = py.Set(initialize=sheet_home.home) +model.alpha = py.Param(initialize=39, doc='time of arrival ') # 6:30 = 39 +model.beta = py.Param(initialize=89, doc='time of departure ') # 7:00 = 89 +model.t_final = py.Param(initialize=96, doc='final time of the day ') +model.second = py.Param(initialize=2, doc='second time of the day ') +model.t = py.Set(initialize=sheet_data.time, ordered=True, doc='time period') +model.EV_Dur = py.Set(initialize=RangeSet(model.alpha.value, model.beta.value), + doc='time interval for EV') +model.EV_Dur1 = py.Set(initialize=RangeSet(model.alpha.value + 1, model.beta.value), + doc='time interval for EV after arrival') +model.t_second = py.Set(initialize=RangeSet(model.second.value, model.t_final.value), + doc='time greater then 1') +model.t_first = py.Param(initialize=1, doc='first time of the day ') +model.ist_interval = py.Set(initialize=RangeSet(model.t_first.value, model.alpha.value - 1), + doc='time period before the arrival time') +model.second_interval = py.Set(initialize=RangeSet(model.beta.value + 1, model.t_final.value), + doc='time interval after departure') +# + +# NOTE:Parameter + +# Base load : uncontrollable loads + + +model.base_load = py.Param(model.h, model.t, + initialize=dict(zip(list(itertools.product(model.h.data(), model.t.data())), + [i for x in sheet_load.columns if "load" in x for i in + sheet_load.loc[:, x].values])), doc="Base Load") + +model.PV = py.Param(model.h, model.t, initialize=dict(zip(list(itertools.product(model.h.data(), model.t.data())), + [i for x in sheet_PV.columns if "PV" in x for i in + sheet_PV.loc[:, x].values])), doc="PV Production") +model.max = py.Param(initialize=100, doc='maximum value for selling and buying power') +# +# todo this is need to be done for single value but tow home has seperate values +model.ChargeRate_EV = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())), + [i for x in sheet_EV.columns if "charging_rate" in x for i + in + sheet_EV.loc[:, x].values])), doc="charging rate of EV") + +model.DischRate_EV = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())), + [i for x in sheet_EV.columns if "rate_of_discharging" in x + for i in + sheet_EV.loc[:, x].values])), doc="Discharging rate of EV") + +model.Ch_Effc_EV = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())), + [i for x in sheet_EV.columns if "charging_efficiency" in x for + i in + sheet_EV.loc[:, x].values])), doc="charging efficency of EV") + +model.DischEffc_EV = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())), + [i for x in sheet_EV.columns if + "efficiency_of_discharging" in x for i + in sheet_EV.loc[:, x].values])), + doc="Discharging efficency of EV") + +model.Cap_EV = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EV.columns if "capacity" in x for i in + sheet_EV.loc[:, x].values])), doc="Capacity of EV") +model.End_percentage_EV = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EV.columns if "end" in x + for i in + sheet_EV.loc[:, x].values])), + doc="Departure energy of EV") + +model.In_Percentage_EV = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EV.columns if "initial" in x + for i in + sheet_EV.loc[:, x].values])), + doc="Initial energy of EV") + +model.Energy_EV_dep = py.Param(model.h, initialize=dict( + zip(model.h.data(), np.array(model.Cap_EV.values()) * model.End_percentage_EV.values())), # just changeing np>py + doc="Departure temperature of EV") + +model.Energy_EV_In = py.Param(model.h, initialize=dict( + zip(model.h.data(), np.array(model.Cap_EV.values()) * model.In_Percentage_EV.values())), + doc="Initial energy of EV") + +# +# + +model.ChargeRate_ESS = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())), + [i for x in sheet_ESS.columns if "charging_rate" in x for i + in + sheet_ESS.loc[:, x].values])), doc="charging rate of ESS") + +model.DischRate_ESS = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())), + [i for x in sheet_ESS.columns if "rate_of_discharging" in x + for + i in + sheet_ESS.loc[:, x].values])), + doc="Discharging rate of ESS") + +# model.ChargeRate_ESS = py.Param(initialize=float(sheet_ESS.charging_rate1), doc='Charging rate of ESS ') +# model.DischRate_ESS = py.Param(initialize=float(sheet_ESS.discharging_rate1), +# doc='Discharging rate of ESS ') +model.Cap_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_ESS.columns if "capacity" in x for i in + sheet_ESS.loc[:, x].values])), doc="Capacity of ESS") + +model.End_percentage_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_ESS.columns if "end" in x + for i in + sheet_ESS.loc[:, x].values])), + doc="Departure energy of ESS") + +model.In_Percentage_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_ESS.columns if + "initial" in x for i in + sheet_ESS.loc[:, x].values])), + doc="Initial energy of ESS") + +model.Ch_Effc_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_ESS.columns if "charging_efficiency" in x + for i in + sheet_ESS.loc[:, x].values])), + doc="charging efficiency of ESS") + +model.DischEffc_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_ESS.columns if + "efficiency_of_dicharging" in x for + i in + sheet_ESS.loc[:, x].values])), + doc="Discharging efficiency of ESS") + +model.Energy_ESS_In = py.Param(model.h, initialize=dict( + zip(model.h.data(), np.array(model.Cap_ESS.values()) * model.In_Percentage_ESS.values())), + doc="Initial energy of ESS") + +model.End_En_ESS = py.Param(model.h, initialize=dict( + zip(model.h.data(), np.array(model.Cap_ESS.values()) * model.End_percentage_ESS.values())), + doc="Departure energy of ESS") + +# +# + +model.tetta_low = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EWH.columns if "tetta_low" in x for i in + sheet_EWH.loc[:, x].values])), + doc="Lower bound of the water temperature") +model.tetta_up = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EWH.columns if "tetta_up" in x for i in + sheet_EWH.loc[:, x].values])), + doc="Upper bound of the water temperature") +model.tetta_amb_int = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EWH.columns if "tetta_amb_init" in x for i + in + sheet_EWH.loc[:, x].values])), + doc="Initial ambient temperature") + + +# todo need to use iteration for this below code + +model.tetta_amb = py.Param(model.h, model.t, + initialize=dict(zip(list(itertools.product(model.h.data(), model.t.data())), + [i for x in sheet_ambient.columns if "celsius" in x for i in + sheet_ambient.loc[:, x].values])), doc="Outdoor temperature") + +model.Q = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EWH.columns if "capacity" in x for i in + sheet_EWH.loc[:, x].values])), doc="Power of the EWH") + +model.R = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EWH.columns if "thermal_resistance" in x for i in + sheet_EWH.loc[:, x].values])), doc="Thermal resistance") +model.C = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EWH.columns if "thermal_capacitance" in x for i in + sheet_EWH.loc[:, x].values])), doc="Thermal resistance") + +model.M = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EWH.columns if "water_cap" in x for i in + sheet_EWH.loc[:, x].values])), doc="Water Capacity (L)") + +model.tetta_EWH_int = py.Param(model.h, initialize=dict(zip((model.h.data()), + [i for x in sheet_EWH.columns if "tetta_wat_init" in x for i + in + sheet_EWH.loc[:, x].values])), doc="Initial temperature") + +model.water_use = py.Param(model.h, model.t, + initialize=dict(zip(list(itertools.product(model.h.data(), model.t.data())), + [i for x in sheet_wateruse.columns if "Litre" in x for i in + sheet_wateruse.loc[:, x].values])), doc="Hot Water usage") + +# +model.Buy_price = py.Param(model.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price)), + doc='Buying Price') +model.Sell_price = py.Param(model.t, initialize=dict(zip(sheet_data.time, sheet_data.Sell_price)), doc='Selling Price') +# Time duration: we took 15 mint granularity so in one hour it will be 1/4 +model.time_d = py.Param(initialize=(1 / 4), doc='time duration ') +# model.DPT = py.Param(initialize=0.069, doc='Daily power tariff ') # .069 +# + +# NOTE:Variable + +# +model.P_Buy_Grid = py.Var(model.h, model.t, bounds=(0, None)) +model.P_Sell_Grid = py.Var(model.h, model.t, bounds=(0, None)) +model.S_P_sell = py.Var(model.h, model.t, within=py.Binary) +model.S_P_buy = py.Var(model.h, model.t, within=py.Binary) +model.peak = py.Var(within=py.NonNegativeReals) + +# + +# +model.P_EV_Charge = py.Var(model.h, model.t, bounds=(0, None)) +model.P_EV_Disch = py.Var(model.h, model.t, bounds=(0, None)) +model.P_EV_Disch_Home = py.Var(model.h, model.t, bounds=(0, None)) +model.P_EV_Disch_Grid = py.Var(model.h, model.t, bounds=(0, None)) +model.S_EV_Charge = py.Var(model.h, model.t, within=py.Binary) +model.S_EV_Disch = py.Var(model.h, model.t, within=py.Binary) +model.Energy_EV = py.Var(model.h, model.t, bounds=(0, None)) # SOC of the EV +# + +# +model.Energy_ESS = py.Var(model.h, model.t, bounds=(0, None)) # SOC of the ESS +model.S_ESS_Charge = py.Var(model.h, model.t, within=py.Binary) +model.P_ESS_Disch = py.Var(model.h, model.t, bounds=(0, None)) +model.P_ESS_Charge = py.Var(model.h, model.t, bounds=(0, None)) +model.P_ESS_Charge_Grid = py.Var(model.h, model.t, bounds=(0, None), doc='ESS charging from the Grid') +model.P_ESS_Disch_Home = py.Var(model.h, model.t, bounds=(0, None)) +model.P_ESS_Disch_Grid = py.Var(model.h, model.t, bounds=(0, None)) +model.S_ESS_Disch = py.Var(model.h, model.t, within=py.Binary) +# + +# +model.PV_Home = py.Var(model.h, model.t, bounds=(0, None)) +model.PV_Grid = py.Var(model.h, model.t, bounds=(0, None)) +model.PV_Battery = py.Var(model.h, model.t, bounds=(0, None)) +# + +# +model.tetta_EWH_wat = py.Var(model.h, model.t) +model.S_EWH = py.Var(model.h, model.t, within=py.Binary) +model.P_EWH = py.Var(model.h, model.t, doc='power of EWH') + + +# + +# NOTE:Constraints + +# +def Power_buy(model, h, i): + return model.P_Buy_Grid[h, i] == model.base_load[h, i] + model.Q[h] * model.S_EWH[h, i] + model.P_EV_Charge[h, i] - \ + model.P_EV_Disch_Home[h, i] + model.P_ESS_Charge[h, i] - \ + model.P_ESS_Disch_Home[h, i] - model.PV_Home[h, i]- model.PV_Battery[h, i] + + +model.Const_1 = py.Constraint(model.h, model.t, rule=Power_buy, doc='Power buy from the Grid') + + +def Power_buy2(model, h, i): + return model.P_Buy_Grid[h, i] <= model.max * model.S_P_buy[h, i] + + +model.Const_1a = py.Constraint(model.h, model.t, rule=Power_buy2, + doc='removing the nonlinearity in the objective fucntion') + + +def Power_sell1(model, h, i): + return model.P_Sell_Grid[h, i] == model.P_EV_Disch_Grid[h, i] + model.PV_Grid[h, i] + model.P_ESS_Disch_Grid[h, i] + + +model.Const_2 = py.Constraint(model.h, model.t, rule=Power_sell1, doc='Power sell to the Grid') + + +def Power_sell2(model, h, i): + return model.P_Sell_Grid[h, i] <= model.max * model.S_P_sell[h, i] + + +model.Const_2a = py.Constraint(model.h, model.t, rule=Power_sell2, doc='Power sell to the Grid') + + +def Status_Power(model, h, i): + return model.S_P_buy[h, i] + model.S_P_sell[h, i] <= 1 + + +model.Const_3 = py.Constraint(model.h, model.t, rule=Status_Power, + doc='Buying and selling power will not occur at same time') + + +# + +# +# model.P_max = py.Var(model.h, model.t, bounds=(0, None), doc=" Updated maximum demand of house") +# model.P_flex = py.Var(model.h, model.t, bounds=(0, None), doc=" Flexibility available in one house") +# model.P_flex_all = py.Var(model.h, model.t, bounds=(0, None), doc=" Sum of the Flexibility in all") +# model.S_device = py.Var(model.h, model.t, within=py.Binary, doc=" Status of the devices") # but i not using this +# model.flex_index = py.Param(doc='flexibility index ') +# +# #Todo need to verify these equation with omer: +# +# # +# def P_h_flex(model, h, i): +# return model.P_flex_all[i] == sum( +# model.base_load[h, i] * model.S_device[h, i] - model.base_load[h, i] * model.S_device[h, i] for h in model.h +# for i in model.t) +# +# +# model.Const_4a = py.Constraint(model.h, model.t, rule=P_h_flex, doc='sum of the flexibility') +# +# +# def P_flexibility(model, h, i): +# return model.P_flex_all[i] == sum([model.P_flex[h, i] for h in model.h for i in model.t]) +# +# +# model.Const_4b = py.Constraint(model.h, model.t, rule=P_flexibility, doc='sum of the flexibility') +# +# +# def Flexibility_index(model, h, i): +# return model.flex_index[h, i] == model.P_flex[h, i] / model.P_flex_all[i] +# +# +# model.Const_4c = py.Constraint(model.h, model.t, rule=Flexibility_index, doc='sum of the flexibility') +# +# +# # this equation make the two way communication +# def P_maximum(model, h, i): +# return model.P_max[h, i] == model.P_max[h, i - 1] - model.P_flex[h, i] +# +# +# model.Const_4d = py.Constraint(model.h, model.t, rule=P_maximum, doc='updated maximum power') + + +# + +# +def Power_EV_Charge_limit1(model, h, i): + return model.P_EV_Charge[h, i] <= model.ChargeRate_EV[h] * model.S_EV_Charge[h, i] + + +model.Const_EV_1 = py.Constraint(model.h, model.t, rule=Power_EV_Charge_limit1, + doc='Charging power of EV (Upper limit)') + + +def Power_EV_Charge_limit2(model, h, i): + return model.P_EV_Charge[h, i] >= 0.2 * (model.ChargeRate_EV[h] * model.S_EV_Charge[h, i]) + + +model.Const_EV_2 = py.Constraint(model.h, model.t, rule=Power_EV_Charge_limit2, + doc='Charging power of EV (lower limit)') + + +def Power_EV_Disch_limit1(model, h, i): + return model.P_EV_Disch[h, i] <= (model.DischRate_EV[h] * model.S_EV_Disch[h, i]) + + +model.Const_EV_3 = py.Constraint(model.h, model.t, rule=Power_EV_Disch_limit1, + doc='Discharging power of EV (Upper limit)') + + +def Power_EV_Disch_limit2(model, h, i): + return model.P_EV_Disch[h, i] >= 0.2 * (model.DischRate_EV[h] * model.S_EV_Disch[h, i]) + + +model.Const_EV_4 = py.Constraint(model.h, model.t, rule=Power_EV_Disch_limit2, + doc='Discharging power of EV (Lower limit)') + + +def Status_EV(model, h, i): + return model.S_EV_Disch[h, i] + model.S_EV_Charge[h, i] <= 1 + + +model.Const_EV_5 = py.Constraint(model.h, model.t, rule=Status_EV, + doc='Charging and discharging will not occur at same time') + + +def Power_EV_Disch(model, h, i): + return model.P_EV_Disch[h, i] == ( + (model.P_EV_Disch_Home[h, i] + model.P_EV_Disch_Grid[h, i]) / model.DischEffc_EV[h]) + + +model.Const_EV_6 = py.Constraint(model.h, model.t, rule=Power_EV_Disch, doc='Discharging power of EV to ' + 'Home and Grid') + + +def SoC_EV1(model, h, i): + return model.Energy_EV[h, i] == model.Energy_EV_In[h] + + +# + ( +# model.P_ESS_Charge[h, i] * model.Ch_Effc_ESS[h] - model.P_ESS_Disch[h, i]) * model.time_d +model.Const_EV_7a = py.Constraint(model.h, [model.alpha.value], rule=SoC_EV1, doc='SoC of the EV at arrival') + + +def SoC_EV2(model, h, i): + return model.Energy_EV[h, i] == model.Energy_EV[h, i - 1] + ( + model.P_EV_Charge[h, i] * model.Ch_Effc_EV[h] - model.P_EV_Disch[h, i]) * model.time_d + + +model.Const_EV_7b = py.Constraint(model.h, model.EV_Dur1, rule=SoC_EV2, doc='SoC of the EV after arrival time') + + +# one equation of the SoC +# def SoC_EV(model, i): +# return model.Energy_EV[i] == model.Energy_EV_In + (model.P_EV_Charge[i] * model.Ch_Effc_EV - model.P_EV_Disch[i]) * model.time_d +# +# +# model.Const_EV_7 = py.Constraint(model.EV_Dur, rule=SoC_EV, doc='SoC of the EV') + +def EV_availability1(model, h, i): + return model.Energy_EV[h, i] == 0 + + +model.Const_EV_8a = py.Constraint(model.h, model.ist_interval, rule=EV_availability1, + doc='SOC available before arrival time') + + +def EV_availability2(model, h, i): + return model.Energy_EV[h, i] == 0 + + +model.Const_EV_8b = py.Constraint(model.h, model.second_interval, rule=EV_availability2, + doc='SOC available after departure time') + + +def EV_status_available1(model, h, i): + return model.S_EV_Disch[h, i] + model.S_EV_Charge[h, i] == 0 + + +model.Const_EV_9a = py.Constraint(model.h, model.ist_interval, rule=EV_status_available1, + doc='EV availability before arrival time') + + +def EV_status_available2(model, h, i): + return model.S_EV_Disch[h, i] + model.S_EV_Charge[h, i] == 0 + + +model.Const_EV_9b = py.Constraint(model.h, model.second_interval, rule=EV_status_available2, + doc='EV availability after arrival time') + + +def EV_SoC_limit1(model, h, i): + return model.Energy_EV[h, i] >= 0.1 * model.Cap_EV[h] + + +model.Const_EV_10 = py.Constraint(model.h, model.EV_Dur, rule=EV_SoC_limit1, doc='Minimum SoC of EV') + + +def EV_SoC_limit2(model, h, i): + return model.Energy_EV[h, i] <= model.Cap_EV[h] + + +model.Const_EV_11 = py.Constraint(model.h, model.t, rule=EV_SoC_limit2, doc='Maximum SoC of EV') + + +def EV_final_SoC(model, h, i): + return model.Energy_EV[h, i] == model.Energy_EV_dep[h] + + +model.Const_EV_12 = py.Constraint(model.h, [model.beta.value], rule=EV_final_SoC, + doc='Final SoC of EV at departure time') + + +# + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +def SoC_ESS1(model, h, i): + return model.Energy_ESS[h, i] == model.Energy_ESS_In[h] + + +# #+ ( +# model.P_ESS_Charge[h, i] * model.Ch_Effc_ESS[h] - model.P_ESS_Disch[h, i]) * model.time_d + + +model.Const_ESS_1a = py.Constraint(model.h, [model.t_first.value], rule=SoC_ESS1, + doc='SoC of ESS') # need to check wather we need to put in the brackets or not + + +def SoC_ESS2(model, h, i): + return model.Energy_ESS[h, i] == model.Energy_ESS[h, i - 1] + ( + model.P_ESS_Charge[h, i] * model.Ch_Effc_ESS[h] - model.P_ESS_Disch[h, i]) * model.time_d + + +model.Const_ESS_1b = py.Constraint(model.h, model.t_second, rule=SoC_ESS2, doc='SoC of ESS') + + +def Power_ESS_Charge(model, h, i): + return model.P_ESS_Charge[h, i] == model.P_ESS_Charge_Grid[h, i] + model.PV_Battery[h, i] + + +model.Const_ESS_2 = py.Constraint(model.h, model.t, rule=Power_ESS_Charge, + doc='Charging power of ESS ') + + +def Power_ESS_Charge_limit1(model, h, i): + return model.P_ESS_Charge[h, i] <= model.ChargeRate_ESS[h] * model.S_ESS_Charge[h, i] + + +model.Const_ESS_2a = py.Constraint(model.h, model.t, rule=Power_ESS_Charge_limit1, + doc='Charging power of ESS (Upper limit)') + + +def Power_ESS_Charge_limit2(model, h, i): + return model.P_ESS_Charge[h, i] >= 0.2 * (model.ChargeRate_ESS[h] * model.S_ESS_Charge[h, i]) + + +model.Const_ESS_2b = py.Constraint(model.h, model.t, rule=Power_ESS_Charge_limit2, + doc='Charging power of ESS (Lower limit)') + + +def Power_ESS_Disch_limit1(model, h, i): + return model.P_ESS_Disch[h, i] <= (model.DischRate_ESS[h] * model.S_ESS_Disch[h, i]) + + +model.Const_ESS_4 = py.Constraint(model.h, model.t, rule=Power_ESS_Disch_limit1, + doc='Discharging power of ESS (Upper limit)') + + +def Power_ESS_Disch_limit2(model, h, i): + return model.P_ESS_Disch[h, i] >= 0.2 * (model.DischRate_ESS[h] * model.S_ESS_Disch[h, i]) + + +model.Const_ESS_5 = py.Constraint(model.h, model.t, rule=Power_ESS_Disch_limit2, + doc='Discharging power of ESS (Lower limit)') + + +def Status_ESS(model, h, i): + return model.S_ESS_Disch[h, i] + model.S_ESS_Charge[h, i] <= 1 + + +model.Const_ESS_6 = py.Constraint(model.h, model.t, rule=Status_ESS, + doc='Charging and discharging will not occur at same time') + + +def Power_ESS_Disch(model, h, i): + return (model.P_ESS_Disch[h, i] * model.DischEffc_ESS[h]) == ( + model.P_ESS_Disch_Home[h, i] + model.P_ESS_Disch_Grid[h, i]) + + +model.Const_ESS_7 = py.Constraint(model.h, model.t, rule=Power_ESS_Disch, + doc='Discharging power of ESS to Home and Grid') + + +def ESS_SoC_limit1(model, h, i): + return model.Energy_ESS[h, i] >= 0.2 * model.Cap_ESS[h] + + +model.Const_ESS_8 = py.Constraint(model.h, model.t, rule=ESS_SoC_limit1, doc='Minimum SoC of ESS') + + +def ESS_SoC_limit2(model, h, i): + return model.Energy_ESS[h, i] <= model.Cap_ESS[h] + + +model.Const_ESS_9 = py.Constraint(model.h, model.t, rule=ESS_SoC_limit2, doc='Maximum SoC of ESS') + + +def ESS_final_SoC(model, h, i): + return model.Energy_ESS[h, i] >= model.End_En_ESS[h] + + +model.Const_ESS_10 = py.Constraint(model.h, [model.t_final.value], rule=ESS_final_SoC, + doc='Final SoC of ESS at departure time') + + +# + +# +def EWH_limit1(model, h, i): + return model.tetta_EWH_wat[h, i] <= model.tetta_up[h] + + +model.Const_EWH_1 = py.Constraint(model.h, model.t, rule=EWH_limit1, doc='Maximum Limit') + + +def EWH_limit2(model, h, i): + return model.tetta_EWH_wat[h, i] >= model.tetta_low[h] + + +model.Const_EWH_2 = py.Constraint(model.h, model.t, rule=EWH_limit2, doc='Minimum Limit') + + +def EWH_power(model, h, i): + return model.P_EWH[h, i] == model.Q[h] * model.S_EWH[h, i] + + +model.Const_EWH_3 = py.Constraint(model.h, model.t, rule=EWH_power, doc='Electric water heater power') + + +# +# # NOTE: I write this function different from the GAMS +def EWH_temp_1(model, h, i): + return model.tetta_EWH_wat[h, i] == model.tetta_amb[h, i] + model.Q[h] * model.S_EWH[h, i] * model.R[ + h] * model.time_d - ( + (model.M[h] - model.water_use[h, i]) / model.M[h]) * ( + model.tetta_amb_int[h] - model.tetta_EWH_int[h]) * py.exp(-model.time_d / (model.R[h] * model.C[h])) + + +model.Const_EWH_4 = py.Constraint(model.h, [model.t_first.value], rule=EWH_temp_1, doc='EWH Model') + + +def EWH_temp_2(model, h, i): + return model.tetta_EWH_wat[h, i] == model.tetta_amb[h, i] + model.Q[h] * model.S_EWH[h, i] * model.R[ + h] * model.time_d - ( + (model.M[h] - model.water_use[h, i]) / model.M[h]) * ( + model.tetta_amb[h, i] - model.tetta_EWH_wat[h, i - 1]) * py.exp( + -model.time_d / (model.R[h] * model.C[h])) + + +model.Const_EWH_5 = py.Constraint(model.h, model.t_second, rule=EWH_temp_2, doc='EWH Model') + + +# + +# # need to include the PV equation of omer paper +# the below equation did not work got an error ( i dont know why!!!) +# def PV_production(model, h, i): +# return model.PV[ h,i] == model.PV_Grid[ h,i] + model.PV_Home[ h,i] +# +# +# model.Const_PV = py.Constraint(model.h, model.t, rule=PV_production, doc='PV Production') + + +def PV_production(model, h, i): + return model.PV_Grid[h, i] + model.PV_Home[h, i] + model.PV_Battery[h, i] == model.PV[h, i] + + +model.Const_PV = py.Constraint(model.h, model.t, rule=PV_production, doc='PV Production') + + +# + +# NOTE: I need to include the third term which is in the paper: daily peak demand * charge of the power +def objective_rule(model): + return sum((model.P_Buy_Grid[h, i] * model.Buy_price[i]) - ( + model.P_Sell_Grid[h, i] * model.Sell_price[i]) for i in + model.t for h in model.h) * model.time_d + + +model.objective = py.Objective(rule=objective_rule, sense=py.minimize, doc='Definition of objective function') + +# +model.write('model2.lp', io_options={'symbolic_solver_labels': True}) +opt = py.SolverFactory('cplex') +opt.options["mipgap"] = 0.09 # 0.155 for load and 0.8 for other load +result = opt.solve(model, tee=True) +# result = opt.solve(model, tee = True) +# model.pprint() +print(result) +# + +# + +if (result.solver.status == py.SolverStatus.ok) and ( + result.solver.termination_condition == py.TerminationCondition.optimal): + print('optimal solution') # Do something when the solution in optimal and feasible +elif result.solver.termination_condition == py.TerminationCondition.infeasible: + print('solution is not feasible') # Do something when model in infeasible +else: + # Something else is wrong + print('Solver Status: ', result.solver.status) +print('Sum of all homes objectives = ', py.value(model.objective)) +# + +for h in model.h: + print( + f"Objective of Home {h} : " f"{py.value((sum((model.P_Buy_Grid[h, i] * model.Buy_price[i]) - (model.P_Sell_Grid[h, i] * model.Sell_price[i]) for i in model.t) * model.time_d))}") +# + + +# Plotting for automatic ploting +time = [i for i in model.t] +pbuy = [[] for j in model.h] +psell = [[] for j in model.h] +PEVCharge = [[] for j in model.h] +PESSCharge = [[] for j in model.h] +PEVDischHome = [[] for j in model.h] +PESSDischHome = [[] for j in model.h] +SEVCharge = [[] for j in model.h] +SEVDisch = [[] for j in model.h] +PEVDischGrid = [[] for j in model.h] +PESSDischGrid = [[] for j in model.h] +PVHome = [[] for j in model.h] +PVGrid = [[] for j in model.h] +PVBattery = [[] for j in model.h] +PPV = [[] for j in model.h] +baseload = [[] for j in model.h] +EnergyESS = [[] for j in model.h] +EnergyEV = [[] for j in model.h] +tettaEWHwat = [[] for j in model.h] +SEWH = [[] for j in model.h] +PEWH = [[] for j in model.h] +Buyprice = [] # cost = [i for i in model.c.values()] <<< check this +Sellprice = [] + +for j in model.h: + for k, v in model.P_Buy_Grid.items(): + if k[0] == j: + pbuy[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.P_Sell_Grid.items(): + if k[0] == j: + psell[j - 1].append(py.value(v)) + +for j in model.h: + for k, v in model.P_EV_Charge.items(): + if k[0] == j: + PEVCharge[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.P_EV_Disch_Home.items(): + if k[0] == j: + PEVDischHome[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.P_EV_Disch_Grid.items(): + if k[0] == j: + PEVDischGrid[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.P_ESS_Charge.items(): + if k[0] == j: + PESSCharge[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.P_ESS_Disch_Home.items(): + if k[0] == j: + PESSDischHome[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.P_ESS_Disch_Grid.items(): + if k[0] == j: + PESSDischGrid[j - 1].append(py.value(v)) + +for j in model.h: + for k, v in model.PV_Home.items(): + if k[0] == j: + PVHome[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.PV_Grid.items(): + if k[0] == j: + PVGrid[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.PV_Battery.items(): + if k[0] == j: + PVBattery[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.PV.items(): + if k[0] == j: + PPV[j - 1].append(py.value(v)) + +for j in model.h: + for k, v in model.base_load.items(): + if k[0] == j: + baseload[j - 1].append(py.value(v)) + +for j in model.h: + for k, v in model.Energy_ESS.items(): + if k[0] == j: + EnergyESS[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.Energy_EV.items(): + if k[0] == j: + EnergyEV[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.P_EWH.items(): + if k[0] == j: + PEWH[j - 1].append(py.value(v)) +for j in model.h: + for k, v in model.tetta_EWH_wat.items(): + if k[0] == j: + tettaEWHwat[j - 1].append(py.value(v)) + +for j in model.h: + for k, v in model.S_EWH.items(): + if k[0] == j: + SEWH[j - 1].append(py.value(v)) + +for i in model.Buy_price: + Buyprice.append(value(model.Buy_price[i])) + +for i in model.Sell_price: + Sellprice.append(value(model.Sell_price[i])) + +# alone +# for k in range(len(model.h)): +# fig, ax = plt.subplots(5, 3, figsize=(10, 10)) +# ax[0, 0].bar(time, pbuy[k], label='Buying power') +# ax[0, 0].bar(time, psell[k], label='Selling power') +# ax[0, 0].legend(loc='best', fontsize='small', ncol=3) +# ax[0, 1].bar(time, baseload[k], label='Base load', color='r') +# ax[0, 1].legend(loc='best', fontsize='small', ncol=3) +# ax[0, 2].plot(time, Buyprice, label='Buyprice') +# ax[0, 2].plot(time, Sellprice, label='Sellprice') +# ax[0, 2].legend(loc='best', fontsize='small', ncol=3) +# ax[1, 0].bar(time, PEWH[k], label='EWH Power', color='r') +# ax[1, 0].legend(loc='best', fontsize='small', ncol=3) +# ax[1, 1].plot(time, tettaEWHwat[k], label='EWH Temp') +# ax[1, 1].legend(loc='best', fontsize='small', ncol=3) +# ax[1, 2].bar(time, SEWH[k], label='Status of EWH') +# ax[1, 2].legend(loc='best', fontsize='small', ncol=3) +# ax[2, 0].bar(time, PESSCharge[k], label='ESS Charging power', color='r') +# ax[2, 0].legend(loc='best', fontsize='small', ncol=3) +# ax[2, 1].bar(time, PESSDischHome[k], label='ESS Disch to home') +# ax[2, 1].bar(time, PESSDischGrid[k], label='ESS Disch to grid') +# ax[2, 1].legend(loc='best', fontsize='small', ncol=3) +# ax[2, 2].plot(time, EnergyESS[k], label='Energy of ESS', color='g') +# ax[2, 2].legend(loc='best', fontsize='small', ncol=3) +# ax[3, 0].bar(time, PVHome[k], label='PV to Home') +# ax[3, 0].legend(loc='best', fontsize='small', ncol=3) +# ax[3, 1].bar(time, PVGrid[k], label='PV to Grid') +# ax[3, 1].bar(time, PVBattery[k], label='PV to battery') +# ax[3, 1].legend(loc='best', fontsize='small', ncol=3) +# ax[3, 2].bar(time, PVHome[k], label='PV to Home') +# ax[3, 2].bar(time, PVGrid[k], label='PV to Grid') +# ax[3, 2].legend(loc='best', fontsize='small', ncol=3) +# ax[4, 0].bar(time, PEVCharge[k], label='EV Charging power', color='r') +# ax[4, 0].legend(loc='best', fontsize='small', ncol=3) +# ax[4, 1].bar(time, PEVDischHome[k], label='EV Disch to home') +# ax[4, 1].bar(time, PEVDischGrid[k], label='EV Disch to grid') +# ax[4, 1].legend(loc='best', fontsize='small', ncol=3) +# ax[4, 2].plot(time, EnergyEV[k], label='Energy of EV', color='r') +# ax[4, 2].legend(loc='best', fontsize='small', ncol=3) +# ax[4, 0].set_xlabel('Time (step)') +# ax[4, 1].set_xlabel('Time (step)') +# ax[4, 2].set_xlabel('Time (step)') +# ax[0, 2].set_ylabel('price ($/W/h)') +# ax[0, 0].set_ylabel('Power (kW)') +# ax[1, 0].set_ylabel('Power (kW)') +# ax[2, 0].set_ylabel('Power (kW)') +# ax[3, 0].set_ylabel('Power (kW)') +# ax[4, 0].set_ylabel('Power (kW)') +# ax[2, 2].set_ylabel('Energy (kWh)') +# ax[4, 2].set_ylabel('Energy (kWh)') +# plt.suptitle( +# f" Home {k + 1} : " f" Cost is : {py.value((sum((model.P_Buy_Grid[k + 1, i] * model.Buy_price[i]) - (model.P_Sell_Grid[k + 1, i] * model.Sell_price[i]) for i in model.t) * model.time_d))}") +# pl.tight_layout() +# plt.show() + +# code for if you want to plot all the home togathoer +fig, ax = plt.subplots(8, len(model.h), figsize=(10,10)) +for k in range(len(model.h)): + ax[0, k].plot(time, pbuy[k], label='Buying power') + ax[0, k].plot(time, psell[k], label='Selling power') + ax[0, k].plot(time, baseload[k], label='Base load') + ax[0, k].legend(loc='best') + ax[1, k].plot(time, Buyprice, label='Buyprice') + ax[1, k].legend(loc='best') + ax[1, k].plot(time, Sellprice, label='Sellprice') + ax[1, k].legend(loc='best') + ax[2, k].plot(time, PEWH[k], label='EWH Power') + ax[2, k].plot(time, tettaEWHwat[k], label='EWH Temp') + ax[2, k].legend(loc='best') + ax[3, k].plot(time, PESSCharge[k],label='ESS Charging power') + ax[3, k].plot(time, PESSDischHome[k], label='ESS Disch to home') + ax[3, k].plot(time, PESSDischGrid[k], label='ESS Disch to grid') + ax[3, k].legend(loc='best') + ax[4, k].plot(time, PVHome[k], label='PV to Home') + ax[4, k].plot(time, PVGrid[k], label='PV to Grid') + ax[4, k].plot(time, PPV[k], label='All PV Production') + ax[4, k].legend(loc='best') + ax[5, k].plot(time, PEVCharge[k], label='EV Charging power') + ax[5, k].plot(time, PEVDischHome[k], label='EV Disch to home') + ax[5, k].plot(time, PEVDischGrid[k], label='EV Disch to grid') + ax[5, k].legend(loc='best') + ax[6, k].plot(time, EnergyESS[k], label='Energy of ESS') + ax[6, k].legend(loc='best') + ax[7, k].plot(time, EnergyEV[k], label='Energy of EV') + ax[7, k].legend(loc='best') + ax[4, k].set_xlabel('Time (step)') + ax[1, k].set_ylabel('Electricity price ($/W/h)') + ax[0, k].set_ylabel('Power (kW)') + ax[2, k].set_ylabel('Power (kW)') + ax[0, k].set_title(f"Home {k+1}") + +pl.tight_layout() + +plt.show() + + +# if you want to plot each home alone in the figure +# for k in range(len(model.h)): +# fig, ax = plt.subplots(2) +# ax[0].plot(time, pbuy[k], 'b-', label='Base') +# ax[0].plot(time, psell[k], 'g--', label='Controllable Load') +# ax[1].plot(time, Buyprice, label='buying price') +# ax[1].plot(time, Sellprice, 'y-+', label='selling price') +# ax[1].legend(loc='upper left') +# ax[1].set_xlabel('Time (hour)') +# ax[1].set_ylabel('Electricity price ($/W/h)') +# plt.suptitle(f"Home {k+1}") +# pl.tight_layout() +# plt.show() + + +# print('this is the end of code') +# #