model = AbstractModel()

model.T = Set()  # Index Set for time steps of optimization horizon
model.T_SoC = Set()  # SoC of the ESSs at the end of optimization horizon are also taken into account

################################## PARAMETERS #################################

model.dT = Param(within=PositiveIntegers)  # Number of seconds in one time step

model.P_PV = Param(model.T, within=NonNegativeReals)  # PV PMPP forecast

model.ESS_Min_SoC = Param(within=PositiveReals)  # Minimum SoC of ESSs
model.ESS_Max_SoC = Param(within=PositiveReals)  # Maximum SoC of ESSs
model.SoC_Value = Param(within=PositiveReals)
model.ESS_Capacity = Param(within=PositiveReals)  # Storage Capacity of ESSs
model.ESS_Max_Charge_Power = Param(within=PositiveReals)  # Max Charge Power of ESSs
model.ESS_Max_Discharge_Power = Param(within=PositiveReals)  # Max Discharge Power of ESSs

model.P_Grid_Max_Export_Power = Param(within=NonNegativeReals)  # Max active power export
model.Q_Grid_Max_Export_Power = Param(within=NonNegativeReals)  # Max reactive power export

model.PV_Inv_Max_Power = Param(within=PositiveReals)  # PV inverter capacity

################################## VARIABLES #################################

model.P_Grid_Output = Var(model.T, within=Reals)
model.Q_Grid_Output = Var(model.T, within=Reals)
model.P_PV_Output = Var(model.T, within=NonNegativeReals, bounds=(0, model.PV_Inv_Max_Power))
model.P_ESS_Output = Var(model.T, within=Reals, bounds=(-model.ESS_Max_Charge_Power, model.ESS_Max_Discharge_Power))
model.SoC_ESS = Var(model.T_SoC, within=NonNegativeReals, bounds=(model.ESS_Min_SoC, model.ESS_Max_SoC))

################################ CONSTRAINTS ##################################

# PV constraints
def con_rule_pv_potential(model, t):
    return model.P_PV_Output[t] <= model.P_PV[t]

# Import/Export constraints for Active power
def con_rule_grid_P_inv(model, t):
    return model.P_Grid_Output[t] >= -model.P_Grid_Max_Export_Power

# Import/Export constraints for Active power
def con_rule_grid_Q_inv(model, t):
    return model.Q_Grid_Output[t] >= -model.Q_Grid_Max_Export_Power

# ESS SoC balance
def con_rule_socBalance(model, t):
    return model.SoC_ESS[t + 1] == model.SoC_ESS[t] - model.P_ESS_Output[t] * model.dT / model.ESS_Capacity / 3600

# Initializing SoC
def con_rule_iniSoC(model):
    return model.SoC_ESS[0] == model.SoC_Value

# Generation-feed in balance
def con_rule_generation_feedin(model, t):
    return model.P_Grid_Output[t] * model.P_Grid_Output[t] + model.Q_Grid_Output[t] * model.Q_Grid_Output[t] == (
            model.P_PV_Output[t] + model.P_ESS_Output[t]) * (model.P_PV_Output[t] + model.P_ESS_Output[t])


# Gnerating constraints
model.con_pv_pmax = Constraint(model.T, rule=con_rule_pv_potential)
model.con_grid_inv_P = Constraint(model.T, rule=con_rule_grid_P_inv)
model.con_grid_inv_Q = Constraint(model.T, rule=con_rule_grid_Q_inv)
model.con_ess_soc = Constraint(model.T, rule=con_rule_socBalance)
model.con_ess_Inisoc = Constraint(rule=con_rule_iniSoC)
model.con_gen_feedin = Constraint(model.T, rule=con_rule_generation_feedin)


############################## OBJECTIVE ##########################################

def obj_rule(model):
    return sum(model.P_PV[t] - model.P_PV_Output[t] for t in model.T)

model.obj = Objective(rule=obj_rule, sense=minimize)