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