diff --git a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/cop_models.py b/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/cop_models.py index 23028697ac..0b5f8f3c23 100644 --- a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/cop_models.py +++ b/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,10 +239,11 @@ 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 @@ -133,6 +251,24 @@ def calc_cop( 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 diff --git a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/AERMEC_NXP_0600_4L_COOL_model_parameters.json b/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/AERMEC_NXP_0600_4L_COOL_model_parameters.json new file mode 100644 index 0000000000..2ca253e3e8 --- /dev/null +++ b/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} diff --git a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/AERMEC_NXP_0600_4L_HEAT_model_parameters.json b/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/AERMEC_NXP_0600_4L_HEAT_model_parameters.json new file mode 100644 index 0000000000..73e28c0f57 --- /dev/null +++ b/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} diff --git a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/Elco_Aerotop_G07-14M_model_parameters.json b/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/Elco_Aerotop_G07-14M_model_parameters.json index d9ce3117b7..bc9569e4ca 100644 --- a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/Elco_Aerotop_G07-14M_model_parameters.json +++ b/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} diff --git a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/Elco_Aerotop_S09M-IR_model_parameters.json b/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/Elco_Aerotop_S09M-IR_model_parameters.json index 50ded42b77..650ff4ff1c 100644 --- a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/Elco_Aerotop_S09M-IR_model_parameters.json +++ b/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} diff --git a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/hoval_UltraSource_B_comfort_C_11_model_parameters.json b/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/hoval_UltraSource_B_comfort_C_11_model_parameters.json index e0a5222016..d327f4d918 100644 --- a/src/gsy_e/models/strategy/energy_parameters/heatpump/cop_models/model_data/hoval_UltraSource_B_comfort_C_11_model_parameters.json +++ b/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} diff --git a/src/gsy_e/models/strategy/energy_parameters/heatpump/heat_pump.py b/src/gsy_e/models/strategy/energy_parameters/heatpump/heat_pump.py index 47eedad618..82a4235b24 100644 --- a/src/gsy_e/models/strategy/energy_parameters/heatpump/heat_pump.py +++ b/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,9 +273,10 @@ 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 @@ -280,7 +284,11 @@ def update_cop_after_dis_charging( 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 diff --git a/src/gsy_e/setup/tekniker_cop_model/GSHP_model_fitter.py b/src/gsy_e/setup/tekniker_cop_model/GSHP_model_fitter.py new file mode 100644 index 0000000000..60858eda6f --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/GSHP_model_fitter.py @@ -0,0 +1,175 @@ +import json +import os + +import numpy as np +import pandas as pd +from scipy.optimize import curve_fit + +from gsy_e.setup.tekniker_cop_model.base_model_fitter import BaseModelFitter + +MODULE_PATH = os.path.dirname(os.path.abspath(__file__)) +INPUT_PATH = os.path.join(MODULE_PATH, "input_data") +EXPORT_PATH = os.path.join(MODULE_PATH, "model_parameters") + + +class GshpModelFitter(BaseModelFitter): + """Fitter for the GSHP model""" + + def __init__(self, filename: str, ref_evap_temp: int = 7, ref_cond_temp: int = 35): + assert filename.endswith(".csv") + self.calibration_data_filename = filename + self.raw_data = pd.read_csv(os.path.join(INPUT_PATH, filename)) + self.ref_evap_temp = ref_evap_temp + self.ref_cond_temp = ref_cond_temp + + def fit_end_export_data(self): + """Fit and export model data""" + export_filename = os.path.join( + EXPORT_PATH, + f"{self.calibration_data_filename.split('.')[0]}" f"_model_parameters.json", + ) + + model_parameters = self._fit_model() + json_model_parameters = {} + for name, value in model_parameters.items(): + json_model_parameters[name] = ( + value.tolist() if isinstance(value, np.ndarray) else value + ) + with open(export_filename, "w") as fp: + json.dump(json_model_parameters, fp) + + @staticmethod + def _poly22(X, a, b, c, d, e, f): + """ + Biquadratic polynomial self._poly22(Tevap, Tcond): + + f(Tev, Tco) = a + b*Tev + c*Tco + d*Tev^2 + e*Tev*Tco + f*Tco^2 + """ + Tev, Tco = X + return a + b * Tev + c * Tco + d * Tev**2 + e * Tev * Tco + f * Tco**2 + + @staticmethod + def _poly2(PLR, a, b, c): + """Quadratic polynomial: a*PLR^2 + b*PLR + c.""" + return a * PLR**2 + b * PLR + c + + def _fit_model(self): + df = self.raw_data.copy() + + # -------------------------------------------------------------- + # 1) Clean data and define reference point + # -------------------------------------------------------------- + df = df[df["Q"] != 0] # remove zero-capacity rows if any + + # Reference conditions for NXP: + df_ref = df[(df["Tevap"] == self.ref_evap_temp) & (df["Tcond"] == self.ref_cond_temp)] + if df_ref.empty: + raise ValueError( + f"Reference point (Tevap={self.ref_evap_temp}, Tcond={self.ref_cond_temp}) " + f"not found in dataset." + ) + + Qref = float(df_ref["Q"].iloc[0]) + Pref = float(df_ref["P"].iloc[0]) + + # -------------------------------------------------------------- + # 2) Compute CAPFT and EIRFT for full-load data + # -------------------------------------------------------------- + df["CAPFT"] = df["Q"] / Qref + df["EIRFT"] = df["P"] / (Pref * df["CAPFT"]) + + # Fit CAPFT(T, T) biquadratic + popt_capft, _ = curve_fit( + self._poly22, + (df["Tevap"].values, df["Tcond"].values), + df["CAPFT"].values, + ) + + # Fit EIRFT(T, T) biquadratic + popt_eirft, _ = curve_fit( + self._poly22, + (df["Tevap"].values, df["Tcond"].values), + df["EIRFT"].values, + ) + + # Evaluate CAPFT/EIRFT with fitted surfaces + df["CAPFT_eval"] = self._poly22( + (df["Tevap"].values, df["Tcond"].values), + *popt_capft, + ) + df["EIRFT_eval"] = self._poly22( + (df["Tevap"].values, df["Tcond"].values), + *popt_eirft, + ) + + # Compute PLR and EIRFPLR for these full-load points + df["PLR"] = df["Q"] / (Qref * df["CAPFT_eval"]) + df["EIRFPLR"] = df["P"] / (Pref * df["CAPFT_eval"] * df["EIRFT_eval"]) + + # -------------------------------------------------------------- + # 3) Generate synthetic part-load points + # -------------------------------------------------------------- + plr_targets = [0.75, 0.50, 0.25] + + tev_ref = float(df_ref["Tevap"].iloc[0]) + tco_ref = float(df_ref["Tcond"].iloc[0]) + + capft_ref = self._poly22((tev_ref, tco_ref), *popt_capft) + eirft_ref = self._poly22((tev_ref, tco_ref), *popt_eirft) + + synthetic_rows = [] + for plr in plr_targets: + Q_syn = Qref * capft_ref * plr + EER_syn = self.eer_fplr(plr) + P_syn = Q_syn / EER_syn + EIR_syn = P_syn / (Pref * capft_ref * eirft_ref) + + synthetic_rows.append( + { + "Tevap": tev_ref, + "Tcond": tco_ref, + "Q": Q_syn, + "P": P_syn, + "CAPFT_eval": capft_ref, + "EIRFT_eval": eirft_ref, + "PLR": plr, + "EIRFPLR": EIR_syn, + } + ) + + df_syn = pd.DataFrame(synthetic_rows) + df_plrfit = pd.concat([df, df_syn], ignore_index=True) + + # Use only reference-temperature rows for the PLR fit + df_ref_plr = df_plrfit[(df_plrfit["Tevap"] == tev_ref) & (df_plrfit["Tcond"] == tco_ref)] + + PLR_fit = df_ref_plr["PLR"].values + EIRFPLR_fit = df_ref_plr["EIRFPLR"].values + + # -------------------------------------------------------------- + # 4) Fit quadratic EIRFPLR(PLR) and inverse the order to fit to gsy-e implementation + # -------------------------------------------------------------- + popt_eirfplr, _ = curve_fit(self._poly2, PLR_fit, EIRFPLR_fit) + # -------------------------------------------------------------- + # 5) Pack output model dictionary + # -------------------------------------------------------------- + model = { + "CAPFT": popt_capft, + "HEIRFT": popt_eirft, + "HEIRFPLR": popt_eirfplr, + "Qref": Qref, + "Pref": Pref, + "Q_min": df["Q"].min(), + "Q_max": df["Q"].max(), + "PLR_min": PLR_fit.min(), + "COP_min": df["EER"].min(), + "COP_max": df["EER"].max(), + "COP_med": df["EER"].median(), + } + + return model + + +if __name__ == "__main__": + GshpModelFitter("AERMEC_NXP_0600_4L_COOL.csv").fit_end_export_data() + GshpModelFitter("AERMEC_NXP_0600_4L_HEAT.csv", ref_evap_temp=8).fit_end_export_data() diff --git a/src/gsy_e/setup/tekniker_cop_model/base_model_fitter.py b/src/gsy_e/setup/tekniker_cop_model/base_model_fitter.py new file mode 100644 index 0000000000..887664cdee --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/base_model_fitter.py @@ -0,0 +1,31 @@ +import numpy as np + + +class BaseModelFitter: + """Collection of methods that are used in all model fitters""" + + @staticmethod + def pack_model_into_dict( + capft: list, + heirft: list, + heirfplr: list, + q_ref: float, + p_ref: float, + Q: list, + cop: np.ndarray, + ): + # pylint: disable=too-many-positional-arguments, too-many-arguments + """Return model parameters.""" + + return { + "CAPFT": capft, + "HEIRFT": heirft, + "HEIRFPLR": heirfplr, + "Qref": q_ref, + "Pref": p_ref, + "Q_min": min(Q), + "Q_max": max(Q), + "COP_min": np.nanmin(cop), + "COP_max": np.nanmax(cop), + "COP_med": np.nanmedian(cop), + } diff --git a/src/gsy_e/setup/tekniker_cop_model/cop_model_fitter.py b/src/gsy_e/setup/tekniker_cop_model/cop_model_fitter.py index 7710535301..ea937e3a7c 100644 --- a/src/gsy_e/setup/tekniker_cop_model/cop_model_fitter.py +++ b/src/gsy_e/setup/tekniker_cop_model/cop_model_fitter.py @@ -7,10 +7,11 @@ import pandas as pd from scipy.optimize import curve_fit -from calibration_data_reader import CalibrationDataReader +from gsy_e.setup.tekniker_cop_model.calibration_data_reader import CalibrationDataReader +from gsy_e.setup.tekniker_cop_model.base_model_fitter import BaseModelFitter -class COPModelFitter: +class COPModelFitter(BaseModelFitter): """class that handles fitting teknikers' COP model""" def __init__(self, calibration_data_filename: str): @@ -63,6 +64,8 @@ def fit_model_q(self): data_fT["Q"][idx] = data_aux["Q_fT"][i, j] data_fT["P"][idx] = data_aux["P_fT"][i, j] + COP = data_fT["Q"] / data_fT["P"] + data_fT_df = pd.DataFrame(data_fT) # Remove NaN values data_fT_df = data_fT_df.dropna(subset=["Q"]) @@ -121,18 +124,19 @@ def poly2(PLR, b0, b1, b2): poly2, data_fTPLR_df["x_PLR"], data_fTPLR_df["x_heirfplr"] ) - # Save curve coefficients in a dictionary (the equivalent of a MATLAB struct) - model = { - "CAPFT": x_capft_params, - "HEIRFT": x_heirft_params, - "HEIRFPLR": x_heirfplr_params, - "Qref": self.calibration_data.q_ref, - "Pref": self.calibration_data.p_ref, - "PLR_min": self.calibration_data.plr_min, - } - return model + # Pack output model dictionary + return self.pack_model_into_dict( + x_capft_params, + x_heirft_params, + x_heirfplr_params, + self.calibration_data.q_ref, + self.calibration_data.p_ref, + data_fT["Q"].tolist(), + COP, + ) if __name__ == "__main__": - for filename in glob.glob("input_data/*.xlsx"): + for filename in glob.glob(os.path.join(os.path.dirname(__file__), "input_data/*.xlsx")): + print(filename) COPModelFitter(filename).export_model_parameters_to_json("model_parameters/") diff --git a/src/gsy_e/setup/tekniker_cop_model/input_data/AERMEC_NXP_0600_4L_COOL.csv b/src/gsy_e/setup/tekniker_cop_model/input_data/AERMEC_NXP_0600_4L_COOL.csv new file mode 100644 index 0000000000..c79c0ee811 --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/input_data/AERMEC_NXP_0600_4L_COOL.csv @@ -0,0 +1,30 @@ +Tcond,Tevap,Q,P,EER,hotWaterFlow,coldWaterFlow,load +25,5,143.1,25.7,5.57,7.9906,6.8339,100 +25,7,151.2,26,5.82,8.3886,7.2267,100 +25,9,159.3,26.3,6.06,8.7858,7.6197,100 +25,11,167.4,26.6,6.29,9.1858,8.0156,100 +25,13,175.6,27,6.51,9.5911,8.4169,100 +30,5,138.9,28,4.96,7.9069,6.6311,100 +30,7,147,28.3,5.2,8.3083,7.0278,100 +30,9,155.2,28.5,5.44,8.7111,7.4261,100 +30,11,163.5,28.9,5.67,9.1181,7.8286,100 +30,13,171.9,29.2,5.89,9.5322,8.2386,100 +30,15,180.5,29.6,6.1,9.9567,8.6583,100 +35,5,133.4,30.7,4.35,7.7833,6.3697,100 +35,7,141.5,30.9,4.58,8.1819,6.7636,100 +35,9,149.7,31.2,4.8,8.5833,7.1606,100 +35,11,158,31.5,5.02,8.9911,7.5639,100 +35,13,166.5,31.8,5.24,9.4078,7.9758,100 +35,15,175.2,32.2,5.44,9.8361,8.3994,100 +40,5,126.7,33.7,3.76,7.6164,6.0511,100 +40,7,134.7,33.9,3.97,8.0056,6.4358,100 +40,9,142.7,34.2,4.18,8.3997,6.8256,100 +40,11,150.9,34.4,4.38,8.8017,7.2228,100 +40,13,159.3,34.7,4.59,9.2142,7.6306,100 +40,15,167.9,35.1,4.79,9.6403,8.0511,100 +45,5,118.9,36.9,3.22,7.4028,5.6772,100 +45,7,126.5,37.1,3.41,7.7764,6.0467,100 +45,9,134.3,37.3,3.6,8.1564,6.4222,100 +45,11,142.2,37.5,3.79,8.5461,6.8072,100 +45,13,150.4,37.8,3.98,8.9481,7.2042,100 +45,15,158.9,38.1,4.17,9.3656,7.6156,100 diff --git a/src/gsy_e/setup/tekniker_cop_model/input_data/AERMEC_NXP_0600_4L_HEAT.csv b/src/gsy_e/setup/tekniker_cop_model/input_data/AERMEC_NXP_0600_4L_HEAT.csv new file mode 100644 index 0000000000..c84bea16b2 --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/input_data/AERMEC_NXP_0600_4L_HEAT.csv @@ -0,0 +1,40 @@ +Tcond,Tevap,Q,P,EER,hotWaterFlow,coldWaterFlow,load +25,5,169.4,26.7,6.33,8.1064,11.5425,100 +30,5,166.4,29,5.74,7.9761,11.1253,100 +35,5,163.1,31.7,5.15,7.8353,10.6569,100 +40,5,159.4,34.7,4.59,7.6706,10.1253,100 +45,5,154.9,37.9,4.09,7.4694,9.5178,100 +50,5,149.5,41.3,3.62,7.2189,8.8211,100 +55,5,142.7,44.5,3.21,6.9058,8.0239,100 +25,8,182.5,27.4,6.66,8.7322,12.5725,100 +30,8,179.3,29.6,6.06,8.5931,12.1408,100 +35,8,175.7,32.2,5.46,8.4375,11.6486,100 +40,8,171.5,35.1,4.88,8.2522,11.0831,100 +45,8,166.5,38.3,4.35,8.0247,10.4322,100 +50,8,160.3,41.6,3.86,7.7419,9.6831,100 +55,8,152.7,44.8,3.41,7.3906,8.8233,100 +25,11,195.9,28.1,6.97,9.3753,13.6333,100 +30,11,192.5,30.2,6.37,9.2281,13.1869,100 +35,11,188.7,32.8,5.75,9.0581,12.6706,100 +40,11,184.1,35.7,5.16,8.8528,12.0717,100 +45,11,178.5,38.8,4.6,8.5992,11.3775,100 +50,11,171.6,41.9,4.09,8.2844,10.5758,100 +55,11,163.2,45.1,3.62,7.8953,9.6536,100 +30,14,206.3,31,6.66,9.8933,14.2842,100 +35,14,202.2,33.5,6.04,9.7094,13.7439,100 +40,14,197.1,36.3,5.43,9.4847,13.1114,100 +45,14,191,39.3,4.85,9.2056,12.3739,100 +50,14,183.5,42.4,4.33,8.8592,11.5192,100 +55,14,174.3,45.5,3.83,8.4322,10.5347,100 +30,17,221.1,31.8,6.95,10.6017,15.4525,100 +35,17,216.6,34.3,6.31,10.4047,14.8883,100 +40,17,211.2,37.1,5.7,10.1606,14.2222,100 +45,17,204.5,40,5.11,9.8561,13.4414,100 +50,17,196.3,43.1,4.56,9.4789,12.5339,100 +55,17,186.4,46.1,4.05,9.015,11.4869,100 +30,18,226.3,32.1,7.04,10.8497,15.8614,100 +35,18,221.7,34.6,6.4,10.6483,15.2889,100 +40,18,216.1,37.4,5.78,10.3978,14.6117,100 +45,18,209.2,40.3,5.19,10.0853,13.8167,100 +50,18,200.8,43.3,4.64,9.6975,12.8917,100 +55,18,190.6,46.3,4.12,9.2214,11.8236,100 diff --git a/src/gsy_e/setup/tekniker_cop_model/input data/Elco_Aerotop_G07-14M.xlsx b/src/gsy_e/setup/tekniker_cop_model/input_data/Elco_Aerotop_G07-14M.xlsx similarity index 100% rename from src/gsy_e/setup/tekniker_cop_model/input data/Elco_Aerotop_G07-14M.xlsx rename to src/gsy_e/setup/tekniker_cop_model/input_data/Elco_Aerotop_G07-14M.xlsx diff --git a/src/gsy_e/setup/tekniker_cop_model/input data/Elco_Aerotop_S09M-IR.xlsx b/src/gsy_e/setup/tekniker_cop_model/input_data/Elco_Aerotop_S09M-IR.xlsx similarity index 100% rename from src/gsy_e/setup/tekniker_cop_model/input data/Elco_Aerotop_S09M-IR.xlsx rename to src/gsy_e/setup/tekniker_cop_model/input_data/Elco_Aerotop_S09M-IR.xlsx diff --git a/src/gsy_e/setup/tekniker_cop_model/input data/hoval_UltraSource_B_comfort_C_11.xlsx b/src/gsy_e/setup/tekniker_cop_model/input_data/hoval_UltraSource_B_comfort_C_11.xlsx similarity index 100% rename from src/gsy_e/setup/tekniker_cop_model/input data/hoval_UltraSource_B_comfort_C_11.xlsx rename to src/gsy_e/setup/tekniker_cop_model/input_data/hoval_UltraSource_B_comfort_C_11.xlsx diff --git a/src/gsy_e/setup/tekniker_cop_model/model_parameters/AERMEC_NXP_0600_4L_COOL_model_parameters.json b/src/gsy_e/setup/tekniker_cop_model/model_parameters/AERMEC_NXP_0600_4L_COOL_model_parameters.json new file mode 100644 index 0000000000..2ca253e3e8 --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/model_parameters/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} diff --git a/src/gsy_e/setup/tekniker_cop_model/model_parameters/AERMEC_NXP_0600_4L_HEAT_model_parameters.json b/src/gsy_e/setup/tekniker_cop_model/model_parameters/AERMEC_NXP_0600_4L_HEAT_model_parameters.json new file mode 100644 index 0000000000..73e28c0f57 --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/model_parameters/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} diff --git a/src/gsy_e/setup/tekniker_cop_model/model_parameters/Elco_Aerotop_G07-14M_model_parameters.json b/src/gsy_e/setup/tekniker_cop_model/model_parameters/Elco_Aerotop_G07-14M_model_parameters.json new file mode 100644 index 0000000000..bc9569e4ca --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/model_parameters/Elco_Aerotop_G07-14M_model_parameters.json @@ -0,0 +1 @@ +{"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} diff --git a/src/gsy_e/setup/tekniker_cop_model/model_parameters/Elco_Aerotop_S09M-IR_model_parameters.json b/src/gsy_e/setup/tekniker_cop_model/model_parameters/Elco_Aerotop_S09M-IR_model_parameters.json new file mode 100644 index 0000000000..650ff4ff1c --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/model_parameters/Elco_Aerotop_S09M-IR_model_parameters.json @@ -0,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} diff --git a/src/gsy_e/setup/tekniker_cop_model/model_parameters/hoval_UltraSource_B_comfort_C_11_model_parameters.json b/src/gsy_e/setup/tekniker_cop_model/model_parameters/hoval_UltraSource_B_comfort_C_11_model_parameters.json new file mode 100644 index 0000000000..d327f4d918 --- /dev/null +++ b/src/gsy_e/setup/tekniker_cop_model/model_parameters/hoval_UltraSource_B_comfort_C_11_model_parameters.json @@ -0,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} diff --git a/tests/strategies/energy_parameters/test_combined_heatpump_state.py b/tests/strategies/energy_parameters/test_combined_heatpump_state.py index 8b0903a856..831ce39d51 100644 --- a/tests/strategies/energy_parameters/test_combined_heatpump_state.py +++ b/tests/strategies/energy_parameters/test_combined_heatpump_state.py @@ -40,7 +40,7 @@ def test_get_energy_to_buy_maximum_kWh_calls_correct_methods(self, combined_stat CURRENT_MARKET_SLOT, 5000 ) combined_state._cop_model.calc_cop.assert_called_with( - source_temp_C=20, condenser_temp_C=30, heat_demand_kW=2.0 + source_temp_C=20, condenser_temp_C=30, heat_demand_kW=2.0, electrical_demand_kW=None ) def test_get_energy_to_buy_maximum_kWh_limits_cop_to_global_setting(self, combined_state): @@ -69,7 +69,7 @@ def test_get_energy_to_buy_minimum_kWh_calls_correct_methods(self, combined_stat CURRENT_MARKET_SLOT, 5000 ) combined_state._cop_model.calc_cop.assert_called_with( - source_temp_C=20, condenser_temp_C=30, heat_demand_kW=2.0 + source_temp_C=20, condenser_temp_C=30, heat_demand_kW=2.0, electrical_demand_kW=None ) @pytest.mark.parametrize( diff --git a/tests/strategies/energy_parameters/test_heat_pump.py b/tests/strategies/energy_parameters/test_heat_pump.py index 3e31201f43..a3f290f8ce 100644 --- a/tests/strategies/energy_parameters/test_heat_pump.py +++ b/tests/strategies/energy_parameters/test_heat_pump.py @@ -222,7 +222,7 @@ def test_cop_model_is_correctly_selected(energy_params_heat_profile): energy_params_heat_profile.event_market_cycle(CURRENT_MARKET_SLOT + duration(minutes=60)) assert isclose( energy_params_heat_profile._state.heatpump.get_cop(CURRENT_MARKET_SLOT), - 4.42, + 3.163, abs_tol=0.001, )