EMS/HEMS.py

# NOTE: for two homes

# <editor-fold desc="Import libraries and packages">
import itertools

import numpy as np
import pandas as pd
import pylab as pl
import pyomo.environ as py
from matplotlib import pyplot as plt

# data = pd.ExcelFile('input_new_cluster.xlsx') #it will read the excel one time you do not need to read again and gain
sheet_data = pd.read_excel("input_new_cluster.xlsx", sheet_name='data')
sheet_EV = pd.read_excel("input_new_cluster.xlsx", sheet_name='EV')
sheet_EWH = pd.read_excel("input_new_cluster.xlsx", sheet_name='EWH')
sheet_wateruse = pd.read_excel("input_new_cluster.xlsx", sheet_name='Water_Use')
sheet_watertemp = pd.read_excel("input_new_cluster.xlsx", sheet_name='Water_Temp')
sheet_home = pd.read_excel("input_new_cluster.xlsx", sheet_name='Home')
sheet_ESS = pd.read_excel("input_new_cluster.xlsx", sheet_name='ESS')
sheet_PV = pd.read_excel("input_new_cluster.xlsx", sheet_name='PV')
sheet_ambient = pd.read_excel("input_new_cluster.xlsx", sheet_name='Ambient_Temp')
sheet_load = pd.read_excel("input_new_cluster.xlsx", sheet_name='Base_Load')

# </editor-fold>


model = py.ConcreteModel()
# <editor-fold desc="Sets: time intervals">
model.h = py.Set(initialize=sheet_home.home)
model.alpha = py.Param(initialize=39, doc='time of arrival ')  # 6:30 = 39
model.beta = py.Param(initialize=89, doc='time of departure ')  # 7:00 = 89
model.t_final = py.Param(initialize=96, doc='final time of the day ')
model.second = py.Param(initialize=2, doc='second time of the day ')
model.t = py.Set(initialize=sheet_data.time, ordered=True, doc='time period')
model.EV_Dur = py.Set(initialize=py.RangeSet(model.alpha.value, model.beta.value),
                      doc='time interval for EV')
model.EV_Dur1 = py.Set(initialize=py.RangeSet(model.alpha.value + 1, model.beta.value),
                       doc='time interval for EV after arrival')
model.t_second = py.Set(initialize=py.RangeSet(model.second.value, model.t_final.value),
                        doc='time greater then 1')
model.t_first = py.Param(initialize=1, doc='first time of the day ')
model.ist_interval = py.Set(initialize=py.RangeSet(model.t_first.value, model.alpha.value - 1),
                            doc='time period before the arrival time')
model.second_interval = py.Set(initialize=py.RangeSet(model.beta.value + 1, model.t_final.value),
                               doc='time interval after departure')
# </editor-fold>

# NOTE:Parameter

# Base load : uncontrollable loads


model.base_load = py.Param(model.h, model.t,
                           initialize=dict(zip(list(itertools.product(model.h.data(), model.t.data())),
                                               [i for x in sheet_load.columns if "load" in x for i in
                                                sheet_load.loc[:, x].values])), doc="Base Load")

model.PV = py.Param(model.h, model.t, initialize=dict(zip(list(itertools.product(model.h.data(), model.t.data())),
                                                          [i for x in sheet_PV.columns if "PV" in x for i in
                                                           sheet_PV.loc[:, x].values])), doc="PV Production")
model.max = py.Param(initialize=100, doc='maximum value for selling and buying power')


# <editor-fold desc="Electric Vehicle Parameters">
model.ChargeRate_EV = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())),
                                                            [i for x in sheet_EV.columns if "charging_rate" in x for i
                                                             in
                                                             sheet_EV.loc[:, x].values])), doc="charging rate of EV")

model.DischRate_EV = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())),
                                                           [i for x in sheet_EV.columns if "rate_of_discharging" in x
                                                            for i in
                                                            sheet_EV.loc[:, x].values])), doc="Discharging rate of EV")

model.Ch_Effc_EV = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())),
                                                         [i for x in sheet_EV.columns if "charging_efficiency" in x for
                                                          i in
                                                          sheet_EV.loc[:, x].values])), doc="charging efficency of EV")

model.DischEffc_EV = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())),
                                                           [i for x in sheet_EV.columns if
                                                            "efficiency_of_discharging" in x for i
                                                            in sheet_EV.loc[:, x].values])),
                              doc="Discharging efficency of EV")

model.Cap_EV = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                     [i for x in sheet_EV.columns if "capacity" in x for i in
                                                      sheet_EV.loc[:, x].values])), doc="Capacity of  EV")
model.End_percentage_EV = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                                [i for x in sheet_EV.columns if "end" in x
                                                                 for i in
                                                                 sheet_EV.loc[:, x].values])),
                                   doc="Departure energy of EV")

model.In_Percentage_EV = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                               [i for x in sheet_EV.columns if "initial" in x
                                                                for i in
                                                                sheet_EV.loc[:, x].values])),
                                  doc="Initial energy of EV")

model.Energy_EV_dep = py.Param(model.h, initialize=dict(
    zip(model.h.data(), np.array(model.Cap_EV.values()) * model.End_percentage_EV.values())),  # just changeing np>py
                               doc="Departure temperature of EV")

model.Energy_EV_In = py.Param(model.h, initialize=dict(
    zip(model.h.data(), np.array(model.Cap_EV.values()) * model.In_Percentage_EV.values())),
                              doc="Initial energy of EV")

# </editor-fold>
# <editor-fold desc="Energy Storage System Parameter">

model.ChargeRate_ESS = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())),
                                                             [i for x in sheet_ESS.columns if "charging_rate" in x for i
                                                              in
                                                              sheet_ESS.loc[:, x].values])), doc="charging rate of ESS")

model.DischRate_ESS = py.Param(model.h, initialize=dict(zip(list(itertools.product(model.h.data())),
                                                            [i for x in sheet_ESS.columns if "rate_of_discharging" in x
                                                             for
                                                             i in
                                                             sheet_ESS.loc[:, x].values])),
                               doc="Discharging rate of ESS")
model.Cap_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                      [i for x in sheet_ESS.columns if "capacity" in x for i in
                                                       sheet_ESS.loc[:, x].values])), doc="Capacity of  ESS")

model.End_percentage_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                                 [i for x in sheet_ESS.columns if "end" in x
                                                                  for i in
                                                                  sheet_ESS.loc[:, x].values])),
                                    doc="Departure energy of ESS")

model.In_Percentage_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                                [i for x in sheet_ESS.columns if
                                                                 "initial" in x for i in
                                                                 sheet_ESS.loc[:, x].values])),
                                   doc="Initial energy of ESS")

model.Ch_Effc_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                          [i for x in sheet_ESS.columns if "charging_efficiency" in x
                                                           for i in
                                                           sheet_ESS.loc[:, x].values])),
                             doc="charging efficiency of ESS")

model.DischEffc_ESS = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                            [i for x in sheet_ESS.columns if
                                                             "efficiency_of_dicharging" in x for
                                                             i in
                                                             sheet_ESS.loc[:, x].values])),
                               doc="Discharging efficiency of ESS")

model.Energy_ESS_In = py.Param(model.h, initialize=dict(
    zip(model.h.data(), np.array(model.Cap_ESS.values()) * model.In_Percentage_ESS.values())),
                               doc="Initial energy of ESS")

model.End_En_ESS = py.Param(model.h, initialize=dict(
    zip(model.h.data(), np.array(model.Cap_ESS.values()) * model.End_percentage_ESS.values())),
                            doc="Departure energy of ESS")

# </editor-fold>
# <editor-fold desc="Electric Water Heater Parameter">

model.tetta_low = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                        [i for x in sheet_EWH.columns if "tetta_low" in x for i in
                                                         sheet_EWH.loc[:, x].values])),
                           doc="Lower bound of the water temperature")
model.tetta_up = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                       [i for x in sheet_EWH.columns if "tetta_up" in x for i in
                                                        sheet_EWH.loc[:, x].values])),
                          doc="Upper bound of the water temperature")
model.tetta_amb_int = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                            [i for x in sheet_EWH.columns if "tetta_amb_init" in x for i
                                                             in
                                                             sheet_EWH.loc[:, x].values])),
                               doc="Initial ambient temperature")


# todo need to use iteration for this below code

model.tetta_amb = py.Param(model.h, model.t,
                           initialize=dict(zip(list(itertools.product(model.h.data(), model.t.data())),
                                               [i for x in sheet_ambient.columns if "celsius" in x for i in
                                                sheet_ambient.loc[:, x].values])), doc="Outdoor temperature")

model.Q = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                [i for x in sheet_EWH.columns if "capacity" in x for i in
                                                 sheet_EWH.loc[:, x].values])), doc="Power of the EWH")

model.R = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                [i for x in sheet_EWH.columns if "thermal_resistance" in x for i in
                                                 sheet_EWH.loc[:, x].values])), doc="Thermal resistance")
model.C = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                [i for x in sheet_EWH.columns if "thermal_capacitance" in x for i in
                                                 sheet_EWH.loc[:, x].values])), doc="Thermal resistance")

model.M = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                [i for x in sheet_EWH.columns if "water_cap" in x for i in
                                                 sheet_EWH.loc[:, x].values])), doc="Water Capacity (L)")

model.tetta_EWH_int = py.Param(model.h, initialize=dict(zip((model.h.data()),
                                                            [i for x in sheet_EWH.columns if "tetta_wat_init" in x for i
                                                             in
                                                             sheet_EWH.loc[:, x].values])), doc="Initial temperature")

model.water_use = py.Param(model.h, model.t,
                           initialize=dict(zip(list(itertools.product(model.h.data(), model.t.data())),
                                               [i for x in sheet_wateruse.columns if "Litre" in x for i in
                                                sheet_wateruse.loc[:, x].values])), doc="Hot Water usage")

# </editor-fold2

# <editor-fold desc="Price, Base load, and time duration parameters">
model.Buy_price = py.Param(model.t, initialize=dict(zip(sheet_data.time, sheet_data.Buy_price)),
                           doc='Buying Price')
model.Sell_price = py.Param(model.t, initialize=dict(zip(sheet_data.time, sheet_data.Sell_price)), doc='Selling Price')
# Time duration: we took 15 mint granularity so in one hour it will be 1/4
model.time_d = py.Param(initialize=(1 / 4), doc='time duration ')
# </editor-fold>

# NOTE:Variable

# <editor-fold desc="buying and selling variables">
model.P_Buy_Grid = py.Var(model.h, model.t, bounds=(0, None))
model.P_Sell_Grid = py.Var(model.h, model.t, bounds=(0, None))
model.S_P_sell = py.Var(model.h, model.t, within=py.Binary)
model.S_P_buy = py.Var(model.h, model.t, within=py.Binary)
model.peak = py.Var(within=py.NonNegativeReals)

# </editor-fold>

# <editor-fold desc="Electric vehicle Variables">
model.P_EV_Charge = py.Var(model.h, model.t, bounds=(0, None))
model.P_EV_Disch = py.Var(model.h, model.t, bounds=(0, None))
model.P_EV_Disch_Home = py.Var(model.h, model.t, bounds=(0, None))
model.P_EV_Disch_Grid = py.Var(model.h, model.t, bounds=(0, None))
model.S_EV_Charge = py.Var(model.h, model.t, within=py.Binary)
model.S_EV_Disch = py.Var(model.h, model.t, within=py.Binary)
model.Energy_EV = py.Var(model.h, model.t, bounds=(0, None))  # SOC of the EV
# </editor-fold>

# <editor-fold desc="Energy storage system Variables">
model.Energy_ESS = py.Var(model.h, model.t, bounds=(0, None))  # SOC of the ESS
model.S_ESS_Charge = py.Var(model.h, model.t, within=py.Binary)
model.P_ESS_Disch = py.Var(model.h, model.t, bounds=(0, None))
model.P_ESS_Charge = py.Var(model.h, model.t, bounds=(0, None))
model.P_ESS_Charge_Grid = py.Var(model.h, model.t, bounds=(0, None), doc='ESS charging from the Grid')
model.P_ESS_Disch_Home = py.Var(model.h, model.t, bounds=(0, None))
model.P_ESS_Disch_Grid = py.Var(model.h, model.t, bounds=(0, None))
model.S_ESS_Disch = py.Var(model.h, model.t, within=py.Binary)
# </editor-fold>

# <editor-fold desc="PV variables">
model.PV_Home = py.Var(model.h, model.t, bounds=(0, None))
model.PV_Grid = py.Var(model.h, model.t, bounds=(0, None))
model.PV_Battery = py.Var(model.h, model.t, bounds=(0, None))
# </editor-fold>

# <editor-fold desc="EWH Variables">
model.tetta_EWH_wat = py.Var(model.h, model.t)
model.S_EWH = py.Var(model.h, model.t, within=py.Binary)
model.P_EWH = py.Var(model.h, model.t, doc='power of EWH')


# </editor-fold>

# NOTE:Constraints

# <editor-fold desc="Buying and Selling Power">
def Power_buy(model, h, i):
    return model.P_Buy_Grid[h, i] == model.base_load[h, i] + model.Q[h] * model.S_EWH[h, i] + model.P_EV_Charge[h, i] - \
           model.P_EV_Disch_Home[h, i] + model.P_ESS_Charge[h, i] - \
           model.P_ESS_Disch_Home[h, i] - model.PV_Home[h, i]- model.PV_Battery[h, i]


model.Const_1 = py.Constraint(model.h, model.t, rule=Power_buy, doc='Power buy from the Grid')


def Power_buy2(model, h, i):
    return model.P_Buy_Grid[h, i] <= model.max * model.S_P_buy[h, i]


model.Const_1a = py.Constraint(model.h, model.t, rule=Power_buy2,
                               doc='removing the nonlinearity in the objective fucntion')


def Power_sell1(model, h, i):
    return model.P_Sell_Grid[h, i] == model.P_EV_Disch_Grid[h, i] + model.PV_Grid[h, i] + model.P_ESS_Disch_Grid[h, i]


model.Const_2 = py.Constraint(model.h, model.t, rule=Power_sell1, doc='Power sell to the Grid')


def Power_sell2(model, h, i):
    return model.P_Sell_Grid[h, i] <= model.max * model.S_P_sell[h, i]


model.Const_2a = py.Constraint(model.h, model.t, rule=Power_sell2, doc='Power sell to the Grid')


def Status_Power(model, h, i):
    return model.S_P_buy[h, i] + model.S_P_sell[h, i] <= 1


model.Const_3 = py.Constraint(model.h, model.t, rule=Status_Power,
                              doc='Buying and selling power will not occur at same time')


# </editor-fold>

# </editor-fold>

# <editor-fold desc="Electric vehicle Constraints">
def Power_EV_Charge_limit1(model, h, i):
    return model.P_EV_Charge[h, i] <= model.ChargeRate_EV[h] * model.S_EV_Charge[h, i]


model.Const_EV_1 = py.Constraint(model.h, model.t, rule=Power_EV_Charge_limit1,
                                 doc='Charging power of EV (Upper limit)')


def Power_EV_Charge_limit2(model, h, i):
    return model.P_EV_Charge[h, i] >= 0.2 * (model.ChargeRate_EV[h] * model.S_EV_Charge[h, i])


model.Const_EV_2 = py.Constraint(model.h, model.t, rule=Power_EV_Charge_limit2,
                                 doc='Charging power of EV (lower limit)')


def Power_EV_Disch_limit1(model, h, i):
    return model.P_EV_Disch[h, i] <= (model.DischRate_EV[h] * model.S_EV_Disch[h, i])


model.Const_EV_3 = py.Constraint(model.h, model.t, rule=Power_EV_Disch_limit1,
                                 doc='Discharging power of EV (Upper limit)')


def Power_EV_Disch_limit2(model, h, i):
    return model.P_EV_Disch[h, i] >= 0.2 * (model.DischRate_EV[h] * model.S_EV_Disch[h, i])


model.Const_EV_4 = py.Constraint(model.h, model.t, rule=Power_EV_Disch_limit2,
                                 doc='Discharging power of EV (Lower limit)')


def Status_EV(model, h, i):
    return model.S_EV_Disch[h, i] + model.S_EV_Charge[h, i] <= 1


model.Const_EV_5 = py.Constraint(model.h, model.t, rule=Status_EV,
                                 doc='Charging and discharging will not occur at same time')


def Power_EV_Disch(model, h, i):
    return model.P_EV_Disch[h, i] == (
            (model.P_EV_Disch_Home[h, i] + model.P_EV_Disch_Grid[h, i]) / model.DischEffc_EV[h])


model.Const_EV_6 = py.Constraint(model.h, model.t, rule=Power_EV_Disch, doc='Discharging power of EV to '
                                                                            'Home and Grid')


def SoC_EV1(model, h, i):
    return model.Energy_EV[h, i] == model.Energy_EV_In[h]
model.Const_EV_7a = py.Constraint(model.h, [model.alpha.value], rule=SoC_EV1, doc='SoC of the EV at arrival')


def SoC_EV2(model, h, i):
    return model.Energy_EV[h, i] == model.Energy_EV[h, i - 1] + (
            model.P_EV_Charge[h, i] * model.Ch_Effc_EV[h] - model.P_EV_Disch[h, i]) * model.time_d


model.Const_EV_7b = py.Constraint(model.h, model.EV_Dur1, rule=SoC_EV2, doc='SoC of the EV after arrival time')

def EV_availability1(model, h, i):
    return model.Energy_EV[h, i] == 0


model.Const_EV_8a = py.Constraint(model.h, model.ist_interval, rule=EV_availability1,
                                  doc='SOC available before arrival time')


def EV_availability2(model, h, i):
    return model.Energy_EV[h, i] == 0


model.Const_EV_8b = py.Constraint(model.h, model.second_interval, rule=EV_availability2,
                                  doc='SOC available after departure time')


def EV_status_available1(model, h, i):
    return model.S_EV_Disch[h, i] + model.S_EV_Charge[h, i] == 0


model.Const_EV_9a = py.Constraint(model.h, model.ist_interval, rule=EV_status_available1,
                                  doc='EV availability before arrival time')


def EV_status_available2(model, h, i):
    return model.S_EV_Disch[h, i] + model.S_EV_Charge[h, i] == 0


model.Const_EV_9b = py.Constraint(model.h, model.second_interval, rule=EV_status_available2,
                                  doc='EV availability after arrival time')


def EV_SoC_limit1(model, h, i):
    return model.Energy_EV[h, i] >= 0.1 * model.Cap_EV[h]


model.Const_EV_10 = py.Constraint(model.h, model.EV_Dur, rule=EV_SoC_limit1, doc='Minimum SoC of EV')


def EV_SoC_limit2(model, h, i):
    return model.Energy_EV[h, i] <= model.Cap_EV[h]


model.Const_EV_11 = py.Constraint(model.h, model.t, rule=EV_SoC_limit2, doc='Maximum SoC of EV')


def EV_final_SoC(model, h, i):
    return model.Energy_EV[h, i] == model.Energy_EV_dep[h]


model.Const_EV_12 = py.Constraint(model.h, [model.beta.value], rule=EV_final_SoC,
                                  doc='Final SoC of EV at departure time')


# </editor-fold>

# <editor-fold desc="Energy Storage System Constraints">
def SoC_ESS1(model, h, i):
    return model.Energy_ESS[h, i] == model.Energy_ESS_In[h]
model.Const_ESS_1a = py.Constraint(model.h, [model.t_first.value], rule=SoC_ESS1,
                                   doc='SoC of ESS')  # need to check wather we need to put in the brackets or not


def SoC_ESS2(model, h, i):
    return model.Energy_ESS[h, i] == model.Energy_ESS[h, i - 1] + (
            model.P_ESS_Charge[h, i] * model.Ch_Effc_ESS[h] - model.P_ESS_Disch[h, i]) * model.time_d
model.Const_ESS_1b = py.Constraint(model.h, model.t_second, rule=SoC_ESS2, doc='SoC of ESS')


def Power_ESS_Charge(model, h, i):
    return model.P_ESS_Charge[h, i] == model.P_ESS_Charge_Grid[h, i] + model.PV_Battery[h, i]
model.Const_ESS_2 = py.Constraint(model.h, model.t, rule=Power_ESS_Charge,
                                   doc='Charging power of ESS ')


def Power_ESS_Charge_limit1(model, h, i):
    return model.P_ESS_Charge[h, i] <= model.ChargeRate_ESS[h] * model.S_ESS_Charge[h, i]


model.Const_ESS_2a = py.Constraint(model.h, model.t, rule=Power_ESS_Charge_limit1,
                                   doc='Charging power of ESS (Upper limit)')


def Power_ESS_Charge_limit2(model, h, i):
    return model.P_ESS_Charge[h, i] >= 0.2 * (model.ChargeRate_ESS[h] * model.S_ESS_Charge[h, i])


model.Const_ESS_2b = py.Constraint(model.h, model.t, rule=Power_ESS_Charge_limit2,
                                   doc='Charging power of ESS (Lower limit)')


def Power_ESS_Disch_limit1(model, h, i):
    return model.P_ESS_Disch[h, i] <= (model.DischRate_ESS[h] * model.S_ESS_Disch[h, i])


model.Const_ESS_4 = py.Constraint(model.h, model.t, rule=Power_ESS_Disch_limit1,
                                  doc='Discharging power of ESS (Upper limit)')


def Power_ESS_Disch_limit2(model, h, i):
    return model.P_ESS_Disch[h, i] >= 0.2 * (model.DischRate_ESS[h] * model.S_ESS_Disch[h, i])


model.Const_ESS_5 = py.Constraint(model.h, model.t, rule=Power_ESS_Disch_limit2,
                                  doc='Discharging power of ESS (Lower limit)')


def Status_ESS(model, h, i):
    return model.S_ESS_Disch[h, i] + model.S_ESS_Charge[h, i] <= 1


model.Const_ESS_6 = py.Constraint(model.h, model.t, rule=Status_ESS,
                                  doc='Charging and discharging will not occur at same time')


def Power_ESS_Disch(model, h, i):
    return (model.P_ESS_Disch[h, i] * model.DischEffc_ESS[h]) == (
            model.P_ESS_Disch_Home[h, i] + model.P_ESS_Disch_Grid[h, i])


model.Const_ESS_7 = py.Constraint(model.h, model.t, rule=Power_ESS_Disch,
                                  doc='Discharging power of ESS to Home and Grid')


def ESS_SoC_limit1(model, h, i):
    return model.Energy_ESS[h, i] >= 0.2 * model.Cap_ESS[h]


model.Const_ESS_8 = py.Constraint(model.h, model.t, rule=ESS_SoC_limit1, doc='Minimum SoC of ESS')


def ESS_SoC_limit2(model, h, i):
    return model.Energy_ESS[h, i] <= model.Cap_ESS[h]


model.Const_ESS_9 = py.Constraint(model.h, model.t, rule=ESS_SoC_limit2, doc='Maximum SoC of ESS')


def ESS_final_SoC(model, h, i):
    return model.Energy_ESS[h, i] >= model.End_En_ESS[h]


model.Const_ESS_10 = py.Constraint(model.h, [model.t_final.value], rule=ESS_final_SoC,
                                   doc='Final SoC of ESS at departure time')


# </editor-fold>

# <editor-fold desc="Electric Water Heater (EWH) Constraints">
def EWH_limit1(model, h, i):
    return model.tetta_EWH_wat[h, i] <= model.tetta_up[h]


model.Const_EWH_1 = py.Constraint(model.h, model.t, rule=EWH_limit1, doc='Maximum Limit')


def EWH_limit2(model, h, i):
    return model.tetta_EWH_wat[h, i] >= model.tetta_low[h]


model.Const_EWH_2 = py.Constraint(model.h, model.t, rule=EWH_limit2, doc='Minimum Limit')


def EWH_power(model, h, i):
    return model.P_EWH[h, i] == model.Q[h] * model.S_EWH[h, i]


model.Const_EWH_3 = py.Constraint(model.h, model.t, rule=EWH_power, doc='Electric water heater power')


def EWH_temp_1(model, h, i):
    return model.tetta_EWH_wat[h, i] == model.tetta_amb[h, i] + model.Q[h] * model.S_EWH[h, i] * model.R[
        h] * model.time_d - (
                   (model.M[h] - model.water_use[h, i]) / model.M[h]) * (
                   model.tetta_amb_int[h] - model.tetta_EWH_int[h]) * py.exp(-model.time_d / (model.R[h] * model.C[h]))


model.Const_EWH_4 = py.Constraint(model.h, [model.t_first.value], rule=EWH_temp_1, doc='EWH Model')


def EWH_temp_2(model, h, i):
    return model.tetta_EWH_wat[h, i] == model.tetta_amb[h, i] + model.Q[h] * model.S_EWH[h, i] * model.R[
        h] * model.time_d - (
                   (model.M[h] - model.water_use[h, i]) / model.M[h]) * (
                   model.tetta_amb[h, i] - model.tetta_EWH_wat[h, i - 1]) * py.exp(
        -model.time_d / (model.R[h] * model.C[h]))


model.Const_EWH_5 = py.Constraint(model.h, model.t_second, rule=EWH_temp_2, doc='EWH Model')


# </editor-fold>

# <editor-fold desc="PV Constraints"> # need to include the PV equation of omer paper
def PV_production(model, h, i):
    return model.PV_Grid[h, i] + model.PV_Home[h, i] + model.PV_Battery[h, i] == model.PV[h, i]

model.Const_PV = py.Constraint(model.h, model.t, rule=PV_production, doc='PV Production')


# </editor-fold>

def objective_rule(model):
    return sum((model.P_Buy_Grid[h, i] * model.Buy_price[i]) - (
            model.P_Sell_Grid[h, i] * model.Sell_price[i]) for i in
               model.t for h in model.h) * model.time_d


model.objective = py.Objective(rule=objective_rule, sense=py.minimize, doc='Definition of objective function')

# <editor-fold desc="Solver">
model.write('model2.lp', io_options={'symbolic_solver_labels': True})
opt = py.SolverFactory('cplex')
opt.options["mipgap"] = 0.09  # 0.155 for load and 0.8 for other load
result = opt.solve(model, tee=True)
print(result)
# </editor-fold>

# <editor-fold desc="Check optimility">
if (result.solver.status == py.SolverStatus.ok) and (
        result.solver.termination_condition == py.TerminationCondition.optimal):
    print('optimal solution')  # Do something when the solution in optimal and feasible
elif result.solver.termination_condition == py.TerminationCondition.infeasible:
    print('solution is not feasible')  # Do something when model in infeasible
else:
    # Something else is wrong
    print('Solver Status: ', result.solver.status)
print('Sum of all homes objectives = ', py.value(model.objective))
# </editor-fold>

for h in model.h:
    print(
        f"Objective of Home {h} : " f"{py.value((sum((model.P_Buy_Grid[h, i] * model.Buy_price[i]) - (model.P_Sell_Grid[h, i] * model.Sell_price[i]) for i in model.t) * model.time_d))}")
# <editor-fold desc="Ploting">


# Plotting for automatic ploting
time = [i for i in model.t]
pbuy = [[] for j in model.h]
psell = [[] for j in model.h]
PEVCharge = [[] for j in model.h]
PESSCharge = [[] for j in model.h]
PEVDischHome = [[] for j in model.h]
PESSDischHome = [[] for j in model.h]
SEVCharge = [[] for j in model.h]
SEVDisch = [[] for j in model.h]
PEVDischGrid = [[] for j in model.h]
PESSDischGrid = [[] for j in model.h]
PVHome = [[] for j in model.h]
PVGrid = [[] for j in model.h]
PVBattery = [[] for j in model.h]
PPV = [[] for j in model.h]
baseload = [[] for j in model.h]
EnergyESS = [[] for j in model.h]
EnergyEV = [[] for j in model.h]
tettaEWHwat = [[] for j in model.h]
SEWH = [[] for j in model.h]
PEWH = [[] for j in model.h]
Buyprice = []  # cost = [i for i in model.c.values()]  <<< check this
Sellprice = []

for j in model.h:
    for k, v in model.P_Buy_Grid.items():
        if k[0] == j:
            pbuy[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.P_Sell_Grid.items():
        if k[0] == j:
            psell[j - 1].append(py.value(v))

for j in model.h:
    for k, v in model.P_EV_Charge.items():
        if k[0] == j:
            PEVCharge[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.P_EV_Disch_Home.items():
        if k[0] == j:
            PEVDischHome[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.P_EV_Disch_Grid.items():
        if k[0] == j:
            PEVDischGrid[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.P_ESS_Charge.items():
        if k[0] == j:
            PESSCharge[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.P_ESS_Disch_Home.items():
        if k[0] == j:
            PESSDischHome[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.P_ESS_Disch_Grid.items():
        if k[0] == j:
            PESSDischGrid[j - 1].append(py.value(v))

for j in model.h:
    for k, v in model.PV_Home.items():
        if k[0] == j:
            PVHome[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.PV_Grid.items():
        if k[0] == j:
            PVGrid[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.PV_Battery.items():
        if k[0] == j:
            PVBattery[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.PV.items():
        if k[0] == j:
            PPV[j - 1].append(py.value(v))

for j in model.h:
    for k, v in model.base_load.items():
        if k[0] == j:
            baseload[j - 1].append(py.value(v))

for j in model.h:
    for k, v in model.Energy_ESS.items():
        if k[0] == j:
            EnergyESS[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.Energy_EV.items():
        if k[0] == j:
            EnergyEV[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.P_EWH.items():
        if k[0] == j:
            PEWH[j - 1].append(py.value(v))
for j in model.h:
    for k, v in model.tetta_EWH_wat.items():
        if k[0] == j:
            tettaEWHwat[j - 1].append(py.value(v))

for j in model.h:
    for k, v in model.S_EWH.items():
        if k[0] == j:
            SEWH[j - 1].append(py.value(v))

for i in model.Buy_price:
    Buyprice.append(py.value(model.Buy_price[i]))

for i in model.Sell_price:
    Sellprice.append(py.value(model.Sell_price[i]))

# alone
# for k in range(len(model.h)):
#     fig, ax = plt.subplots(5, 3, figsize=(10, 10))
#     ax[0, 0].bar(time, pbuy[k], label='Buying power')
#     ax[0, 0].bar(time, psell[k], label='Selling power')
#     ax[0, 0].legend(loc='best', fontsize='small', ncol=3)
#     ax[0, 1].bar(time, baseload[k], label='Base load', color='r')
#     ax[0, 1].legend(loc='best', fontsize='small', ncol=3)
#     ax[0, 2].plot(time, Buyprice, label='Buyprice')
#     ax[0, 2].plot(time, Sellprice, label='Sellprice')
#     ax[0, 2].legend(loc='best', fontsize='small', ncol=3)
#     ax[1, 0].bar(time, PEWH[k], label='EWH Power', color='r')
#     ax[1, 0].legend(loc='best', fontsize='small', ncol=3)
#     ax[1, 1].plot(time, tettaEWHwat[k], label='EWH Temp')
#     ax[1, 1].legend(loc='best', fontsize='small', ncol=3)
#     ax[1, 2].bar(time, SEWH[k], label='Status of EWH')
#     ax[1, 2].legend(loc='best', fontsize='small', ncol=3)
#     ax[2, 0].bar(time, PESSCharge[k], label='ESS Charging power', color='r')
#     ax[2, 0].legend(loc='best', fontsize='small', ncol=3)
#     ax[2, 1].bar(time, PESSDischHome[k], label='ESS Disch to home')
#     ax[2, 1].bar(time, PESSDischGrid[k], label='ESS Disch to grid')
#     ax[2, 1].legend(loc='best', fontsize='small', ncol=3)
#     ax[2, 2].plot(time, EnergyESS[k], label='Energy of ESS', color='g')
#     ax[2, 2].legend(loc='best', fontsize='small', ncol=3)
#     ax[3, 0].bar(time, PVHome[k], label='PV to Home')
#     ax[3, 0].legend(loc='best', fontsize='small', ncol=3)
#     ax[3, 1].bar(time, PVGrid[k], label='PV to Grid')
#     ax[3, 1].bar(time, PVBattery[k], label='PV to battery')
#     ax[3, 1].legend(loc='best', fontsize='small', ncol=3)
#     ax[3, 2].bar(time, PVHome[k], label='PV to Home')
#     ax[3, 2].bar(time, PVGrid[k], label='PV to Grid')
#     ax[3, 2].legend(loc='best', fontsize='small', ncol=3)
#     ax[4, 0].bar(time, PEVCharge[k], label='EV Charging power', color='r')
#     ax[4, 0].legend(loc='best', fontsize='small', ncol=3)
#     ax[4, 1].bar(time, PEVDischHome[k], label='EV Disch to home')
#     ax[4, 1].bar(time, PEVDischGrid[k], label='EV Disch to grid')
#     ax[4, 1].legend(loc='best', fontsize='small', ncol=3)
#     ax[4, 2].plot(time, EnergyEV[k], label='Energy of EV', color='r')
#     ax[4, 2].legend(loc='best', fontsize='small', ncol=3)
#     ax[4, 0].set_xlabel('Time (step)')
#     ax[4, 1].set_xlabel('Time (step)')
#     ax[4, 2].set_xlabel('Time (step)')
#     ax[0, 2].set_ylabel('price ($/W/h)')
#     ax[0, 0].set_ylabel('Power (kW)')
#     ax[1, 0].set_ylabel('Power (kW)')
#     ax[2, 0].set_ylabel('Power (kW)')
#     ax[3, 0].set_ylabel('Power (kW)')
#     ax[4, 0].set_ylabel('Power (kW)')
#     ax[2, 2].set_ylabel('Energy (kWh)')
#     ax[4, 2].set_ylabel('Energy (kWh)')
#     plt.suptitle(
#         f" Home {k + 1} : " f" Cost is :  {py.value((sum((model.P_Buy_Grid[k + 1, i] * model.Buy_price[i]) - (model.P_Sell_Grid[k + 1, i] * model.Sell_price[i]) for i in model.t) * model.time_d))}")
#     pl.tight_layout()
#     plt.show()

# code for if you want to plot all the home togathoer
# fig, ax = plt.subplots(6, len(model.h), figsize=(10,10))
# for k in range(len(model.h)):
#     ax[0, k].plot(time, pbuy[k], label='Buying power')
#     ax[0, k].plot(time, psell[k], label='Selling power')
#     # ax[0, k].plot(time, baseload[k], label='Base load')
#     ax[0, 0].legend(loc='best')
#     ax[1, k].plot(time, Buyprice,  label='Buyprice')
#     ax[1, 0].legend(loc='best')
#     ax[1, k].plot(time, Sellprice, label='Sellprice')
#     ax[1, 0].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, 0].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, 0].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, 0].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, 0].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, 0].legend(loc='best')
#     ax[4, 5].set_xlabel('Time (step)')
#     ax[1, 0].set_ylabel('Electricity price ($/W/h)')
#     ax[0, 0].set_ylabel('Power (kW)')
#     ax[2,0].set_ylabel('Power (kW)')
#     ax[0, k].set_title(f"Home {k+1}")
#

# plt.show()


# if you want to plot each home alone in the figure
bar_width =0.4
time1 = np.array(time)
for k in range(len(model.h)):
    fig, ax = plt.subplots(2,constrained_layout=True, )
    ax[0].bar(time1 - bar_width/2 , pbuy[k], width=bar_width, label='Base')
    ax[0].bar(time1 + bar_width/2 , psell[k], width=bar_width, 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}")
    plt.show()
# # </editor-fold>