Feature/gsye 901

src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/cop_models.py

@@ -4,6 +4,7 @@
 from enum import Enum
 from typing import Optional
 from logging import getLogger
+import sympy as sp
 
 from gsy_framework.enums import HeatPumpSourceType
 
@@ -17,6 +18,7 @@ class COPModelType(Enum):
     ELCO_AEROTOP_S09_IR = 1
     ELCO_AEROTOP_G07_14M = 2
     HOVAL_ULTRASOURCE_B_COMFORT_C11 = 3
+    AERMEC_NXP_0600_4L_HEAT = 4
 
 
 MODEL_FILE_DIR = os.path.join(os.path.dirname(__file__), "model_data")
@@ -26,6 +28,7 @@ class COPModelType(Enum):
     COPModelType.ELCO_AEROTOP_G07_14M: "Elco_Aerotop_G07-14M_model_parameters.json",
     COPModelType.HOVAL_ULTRASOURCE_B_COMFORT_C11: "hoval_UltraSource_B_comfort_C_11_model_"
     "parameters.json",
+    COPModelType.AERMEC_NXP_0600_4L_HEAT: "AERMEC_NXP_0600_4L_HEAT_model_parameters.json",
 }
 
 
@@ -34,10 +37,23 @@ class BaseCOPModel:
 
     @abstractmethod
     def calc_cop(
-        self, source_temp_C: float, condenser_temp_C: float, heat_demand_kW: Optional[float]
+        self,
+        source_temp_C: float,
+        condenser_temp_C: float,
+        heat_demand_kW: Optional[float],
+        electrical_demand_kW: Optional[float] = None,
     ):
         """Return COP value for provided inputs"""
 
+    @abstractmethod
+    def calc_q_from_p_kW(
+        self,
+        source_temp_C: float,
+        condenser_temp_C: float,
+        electrical_demand_kW: Optional[float] = None,
+    ):
+        """Calculate heat energy from provided inputs."""
+
 
 class IndividualCOPModel(BaseCOPModel):
     """Handles cop models for specific heat pump models"""
@@ -54,8 +70,8 @@ def __init__(self, model_type: COPModelType):
     def _calc_power(self, source_temp_C: float, condenser_temp_C: float, heat_demand_kW: float):
         CAPFT = (
             self._model["CAPFT"][0]
-            + self._model["CAPFT"][1] * condenser_temp_C
-            + self._model["CAPFT"][3] * condenser_temp_C**2
+            + self._model["CAPFT"][1] * source_temp_C
+            + self._model["CAPFT"][3] * source_temp_C**2
             + self._model["CAPFT"][2] * condenser_temp_C
             + self._model["CAPFT"][5] * condenser_temp_C**2
             + self._model["CAPFT"][4] * source_temp_C * condenser_temp_C
@@ -83,8 +99,109 @@ def _calc_power(self, source_temp_C: float, condenser_temp_C: float, heat_demand
         # Power consumption (P) calculation
         return self._model["Pref"] * CAPFT * HEIRFT * HEIRFPLR
 
-    def calc_cop(self, source_temp_C: float, condenser_temp_C: float, heat_demand_kW: float):
+    def _resolve_heat(
+        self, source_temp_C: float, condenser_temp_C: float, electricity_demand_kW: float
+    ):
+        CAPFT = (
+            self._model["CAPFT"][0]
+            + self._model["CAPFT"][1] * source_temp_C
+            + self._model["CAPFT"][3] * source_temp_C**2
+            + self._model["CAPFT"][2] * condenser_temp_C
+            + self._model["CAPFT"][5] * condenser_temp_C**2
+            + self._model["CAPFT"][4] * source_temp_C * condenser_temp_C
+        )
+
+        HEIRFT = (
+            self._model["HEIRFT"][0]
+            + self._model["HEIRFT"][1] * source_temp_C
+            + self._model["HEIRFT"][3] * source_temp_C**2
+            + self._model["HEIRFT"][2] * condenser_temp_C
+            + self._model["HEIRFT"][5] * condenser_temp_C**2
+            + self._model["HEIRFT"][4] * source_temp_C * condenser_temp_C
+        )
+
+        # Partial Load Ratio (PLR)
+        Q = sp.symbols("Q")
+        PLR = Q / (self._model["Qref"] * CAPFT)
+
+        # HEIRFPLR calculation
+        HEIRFPLR = (
+            self._model["HEIRFPLR"][0]
+            + self._model["HEIRFPLR"][1] * PLR
+            + self._model["HEIRFPLR"][2] * PLR**2
+        )
+
+        solutions = sp.solve(
+            sp.Eq(electricity_demand_kW, self._model["Pref"] * CAPFT * HEIRFT * HEIRFPLR), Q
+        )
+        Q = self._select_Q_solution(solutions, CAPFT)
+        if Q is None:
+            # fallback: use median COP of training dataset to calculate Q
+            Q = self._model["COP_med"] * electricity_demand_kW
+        return Q
+
+    def _select_Q_solution(self, Q_solutions, CAPFT) -> Optional[float]:
+        """
+        Selects the correct Q and PLR solution based on:
+        - PLR = Q / (Qref * fCAPFT)
+        - both PLRs must be between 0 and 1
+        - the correct branch is the one with the LARGER PLR
+          (as indicated by the training dataset)
+        """
+
+        PLR_dict = {
+            q / (self._model["Qref"] * CAPFT): q
+            for q in Q_solutions
+            if 0 <= q / (self._model["Qref"] * CAPFT) <= 1
+        }
+        if not PLR_dict:
+            PLR_list = [q / (self._model["Qref"] * CAPFT) for q in Q_solutions]
+            log.error("IndividualCOPModel: No physically feasible PLR solutions. %s", PLR_list)
+            return None
+        return float(PLR_dict[max(PLR_dict)])
+
+    def _limit_heat_demand_kW(self, heat_demand_kW: float) -> float:
         assert heat_demand_kW is not None, "heat demand should be provided"
+        if heat_demand_kW > self._model["Q_max"]:
+            log.debug(
+                "calc_cop: heat demand exceeds maximum heat_demand_kW: %s", self._model["Q_max"]
+            )
+            return self._model["Q_max"]
+        if heat_demand_kW < self._model["Q_min"]:
+            log.debug(
+                "calc_cop: heat demand exceeds minimum heat_demand_kW: %s", self._model["Q_min"]
+            )
+            return self._model["Q_min"]
+        return heat_demand_kW
+
+    def calc_q_from_p_kW(
+        self,
+        source_temp_C: float,
+        condenser_temp_C: float,
+        electrical_demand_kW: Optional[float] = None,
+    ):
+        return self._resolve_heat(
+            source_temp_C=source_temp_C,
+            condenser_temp_C=condenser_temp_C,
+            electricity_demand_kW=electrical_demand_kW,
+        )
+
+    def calc_cop(
+        self,
+        source_temp_C: float,
+        condenser_temp_C: float,
+        heat_demand_kW: float,
+        electrical_demand_kW: Optional[float] = None,
+    ):
+
+        if (electrical_demand_kW is None) == (heat_demand_kW is None):
+            assert False, "heat_demand_kW and electrical_demand_kW can only be set exclusively"
+
+        if electrical_demand_kW:
+            # estimate the heat demand by using the median COP of the model fitting data
+            heat_demand_kW = electrical_demand_kW * self._model["COP_med"]
+
+        heat_demand_kW = self._limit_heat_demand_kW(heat_demand_kW)
         if heat_demand_kW == 0:
             return 0
         electrical_power_kW = self._calc_power(source_temp_C, condenser_temp_C, heat_demand_kW)
@@ -101,7 +218,7 @@ def calc_cop(self, source_temp_C: float, condenser_temp_C: float, heat_demand_kW
             )
             return 0
         cop = heat_demand_kW / electrical_power_kW
-        if cop > 6:
+        if cop > self._model["COP_max"] or cop < self._model["COP_min"]:
             log.error(
                 "calculated COP (%s) is unrealistic: "
                 "hp model: %s  source_temp: %s, "
@@ -122,17 +239,36 @@ class UniversalCOPModel(BaseCOPModel):
     def __init__(self, source_type: int = HeatPumpSourceType.AIR.value):
         self._source_type = source_type
 
-    def calc_cop(
-        self, source_temp_C: float, condenser_temp_C: float, heat_demand_kW: Optional[float]
-    ) -> float:
-        """COP model following https://www.nature.com/articles/s41597-019-0199-y"""
+    def _calc_cop_from_temps(
+        self,
+        source_temp_C: float,
+        condenser_temp_C: float,
+    ):
         delta_temp = condenser_temp_C - source_temp_C
         if self._source_type == HeatPumpSourceType.AIR.value:
             return 6.08 - 0.09 * delta_temp + 0.0005 * delta_temp**2
         if self._source_type == HeatPumpSourceType.GROUND.value:
             return 10.29 - 0.21 * delta_temp + 0.0012 * delta_temp**2
         assert False, "Source type not supported"
 
+    def calc_cop(
+        self,
+        source_temp_C: float,
+        condenser_temp_C: float,
+        heat_demand_kW: Optional[float],
+        electrical_demand_kW: Optional[float] = None,
+    ) -> float:
+        """COP model following https://www.nature.com/articles/s41597-019-0199-y"""
+        return self._calc_cop_from_temps(source_temp_C, condenser_temp_C)
+
+    def calc_q_from_p_kW(
+        self,
+        source_temp_C: float,
+        condenser_temp_C: float,
+        electrical_demand_kW: Optional[float] = None,
+    ):
+        return electrical_demand_kW * self._calc_cop_from_temps(source_temp_C, condenser_temp_C)
+
 
 def cop_model_factory(
     model_type: COPModelType, source_type: int = HeatPumpSourceType.AIR.value

src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/AERMEC_NXP_0600_4L_COOL_model_parameters.json

@@ -0,0 +1 @@
+{"CAPFT": [0.8454605267956957, 0.027689073804528678, 0.006404959469341046, 0.0001496419458169494, -4.569496057765451e-05, -0.00021118596972354869], "HEIRFT": [0.7418220836760441, -0.0022664853115584786, -0.007181589206940325, 0.000674105489882314, -0.000946546960846392, 0.0005940100532001941], "HEIRFPLR": [-0.5215602319631321, 1.6972347812119102, -0.1618280142560721], "Qref": 141.5, "Pref": 30.9, "Q_min": 118.9, "Q_max": 180.5, "PLR_min": 0.25, "COP_min": 3.22, "COP_max": 6.51, "COP_med": 4.8}

src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/AERMEC_NXP_0600_4L_HEAT_model_parameters.json

@@ -0,0 +1 @@
+{"CAPFT": [0.7841944915827704, 0.027025683620693207, 0.005222590512291558, 0.00021692430116427368, -0.00019218411177690167, -0.00011531764270024707], "HEIRFT": [0.7162755985430916, -0.004161089855566163, -0.0005729077062670702, 0.00041800740384663796, -0.0006154588026100481, 0.00039346652665741823], "HEIRFPLR": [-1.3954158938415628, 2.7629166672566634, -0.33808079196917734], "Qref": 175.7, "Pref": 32.2, "Q_min": 142.7, "Q_max": 226.3, "PLR_min": 0.25, "COP_min": 3.21, "COP_max": 7.04, "COP_med": 5.15}

src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/Elco_Aerotop_G07-14M_model_parameters.json

@@ -1 +1 @@
-{"CAPFT": [1.3674619178648815, 0.03804659465962401, -0.019681855839990714, 0.00035564628176876223, -0.0003621972557350528, 0.00020130672570195518], "HEIRFT": [0.0013367484492077253, -0.005549193693871857, 0.020206838944944794, -6.42720175278999e-05, -0.00023069214578574915, -7.601014281322094e-06], "HEIRFPLR": [-0.194398846570212, 1.468913188968992, -0.2833766224760481], "Qref": 16.18, "Pref": 7.6, "PLR_min": 0.2}
+{"CAPFT": [1.367461917867454, 0.0380465946617361, -0.019681855837919038, 0.0003556462839485741, -0.0003621972535536866, 0.00020130672788409854], "HEIRFT": [0.001336748450156744, -0.005549193691562149, 0.02020683894711295, -6.427201534764393e-05, -0.0002306921436066034, -7.601012100844073e-06], "HEIRFPLR": [-0.19439885055755093, 1.4689131623547622, -0.2833766040414506], "Qref": 16.18, "Pref": 7.6, "Q_min": 11.14, "Q_max": 25.99, "COP_min": 1.6940509915014166, "COP_max": 5.208416833667334, "COP_med": 2.6843220338983054}

src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/Elco_Aerotop_S09M-IR_model_parameters.json

@@ -1 +1 @@
-{"CAPFT": [0.9300492302752958, 0.03526591169765436, -0.004968669510962087, 0.00041747226581356767, -0.0003022915023160877, 9.773660370027137e-06], "HEIRFT": [0.3556666343811108, -0.02211409627372074, 0.0089215929159423, 5.995243647605175e-05, -3.4550553532408657e-05, 0.00014357037230472436], "HEIRFPLR": [-0.049758215126849414, 1.5407231145753069, -0.4967040777633469], "Qref": 16.2, "Pref": 6.47, "PLR_min": 0.1}
+{"CAPFT": [0.930049230273987, 0.03526591169977844, -0.00496866950870567, 0.0004174722679944898, -0.00030229150013494355, 9.77366255017209e-06], "HEIRFT": [0.3556666343824486, -0.022114096271478978, 0.008921592918107346, 5.995243865708488e-05, -3.455055135148655e-05, 0.00014357037448553545], "HEIRFPLR": [-0.04975819746928794, 1.540722949813866, -0.4967039408074606], "Qref": 16.2, "Pref": 6.47, "Q_min": 7.76, "Q_max": 23.15, "COP_min": 1.533596837944664, "COP_max": 5.701970443349754, "COP_med": 2.563652692071984}

src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/hoval_UltraSource_B_comfort_C_11_model_parameters.json

@@ -1 +1 @@
-{"CAPFT": [0.5152305055168107, 0.016858287502234393, 0.04002915190886247, -0.00029916924915052157, 6.765446186174362e-05, -0.000496776482422856], "HEIRFT": [0.8811301356647324, -0.0019346380513800977, -0.024648245650783343, -0.0003083507655603219, -0.0001586422096891705, 0.00039208468135561403], "HEIRFPLR": [-0.3720222858742652, 2.375345746851973, -1.0047114539738535], "Qref": 8.3, "Pref": 5.7, "PLR_min": 0.1}
+{"CAPFT": [0.515230505521292, 0.016858287504343594, 0.040029151910794813, -0.0002991692469687113, 6.765446404344289e-05, -0.0004967764802388253], "HEIRFT": [0.8811300289423628, -0.001934636544902224, -0.02464824123871967, -0.0003083507728460777, -0.00015864223856865145, 0.00039208463749784705], "HEIRFPLR": [-0.372022346462102, 2.3753458974617483, -1.0047115475588355], "Qref": 8.3, "Pref": 5.7, "Q_min": 8.3, "Q_max": 12.8, "COP_min": 1.456140350877193, "COP_max": 5.565217391304349, "COP_med": 2.7027027027027026}

src/gsy_e/models/strategy/energy_parameters/heatpump/heat_pump.py

@@ -10,6 +10,9 @@
     convert_kJ_to_kWh,
     convert_kWh_to_kJ,
     convert_kJ_to_kW,
+    convert_kWh_to_W,
+    convert_kWh_to_kW,
+    convert_kW_to_kWh,
 )
 from pendulum import DateTime
 
@@ -187,7 +190,7 @@ def get_energy_to_buy_maximum_kWh(self, time_slot: DateTime, source_temp_C: floa
             time_slot, self.heatpump.get_heat_demand_kJ(time_slot)
         )
         cop = self._calc_cop(
-            produced_heat_kJ=max_heat_demand_kJ, source_temp_C=source_temp_C, time_slot=time_slot
+            heat_energy_kJ=max_heat_demand_kJ, source_temp_C=source_temp_C, time_slot=time_slot
         )
         if cop == 0:
             return 0
@@ -207,7 +210,7 @@ def get_energy_to_buy_minimum_kWh(self, time_slot: DateTime, source_temp_C: floa
             time_slot, self.heatpump.get_heat_demand_kJ(time_slot)
         )
         cop = self._calc_cop(
-            produced_heat_kJ=min_heat_demand_kJ, source_temp_C=source_temp_C, time_slot=time_slot
+            heat_energy_kJ=min_heat_demand_kJ, source_temp_C=source_temp_C, time_slot=time_slot
         )
         if cop == 0:
             return 0
@@ -220,7 +223,7 @@ def get_energy_demand_kWh(self, time_slot: DateTime, source_temp_C: float) -> fl
         """Return the energy demand in kWh that is needed to be traded by the heat pump."""
         heat_demand_kJ = self.heatpump.get_heat_demand_kJ(time_slot)
         cop = self._calc_cop(
-            produced_heat_kJ=heat_demand_kJ, source_temp_C=source_temp_C, time_slot=time_slot
+            heat_energy_kJ=heat_demand_kJ, source_temp_C=source_temp_C, time_slot=time_slot
         )
         if cop == 0:
             return 0
@@ -270,17 +273,22 @@ def update_cop_after_dis_charging(
             cop = self._hp_state.get_cop(last_time_slot)
         else:
             cop = self._calc_cop(
-                produced_heat_kJ=convert_kWh_to_kJ(bought_energy_kWh),
+                heat_energy_kJ=None,
                 source_temp_C=source_temp_C,
                 time_slot=last_time_slot,
+                electrical_energy_kWh=bought_energy_kWh,
             )
 
         # Set the calculated COP on both the last and the current time slot to use in calculations
         self._hp_state.set_cop(last_time_slot, cop)
         self._hp_state.set_cop(time_slot, cop)
 
     def _calc_cop(
-        self, produced_heat_kJ: float, source_temp_C: float, time_slot: DateTime
+        self,
+        heat_energy_kJ: float,
+        source_temp_C: float,
+        time_slot: DateTime,
+        electrical_energy_kWh: Optional[float] = None,
     ) -> float:
         """
         Return the coefficient of performance (COP) for a given ambient and storage temperature.
@@ -289,26 +297,34 @@ def _calc_cop(
         Generally, the higher the temperature difference between the source and the sink,
         the lower the efficiency of the heat pump (the lower COP).
         """
-        if produced_heat_kJ < FLOATING_POINT_TOLERANCE:
-            return 0
-        heat_demand_kW = convert_kJ_to_kW(produced_heat_kJ, GlobalConfig.slot_length)
+        if electrical_energy_kWh is None:
+            if heat_energy_kJ < FLOATING_POINT_TOLERANCE:
+                return 0
+            heat_demand_kW = convert_kJ_to_kW(heat_energy_kJ, GlobalConfig.slot_length)
+            electrical_energy_kW = None
+        else:
+            heat_demand_kW = None
+            electrical_energy_kW = (
+                convert_kWh_to_W(electrical_energy_kWh, GlobalConfig.slot_length) / 1000
+            )
+
         return self._cop_model.calc_cop(
             source_temp_C=source_temp_C,
             condenser_temp_C=self._charger.get_average_inlet_temperature_C(time_slot),
             heat_demand_kW=heat_demand_kW,
+            electrical_demand_kW=electrical_energy_kW,
         )
 
     def calc_Q_kJ_from_energy_kWh(
         self, time_slot: DateTime, energy_kWh: float, source_temp_C: float
     ) -> float:
         """Calculate heat in kJ from energy in kWh."""
-        energy_kJ = convert_kWh_to_kJ(energy_kWh)
-        cop = self._calc_cop(
-            produced_heat_kJ=energy_kJ,
+        heat_energy_kW = self._cop_model.calc_q_from_p_kW(
             source_temp_C=source_temp_C,
-            time_slot=time_slot,
+            condenser_temp_C=self._charger.get_average_inlet_temperature_C(time_slot),
+            electrical_demand_kW=convert_kWh_to_kW(energy_kWh, GlobalConfig.slot_length),
         )
-        return energy_kJ * cop
+        return convert_kWh_to_kJ(convert_kW_to_kWh(heat_energy_kW, GlobalConfig.slot_length))
 
     def calc_energy_kWh_from_Q_kJ(self, time_slot: DateTime, Q_energy_kJ: float) -> float:
         """Calculate energy in kWh from heat in kJ."""
@@ -558,9 +574,9 @@ def _populate_state(self, time_slot: DateTime):
 
         if not self._heat_demand_Q_J:
             produced_heat_energy_kJ = self.combined_state.calc_Q_kJ_from_energy_kWh(
-                time_slot,
-                self._consumption_kWh.get_value(time_slot),
-                self._source_temp_C.get_value(time_slot),
+                time_slot=time_slot,
+                energy_kWh=self._consumption_kWh.get_value(time_slot),
+                source_temp_C=self._source_temp_C.get_value(time_slot),
             )
         else:
             produced_heat_energy_kJ = self._heat_demand_Q_J.get_value(time_slot) / 1000.0