From 254a89ed1b2cec932097e86eee1bd9536ae7752e Mon Sep 17 00:00:00 2001 From: Sadam Hussain Date: Fri, 3 May 2024 09:54:38 -0400 Subject: [PATCH] Upload files to "/" --- HEMS.txt | 2823 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 2823 insertions(+) create mode 100644 HEMS.txt diff --git a/HEMS.txt b/HEMS.txt new file mode 100644 index 0000000..2666e54 --- /dev/null +++ b/HEMS.txt @@ -0,0 +1,2823 @@ +# HERE WE ARE CONSIDERING DIFFEREN TIME of arival and depature of EV, + +import itertools +import pyomo.environ as py +import pandas as pd +import numpy as np +from matplotlib import pyplot as plt +from matplotlib import rc + +# to ignore the warning for dividing by 0 +np.seterr(divide='ignore', invalid='ignore') +m1 = py.ConcreteModel() + +# +sheet_data = pd.read_excel("input_data.xlsx", sheet_name='data') +sheet_tr = pd.read_excel("input_data.xlsx", sheet_name='Ambient_Temp') +sheet_EV = pd.read_excel("input_data.xlsx", sheet_name='EV_100') +# sheet_EV = pd.read_excel("input_data.xlsx", sheet_name='EV_10') +# sheet_EV = pd.read_excel("input_data.xlsx", sheet_name='EV_50') +# sheet_EV = pd.read_excel("input_data.xlsx", sheet_name='EV_80') +sheet_EWH = pd.read_excel("input_data.xlsx", sheet_name='EWH') +sheet_wateruse = pd.read_excel("input_data.xlsx", sheet_name='Water_Use') +sheet_home = pd.read_excel("input_data.xlsx", sheet_name='Home') +sheet_ESS = pd.read_excel("input_data.xlsx", sheet_name='ESS_100') +# sheet_ESS = pd.read_excel("input_data.xlsx", sheet_name='ESS_10') +# sheet_ESS = pd.read_excel("input_data.xlsx", sheet_name='ESS_50') +# sheet_ESS = pd.read_excel("input_data.xlsx", sheet_name='ESS_80') +sheet_PV = pd.read_excel("input_data.xlsx", sheet_name='PV_100') +# sheet_PV = pd.read_excel("input_data.xlsx", sheet_name='PV_10') +# sheet_PV = pd.read_excel("input_data.xlsx", sheet_name='PV_50') +# sheet_PV = pd.read_excel("input_data.xlsx", sheet_name='PV_80') +sheet_load = pd.read_excel("input_data.xlsx", sheet_name='Base_Load') +sheet_bus = pd.read_excel("input_data.xlsx", sheet_name='Bus') +sheet_bus1 = pd.read_excel("input_data.xlsx", sheet_name='bus1') +# sheet_bus1 = pd.read_excel("input_data.xlsx", sheet_name='bus1') +sheet_sumdict = pd.read_excel("input_data.xlsx", sheet_name='sum_dict') +sheet_ambient = pd.read_excel("input_data.xlsx", sheet_name='Ambient_Temp') + +# + + +# + +# +m1.t = py.Set(initialize=sheet_data.time, ordered=True, doc='time period') +m1.b = py.Set(initialize=sheet_bus.bus, doc='bus') +m1.b1 = py.Set(initialize=sheet_bus1.bus1, doc='bus') +# 6:30 = 39 +m1.h = py.Set(initialize=sheet_home.home) + +m1.cost = py.Param(m1.h, mutable=True) + +# + +# +# Base load : uncontrollable loads +m1.base_load = py.Param(m1.h, m1.t, + initialize=dict(zip(list(itertools.product(m1.h.data(), m1.t.data())), + [i for x in sheet_load.columns if "load" in x for i in + sheet_load.loc[:, x].values])), doc="Base Load") + +m1.PV = py.Param(m1.h, m1.t, initialize=dict(zip(list(itertools.product(m1.h.data(), m1.t.data())), + [i for x in sheet_PV.columns if "PV" in x for i in + sheet_PV.loc[:, x].values])), doc="PV Production") +m1.max = py.Param(initialize=1000000, doc='maximum value for selling and buying power') +# m1.Transformer_Limit = py.Param(initialize=75, doc='Power limit of transformer') + +m1.UDC = py.Param(m1.t, mutable=True, initialize=80, doc='Uniform Distribution Capacity of transformer') +m1.Buy_Price_p2p = py.Param(m1.h, m1.t, initialize=0, mutable=True) +m1.combine = py.Param(m1.t, initialize=0, mutable=True) +m1.Sell_Price_p2p = py.Param(m1.h, m1.t, initialize=0, mutable=True) +m1.P_Buy_extra = py.Param(m1.h, m1.t, initialize=0) +m1.P_Buy_extra_1 = py.Param(m1.h, m1.t, initialize=0) +m1.P_sell_extra = py.Param(m1.h, m1.t, initialize=0) +m1.P_sell_extra_1 = py.Param(m1.h, m1.t, initialize=0) +m1.Sell_price_p2p = py.Param(m1.h, m1.t, initialize=0) +# Calculating the uniform distribution capacity +# m1.Transformer_Limit = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Transformer_limit)), +# doc='Maximum power limit on transformer') + +# for i in m1.t: +# m1.UDC[i] = m1.Transformer_Limit[i] / m1.h[-1] + +# + +# +# m1.alpha = py.Param(initialize=37, doc='time of arrival ') # (40 = 7:00 pm) 37 = 6 pm +m1.alpha = py.Param(m1.h, initialize=dict(zip(list(itertools.product(m1.h.data())), + [i for x in sheet_EV.columns if "arrival" in x for i + in + sheet_EV.loc[:, x].values])), doc="Arrival time of EV") + +m1.beta = py.Param(m1.h, initialize=dict(zip(list(itertools.product(m1.h.data())), + [i for x in sheet_EV.columns if "departure" in x for i + in + sheet_EV.loc[:, x].values])), doc="departure time of EV") +# m1.beta = py.Param(initialize=92, doc='time of departure ') # (92 = 8:00 am) +m1.t_final = py.Param(initialize=96, doc='final time of the day ') # 96 +m1.second = py.Param(initialize=2, doc='second time of the day ') +m1.EV_Dur = py.Set(m1.h, initialize=lambda m, h: range(m.alpha[h], m.beta[h] + 1), doc='time interval for EV') +# m1.EV_Dur = py.Set(m1.h, initialize=py.RangeSet((m1.alpha[h].value, m1.beta[h].value) for h in m1.h), +# doc='time interval for EV') +# m1.EV_Dur = py.Set(initialize=((h, t) for h in m1.h for t in range(m1.alpha[h], m1.beta[h]+1)), doc='time interval for EV') + +# m1.EV_Dur = py.Set(initialize={(i,j) for i in m1.h for j in range(m1.alpha[i].value, m1.beta[i].value + 1)}, +# dimen=2, +# doc='time interval for EV') + + +m1.EV_Dur1 = py.Set(m1.h, initialize=lambda m, h: range(m.alpha[h] + 1, m.beta[h] + 1), doc='time interval for EV') +# +# +# def init_EV_Dur1(m, h): +# return range(m.alpha[h] + 1, m.beta[h] + 1) +# +# m1.EV_Dur1 = py.Set(m1.h, initialize=init_EV_Dur1, ordered=True, dimen=1, within=OrderedSimpleSet, doc='time interval for EV') + + +# m1.EV_Dur1 = py.Set(m1.h,initialize=py.RangeSet(m1.alpha[h].value + 1, m1.beta[h].value), +# doc='time interval for EV after arrival') +m1.t_second = py.Set(initialize=py.RangeSet(m1.second.value, m1.t_final.value), + doc='time greater then 1') +m1.t_first = py.Param(initialize=1, doc='first time of the day ') +# m1.ist_interval = py.Set(initialize=py.RangeSet(m1.t_first.value, m1.alpha - 1), +# doc='time period before the arrival time') + +m1.ist_interval = py.Set(m1.h, initialize=lambda m, h: range(m1.t_first.value, m.alpha[h]), doc='time interval for EV') +# m1.second_interval = py.Set(initialize=py.RangeSet(m1.beta.value + 1, m1.t_final.value), +# doc='time interval after departure') + +m1.second_interval = py.Set(m1.h, initialize=lambda m, h: range(m.beta[h] + 1, m1.t_final.value + 1), + doc='time interval after departure') + +m1.ChargeRate_EV = py.Param(m1.h, initialize=dict(zip(list(itertools.product(m1.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") + +m1.DischRate_EV = py.Param(m1.h, initialize=dict(zip(list(itertools.product(m1.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") + +m1.Ch_Effc_EV = py.Param(m1.h, initialize=dict(zip(list(itertools.product(m1.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") + +m1.DischEffc_EV = py.Param(m1.h, initialize=dict(zip(list(itertools.product(m1.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") + +m1.Cap_EV = py.Param(m1.h, initialize=dict(zip((m1.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") +m1.End_percentage_EV = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.In_Percentage_EV = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.Energy_EV_dep = py.Param(m1.h, initialize=dict( + zip(m1.h.data(), np.array(m1.Cap_EV.values()) * m1.End_percentage_EV.values())), # just changing np>py + doc="Departure temperature of EV") + +m1.Energy_EV_In = py.Param(m1.h, initialize=dict( + zip(m1.h.data(), np.array(m1.Cap_EV.values()) * m1.In_Percentage_EV.values())), + doc="Initial energy of EV") + +# + +# + +m1.ChargeRate_ESS = py.Param(m1.h, initialize=dict(zip(list(itertools.product(m1.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") + +m1.DischRate_ESS = py.Param(m1.h, initialize=dict(zip(list(itertools.product(m1.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") + +# m1.ChargeRate_ESS = py.Param(initialize=float(sheet_ESS.charging_rate1), doc='Charging rate of ESS ') +# m1.DischRate_ESS = py.Param(initialize=float(sheet_ESS.discharging_rate1), +# doc='Discharging rate of ESS ') +m1.Cap_ESS = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.End_percentage_ESS = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.In_Percentage_ESS = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.Ch_Effc_ESS = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.DischEffc_ESS = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.Energy_ESS_In = py.Param(m1.h, initialize=dict( + zip(m1.h.data(), np.array(m1.Cap_ESS.values()) * m1.In_Percentage_ESS.values())), + doc="Initial energy of ESS") + +m1.End_En_ESS = py.Param(m1.h, initialize=dict( + zip(m1.h.data(), np.array(m1.Cap_ESS.values()) * m1.End_percentage_ESS.values())), + doc="Departure energy of ESS") + +# + +# + +m1.tetta_low = py.Param(m1.h, initialize=dict(zip((m1.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") +m1.tetta_up = py.Param(m1.h, initialize=dict(zip((m1.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") +m1.tetta_amb_int = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.tetta_amb = py.Param(m1.h, m1.t, + initialize=dict(zip(list(itertools.product(m1.h.data(), m1.t.data())), + [i for x in sheet_ambient.columns if "celsius" in x for i in + sheet_ambient.loc[:, x].values])), doc="Outdoor temperature") + +m1.Q = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.R = py.Param(m1.h, initialize=dict(zip((m1.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") +m1.C = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.M = py.Param(m1.h, initialize=dict(zip((m1.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)") + +m1.tetta_EWH_int = py.Param(m1.h, initialize=dict(zip((m1.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") + +m1.water_use = py.Param(m1.h, m1.t, + initialize=dict(zip(list(itertools.product(m1.h.data(), m1.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") + +# +# m1.Buy_price = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_RTP)), doc='Buying Price') +# m1.Buy_price = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_TOU)), doc='Buying Price') +# m1.Buy_price = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_TOU_winter)), doc='Buying Price') +# m1.Buy_price = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_flex_D)), doc='Buying Price') +m1.Buy_price = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.TOU_PGE)), doc='Buying Price') +# m1.Buy_price = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.TOU_Texas_RTP)), doc='Buying Price') +# m1.Buy_price_Peer = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_TOU)), doc='Buying Price') +# m1.Buy_price_Peer = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_Peer)), +# doc='Buying Price of peer') +# m1.Sell_price_Peer = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Sell_price_new)), +# doc='Selling Price of peer') +m1.Sell_price = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.FiT)), + doc='Selling Price') # Fit is TX US +# m1.Sell_price_Peer = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.fixed_incentive)), doc='Selling Price') +# Time duration: we took 15 mint granularity so in one hour it will be 1/4 +m1.time_d = py.Param(initialize=(1 / 4), doc='time duration ') # 1/4 +m1.Export_Percent = py.Param(initialize=1, doc='pecentage export ') # 0.814 + +m1.cost = py.Param(m1.h, mutable=True) +# + +for h in m1.h: + for t in m1.t: + m1.Buy_Price_p2p[h, t] = m1.Buy_price[t] +for h in m1.h: + for t in m1.t: + m1.Sell_Price_p2p[h, t] = m1.Sell_price[t] +# m1.max_limit = py.Param(m1.t, initialize=75, doc='Maximum power limit on transformer ') +# m1.DPT = py.Param(initialize=0.069, doc='Daily power tariff ') # .069 +# + +print('Code Starts for HEMS for cost') +# Variable + +# +m1.P_Buy_Grid = py.Var(m1.h, m1.t, bounds=(0, None)) + +m1.E_Buy_Grid = py.Param(m1.h, mutable=True) +m1.P_Buy_Total = py.Var(m1.h, m1.t, bounds=(0, None)) +# m1.P_Buy_Peer = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.P_Sell_Grid = py.Var(m1.h, m1.t, bounds=(0, 1000000)) +m1.P_Sell_Total = py.Var(m1.h, m1.t, bounds=(0, None)) +# m1.P_Sell_Peer = py.Var(m1.h, m1.t, bounds=(0, None)) + +m1.S_P_sell = py.Var(m1.h, m1.t, within=py.Binary) +m1.S_P_buy = py.Var(m1.h, m1.t, within=py.Binary) + +# m1.pflex = py.Var(m1.b, m1.t, doc='active Flexibility') +# m1.qflex = py.Var(m1.b, m1.t, +# doc='reactive Flexibility') # i comment this because i use m1.pflex variable * 0.48 + +# Tranformal variable +m1.Power_T_P = py.Var(m1.t, bounds=(0, None), doc="power from transformer to Peer") +m1.Power_Total_power = py.Param(m1.t, mutable=True, doc="power from transformer to Peer") +m1.HST = py.Var(m1.t, bounds=(0, None), doc="HST") +m1.DTOT = py.Var(m1.t, bounds=(0, None), doc="HST") +m1.DHST = py.Var(m1.t, bounds=(0, None), doc="HST") +m1.DTOTU = py.Var(m1.t, bounds=(0, None), doc="HST") +m1.DHSTU = py.Var(m1.t, bounds=(0, None), doc="HST") +m1.Power_P_T = py.Var(m1.t, bounds=(0, None), doc="power from Peer to Transformer") +m1.S_P_T = py.Var(m1.h, m1.t, within=py.Binary) +# + +# +m1.P_EV_Charge = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.P_EV_Disch = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.P_EV_Disch_Home = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.P_EV_Disch_Grid = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.S_EV_Charge = py.Var(m1.h, m1.t, within=py.Binary) +m1.S_EV_Disch = py.Var(m1.h, m1.t, within=py.Binary) +m1.Energy_EV = py.Var(m1.h, m1.t, bounds=(0, None)) # SOC of the EV +# + +# +m1.Energy_ESS = py.Var(m1.h, m1.t, bounds=(0, None)) # SOC of the ESS +m1.S_ESS_Charge = py.Var(m1.h, m1.t, within=py.Binary) +m1.P_ESS_Disch = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.P_ESS_Charge = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.P_ESS_Charge_Grid = py.Var(m1.h, m1.t, bounds=(0, None), doc='ESS charging from the Grid') +m1.P_ESS_Disch_Home = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.P_ESS_Disch_Grid = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.S_ESS_Disch = py.Var(m1.h, m1.t, within=py.Binary) +# + +# +m1.PV_Home = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.PV_Grid = py.Var(m1.h, m1.t, bounds=(0, None)) +m1.PV_Battery = py.Var(m1.h, m1.t, bounds=(0, None)) +# + +# +m1.tetta_EWH_wat = py.Var(m1.h, m1.t) +m1.S_EWH = py.Var(m1.h, m1.t, within=py.Binary) +m1.P_EWH = py.Var(m1.h, m1.t, doc='power of EWH') + + +# +# Constraints + + +# +def Power_buy(m1, h, i): + return m1.P_Buy_Grid[h, i] == m1.base_load[h, i] + m1.Q[h] * m1.S_EWH[h, i] + m1.P_EV_Charge[h, i] - \ + m1.P_EV_Disch_Home[h, i] + m1.P_ESS_Charge[h, i] - \ + m1.P_ESS_Disch_Home[h, i] - m1.PV_Home[h, i] - m1.PV_Battery[h, i] # + + +m1.Const_1 = py.Constraint(m1.h, m1.t, rule=Power_buy, doc='Power buy from the Grid') + + +def Power_buy2(m1, h, i): + return m1.P_Buy_Grid[h, i] <= m1.max * m1.S_P_buy[h, i] + + +m1.Const_1a = py.Constraint(m1.h, m1.t, rule=Power_buy2, + doc='removing the nonlinearity in the objective fucntion') + + +def Power_sell1(m1, h, i): + return m1.P_Sell_Grid[h, i] == m1.P_EV_Disch_Grid[h, i] + m1.PV_Grid[h, i] + m1.P_ESS_Disch_Grid[h, i] + + +m1.Const_2 = py.Constraint(m1.h, m1.t, rule=Power_sell1, doc='Power sell to the Grid') + + +def Power_sell2(m1, h, i): + # return m1.P_Sell_Grid[h, i] <= m1.max.value * (1 - m1.S_P_buy[h, i]) + return m1.P_Sell_Grid[h, i] <= m1.max.value * (1 - m1.S_P_buy[h, i]) + + +m1.Const_2a = py.Constraint(m1.h, m1.t, rule=Power_sell2, doc='Power sell to the Grid') + + +# put some percentage to export power the grid +def percentage(m1, h, i): + return sum(m1.P_Sell_Grid[h, i] for i in m1.t for h in m1.h) >= m1.Export_Percent.value * sum( + m1.P_EV_Disch_Home[h, i] + m1.P_ESS_Disch_Home[h, i] + + m1.PV_Battery[h, i] + m1.PV_Home[h, i] for i in m1.t for h in m1.h) + + # return sum(m1.P_Sell_Grid[y, i] for i in m1.t for y in m1.h) >= m1.Export_Percent.value * sum( + # m1.P_EV_Disch_Home[y, i] + m1.P_ESS_Disch_Home[y, i] + # + m1.PV_Battery[y, i] + m1.PV_Home[y, i] for i in m1.t for y in m1.h) + + +m1.connn1 = py.Constraint(m1.h, m1.t, rule=percentage, doc='Export constraint') + + +# +# def Status_Power(m1, h, i): +# return m1.S_P_buy[h, i] + m1.S_P_sell[h, i] <= 1 +# +# +# m1.Const_3 = py.Constraint(m1.h, m1.t, rule=Status_Power, +# doc='Buying and selling power will not occur at same time') +# + +# def Power_total_buy(m1, h, i): +# return m1.P_Buy_Total[h, i] == m1.P_Buy_Grid[h, i] +# # return m1.P_Buy_Total[h, i] == m1.P_Buy_Grid[h, i] + m1.P_Buy_Peer[h, i] +# +# +# m1.Const_3a = py.Constraint(m1.h, m1.t, rule=Power_total_buy, +# doc='Total Buying power equal to power buy from grid and form ' +# 'other peers') + + +# +# def Power_total_Sell(m1, h, i): +# return m1.P_Sell_Total[h, i] == m1.P_Sell_Grid[h, i] +# # return m1.P_Sell_Total[h, i] == m1.P_Sell_Grid[h, i] + m1.P_Sell_Peer[h, i] +# +# +# m1.Const_3b = py.Constraint(m1.h, m1.t, rule=Power_total_Sell, +# doc='Total Selling power equal to power sell to grid and to ' +# 'other peers') + + +# +# def Power_total(m1, h, i): +# return sum(m1.P_Buy_Peer[h, i] for h in m1.h) == sum(m1.P_Sell_Peer[h, i] for h in m1.h) +# +# +# m1.Const_3c = py.Constraint(m1.h, m1.t, rule=Power_total, +# doc='Total Buying power equal to total Selling power') + + +def Power_T_P(m1, h, i): + return m1.Power_T_P[i] == sum(m1.P_Buy_Grid[h, i] for h in m1.h) + + +# +m1.Const_3d = py.Constraint(m1.h, m1.t, rule=Power_T_P, + doc='Power bought from the grid is equal to the Power of transformer to grid direction') + + +def Power_P_T(m1, h, i): + return m1.Power_P_T[i] == sum(m1.P_Sell_Grid[h, i] for h in m1.h) + + +m1.Const_3e = py.Constraint(m1.h, m1.t, rule=Power_P_T, + doc='Power sell to the grid is equal to the Power on the transformer from grid direction') + + +# Transformer limits + +# def Power_limit_T_P(m1, h, i): +# return m1.Power_T_P[i] <= m1.max.value * (m1.S_P_T[h, i]) +# +# +# m1.Const_3f = py.Constraint(m1.h, m1.t, rule=Power_limit_T_P, doc='Power limit on the transformer') +# +# +# def Power_limit_P_T(m1, h, i): +# return m1.Power_P_T[i] <= m1.max.value * (1 - m1.S_P_T[h, i]) +# +# +# m1.Const_3g = py.Constraint(m1.h, m1.t, rule=Power_limit_P_T, doc='Power limit on the transformer') + + +# def Power_UDC_Limit(m1, h, i): +# return m1.P_Buy_Grid[h, i] <= m1.UDC[i] +# +# +# m1.Const_3h = py.Constraint(m1.h, m1.t, rule=Power_UDC_Limit, doc='UDC limit') + + +# + +# +# +def Power_EV_Charge_limit1(m1, h, i): + return m1.P_EV_Charge[h, i] <= m1.ChargeRate_EV[h] * m1.S_EV_Charge[h, i] + + +m1.Const_EV_1 = py.Constraint(m1.h, m1.t, rule=Power_EV_Charge_limit1, + doc='Charging power of EV (Upper limit)') + + +def Power_EV_Charge_limit2(m1, h, i): + return m1.P_EV_Charge[h, i] >= 0 + + +m1.Const_EV_2 = py.Constraint(m1.h, m1.t, rule=Power_EV_Charge_limit2, + doc='Charging power of EV (lower limit)') + + +def Power_EV_Disch_limit1(m1, h, i): + return m1.P_EV_Disch[h, i] <= (m1.DischRate_EV[h] * m1.S_EV_Disch[h, i]) + + +m1.Const_EV_3 = py.Constraint(m1.h, m1.t, rule=Power_EV_Disch_limit1, + doc='Discharging power of EV (Upper limit)') + + +def Power_EV_Disch_limit2(m1, h, i): + return m1.P_EV_Disch[h, i] >= 0 + + +m1.Const_EV_4 = py.Constraint(m1.h, m1.t, rule=Power_EV_Disch_limit2, + doc='Discharging power of EV (Lower limit)') + + +def Status_EV(m1, h, i): + return m1.S_EV_Disch[h, i] + m1.S_EV_Charge[h, i] <= 1 + + +m1.Const_EV_5 = py.Constraint(m1.h, m1.t, rule=Status_EV, + doc='Charging and discharging will not occur at same time') + + +def Power_EV_Disch(m1, h, i): + return m1.P_EV_Disch[h, i] * m1.DischEffc_EV[h] == ( + (m1.P_EV_Disch_Home[h, i] + m1.P_EV_Disch_Grid[h, i])) + + +m1.Const_EV_6 = py.Constraint(m1.h, m1.t, rule=Power_EV_Disch, doc='Discharging power of EV to ' + 'Home and Grid') + + +def SoC_EV1(m1, h, i): + return m1.Energy_EV[h, i] == m1.Energy_EV_In[h] + + +m1.Const_EV_7a = py.Constraint([(h, py.value(m1.alpha[h])) for h in m1.h], rule=SoC_EV1, doc='SoC of the EV at arrival') + + +# def SoC_EV1(m1, h): +# alpha_h = m1.alpha[h] +# return m1.Energy_EV[h, alpha_h] == m1.Energy_EV_In[h] +# +# m1.Const_EV_7a = py.Constraint(m1.h, rule=SoC_EV1, doc='SoC of the EV at arrival') + +def SoC_EV2(m1, h, i): + return m1.Energy_EV[h, i] == m1.Energy_EV[h, i - 1] + ( + m1.P_EV_Charge[h, i] * m1.Ch_Effc_EV[h] * m1.DischEffc_EV[h] - m1.P_EV_Disch[h, i]) * m1.time_d + + +m1.Const_EV_7b = py.ConstraintList() +for h in m1.h: + for i in m1.EV_Dur1[h].value: + m1.Const_EV_7b.add(SoC_EV2(m1, h, i)) + + +def EV_availability1(m1, h, i): + return m1.Energy_EV[h, i] == 0 + + +# m1.Const_EV_8a = py.Constraint(m1.h, m1.ist_interval[h], rule=EV_availability1,doc='SOC available before arrival time') +# m1.Const_EV_8a = py.Constraint([(h, m1.ist_interval[h].value) for h in m1.h], rule=EV_availability1, doc='SOC available before arrival time') +m1.Const_EV_8a = py.ConstraintList() +for h in m1.h: + for i in m1.ist_interval[h].value: + m1.Const_EV_8a.add(EV_availability1(m1, h, i)) + + +def EV_availability2(m1, h, i): + return m1.Energy_EV[h, i] == 0 + + +# m1.Const_EV_8b = py.Constraint([(h, m1.second_interval[h].value) for h in m1.h], rule=EV_availability2, doc='SOC available after departure time') +# m1.Const_EV_8b = py.Constraint(m1.h, m1.second_interval[h], rule=EV_availability2,doc='SOC available after departure time') +m1.Const_EV_8b = py.ConstraintList() +for h in m1.h: + for i in m1.second_interval[h].value: + m1.Const_EV_8b.add(EV_availability2(m1, h, i)) + + +def EV_status_available1(m1, h, i): + return m1.S_EV_Disch[h, i] + m1.S_EV_Charge[h, i] == 0 + + +# m1.Const_EV_9a = py.Constraint(m1.h, m1.ist_interval[h], rule=EV_status_available1, doc='EV availability before arrival time') +# m1.Const_EV_9a = py.Constraint([(h, m1.ist_interval[h].value) for h in m1.h], rule=EV_status_available1, doc='EV availability before arrival time') +m1.Const_EV_9a = py.ConstraintList() +for h in m1.h: + for i in m1.ist_interval[h].value: + m1.Const_EV_9a.add(EV_status_available1(m1, h, i)) + + +def EV_status_available2(m1, h, i): + return m1.S_EV_Disch[h, i] + m1.S_EV_Charge[h, i] == 0 + + +# m1.Const_EV_9b = py.Constraint(m1.h, m1.second_interval[h], rule=EV_status_available2,doc='EV availability after arrival time') +# m1.Const_EV_9b = py.Constraint([(h, m1.second_interval[h].value) for h in m1.h], rule=EV_status_available2, doc='EV availability after arrival time') +m1.Const_EV_9b = py.ConstraintList() +for h in m1.h: + for i in m1.second_interval[h].value: + m1.Const_EV_9b.add(EV_status_available2(m1, h, i)) + + +def EV_SoC_limit1(m1, h, i): + return m1.Energy_EV[h, i] >= 0.2 * m1.Cap_EV[h] + + +# m1.Const_EV_10 = py.Constraint(m1.h, m1.EV_Dur[h], rule=EV_SoC_limit1, doc='Minimum SoC of EV') +# m1.Const_EV_10 = py.Constraint([(h, m1.EV_Dur[h].value) for h in m1.h], rule=EV_SoC_limit1, doc='Minimum SoC of EV') +m1.Const_EV_10 = py.ConstraintList() +for h in m1.h: + for i in m1.EV_Dur[h].value: + m1.Const_EV_10.add(EV_SoC_limit1(m1, h, i)) + + +def EV_SoC_limit2(m1, h, i): + return m1.Energy_EV[h, i] <= m1.Cap_EV[h] + + +# m1.Const_EV_11 = py.Constraint(m1.h, m1.t, rule=EV_SoC_limit2, doc='Maximum SoC of EV') +m1.Const_EV_11 = py.ConstraintList() +for h in m1.h: + for i in m1.EV_Dur[h].value: + m1.Const_EV_11.add(EV_SoC_limit2(m1, h, i)) + + +def EV_final_SoC(m1, h, i): + return m1.Energy_EV[h, i] == m1.Energy_EV_dep[h] + + +# m1.Const_EV_12 = py.Constraint(m1.h, [m1.beta[h]], rule=EV_final_SoC,doc='Final SoC of EV at departure time') +# m1.Const_EV_12 = py.Constraint([(h,py.value([m1.beta[h]])) for h in m1.h], rule=EV_final_SoC,doc='Final SoC of EV at departure time') +m1.Const_EV_12 = py.Constraint([(h, py.value(m1.beta[h])) for h in m1.h], rule=EV_final_SoC, + doc='Final SoC of EV at departure time') + + +# m1.Const_EV_12 = py.ConstraintList() +# for h in m1.h: +# for i in py.value(m1.beta[h]): +# m1.Const_EV_12.add(EV_final_SoC(m1, h, i)) + +# + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +def SoC_ESS1(m1, h, i): + return m1.Energy_ESS[h, i] == m1.Energy_ESS_In[h] + + +m1.Const_ESS_1a = py.Constraint(m1.h, [m1.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(m1, h, i): + return m1.Energy_ESS[h, i] == m1.Energy_ESS[h, i - 1] + ( + m1.P_ESS_Charge[h, i] * m1.Ch_Effc_ESS[h] - m1.P_ESS_Disch[h, i]) * m1.time_d + + +m1.Const_ESS_1b = py.Constraint(m1.h, m1.t_second, rule=SoC_ESS2, doc='SoC of ESS') + + +def Power_ESS_Charge(m1, h, i): + return m1.P_ESS_Charge[h, i] == m1.P_ESS_Charge_Grid[h, i] + m1.PV_Battery[h, i] + + +m1.Const_ESS_2 = py.Constraint(m1.h, m1.t, rule=Power_ESS_Charge, + doc='Charging power of ESS ') + + +def Power_ESS_Charge_limit1(m1, h, i): + return (m1.P_ESS_Charge_Grid[h, i] + m1.PV_Battery[h, i]) <= m1.ChargeRate_ESS[h] * m1.S_ESS_Charge[h, i] + + +m1.Const_ESS_2a = py.Constraint(m1.h, m1.t, rule=Power_ESS_Charge_limit1, + doc='Charging power of ESS (Upper limit)') + + +def Power_ESS_Charge_limit2(m1, h, i): + return (m1.P_ESS_Charge_Grid[h, i] + m1.PV_Battery[h, i]) >= 0 + # return (m1.P_ESS_Charge_Grid[h, i] + m1.PV_Battery[h, i]) >= 0.2 * (m1.ChargeRate_ESS[h] * m1.S_ESS_Charge[h, i]) + + +m1.Const_ESS_2b = py.Constraint(m1.h, m1.t, rule=Power_ESS_Charge_limit2, + doc='Charging power of ESS (Lower limit)') + + +def Power_ESS_Disch_limit1(m1, h, i): + return (m1.P_ESS_Disch_Home[h, i] + m1.P_ESS_Disch_Grid[h, i]) <= (m1.DischRate_ESS[h] * m1.S_ESS_Disch[h, i]) + + +m1.Const_ESS_4 = py.Constraint(m1.h, m1.t, rule=Power_ESS_Disch_limit1, + doc='Discharging power of ESS (Upper limit)') + + +def Power_ESS_Disch_limit2(m1, h, i): + # return (m1.P_ESS_Disch[h, i]) >= 0.2 * (m1.DischRate_ESS[h] * m1.S_ESS_Disch[h, i]) + return (m1.P_ESS_Disch_Home[h, i] + m1.P_ESS_Disch_Grid[h, i]) >= 0 + + +m1.Const_ESS_5 = py.Constraint(m1.h, m1.t, rule=Power_ESS_Disch_limit2, + doc='Discharging power of ESS (Lower limit)') + + +def Status_ESS(m1, h, i): + return m1.S_ESS_Disch[h, i] + m1.S_ESS_Charge[h, i] <= 1 + + +m1.Const_ESS_6 = py.Constraint(m1.h, m1.t, rule=Status_ESS, + doc='Charging and discharging will not occur at same time') + + +def Power_ESS_Disch(m1, h, i): + return (m1.P_ESS_Disch[h, i]) * m1.DischEffc_ESS[h] == (m1.P_ESS_Disch_Home[h, i] + m1.P_ESS_Disch_Grid[h, i]) + + +m1.Const_ESS_7 = py.Constraint(m1.h, m1.t, rule=Power_ESS_Disch, + doc='Discharging power of ESS to Home and Grid') + + +def ESS_SoC_limit1(m1, h, i): + return m1.Energy_ESS[h, i] >= 0.1 * m1.Cap_ESS[h] + + +m1.Const_ESS_8 = py.Constraint(m1.h, m1.t, rule=ESS_SoC_limit1, doc='Minimum SoC of ESS') + + +def ESS_SoC_limit2(m1, h, i): + return m1.Energy_ESS[h, i] <= m1.Cap_ESS[h] + + +m1.Const_ESS_9 = py.Constraint(m1.h, m1.t, rule=ESS_SoC_limit2, doc='Maximum SoC of ESS') + + +def ESS_final_SoC(m1, h, i): + return m1.Energy_ESS[h, i] >= m1.End_En_ESS[h] + + +m1.Const_ESS_10 = py.Constraint(m1.h, [m1.t_final.value], rule=ESS_final_SoC, + doc='Final SoC of ESS at departure time') + + +# + +# +def EWH_limit1(m1, h, i): + return m1.tetta_EWH_wat[h, i] <= m1.tetta_up[h] + + +m1.Const_EWH_1 = py.Constraint(m1.h, m1.t, rule=EWH_limit1, doc='Maximum Limit') + + +def EWH_limit2(m1, h, i): + return m1.tetta_EWH_wat[h, i] >= m1.tetta_low[h] + + +m1.Const_EWH_2 = py.Constraint(m1.h, m1.t, rule=EWH_limit2, doc='Minimum Limit') + + +def EWH_power(m1, h, i): + return m1.P_EWH[h, i] == m1.Q[h] * m1.S_EWH[h, i] + + +m1.Const_EWH_3 = py.Constraint(m1.h, m1.t, rule=EWH_power, doc='Electric water heater power') + + +def EWH_temp_1(m1, h, i): + return m1.tetta_EWH_wat[h, i] == m1.tetta_amb[h, i] + m1.Q[h] * m1.S_EWH[h, i] * m1.R[ + h] * m1.time_d - ( + (m1.M[h] - m1.water_use[h, i]) / m1.M[h]) * ( + m1.tetta_amb_int[h] - m1.tetta_EWH_int[h]) * py.exp(-m1.time_d / (m1.R[h] * m1.C[h])) + + +m1.Const_EWH_4 = py.Constraint(m1.h, [m1.t_first.value], rule=EWH_temp_1, doc='EWH m1') + + +def EWH_temp_2(m1, h, i): + return m1.tetta_EWH_wat[h, i] == m1.tetta_amb[h, i] + m1.Q[h] * m1.S_EWH[h, i] * m1.R[ + h] * m1.time_d - ( + (m1.M[h] - m1.water_use[h, i]) / m1.M[h]) * ( + m1.tetta_amb[h, i] - m1.tetta_EWH_wat[h, i - 1]) * py.exp( + -m1.time_d / (m1.R[h] * m1.C[h])) + + +m1.Const_EWH_5 = py.Constraint(m1.h, m1.t_second, rule=EWH_temp_2, doc='EWH m1') + + +# + +# # need to include the PV equation of omer paper + +def PV_production(m1, h, i): + return m1.PV_Grid[h, i] + m1.PV_Home[h, i] + m1.PV_Battery[h, i] == m1.PV[h, i] + # return m1.PV_Grid[h, i] + m1.PV_Battery[h, i] == m1.PV[h, i] + + +m1.Const_PV = py.Constraint(m1.h, m1.t, rule=PV_production, doc='PV Production') + + +# + + +# + + +def objective_rule(m1): + return sum( + m1.P_Buy_Grid[y, i] * m1.Buy_price[i] - m1.P_Sell_Grid[y, i] * m1.Sell_price[i] for i in + m1.t for y in m1.h) * m1.time_d.value + + # return sum( + # (m1.P_Buy_Grid[y, i] * m1.Buy_price[i]) * m1.S_P_buy[y, i] - (m1.P_Sell_Grid[y, i] * m1.Sell_price[i])* m1.S_P_sell[y, i] for i in + # m1.t for y in m1.h) * m1.time_d.value + + +m1.obj1 = py.Objective(rule=objective_rule, sense=py.minimize, doc='Definition of objective function') +# +# m1.write('m12.lp', io_options={'symbolic_solver_labels': True}) +opt = py.SolverFactory('gurobi') +# opt = py.SolverFactory('mindtpy').solve(m1, mip_solver='cplex', nlp_solver='ipopt', tee=True) +# opt = py.SolverFactory('couenne ') +opt.options["mipgap"] = 0.03 # 0.155 for load and 0.8 for other load +result = opt.solve(m1, tee=True) +# result = opt.solve(m1) + +# this is for calling the neos +# os.environ['NEOS_EMAIL'] = 'sadamengr15@gmail.com' +# opt = py.SolverManagerFactory('neos') +# result = opt.solve(m1, opt='cplex') + +# m1.pprint() +print(result) +for y in m1.h: + m1.cost[y] = py.value(sum(m1.P_Buy_Grid[y, i] * m1.Buy_price[i] - m1.P_Sell_Grid[y, i] * + m1.Sell_price[i] for i in m1.t) * m1.time_d.value) +for h in m1.h: + print( + f"Before Cost Objective of Home {h} : " f"{py.value((sum((m1.P_Buy_Grid[h, i] * m1.Buy_price[i]) - (m1.P_Sell_Grid[h, i] * m1.Sell_price[i]) for i in m1.t) * m1.time_d.value))}") +# + +# total power on the transformer: +for i in m1.t: + m1.Power_Total_power[i] = round(m1.Power_P_T[i].value + m1.Power_T_P[i].value, 3) + +# print('Sum of all homes objectives = ', py.value(m1.obj1)) +# for h in m1.h: +# print( +# f"Before Cost Objective of Hom {h} : " f"{py.value((sum((m1.P_Buy_Grid[h, i] * m1.Buy_price[i]) - (m1.P_Sell_Grid[h, i] * m1.Sell_price[i]) for i in m1.t) * m1.time_d))}") +# < / editor - fold > +# for h in m1.h: +# print( +# f"Energy of Home {h} : " f"{py.value((sum((m1.P_Buy_Grid[h, i] * m1.time_d) for i in m1.t)))}") +# for h in m1.h: +# m1.E_Buy_Grid[h] = ((sum((m1.P_Buy_Grid[h, i] * m1.time_d) for i in m1.t) )) +print('power export ', py.value(sum(m1.P_Sell_Grid[y, i] for i in m1.t for y in m1.h))) +print('power import ', py.value(sum(m1.P_Buy_Grid[y, i] for i in m1.t for y in m1.h))) + +# # # +## +# # # Plotting for automatic ploting +time = [i for i in m1.t] + +pl1_pbuy = [[] for j in m1.h] +pl1_pbuy_extra = [[] for j in m1.h] +pl1_pbuy_extra_1 = [[] for j in m1.h] +pl1_psell_extra = [[] for j in m1.h] +pl1_psell_extra_1 = [[] for j in m1.h] +pl_Buy_Price_p2p = [[] for j in m1.h] +pl_Sell_Price_p2p = [[] for j in m1.h] +pl1_pbuy_peer = [[] for j in m1.h] +pl1_psell = [[] for j in m1.h] +pl1_psell_peer = [[] for j in m1.h] + +PEVCharge = [[] for j in m1.h] +PESSChargeGrid = [[] for j in m1.h] +PEVDischHome = [[] for j in m1.h] +PESSDischHome = [[] for j in m1.h] +PESSCharge = [[] for j in m1.h] +SEVCharge = [[] for j in m1.h] +SEVDisch = [[] for j in m1.h] +PEVDischGrid = [[] for j in m1.h] +PESSDischGrid = [[] for j in m1.h] +PVHome = [[] for j in m1.h] +PVGrid = [[] for j in m1.h] +PVBattery = [[] for j in m1.h] +PPV = [[] for j in m1.h] +baseload = [[] for j in m1.h] +EnergyESS = [[] for j in m1.h] +EnergyEV = [[] for j in m1.h] +tettaEWHwat = [[] for j in m1.h] +SEWH = [[] for j in m1.h] +PEWH = [[] for j in m1.h] +Buyprice = [] # cost = [i for i in m1.c.values()] <<< check this +HST = [] # cost = [i for i in m1.c.values()] <<< check this +DHST = [] # cost = [i for i in m1.c.values()] <<< check this +DTOT = [] # cost = [i for i in m1.c.values()] <<< check this +Transformer_Limit = [] # cost = [i for i in m1.c.values()] <<< check this +Buyprice_peer = [] # cost = [i for i in m1.c.values()] <<< check this +Power_P_T = [] # cost = [i for i in m1.c.values()] <<< check this +Power_T_P = [] # cost = [i for i in m1.c.values()] <<< check this +Sellprice = [] +combine = [] +cost_1 = [[] for j in m1.h] +Sellprice_peer = [] + +for j in m1.h: + for k, v in m1.P_Buy_Grid.items(): + if k[0] == j: + pl1_pbuy[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_Buy_extra.items(): + if k[0] == j: + pl1_pbuy_extra[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_Buy_extra_1.items(): + if k[0] == j: + pl1_pbuy_extra_1[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_sell_extra.items(): + if k[0] == j: + pl1_psell_extra[j - 1].append(py.value(v)) + +for j in m1.h: + for k, v in m1.P_sell_extra_1.items(): + if k[0] == j: + pl1_psell_extra_1[j - 1].append(py.value(v)) + +for j in m1.h: + for k, v in m1.Buy_Price_p2p.items(): + if k[0] == j: + pl_Buy_Price_p2p[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.Sell_Price_p2p.items(): + if k[0] == j: + pl_Sell_Price_p2p[j - 1].append(py.value(v)) + +# for j in m1.h: +# for k, v in m1.P_Sell_Peer.items(): +# if k[0] == j: +# pl1_psell_peer[j - 1].append(py.value(v)) +# for j in m1.h: +# for k, v in m1.P_Buy_Peer.items(): +# if k[0] == j: +# pl1_pbuy_peer[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_Sell_Grid.items(): + if k[0] == j: + pl1_psell[j - 1].append(py.value(v)) + +for j in m1.h: + for k, v in m1.P_EV_Charge.items(): + if k[0] == j: + PEVCharge[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_EV_Disch_Home.items(): + if k[0] == j: + PEVDischHome[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_EV_Disch_Grid.items(): + if k[0] == j: + PEVDischGrid[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_ESS_Charge_Grid.items(): + if k[0] == j: + PESSChargeGrid[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_ESS_Disch_Home.items(): + if k[0] == j: + PESSDischHome[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_ESS_Charge.items(): + if k[0] == j: + PESSCharge[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_ESS_Disch_Grid.items(): + if k[0] == j: + PESSDischGrid[j - 1].append(py.value(v)) + +for j in m1.h: + for k, v in m1.PV_Home.items(): + if k[0] == j: + PVHome[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.PV_Grid.items(): + if k[0] == j: + PVGrid[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.PV_Battery.items(): + if k[0] == j: + PVBattery[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.PV.items(): + if k[0] == j: + PPV[j - 1].append(py.value(v)) + +for j in m1.h: + for k, v in m1.base_load.items(): + if k[0] == j: + baseload[j - 1].append(py.value(v)) + +for j in m1.h: + for k, v in m1.Energy_ESS.items(): + if k[0] == j: + EnergyESS[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.Energy_EV.items(): + if k[0] == j: + EnergyEV[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.P_EWH.items(): + if k[0] == j: + PEWH[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.tetta_EWH_wat.items(): + if k[0] == j: + tettaEWHwat[j - 1].append(py.value(v)) + +for j in m1.h: + for k, v in m1.S_EWH.items(): + if k[0] == j: + SEWH[j - 1].append(py.value(v)) +for j in m1.h: + for k, v in m1.cost.items(): + if k == j: + cost_1[j - 1].append(py.value(v)) +# for i in m1.Buy_price_Peer: +# Buyprice_peer.append(py.value(m1.Buy_price_Peer[i])) +for i in m1.Buy_price: + Buyprice.append(py.value(m1.Buy_price[i])) + +for i in m1.Sell_price: + Sellprice.append(py.value(m1.Sell_price[i])) +# for i in m1.Sell_price_Peer: +# Sellprice_peer.append(py.value(m1.Sell_price_Peer[i])) + +# # alone +# for k in range(len(m1.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, PESSChargeGrid[k], label='ESS Charging from Grid', color='r') +# ax[2, 0].bar(time, PVBattery[k], label='ESS Charging from PV', color='g') +# 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((m1.P_Buy_Grid[k + 1, i] * m1.Buy_price[i]) - (m1.P_Sell_Grid[k + 1, i] * m1.Sell_price[i]) for i in m1.t) * m1.time_d))}") +# pl.tight_layout() +# plt.show() +for k in m1.t: + m1.combine[k] = (PEVCharge[1][k - 1] + PEWH[1][k - 1] + baseload[1][k - 1] + PESSChargeGrid[1][k - 1]) +# # code for if you want to plot all the home togathoer +for i in m1.combine: + combine.append(py.value(m1.combine[i])) +# fig, ax = plt.subplots(6, len(m1.h), constrained_layout=True, sharex='col',sharey='row') # , sharey='row' +# gs = fig.add_gridspec(hspace=0, wspace=0) +# for k in range(len(m1.h)): +# ax[0, k].bar(time, pl1_pbuy[k], color="blue") +# ax[1, k].bar(time, pl1_psell[k], color="black") +# # ax[2, k].bar(time, pl1_pbuy_peer[k], color="green") +# # ax[3, k].bar(time, pl1_psell_peer[k], color="red") +# # ax[4, 0].plot(time, Transformer_Limit, color="red") +# # ax[4, 0].bar(time, Power_P_T, color="blue") +# # ax[4, 0].legend(['Tranformer_limit']) +# # ax[4, 1].plot(time, Transformer_Limit, color="red") +# # ax[4, 1].bar(time, Power_T_P, color="blue") +# # ax[4, 1].legend(['Power_T_P']) +# # # ax[0, k].plot(time, baseload[k], label='Base load') +# +# ax[2, k].plot(time, Buyprice, label='Buyprice_G', color="blue") +# # ax[4, 2].plot(time, Buyprice_peer, label='Buyprice_P', color="red") +# ax[2, k].plot(time, Sellprice, label='Sellprice_G', color="green") +# # ax[4, 2].plot(time, Sellprice_peer, label='Sellprice_P', color="black") +# ax[2, 0].legend(['Buyprice_G', 'Sellprice_G'],loc='upper left', fontsize='xx-small') +# # ax[4, 3].plot(time, HST, label='HST', color="red") +# # ax[4, 3].plot(time, DHST, label='DHST', color="blue") +# # ax[4, 3].plot(time, DTOT, label='DTOT', color="green") +# # ax[4, 3].legend(['HST', 'DHST', 'DTOT']) +# +# # ax[9, k].plot(time, PEWH[k], label='EWH Power') +# ax[3, k].plot(time, tettaEWHwat[k], label='EWH Temp') +# ax[3, 0].legend(ncol=3,loc='upper left', fontsize='xx-small') +# # ax[5, 0].legend(ncol=3,loc='best', fontsize='xx-small') +# # 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, 0].legend(ncol=3,loc='best', fontsize='xx-small') +# # 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, 0].legend(ncol=3,loc='best', fontsize='xx-small') +# # 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, 0].legend(ncol=3,loc='best', fontsize='xx-small') +# ax[4, k].plot(time, EnergyESS[k], label='Energy of ESS') +# ax[4, 0].legend(ncol=3,loc='upper left', fontsize='xx-small') +# ax[5, k].plot(time, EnergyEV[k], label='Energy of EV') +# ax[5, 0].legend(ncol=3,loc='upper left', fontsize='xx-small') +# ax[5, k].set_xlabel('Time (step)') +# # # ax[1, k].set_ylabel('Electricity price ($/W/h)') +# ax[0, 0].set_ylabel('Dem', fontsize='x-small') +# ax[1, 0].set_ylabel('Sel', fontsize='x-small') +# ax[0, k].set_title(f"Home {k + 1}", fontsize='xx-small') +# fig.suptitle('Home Profile before limit', fontsize=16) +# # # +# # pl.tight_layout() +# +# plt.show() + +# # if you want to plot each home alone in the figure +# # for k in range(len(m1.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') +# # # # + +# +# # + +# +# # + +# +# # # # +# # # +# # # # Plotting for automatic ploting +# # time = [i for i in m1.t] +# pl2_pbuy = [[] for j in m1.h] +# # pl2_psell = [[] for j in m1.h] +# # # # PEVCharge = [[] for j in m1.h] +# # # # PESSChargeGrid = [[] for j in m1.h] +# # # # PEVDischHome = [[] for j in m1.h] +# # # # PESSDischHome = [[] for j in m1.h] +# # # # SEVCharge = [[] for j in m1.h] +# # # # SEVDisch = [[] for j in m1.h] +# # # # PEVDischGrid = [[] for j in m1.h] +# # # # PESSDischGrid = [[] for j in m1.h] +# # # # PVHome = [[] for j in m1.h] +# # # # PVGrid = [[] for j in m1.h] +# # # # PVBattery = [[] for j in m1.h] +# # # # PPV = [[] for j in m1.h] +# # # # baseload = [[] for j in m1.h] +# # # # EnergyESS = [[] for j in m1.h] +# # # # EnergyEV = [[] for j in m1.h] +# # # # tettaEWHwat = [[] for j in m1.h] +# # # # SEWH = [[] for j in m1.h] +# # # # PEWH = [[] for j in m1.h] +# # # # Buyprice = [] # cost = [i for i in m1.c.values()] <<< check this +# # # # Sellprice = [] +# +# for j in m1.h: +# for k, v in m1.P_Buy_Grid.items(): +# if k[0] == j: +# pl2_pbuy[j - 1].append(py.value(v)) +# # for j in m1.h: +# # for k, v in m1.P_Sell_Grid.items(): +# # if k[0] == j: +# # pl2_psell[j - 1].append(py.value(v)) +# # # # +# # # # for j in m1.h: +# # # # for k, v in m1.P_EV_Charge.items(): +# # # # if k[0] == j: +# # # # PEVCharge[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.P_EV_Disch_Home.items(): +# # # # if k[0] == j: +# # # # PEVDischHome[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.P_EV_Disch_Grid.items(): +# # # # if k[0] == j: +# # # # PEVDischGrid[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.P_ESS_Charge_Grid.items(): +# # # # if k[0] == j: +# # # # PESSChargeGrid[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.P_ESS_Disch_Home.items(): +# # # # if k[0] == j: +# # # # PESSDischHome[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.P_ESS_Disch_Grid.items(): +# # # # if k[0] == j: +# # # # PESSDischGrid[j - 1].append(py.value(v)) +# # # # +# # # # for j in m1.h: +# # # # for k, v in m1.PV_Home.items(): +# # # # if k[0] == j: +# # # # PVHome[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.PV_Grid.items(): +# # # # if k[0] == j: +# # # # PVGrid[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.PV_Battery.items(): +# # # # if k[0] == j: +# # # # PVBattery[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.PV.items(): +# # # # if k[0] == j: +# # # # PPV[j - 1].append(py.value(v)) +# # # # +# # # # for j in m1.h: +# # # # for k, v in m1.base_load.items(): +# # # # if k[0] == j: +# # # # baseload[j - 1].append(py.value(v)) +# # # # +# # # # for j in m1.h: +# # # # for k, v in m1.Energy_ESS.items(): +# # # # if k[0] == j: +# # # # EnergyESS[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.Energy_EV.items(): +# # # # if k[0] == j: +# # # # EnergyEV[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.P_EWH.items(): +# # # # if k[0] == j: +# # # # PEWH[j - 1].append(py.value(v)) +# # # # for j in m1.h: +# # # # for k, v in m1.tetta_EWH_wat.items(): +# # # # if k[0] == j: +# # # # tettaEWHwat[j - 1].append(py.value(v)) +# # # # +# # # # for j in m1.h: +# # # # for k, v in m1.S_EWH.items(): +# # # # if k[0] == j: +# # # # SEWH[j - 1].append(py.value(v)) +# # # # +# # # # for i in m1.Buy_price: +# # # # Buyprice.append(value(m1.Buy_price[i])) +# # # # +# # # # for i in m1.Sell_price: +# # # # Sellprice.append(value(m1.Sell_price[i])) +# # # # +# # # # # alone +# # # for k in range(len(m1.h)): +# # # fig, ax = plt.subplots(5, 2, figsize=(10, 10)) +# # # ax[0, 0].bar(time, pl2_pbuy[k], label='Buying power-obj2') +# # # ax[0, 0].bar(time, pl1_pbuy[k], label='Buying power-obj1') +# # # ax[1, 0].bar(time, pl2_psell[k], label='Selling power-obj2') +# # # ax[1, 0].bar(time, pl1_psell[k], label='Selling power-obj1') +# # # ax[0, 0].legend(loc='best', fontsize='small', ncol=3) +# # # ax[1, 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, PESSChargeGrid[k], label='ESS Charging from Grid', color='r') +# # # # # ax[2, 0].bar(time, PVBattery[k], label='ESS Charging from PV', color='g') +# # # # # 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('Power1 (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(m1.P_Buy_Grid[k + 1, i] for i in m1.t) * m1.time_d))}") +# # # pl.tight_layout() +# # # plt.show() +# # # # +# # # # # code for if you want to plot all the home togathoer +# # fig, ax = plt.subplots(4, len(m1.h), figsize=(10, 10), constrained_layout=True, sharex='col', sharey='row') +# # gs = fig.add_gridspec(hspace=0, wspace=0) +# # for k in range(len(m1.h)): +# # ax[0, k].bar(time, pl1_pbuy[k], color='green') +# # ax[1, k].bar(time, pl2_pbuy[k], color='blue') +# # ax[2, k].bar(time, pl1_psell[k], color='green') +# # ax[3, k].bar(time, pl2_psell[k], label='blue') +# # # ax[0, k].legend(loc='best') +# # # ax[1, k].legend(loc='best') +# # # ax[2, k].legend(loc='best') +# # # ax[3, k].legend(loc='best') +# # # # ax[0, k].plot(time, baseload[k], label='Base load') +# # # # # 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, 0].set_ylabel( +# # 'Buy-Power1 (kW)') # showing only on the side of the graph if you want to put it on all need to use [0,k] +# # ax[1, 0].set_ylabel('Buy-Power2 (kW)') +# # ax[2, 0].set_ylabel('Sell-Power1 (kW)') +# # ax[3, 0].set_ylabel('Sell-Power2 (kW)') +# # +# # ax[0, k].set_title(f"Home {k + 1}", fontsize='xx-small') +# # fig.suptitle('Cost and Energy minimization', fontsize=16) +# # +# # plt.show() +# # # +# # # # if you want to plot each home alone in the figure +# # # # for k in range(len(m1.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') +# # # +# +# writer = pd.ExcelWriter('Result_HEMS.xlsx', engine='xlsxwriter') +# time_excel = pd.DataFrame({'Time_Step': time}) +# time_excel.to_excel(writer, sheet_name='Power', startcol=0, index=False, header=True) +# for h in range(len(m1.h)): +# pload_excel = pd.DataFrame({f"Power_Demand {h + 1}": pl1_pbuy[h]}) +# pload_excel.to_excel(writer, sheet_name='Power', startcol=h + 1, index=False, header=True) +# +# pload_excel = pd.DataFrame({f"Power_Sell {h + 1}": pl1_psell[h]}) +# pload_excel.to_excel(writer, sheet_name='Power', startcol=h + 4, index=False, header=True) +# writer.save() + + +# fig_HEMS, ax = plt.subplots(5, len(m1.h), constrained_layout=True, sharex='col', sharey='row') # , sharey='row' +# gs = fig_HEMS.add_gridspec(hspace=0, wspace=0) +# for k in range(len(m1.h)): +# ax[0, k].bar(time, pl1_pbuy[k], color="blue") +# ax[1, k].bar(time, pl1_psell[k], color="black") +# # ax[2, k].bar(time, pl1_pbuy_peer[k], color="green") +# # ax[3, k].bar(time, pl1_psell_peer[k], color="red") +# # ax[4, 0].plot(time, Transformer_Limit, color="red") +# # ax[4, 0].bar(time, Power_P_T, color="blue") +# # ax[4, 0].legend(['Tranformer_limit']) +# # ax[4, 1].plot(time, Transformer_Limit, color="red") +# # ax[4, 1].bar(time, Power_T_P, color="blue") +# # ax[4, 1].legend(['Power_T_P']) +# # # ax[0, k].plot(time, baseload[k], label='Base load') +# +# # ax[2, k].plot(time, Buyprice, label='Buyprice_G', color="blue") +# # ax[4, 2].plot(time, Buyprice_peer, label='Buyprice_P', color="red") +# # ax[2, k].plot(time, Sellprice, label='Sellprice_G', color="green") +# # ax[4, 2].plot(time, Sellprice_peer, label='Sellprice_P', color="black") +# # ax[2, 0].legend(['Buyprice_G', 'Sellprice_G'],loc='upper left', fontsize='xx-small') +# # ax[4, 3].plot(time, HST, label='HST', color="red") +# # ax[4, 3].plot(time, DHST, label='DHST', color="blue") +# # ax[4, 3].plot(time, DTOT, label='DTOT', color="green") +# # ax[4, 3].legend(['HST', 'DHST', 'DTOT']) +# +# # ax[9, k].plot(time, PEWH[k], label='EWH Power') +# ax[2, k].plot(time, tettaEWHwat[k], label='EWH Temp') +# ax[2, 0].legend(ncol=3, loc='upper left', fontsize='xx-small') +# # ax[5, 0].legend(ncol=3,loc='best', fontsize='xx-small') +# # 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, 0].legend(ncol=3,loc='best', fontsize='xx-small') +# # 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, 0].legend(ncol=3,loc='best', fontsize='xx-small') +# # 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, 0].legend(ncol=3,loc='best', fontsize='xx-small') +# ax[3, k].plot(time, EnergyESS[k], label='Energy of ESS') +# ax[3, 0].legend(ncol=3, loc='upper left', fontsize='xx-small') +# ax[4, k].bar(time, EnergyEV[k], label='Energy of EV') +# ax[4, 0].legend(ncol=3, loc='upper left', fontsize='xx-small') +# ax[4, k].set_xlabel('Time (step)') +# # # ax[1, k].set_ylabel('Electricity price ($/W/h)') +# ax[0, 0].set_ylabel('Dem', fontsize='x-small') +# ax[1, 0].set_ylabel('Sel', fontsize='x-small') +# ax[0, k].set_title(f"Home {k + 1}", fontsize='xx-small') +# fig_HEMS.suptitle('Home Profile', fontsize=16) +# # # +# # pl.tight_layout() +# +# plt.show() + +# stack plot of power demand of one house : +# Generate some example data +# Plot stacked bar chart +fig_all, axs = plt.subplots(5, 2, figsize=(8, 8),constrained_layout=True, sharex='col', sharey='row') +labels = ('EV_Charging', 'Power_EWH', 'Baseload', 'Battery_Charging') +data_lists = [(PEVCharge[i], PEWH[i], baseload[i], PESSChargeGrid[i]) for i in range(10)] +rc('font', family='Times New Roman') +for i, ax in enumerate(axs.flat): + bottom = np.zeros(len(time)) + ax.set_xlim([0, 98]) + ax.set_ylim([0, 18]) + ax.set_xticks(np.arange(0, 98, 6)) + # ax.set_yticks(np.arange(0, 24, 4)) + for label, data in zip(labels, data_lists[i]): + ax.bar(time, data, bottom=bottom, label=label) + bottom += data + # if i == 3 or i == 6 or i ==9: + # ax.set_title(f'Consumer-{i+1}', fontname='Times New Roman', fontsize=12) + # else: + ax.set_title(f'Prosumer-{i + 1}', fontname='Times New Roman', fontsize=12) + ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) + if i == 0 or i == 2 or i == 4 or i == 6 or i == 8: + ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) + ax.set_yticks(np.arange(0, 18, 4)) + else: + ax.set_ylabel('') + if i == 8 or i == 9: + ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman', fontsize=12) + ax.tick_params(axis='x', which='both', labelsize=12) + else: + ax.tick_params(axis='x', which='both', bottom=True, labelbottom=False) # modified + ax.set_xticks([]) # added +handles, labels = ax.get_legend_handles_labels() + +fig_all.legend(handles, labels, loc='lower center', ncol=4, fontsize=12) +plt.tight_layout() +fig_all.subplots_adjust(bottom=0.129) +plt.show() +# +# fig_h1, ax = plt.subplots(constrained_layout=True,figsize=(8, 4)) +# labels_h1 = ('EV_Charging', 'Power_EWH', 'Baseload', 'ESS_Charging') +# data_list_h1 = [PEVCharge[0], PEWH[0], baseload[0], PESSChargeGrid[0]] +# rc('font', family='Times New Roman') +# bottom_h1 = np.zeros(len(time)) +# for label_h1, data_h1 in zip(labels_h1, data_list_h1): +# ax.bar(time, data_h1, bottom=bottom_h1, label=label_h1) +# bottom_h1 += data_h1 +# # Set plot labels and title +# ax.set_ylabel('Power demand (kW)',fontname='Times New Roman',fontsize=12) +# # ax.set_title('Power demand (Prosumer-1)', fontname='Times New Roman') +# ax.legend(loc='upper left',fontsize=12,ncol=2) +# ax.set_xlabel('Time (15 min interval)',fontname='Times New Roman',fontsize=12) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# # Show plot +# plt.show() +# +# +# fig_h2, ax = plt.subplots(constrained_layout=True,figsize=(8, 4)) +# labels_h2 = ('EV_Charging', 'Power_EWH', 'Baseload', 'ESS_Charging') +# data_list_h2 = [PEVCharge[1], PEWH[1], baseload[1], PESSChargeGrid[1]] +# rc('font', family='Times New Roman') +# bottom_h2 = np.zeros(len(time)) +# for label_h2, data_h2 in zip(labels_h2, data_list_h2): +# ax.bar(time, data_h2, bottom=bottom_h2, label=label_h2) +# bottom_h2 += data_h2 +# # ax.plot(time,combine) +# # Set plot labels and title +# ax.set_ylabel('Power demand (kW)', fontname='Times New Roman',fontsize=12) +# # ax.set_title('Power demand (Prosumer-2)', fontname='Times New Roman',fontsize=12) +# ax.legend(loc='upper left',fontsize=12,ncol=2) +# ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman',fontsize=12) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# # Show plot +# plt.show() +# +# +# fig_h3, ax = plt.subplots(constrained_layout=True,figsize=(8, 4)) +# labels_h3 = ( 'Power_EWH', 'Baseload') +# data_list_h3 = [PEWH[2], baseload[2]] +# rc('font', family='Times New Roman') +# colors = ['orange', 'green'] +# bottom_h3 = np.zeros(len(time)) +# # for i,(label_h3, data_h3) in zip(labels_h3, data_list_h3): +# # ax.bar(time, data_h3, bottom=bottom_h3, label=label_h3, color=colors[label_h3]) +# # bottom_h3 += data_h3 +# for i, (label_h3, data_h3) in enumerate(zip(labels_h3, data_list_h3)): +# ax.bar(time, data_h3, bottom=bottom_h3, label=label_h3, color=colors[i]) +# bottom_h3 += data_h3 +# +# # ax.plot(time,combine) +# # Set plot labels and title +# ax.set_ylabel('Power demand (kW)', fontname='Times New Roman',fontsize=12) +# # ax.set_title('Power demand (Consumer-3)', fontname='Times New Roman',fontsize=12) +# ax.legend(loc='upper left',ncol=2,fontsize=12) +# ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman',fontsize=12) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# # Show plot +# plt.show() + + +# stack plot for the selling energy to the grid +# Plot stacked bar chart + +fig_all_sell, axs = plt.subplots(5, 2, sharex='col',constrained_layout=True, figsize=(8, 8), sharey='row') +labels = ('EV-Discharging_Grid', 'EV-Discharging_Home', 'PV_Grid', 'PV_Home', 'PV_Battery', 'Battery-Discharging_Grid', + 'Battery-Discharging_Home') +data_lists = [(PEVDischGrid[k], PEVDischHome[k], PVGrid[k], PVHome[k], PVBattery[k], PESSDischGrid[k], PESSDischHome[k]) + for k in range(10)] +rc('font', family='Times New Roman') +for i, ax in enumerate(axs.flat): + bottom = np.zeros(len(time)) + ax.set_xlim([0, 98]) + ax.set_ylim([0, 14]) + ax.set_xticks(np.arange(0, 98, 6)) + # ax.set_xticks( + # ['9am', '10am', '11am', '12pm', '1pm', '2pm', '3pm', '4pm', '5pm', '6pm', '7pm', '8pm', '9pm', '10pm', '11pm', + # '12am', '1am','2am', '3am', '4am', '5am', '6am', '7am', '8am']) + + for label, data in zip(labels, data_lists[i]): + ax.bar(time, data, bottom=bottom, label=label) + bottom += data + # if i == 3 or i == 6 or i ==9: + # ax.set_title(f'Consumer-{i+1}', fontname='Times New Roman', fontsize=12) + # else: + ax.set_title(f'Prosumer-{i + 1}', fontname='Times New Roman', fontsize=12) + ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) + if i == 0 or i == 2 or i == 4 or i == 6 or i == 8: + ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) + ax.set_yticks(np.arange(0, 14, 2)) + else: + ax.set_ylabel('') + if i == 8 or i == 9: + ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman', fontsize=12) + ax.tick_params(axis='x', which='both', labelsize=12) + else: + ax.tick_params(axis='x', which='both', bottom=True, labelbottom=False) # modified + ax.set_xticks([]) # added +handles, labels = ax.get_legend_handles_labels() +fig_all_sell.legend(handles, labels, loc='lower center', ncol=4, fontsize=10) +plt.tight_layout() +fig_all_sell.subplots_adjust(bottom=0.14) +plt.show() + +# fig2_h1, ax = plt.subplots(constrained_layout=True,figsize=(8, 4)) +# labels2_h1 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h1 = [PEVDischGrid[0],PEVDischHome[0], PVGrid[0],PVHome[0],PVBattery[0], PESSDischGrid[0],PESSDischHome[0]] +# bottom2_h1 = np.zeros(len(time)) +# for label2_h1, data2_h1 in zip(labels2_h1, data_list2_h1): +# ax.bar(time, data2_h1, bottom=bottom2_h1, label=label2_h1) +# bottom2_h1 += data2_h1 +# # Set plot labels and title +# ax.set_ylabel('Power Export (kW)', fontname='Times New Roman',fontsize=12) +# ax.set_xlabel('Time (15 min interval)',fontname='Times New Roman',fontsize=12) +# # ax.set_title('Power Export (Prosumer-1)', fontname='Times New Roman',fontsize=12) +# ax.legend(loc='upper left',ncol=3,fontsize=12) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 96, 4)) +# # Show plot +# plt.show() + +# fig_one, ax = plt.subplots(constrained_layout=True, figsize=(9, 6)) +# +# labels_h1 = ('EV_Charging', 'Power_EWH', 'Baseload', 'ESS_Charging') +# data_list_h1 = [PEVCharge[0], PEWH[0], baseload[0], PESSChargeGrid[0]] +# bottom_h1 = np.zeros(len(time)) +# colors_h1 = ['red', 'blue', 'green', 'purple'] +# for label_h1, data_h1, color_h1 in zip(labels_h1, data_list_h1, colors_h1): +# ax.bar(time, data_h1, bottom=bottom_h1, label=label_h1, color=color_h1) +# bottom_h1 += data_h1 +# +# labels2_h1 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h1 = [PEVDischGrid[0], PEVDischHome[0], PVGrid[0], PVHome[0], PVBattery[0], PESSDischGrid[0], PESSDischHome[0]] +# bottom2_h1 = np.zeros(len(time)) +# colors2_h1 = ['orange', 'pink', 'brown', 'gray', 'cyan', 'magenta', 'skyblue'] +# for label2_h1, data2_h1, color2_h1 in zip(labels2_h1, data_list2_h1, colors2_h1): +# ax.bar(time, np.negative(data2_h1), bottom=bottom2_h1[-len(data2_h1):], label=label2_h1, color=color2_h1) # invert the y-axis by negating the data and bottom values +# bottom2_h1[-len(data2_h1):] -= data2_h1 +# +# ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) +# ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman', fontsize=12) +# ax.legend(loc='upper left', fontsize=12, ncol=2) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# +# plt.show() +# +# my_palette = ['#FF5733', '#C70039', '#900C3F', '#581845', '#2E86C1', '#28B463', '#F7DC6F', '#F39C12', '#D35400', '#8E44AD', '#5D6D7E', '#34495E', '#2C3E50', '#ECF0F1'] +# +# palette = sns.color_palette(my_palette) +# fig_one, ax = plt.subplots(constrained_layout=True, figsize=(9, 6)) +# labels_h1 = ('EV_Charging', 'Power_EWH', 'Baseload', 'ESS_Charging') +# data_list_h1 = [PEVCharge[0], PEWH[0], baseload[0], PESSChargeGrid[0]] +# bottom_h1 = np.zeros(len(time)) +# for label_h1, data_h1, color_h1 in zip(labels_h1, data_list_h1, palette[:4]): +# ax.bar(time, data_h1, bottom=bottom_h1, label=label_h1, color=color_h1) +# bottom_h1 += data_h1 +# labels2_h1 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h1 = [PEVDischGrid[0], PEVDischHome[0], PVGrid[0], PVHome[0], PVBattery[0], PESSDischGrid[0], PESSDischHome[0]] +# bottom2_h1 = np.zeros(len(time)) +# for label2_h1, data2_h1, color2_h1 in zip(labels2_h1, data_list2_h1, palette[4:]): +# ax.bar(time, np.negative(data2_h1), bottom=bottom2_h1[-len(data2_h1):], label=label2_h1, color=color2_h1) # invert the y-axis by negating the data and bottom values +# bottom2_h1[-len(data2_h1):] -= data2_h1 +# ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) +# ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman', fontsize=12) +# ax.legend(loc='upper left', fontsize=12, ncol=2) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# def pos_int(x, pos): +# return abs(int(x)) +# +# # Format y-axis ticks and labels +# y_ticks = ticker.FuncFormatter(pos_int) +# ax.yaxis.set_major_formatter(y_ticks) +# plt.show() + + +# fig_two, ax = plt.subplots(constrained_layout=True, figsize=(9, 6)) +# labels_h1 = ('EV_Charging', 'Power_EWH', 'Baseload', 'ESS_Charging') +# data_list_h1 = [PEVCharge[1], PEWH[1], baseload[1], PESSChargeGrid[1]] +# bottom_h1 = np.zeros(len(time)) +# for label_h1, data_h1, color_h1 in zip(labels_h1, data_list_h1, palette[:4]): +# ax.bar(time, data_h1, bottom=bottom_h1, label=label_h1, color=color_h1) +# bottom_h1 += data_h1 +# labels2_h1 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h1 = [PEVDischGrid[1], PEVDischHome[1], PVGrid[1], PVHome[1], PVBattery[1], PESSDischGrid[1], PESSDischHome[1]] +# bottom2_h1 = np.zeros(len(time)) +# for label2_h1, data2_h1, color2_h1 in zip(labels2_h1, data_list2_h1, palette[4:]): +# ax.bar(time, np.negative(data2_h1), bottom=bottom2_h1[-len(data2_h1):], label=label2_h1, color=color2_h1) # invert the y-axis by negating the data and bottom values +# bottom2_h1[-len(data2_h1):] -= data2_h1 +# ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) +# ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman', fontsize=12) +# ax.legend(loc='upper left', fontsize=12, ncol=2) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# def pos_int(x, pos): +# return abs(int(x)) +# +# # Format y-axis ticks and labels +# y_ticks = ticker.FuncFormatter(pos_int) +# ax.yaxis.set_major_formatter(y_ticks) +# plt.show() + + +# fig_three, ax = plt.subplots(constrained_layout=True, figsize=(9, 6)) +# labels_h1 = ('EV_Charging', 'Power_EWH', 'Baseload', 'ESS_Charging') +# data_list_h1 = [PEVCharge[2], PEWH[2], baseload[2], PESSChargeGrid[2]] +# bottom_h1 = np.zeros(len(time)) +# for label_h1, data_h1, color_h1 in zip(labels_h1, data_list_h1, palette[:4]): +# ax.bar(time, data_h1, bottom=bottom_h1, label=label_h1, color=color_h1) +# bottom_h1 += data_h1 +# labels2_h1 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h1 = [PEVDischGrid[2], PEVDischHome[2], PVGrid[2], PVHome[2], PVBattery[2], PESSDischGrid[2], PESSDischHome[2]] +# bottom2_h1 = np.zeros(len(time)) +# for label2_h1, data2_h1, color2_h1 in zip(labels2_h1, data_list2_h1, palette[4:]): +# ax.bar(time, np.negative(data2_h1), bottom=bottom2_h1[-len(data2_h1):], label=label2_h1, color=color2_h1) # invert the y-axis by negating the data and bottom values +# bottom2_h1[-len(data2_h1):] -= data2_h1 +# ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) +# ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman', fontsize=12) +# ax.legend(loc='upper left', fontsize=12, ncol=2) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# def pos_int(x, pos): +# return abs(int(x)) +# +# # Format y-axis ticks and labels +# y_ticks = ticker.FuncFormatter(pos_int) +# ax.yaxis.set_major_formatter(y_ticks) +# plt.show() + +# fig_four, ax = plt.subplots(constrained_layout=True, figsize=(9, 6)) +# labels_h1 = ('EV_Charging', 'Power_EWH', 'Baseload', 'ESS_Charging') +# data_list_h1 = [PEVCharge[3], PEWH[3], baseload[3], PESSChargeGrid[3]] +# bottom_h1 = np.zeros(len(time)) +# for label_h1, data_h1, color_h1 in zip(labels_h1, data_list_h1, palette[:4]): +# ax.bar(time, data_h1, bottom=bottom_h1, label=label_h1, color=color_h1) +# bottom_h1 += data_h1 +# labels2_h1 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h1 = [PEVDischGrid[3], PEVDischHome[3], PVGrid[3], PVHome[3], PVBattery[3], PESSDischGrid[3], PESSDischHome[3]] +# bottom2_h1 = np.zeros(len(time)) +# for label2_h1, data2_h1, color2_h1 in zip(labels2_h1, data_list2_h1, palette[4:]): +# ax.bar(time, np.negative(data2_h1), bottom=bottom2_h1[-len(data2_h1):], label=label2_h1, color=color2_h1) # invert the y-axis by negating the data and bottom values +# bottom2_h1[-len(data2_h1):] -= data2_h1 +# ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) +# ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman', fontsize=12) +# ax.legend(loc='upper left', fontsize=12, ncol=2) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# def pos_int(x, pos): +# return abs(int(x)) +# +# # Format y-axis ticks and labels +# y_ticks = ticker.FuncFormatter(pos_int) +# ax.yaxis.set_major_formatter(y_ticks) +# plt.show() + + +# +# fig_one, ax = plt.subplots(constrained_layout=True, figsize=(9, 6)) +# labels_h1 = ('EV_Charging', 'Power_EWH', 'Baseload', 'ESS_Charging') +# data_list_h1 = [PEVCharge[0], PEWH[0], baseload[0], PESSChargeGrid[0]] +# bottom_h1 = np.zeros(len(time)) +# colors_h1 = ['#FF5733', '#2471A3', '#4CAF50', '#9B59B6'] +# for label_h1, data_h1, color_h1 in zip(labels_h1, data_list_h1, colors_h1): +# ax.bar(time, data_h1, bottom=bottom_h1, label=label_h1, color=color_h1) +# bottom_h1 += data_h1 +# labels2_h1 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h1 = [PEVDischGrid[0], PEVDischHome[0], PVGrid[0], PVHome[0], PVBattery[0], PESSDischGrid[0], PESSDischHome[0]] +# bottom2_h1 = np.zeros(len(time)) +# colors2_h1 = ['#FF8C42', '#3498DB', '#2ECC71', '#8E44AD', '#FFC300', '#E74C3C', '#85C1E9'] +# for label2_h1, data2_h1, color2_h1 in zip(labels2_h1, data_list2_h1, colors2_h1): +# ax.bar(time, np.negative(data2_h1), bottom=bottom2_h1[-len(data2_h1):], label=label2_h1, color=color2_h1) +# bottom2_h1[-len(data2_h1):] -= data2_h1 +# ax.set_ylabel('Power (kW)', fontname='Times New Roman', fontsize=12) +# ax.set_xlabel('Time (15 min interval)', fontname='Times New Roman', fontsize=12) +# ax.legend(loc='upper left', fontsize=12, ncol=2) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 98, 4)) +# plt.show() +# + +# +# +# +# +# fig2_h2, ax = plt.subplots(constrained_layout=True,figsize=(8, 4)) +# labels2_h2 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h2 = [PEVDischGrid[1],PEVDischHome[1], PVGrid[1],PVHome[1],PVBattery[1], PESSDischGrid[1],PESSDischHome[1]] +# bottom2_h2 = np.zeros(len(time)) +# for label2_h2, data2_h2 in zip(labels2_h2, data_list2_h2): +# ax.bar(time, data2_h2, bottom=bottom2_h2, label=label2_h2) +# bottom2_h2 += data2_h2 +# # Set plot labels and title +# ax.set_ylabel('Power Export (kW)', fontname='Times New Roman',fontsize=12) +# ax.set_xlabel('Time (15 min interval)',fontname='Times New Roman',fontsize=12) +# # ax.set_title('Power Export (Prosumer-2)', fontname='Times New Roman',fontsize=12) +# ax.legend(loc='upper left',ncol=3,fontsize=12) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 96, 4)) +# # Show plot +# plt.show() +# +# +# fig2_h3, ax = plt.subplots(constrained_layout=True,figsize=(8, 4)) +# labels2_h3 = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list2_h3 = [PEVDischGrid[2],PEVDischHome[2], PVGrid[2],PVHome[2],PVBattery[2], PESSDischGrid[2],PESSDischHome[2]] +# bottom2_h3 = np.zeros(len(time)) +# for label2_h3, data2_h3 in zip(labels2_h3, data_list2_h3): +# ax.bar(time, data2_h3, bottom=bottom2_h3, label=label2_h3) +# bottom2_h3 += data2_h3 +# # Set plot labels and title +# ax.set_ylabel('Power Export (kW)', fontname='Times New Roman') +# ax.set_xlabel('Time (15 min interval)',fontsize=13, fontname='Times New Roman') +# # ax.set_title('Power Export (Consumer-3)', fontname='Times New Roman') +# ax.legend(loc='upper left',ncol=3) +# ax.set_xlim([0, 97]) +# ax.set_xticks(np.arange(0, 96, 4)) +# # Show plot +# plt.show() + + +# +m2 = py.ConcreteModel() + +# +# m2.h1 = py.Set(initialize=sheet_home.home_bus) +# +# m2.h2 = py.Set(initialize=sheet_home.home) +# m2.h = py.Set(m2.h1, m2.h2) +# +# +# def homeset(m): +# return +# # ({h1: h2} for h1 in m2.h1 for h2 in m2.h2) +# # create sets with NO indexes +# m2.setHomes = py.Set(initialize=sheet_home.home.unique()) +# m2.setBuses = py.Set(initialize=sheet_home.home_bus.unique()) +# +# # create sets WITH indexes +# dict_home_givenBus = {j: sheet_home.home[sheet_home.home_bus == j].unique() for j in m2.setBuses} +# m2.setHome_givenBus = py.Set(m2.setBuses, initialize=dict_home_givenBus) +# +# dict_bus_givenHome = {h:int(sheet_home.home_bus[sheet_home.home==h]) for h in m2.setHomes} +# m2.setBuses_givenHome = py.Param(m2.setHomes, initialize=dict_bus_givenHome) +m2.b = py.Set(initialize=sheet_bus.bus, doc='bus') + +m2.h = py.Set(initialize=sheet_home.home) +m2.alpha = py.Param(m2.h, initialize=dict(zip(list(itertools.product(m2.h.data())), + [i for x in sheet_EV.columns if "arrival" in x for i + in + sheet_EV.loc[:, x].values])), doc="Arrival time of EV") + +m2.beta = py.Param(m2.h, initialize=dict(zip(list(itertools.product(m2.h.data())), + [i for x in sheet_EV.columns if "departure" in x for i + in + sheet_EV.loc[:, x].values])), doc="departure time of EV") +m2.t_final = py.Param(initialize=96, doc='final time of the day ') # 96 +m2.second = py.Param(initialize=2, doc='second time of the day ') +m2.t = py.Set(initialize=sheet_data.time, ordered=True, doc='time period') +# m2.EV_Dur = py.Set(initialize=py.RangeSet(m2.alpha.value, m2.beta.value), +# doc='time interval for EV') +# m2.EV_Dur1 = py.Set(initialize=py.RangeSet(m2.alpha.value + 1, m2.beta.value), +# doc='time interval for EV after arrival') +m2.t_second = py.Set(initialize=py.RangeSet(m2.second.value, m2.t_final.value), + doc='time greater then 1') +m2.t_first = py.Param(initialize=1, doc='first time of the day ') +# m2.ist_interval = py.Set(initialize=py.RangeSet(m2.t_first.value, m2.alpha.value - 1), +# doc='time period before the arrival time') +# m2.second_interval = py.Set(initialize=py.RangeSet(m2.beta.value + 1, m2.t_final.value), +# doc='time interval after departure') +m2.EV_Dur1 = py.Set(m2.h, initialize=lambda m, h: range(m.alpha[h] + 1, m.beta[h] + 1), doc='time interval for EV') +m2.EV_Dur = py.Set(m2.h, initialize=lambda m, h: range(m.alpha[h], m.beta[h] + 1), doc='time interval for EV') +m2.ist_interval = py.Set(m2.h, initialize=lambda m, h: range(m2.t_first.value, m.alpha[h]), doc='time interval for EV') +m2.second_interval = py.Set(m2.h, initialize=lambda m, h: range(m.beta[h] + 1, m2.t_final.value + 1), + doc='time interval after departure') + +# + +# +# Base load : uncontrollable loads +m2.base_load = py.Param(m2.h, m2.t, + initialize=dict(zip(list(itertools.product(m2.h.data(), m2.t.data())), + [i for x in sheet_load.columns if "load" in x for i in + sheet_load.loc[:, x].values])), doc="Base Load") + +m2.PV = py.Param(m2.h, m2.t, initialize=dict(zip(list(itertools.product(m2.h.data(), m2.t.data())), + [i for x in sheet_PV.columns if "PV" in x for i in + sheet_PV.loc[:, x].values])), doc="PV Production") +m2.max = py.Param(initialize=10000, doc='maximum value for selling and buying power') + +# + +# + +m2.ChargeRate_EV = py.Param(m2.h, initialize=dict(zip(list(itertools.product(m2.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") + +m2.DischRate_EV = py.Param(m2.h, initialize=dict(zip(list(itertools.product(m2.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") + +m2.Ch_Effc_EV = py.Param(m2.h, initialize=dict(zip(list(itertools.product(m2.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") + +m2.DischEffc_EV = py.Param(m2.h, initialize=dict(zip(list(itertools.product(m2.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") + +m2.Cap_EV = py.Param(m2.h, initialize=dict(zip((m2.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") +m2.End_percentage_EV = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.In_Percentage_EV = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.Energy_EV_dep = py.Param(m2.h, initialize=dict( + zip(m2.h.data(), np.array(m2.Cap_EV.values()) * m2.End_percentage_EV.values())), # just changing np>py + doc="Departure temperature of EV") + +m2.Energy_EV_In = py.Param(m2.h, initialize=dict( + zip(m2.h.data(), np.array(m2.Cap_EV.values()) * m2.In_Percentage_EV.values())), + doc="Initial energy of EV") + +# + +# + +m2.ChargeRate_ESS = py.Param(m2.h, initialize=dict(zip(list(itertools.product(m2.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") + +m2.DischRate_ESS = py.Param(m2.h, initialize=dict(zip(list(itertools.product(m2.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") + +# m2.ChargeRate_ESS = py.Param(initialize=float(sheet_ESS.charging_rate1), doc='Charging rate of ESS ') +# m2.DischRate_ESS = py.Param(initialize=float(sheet_ESS.discharging_rate1), +# doc='Discharging rate of ESS ') +m2.Cap_ESS = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.End_percentage_ESS = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.In_Percentage_ESS = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.Ch_Effc_ESS = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.DischEffc_ESS = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.Energy_ESS_In = py.Param(m2.h, initialize=dict( + zip(m2.h.data(), np.array(m2.Cap_ESS.values()) * m2.In_Percentage_ESS.values())), + doc="Initial energy of ESS") + +m2.End_En_ESS = py.Param(m2.h, initialize=dict( + zip(m2.h.data(), np.array(m2.Cap_ESS.values()) * m2.End_percentage_ESS.values())), + doc="Departure energy of ESS") + +# + +# + +m2.tetta_low = py.Param(m2.h, initialize=dict(zip((m2.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") +m2.tetta_up = py.Param(m2.h, initialize=dict(zip((m2.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") +m2.tetta_amb_int = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.tetta_amb = py.Param(m2.h, m2.t, + initialize=dict(zip(list(itertools.product(m2.h.data(), m2.t.data())), + [i for x in sheet_ambient.columns if "celsius" in x for i in + sheet_ambient.loc[:, x].values])), doc="Outdoor temperature") + +m2.Q = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.R = py.Param(m2.h, initialize=dict(zip((m2.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") +m2.C = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.M = py.Param(m2.h, initialize=dict(zip((m2.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)") + +m2.tetta_EWH_int = py.Param(m2.h, initialize=dict(zip((m2.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") + +m2.water_use = py.Param(m2.h, m2.t, + initialize=dict(zip(list(itertools.product(m2.h.data(), m2.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") + +# +# m2.Buy_price = py.Param(m2.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_RTP)), doc='Buying Price') +# m2.Buy_price = py.Param(m2.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_TOU_winter)), doc='Buying Price') +m2.Buy_price = py.Param(m2.t, initialize=dict(zip(sheet_data.time, sheet_data.TOU_PGE)), doc='Buying Price') +# m2.Buy_price = py.Param(m2.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_flex_D)), doc='Buying Price') +# m2.Buy_price = py.Param(m2.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price_fixed)), doc='Buying Price') +m2.Sell_price = py.Param(m2.t, initialize=dict(zip(sheet_data.time, sheet_data.FiT)), doc='Selling Price') +# Time duration: we took 15 mint granularity so in one hour it will be 1/4 +m2.time_d = py.Param(initialize=(1 / 4), doc='time duration ') +# m2.DPT = py.Param(initialize=0.069, doc='Daily power tariff ') # .069 +# + +print('Code Starts for HEMS Energy Optimization') +# Variable + +# +m2.P_Buy_Grid = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.P_Sell_Grid = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.S_P_sell = py.Var(m2.h, m2.t, within=py.Binary) +m2.S_P_buy = py.Var(m2.h, m2.t, within=py.Binary) +m2.peak = py.Var(within=py.NonNegativeReals) +# m2.pflex = py.Var(m2.b, m2.t, doc='active Flexibility') +# m2.qflex = py.Var(m2.b, m2.t, +# doc='reactive Flexibility') # i comment this because i use m2.pflex variable * 0.48 +m2.pmax = py.Var(m2.h, m2.t, bounds=(0, None), doc='Maximum power') +m2.pmax_int = py.Var(m2.h, m2.t, bounds=(0, None), doc='Initial Power') + +# + +# +m2.P_EV_Charge = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.P_EV_Disch = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.P_EV_Disch_Home = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.P_EV_Disch_Grid = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.S_EV_Charge = py.Var(m2.h, m2.t, within=py.Binary) +m2.S_EV_Disch = py.Var(m2.h, m2.t, within=py.Binary) +m2.Energy_EV = py.Var(m2.h, m2.t, bounds=(0, None)) # SOC of the EV +# + +# +m2.Energy_ESS = py.Var(m2.h, m2.t, bounds=(0, None)) # SOC of the ESS +m2.S_ESS_Charge = py.Var(m2.h, m2.t, within=py.Binary) +m2.P_ESS_Disch = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.P_ESS_Charge = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.P_ESS_Charge_Grid = py.Var(m2.h, m2.t, bounds=(0, None), doc='ESS charging from the Grid') +m2.P_ESS_Disch_Home = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.P_ESS_Disch_Grid = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.S_ESS_Disch = py.Var(m2.h, m2.t, within=py.Binary) +# + +# +m2.PV_Home = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.PV_Grid = py.Var(m2.h, m2.t, bounds=(0, None)) +m2.PV_Battery = py.Var(m2.h, m2.t, bounds=(0, None)) +# + +# +m2.tetta_EWH_wat = py.Var(m2.h, m2.t) +m2.S_EWH = py.Var(m2.h, m2.t, within=py.Binary) +m2.P_EWH = py.Var(m2.h, m2.t, doc='power of EWH') + + +# +# Constraints + + +# +def Power_buy(m2, h, i): + return m2.P_Buy_Grid[h, i] == m2.base_load[h, i] + m2.Q[h] * m2.S_EWH[h, i] + m2.P_EV_Charge[h, i] - \ + m2.P_EV_Disch_Home[h, i] + m2.P_ESS_Charge[h, i] - \ + m2.P_ESS_Disch_Home[h, i] - m2.PV_Home[h, i] - m2.PV_Battery[h, i] + + +m2.Const_1 = py.Constraint(m2.h, m2.t, rule=Power_buy, doc='Power buy from the Grid') + + +def Power_buy2(m2, h, i): + return m2.P_Buy_Grid[h, i] <= m2.max * m2.S_P_buy[h, i] + + +m2.Const_1a = py.Constraint(m2.h, m2.t, rule=Power_buy2, + doc='removing the nonlinearity in the objective fucntion') + + +def Power_sell1(m2, h, i): + return m2.P_Sell_Grid[h, i] == m2.P_EV_Disch_Grid[h, i] + m2.PV_Grid[h, i] + m2.P_ESS_Disch_Grid[h, i] + + +m2.Const_2 = py.Constraint(m2.h, m2.t, rule=Power_sell1, doc='Power sell to the Grid') + + +# def Power_sell2(m2, h, i): +# return m2.P_Sell_Grid[h, i] <= m2.max * m2.S_P_sell[h, i] +# +# +# m2.Const_2a = py.Constraint(m2.h, m2.t, rule=Power_sell2, doc='Power sell to the Grid') +def Power_sell2(m2, h, i): + return m2.P_Sell_Grid[h, i] <= m2.max.value * (1 - m2.S_P_buy[h, i]) + + +m2.Const_2a = py.Constraint(m2.h, m2.t, rule=Power_sell2, doc='Power sell to the Grid') + + +# def Status_Power(m2, h, i): +# return m2.S_P_buy[h, i] + m2.S_P_sell[h, i] <= 1 +# +# +# m2.Const_3 = py.Constraint(m2.h, m2.t, rule=Status_Power, +# doc='Buying and selling power will not occur at same time') + + +# + + +# +def Power_EV_Charge_limit1(m2, h, i): + return m2.P_EV_Charge[h, i] <= m2.ChargeRate_EV[h] * m2.S_EV_Charge[h, i] + + +m2.Const_EV_1 = py.Constraint(m2.h, m2.t, rule=Power_EV_Charge_limit1, + doc='Charging power of EV (Upper limit)') + + +def Power_EV_Charge_limit2(m2, h, i): + return m2.P_EV_Charge[h, i] >= 0 + + +m2.Const_EV_2 = py.Constraint(m2.h, m2.t, rule=Power_EV_Charge_limit2, + doc='Charging power of EV (lower limit)') + + +def Power_EV_Disch_limit1(m2, h, i): + return m2.P_EV_Disch[h, i] <= (m2.DischRate_EV[h] * m2.S_EV_Disch[h, i]) + + +m2.Const_EV_3 = py.Constraint(m2.h, m2.t, rule=Power_EV_Disch_limit1, + doc='Discharging power of EV (Upper limit)') + + +def Power_EV_Disch_limit2(m2, h, i): + return m2.P_EV_Disch[h, i] >= 0 + + +m2.Const_EV_4 = py.Constraint(m2.h, m2.t, rule=Power_EV_Disch_limit2, + doc='Discharging power of EV (Lower limit)') + + +def Status_EV(m2, h, i): + return m2.S_EV_Disch[h, i] + m2.S_EV_Charge[h, i] <= 1 + + +m2.Const_EV_5 = py.Constraint(m2.h, m2.t, rule=Status_EV, + doc='Charging and discharging will not occur at same time') + + +def Power_EV_Disch(m2, h, i): + return m2.P_EV_Disch[h, i] * m2.DischEffc_EV[h] == ( + (m2.P_EV_Disch_Home[h, i] + m2.P_EV_Disch_Grid[h, i])) + + +m2.Const_EV_6 = py.Constraint(m2.h, m2.t, rule=Power_EV_Disch, doc='Discharging power of EV to ' + 'Home and Grid') + + +def SoC_EV1(m2, h, i): + return m2.Energy_EV[h, i] == m2.Energy_EV_In[h] + + +# m2.Const_EV_7a = py.Constraint(m2.h, [m2.alpha.value], rule=SoC_EV1, doc='SoC of the EV at arrival') +m2.Const_EV_7a = py.Constraint([(h, py.value(m2.alpha[h])) for h in m2.h], rule=SoC_EV1, doc='SoC of the EV at arrival') + + +def SoC_EV2(m2, h, i): + return m2.Energy_EV[h, i] == m2.Energy_EV[h, i - 1] + ( + m2.P_EV_Charge[h, i] * m2.Ch_Effc_EV[h] - m2.P_EV_Disch[h, i]) * m2.time_d + + +# m2.Const_EV_7b = py.Constraint(m2.h, m2.EV_Dur1, rule=SoC_EV2, doc='SoC of the EV after arrival time') +m2.Const_EV_7b = py.ConstraintList() +for h in m2.h: + for i in m2.EV_Dur1[h].value: + m2.Const_EV_7b.add(SoC_EV2(m2, h, i)) + + +def EV_availability1(m2, h, i): + return m2.Energy_EV[h, i] == 0 + + +# m2.Const_EV_8a = py.Constraint(m2.h, m2.ist_interval, rule=EV_availability1, +# doc='SOC available before arrival time') +m2.Const_EV_8a = py.ConstraintList() +for h in m2.h: + for i in m2.ist_interval[h].value: + m2.Const_EV_8a.add(EV_availability1(m2, h, i)) + + +def EV_availability2(m2, h, i): + return m2.Energy_EV[h, i] == 0 + + +# m2.Const_EV_8b = py.Constraint(m2.h, m2.second_interval, rule=EV_availability2, +# doc='SOC available after departure time') +m2.Const_EV_8b = py.ConstraintList() +for h in m2.h: + for i in m2.second_interval[h].value: + m2.Const_EV_8b.add(EV_availability2(m2, h, i)) + + +def EV_status_available1(m2, h, i): + return m2.S_EV_Disch[h, i] + m2.S_EV_Charge[h, i] == 0 + + +# m2.Const_EV_9a = py.Constraint(m2.h, m2.ist_interval, rule=EV_status_available1, +# doc='EV availability before arrival time') +m2.Const_EV_9a = py.ConstraintList() +for h in m2.h: + for i in m2.ist_interval[h].value: + m2.Const_EV_9a.add(EV_status_available1(m2, h, i)) + + +def EV_status_available2(m2, h, i): + return m2.S_EV_Disch[h, i] + m2.S_EV_Charge[h, i] == 0 + + +# m2.Const_EV_9b = py.Constraint(m2.h, m2.second_interval, rule=EV_status_available2, +# doc='EV availability after arrival time') +m2.Const_EV_9b = py.ConstraintList() +for h in m2.h: + for i in m2.second_interval[h].value: + m2.Const_EV_9b.add(EV_status_available2(m2, h, i)) + + +def EV_SoC_limit1(m2, h, i): + return m2.Energy_EV[h, i] >= 0.2 * m2.Cap_EV[h] + + +# m2.Const_EV_10 = py.Constraint(m2.h, m2.EV_Dur, rule=EV_SoC_limit1, doc='Minimum SoC of EV') +m2.Const_EV_10 = py.ConstraintList() +for h in m2.h: + for i in m2.EV_Dur[h].value: + m2.Const_EV_10.add(EV_SoC_limit1(m2, h, i)) + + +def EV_SoC_limit2(m2, h, i): + return m2.Energy_EV[h, i] <= m2.Cap_EV[h] + + +# m2.Const_EV_11 = py.Constraint(m2.h, m2.t, rule=EV_SoC_limit2, doc='Maximum SoC of EV') +m2.Const_EV_11 = py.ConstraintList() +for h in m2.h: + for i in m2.EV_Dur[h].value: + m2.Const_EV_11.add(EV_SoC_limit2(m2, h, i)) + + +def EV_final_SoC(m2, h, i): + return m2.Energy_EV[h, i] == m2.Energy_EV_dep[h] + + +# m2.Const_EV_12 = py.Constraint(m2.h, [m2.beta.value], rule=EV_final_SoC, +# doc='Final SoC of EV at departure time') +m2.Const_EV_12 = py.Constraint([(h, py.value(m2.beta[h])) for h in m2.h], rule=EV_final_SoC, + doc='Final SoC of EV at departure time') + + +# + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +def SoC_ESS1(m2, h, i): + return m2.Energy_ESS[h, i] == m2.Energy_ESS_In[h] + + +m2.Const_ESS_1a = py.Constraint(m2.h, [m2.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(m2, h, i): + return m2.Energy_ESS[h, i] == m2.Energy_ESS[h, i - 1] + ( + m2.P_ESS_Charge[h, i] * m2.Ch_Effc_ESS[h] - m2.P_ESS_Disch[h, i]) * m2.time_d + + +m2.Const_ESS_1b = py.Constraint(m2.h, m2.t_second, rule=SoC_ESS2, doc='SoC of ESS') + + +def Power_ESS_Charge(m2, h, i): + return m2.P_ESS_Charge[h, i] == m2.P_ESS_Charge_Grid[h, i] + m2.PV_Battery[h, i] + + +m2.Const_ESS_2 = py.Constraint(m2.h, m2.t, rule=Power_ESS_Charge, + doc='Charging power of ESS ') + + +def Power_ESS_Charge_limit1(m2, h, i): + return (m2.P_ESS_Charge_Grid[h, i] + m2.PV_Battery[h, i]) <= m2.ChargeRate_ESS[h] * m2.S_ESS_Charge[ + h, i] + + +m2.Const_ESS_2a = py.Constraint(m2.h, m2.t, rule=Power_ESS_Charge_limit1, + doc='Charging power of ESS (Upper limit)') + + +def Power_ESS_Charge_limit2(m2, h, i): + return (m2.P_ESS_Charge_Grid[h, i] + m2.PV_Battery[h, i]) >= 0 + + +m2.Const_ESS_2b = py.Constraint(m2.h, m2.t, rule=Power_ESS_Charge_limit2, + doc='Charging power of ESS (Lower limit)') + + +def Power_ESS_Disch_limit1(m2, h, i): + return (m2.P_ESS_Disch_Home[h, i] + m2.P_ESS_Disch_Grid[h, i]) <= ( + m2.DischRate_ESS[h] * m2.S_ESS_Disch[h, i]) + + +m2.Const_ESS_4 = py.Constraint(m2.h, m2.t, rule=Power_ESS_Disch_limit1, + doc='Discharging power of ESS (Upper limit)') + + +def Power_ESS_Disch_limit2(m2, h, i): + return (m2.P_ESS_Disch_Home[h, i] + m2.P_ESS_Disch_Grid[h, i]) >= 0 + + +m2.Const_ESS_5 = py.Constraint(m2.h, m2.t, rule=Power_ESS_Disch_limit2, + doc='Discharging power of ESS (Lower limit)') + + +def Status_ESS(m2, h, i): + return m2.S_ESS_Disch[h, i] + m2.S_ESS_Charge[h, i] <= 1 + + +m2.Const_ESS_6 = py.Constraint(m2.h, m2.t, rule=Status_ESS, + doc='Charging and discharging will not occur at same time') + + +def Power_ESS_Disch(m2, h, i): + return (m2.P_ESS_Disch[h, i]) * m2.DischEffc_ESS[h] == (m2.P_ESS_Disch_Home[h, i] + m2.P_ESS_Disch_Grid[h, i]) + + +m2.Const_ESS_7 = py.Constraint(m2.h, m2.t, rule=Power_ESS_Disch, + doc='Discharging power of ESS to Home and Grid') + + +def ESS_SoC_limit1(m2, h, i): + return m2.Energy_ESS[h, i] >= 0.1 * m2.Cap_ESS[h] + + +m2.Const_ESS_8 = py.Constraint(m2.h, m2.t, rule=ESS_SoC_limit1, doc='Minimum SoC of ESS') + + +def ESS_SoC_limit2(m2, h, i): + return m2.Energy_ESS[h, i] <= m2.Cap_ESS[h] + + +m2.Const_ESS_9 = py.Constraint(m2.h, m2.t, rule=ESS_SoC_limit2, doc='Maximum SoC of ESS') + + +def ESS_final_SoC(m2, h, i): + return m2.Energy_ESS[h, i] >= m2.End_En_ESS[h] + + +m2.Const_ESS_10 = py.Constraint(m2.h, [m2.t_final.value], rule=ESS_final_SoC, + doc='Final SoC of ESS at departure time') + + +# + +# +def EWH_limit1(m2, h, i): + return m2.tetta_EWH_wat[h, i] <= m2.tetta_up[h] + + +m2.Const_EWH_1 = py.Constraint(m2.h, m2.t, rule=EWH_limit1, doc='Maximum Limit') + + +def EWH_limit2(m2, h, i): + return m2.tetta_EWH_wat[h, i] >= m2.tetta_low[h] + + +m2.Const_EWH_2 = py.Constraint(m2.h, m2.t, rule=EWH_limit2, doc='Minimum Limit') + + +def EWH_power(m2, h, i): + return m2.P_EWH[h, i] == m2.Q[h] * m2.S_EWH[h, i] + + +m2.Const_EWH_3 = py.Constraint(m2.h, m2.t, rule=EWH_power, doc='Electric water heater power') + + +def EWH_temp_1(m2, h, i): + return m2.tetta_EWH_wat[h, i] == m2.tetta_amb[h, i] + m2.Q[h] * m2.S_EWH[h, i] * m2.R[ + h] * m2.time_d - ( + (m2.M[h] - m2.water_use[h, i]) / m2.M[h]) * ( + m2.tetta_amb_int[h] - m2.tetta_EWH_int[h]) * py.exp(-m2.time_d / (m2.R[h] * m2.C[h])) + + +m2.Const_EWH_4 = py.Constraint(m2.h, [m2.t_first.value], rule=EWH_temp_1, doc='EWH m2') + + +def EWH_temp_2(m2, h, i): + return m2.tetta_EWH_wat[h, i] == m2.tetta_amb[h, i] + m2.Q[h] * m2.S_EWH[h, i] * m2.R[ + h] * m2.time_d - ( + (m2.M[h] - m2.water_use[h, i]) / m2.M[h]) * ( + m2.tetta_amb[h, i] - m2.tetta_EWH_wat[h, i - 1]) * py.exp( + -m2.time_d / (m2.R[h] * m2.C[h])) + + +m2.Const_EWH_5 = py.Constraint(m2.h, m2.t_second, rule=EWH_temp_2, doc='EWH m2') + + +# + +# # need to include the PV equation of omer paper + +def PV_production(m2, h, i): + return m2.PV_Grid[h, i] + m2.PV_Home[h, i] + m2.PV_Battery[h, i] == m2.PV[h, i] + + +m2.Const_PV = py.Constraint(m2.h, m2.t, rule=PV_production, doc='PV Production') + + +# +# + + +def objective_rule(m2): + return sum(m2.P_Buy_Grid[h, i] for i in m2.t for h in m2.h) * m2.time_d.value + + +m2.obj2 = py.Objective(rule=objective_rule, sense=py.minimize, doc='Definition of objective function') +# +m2.write('m22.lp', io_options={'symbolic_solver_labels': True}) +opt = py.SolverFactory('gurobi') +opt.options["mipgap"] = 0.8 # 0.155 for load and 0.8 for other load +result = opt.solve(m2, tee=True) # report_timing=True + +print(result) +# print('Sum of all homes objectives = ', py.value(m2.obj2)) +# for h in m2.h: +# print(f"Objective of Home {h} : " f"{(sum(m2.P_Buy_Grid[h, i].value for i in m2.t) * m2.time_d.value)}") +# print("second objective Done ") + + +# +# # + +# +# # # +# # +# # # Plotting for automatic ploting +time = [i for i in m2.t] +pl2_pbuy = [[] for j in m2.h] +# pl2_psell = [[] for j in m2.h] +PEVCharge2 = [[] for j in m2.h] +PESSChargeGrid2 = [[] for j in m2.h] +PEVDischHome2 = [[] for j in m2.h] +PESSDischHome2 = [[] for j in m2.h] +# # # SEVCharge = [[] for j in m2.h] +# # # SEVDisch = [[] for j in m2.h] +PEVDischGrid2 = [[] for j in m2.h] +PESSDischGrid2 = [[] for j in m2.h] +PVHome2 = [[] for j in m2.h] +PVGrid2 = [[] for j in m2.h] +PVBattery2 = [[] for j in m2.h] +# # # PPV = [[] for j in m2.h] +baseload2 = [[] for j in m2.h] +# # # EnergyESS = [[] for j in m2.h] +# # # EnergyEV = [[] for j in m2.h] +# # # tettaEWHwat = [[] for j in m2.h] +# # # SEWH = [[] for j in m2.h] +PEWH2 = [[] for j in m2.h] +# # # Buyprice = [] # cost = [i for i in m2.c.values()] <<< check this +# # # Sellprice = [] +# # # + +for j in m2.h: + for k, v in m2.P_Buy_Grid.items(): + if k[0] == j: + pl2_pbuy[j - 1].append(py.value(v)) +# for j in m2.h: +# for k, v in m2.P_Sell_Grid.items(): +# if k[0] == j: +# pl2_psell[j - 1].append(py.value(v)) +# # # +for j in m2.h: + for k, v in m2.P_EV_Charge.items(): + if k[0] == j: + PEVCharge2[j - 1].append(py.value(v)) +for j in m2.h: + for k, v in m2.P_EV_Disch_Home.items(): + if k[0] == j: + PEVDischHome2[j - 1].append(py.value(v)) +for j in m2.h: + for k, v in m2.P_EV_Disch_Grid.items(): + if k[0] == j: + PEVDischGrid2[j - 1].append(py.value(v)) +for j in m2.h: + for k, v in m2.P_ESS_Charge_Grid.items(): + if k[0] == j: + PESSChargeGrid2[j - 1].append(py.value(v)) +for j in m2.h: + for k, v in m2.P_ESS_Disch_Home.items(): + if k[0] == j: + PESSDischHome2[j - 1].append(py.value(v)) +for j in m2.h: + for k, v in m2.P_ESS_Disch_Grid.items(): + if k[0] == j: + PESSDischGrid2[j - 1].append(py.value(v)) +# # # +for j in m2.h: + for k, v in m2.PV_Home.items(): + if k[0] == j: + PVHome2[j - 1].append(py.value(v)) +for j in m2.h: + for k, v in m2.PV_Grid.items(): + if k[0] == j: + PVGrid2[j - 1].append(py.value(v)) +for j in m2.h: + for k, v in m2.PV_Battery.items(): + if k[0] == j: + PVBattery2[j - 1].append(py.value(v)) +# # # for j in m2.h: +# # # for k, v in m2.PV.items(): +# # # if k[0] == j: +# # # PPV[j - 1].append(py.value(v)) +# # # +for j in m2.h: + for k, v in m2.base_load.items(): + if k[0] == j: + baseload2[j - 1].append(py.value(v)) +# # # +# # # for j in m2.h: +# # # for k, v in m2.Energy_ESS.items(): +# # # if k[0] == j: +# # # EnergyESS[j - 1].append(py.value(v)) +# # # for j in m2.h: +# # # for k, v in m2.Energy_EV.items(): +# # # if k[0] == j: +# # # EnergyEV[j - 1].append(py.value(v)) +for j in m2.h: + for k, v in m2.P_EWH.items(): + if k[0] == j: + PEWH2[j - 1].append(py.value(v)) +# # # for j in m2.h: +# # # for k, v in m2.tetta_EWH_wat.items(): +# # # if k[0] == j: +# # # tettaEWHwat[j - 1].append(py.value(v)) +# # # +# # # for j in m2.h: +# # # for k, v in m2.S_EWH.items(): +# # # if k[0] == j: +# # # SEWH[j - 1].append(py.value(v)) +# # # + +# # # +# # # # alone +# # for k in range(len(m2.h)): +# # fig, ax = plt.subplots(5, 2, figsize=(10, 10)) +# # ax[0, 0].bar(time, pl2_pbuy[k], label='Buying power-obj2') +# # ax[0, 0].bar(time, pl1_pbuy[k], label='Buying power-obj1') +# # ax[1, 0].bar(time, pl2_psell[k], label='Selling power-obj2') +# # ax[1, 0].bar(time, pl1_psell[k], label='Selling power-obj1') +# # ax[0, 0].legend(loc='best', fontsize='small', ncol=3) +# # ax[1, 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, PESSChargeGrid[k], label='ESS Charging from Grid', color='r') +# # # # ax[2, 0].bar(time, PVBattery[k], label='ESS Charging from PV', color='g') +# # # # 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('Power1 (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(m2.P_Buy_Grid[k + 1, i] for i in m2.t) * m2.time_d))}") +# # pl.tight_layout() +# # plt.show() +# # # +# # # # code for if you want to plot all the home togathoer +# fig, ax = plt.subplots(4, len(m2.h), figsize=(10, 10), constrained_layout=True, sharex='col', sharey='row') +# gs = fig.add_gridspec(hspace=0, wspace=0) +# for k in range(len(m2.h)): +# ax[0, k].bar(time, pl1_pbuy[k], color='green') +# ax[1, k].bar(time, pl2_pbuy[k], color='blue') +# ax[2, k].bar(time, pl1_psell[k], color='green') +# ax[3, k].bar(time, pl2_psell[k], label='blue') +# # ax[0, k].legend(loc='best') +# # ax[1, k].legend(loc='best') +# # ax[2, k].legend(loc='best') +# # ax[3, k].legend(loc='best') +# # # ax[0, k].plot(time, baseload[k], label='Base load') +# # # # 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, 0].set_ylabel( +# 'Buy-Power1 (kW)') # showing only on the side of the graph if you want to put it on all need to use [0,k] +# ax[1, 0].set_ylabel('Buy-Power2 (kW)') +# ax[2, 0].set_ylabel('Sell-Power1 (kW)') +# ax[3, 0].set_ylabel('Sell-Power2 (kW)') +# +# ax[0, k].set_title(f"Home {k + 1}", fontsize='xx-small') +# fig.suptitle('Cost and Energy minimization', fontsize=16) +# +# plt.show() +# # +# # # if you want to plot each home alone in the figure +# # # for k in range(len(m2.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') +# # +# +# # +# + +# + + +# Plot stacked bar chart +# fig11, ax = plt.subplots() +# labels = ('EV_Charging', 'EWH', 'Baseload', 'ESS_Charging') +# data_list = [PEVCharge2[1], PEWH2[1], baseload2[1], PESSChargeGrid2[1]] +# bottom = np.zeros(len(time)) +# for label, data in zip(labels, data_list): +# ax.bar(time, data, bottom=bottom, label=label) +# bottom += data +# # Set plot labels and title +# ax.set_ylabel('Power (kW)') +# ax.set_xlabel('Time Step') +# ax.set_title('House Appliances (HEMS)_Energy Min') +# ax.legend(loc='upper left') +# +# # Show plot +# plt.show() +# +# +# fig111_en, ax = plt.subplots() +# labels_en = ('EV-Discharging_Grid','EV-Discharging_Home', 'PV_Grid', 'PV_Home','PV_Battery','ESS-Discharging_Grid','ESS-Discharging_Home') +# data_list_en = [PEVDischGrid[1],PEVDischHome[1], PVGrid[1],PVHome[1],PVBattery[1], PESSDischGrid[1],PESSDischHome[1]] +# bottom_en = np.zeros(len(time)) +# for label_en, data_en in zip(labels_en, data_list_en): +# ax.bar(time, data_en, bottom=bottom_en, label=label_en) +# bottom2+= data_en +# # Set plot labels and title +# ax.set_ylabel('Power (kW)') +# ax.set_xlabel('Time Step') +# ax.set_title('Power of House Appliances export to Home/Grid (HEMS)-Energ Min') +# ax.legend(loc='upper left') +# +# # Show plot +# plt.show() + +writer = pd.ExcelWriter('Result_HEMS.xlsx', engine='xlsxwriter') +time_excel = pd.DataFrame({'Time_Step': time}) +time_excel.to_excel(writer, sheet_name='Power', startcol=0, index=False, header=True) +for h in range(len(m1.h)): + pload_excel = pd.DataFrame({f"Power_Demand {h + 1}": pl1_pbuy[h]}) + pload_excel.to_excel(writer, sheet_name='Power', startcol=h + 1, index=False, header=True) + + pload_excel = pd.DataFrame({f"Power_Sell {h + 1}": pl1_psell[h]}) + pload_excel.to_excel(writer, sheet_name='Power', startcol=h + 4, index=False, header=True) +writer.save()