pypsa-eur/scripts/build_cop_profiles.py

# -*- coding: utf-8 -*-
# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
#
# SPDX-License-Identifier: MIT
"""
Build coefficient of performance (COP) time series for air- or ground-sourced
heat pumps.

The COP is approximated as a quatratic function of the temperature difference between source and
sink, based on Staffell et al. 2012.

This rule is executed in ``build_sector.smk``.

Relevant Settings
-----------------

.. code:: yaml
    heat_pump_sink_T:


Inputs:
-------
- ``resources/<run_name>/temp_soil_total_elec_s<simpl>_<clusters>.nc``: Soil temperature (total) time series.
- ``resources/<run_name>/temp_soil_rural_elec_s<simpl>_<clusters>.nc``: Soil temperature (rural) time series.
- ``resources/<run_name>/temp_soil_urban_elec_s<simpl>_<clusters>.nc``: Soil temperature (urban) time series.
- ``resources/<run_name>/temp_air_total_elec_s<simpl>_<clusters>.nc``: Ambient air temperature (total) time series.
- ``resources/<run_name>/temp_air_rural_elec_s<simpl>_<clusters>.nc``: Ambient air temperature (rural) time series.
- ``resources/<run_name>/temp_air_urban_elec_s<simpl>_<clusters>.nc``: Ambient air temperature (urban) time series.

Outputs:
--------
- ``resources/cop_soil_total_elec_s<simpl>_<clusters>.nc``: COP (ground-sourced) time series (total).
- ``resources/cop_soil_rural_elec_s<simpl>_<clusters>.nc``: COP (ground-sourced) time series (rural).
- ``resources/cop_soil_urban_elec_s<simpl>_<clusters>.nc``: COP (ground-sourced) time series (urban).
- ``resources/cop_air_total_elec_s<simpl>_<clusters>.nc``: COP (air-sourced) time series (total).
- ``resources/cop_air_rural_elec_s<simpl>_<clusters>.nc``: COP (air-sourced) time series (rural).
- ``resources/cop_air_urban_elec_s<simpl>_<clusters>.nc``: COP (air-sourced) time series (urban).


References
----------
[1] Staffell et al., Energy & Environmental Science 11 (2012): A review of domestic heat pumps, https://doi.org/10.1039/C2EE22653G.
"""

from typing import Union
from enum import Enum
import xarray as xr
import numpy as np
from _helpers import set_scenario_config


def coefficient_of_performance_individual_heating(delta_T, source="air"):
    if source == "air":
        return 6.81 - 0.121 * delta_T + 0.000630 * delta_T**2
    elif source == "soil":
        return 8.77 - 0.150 * delta_T + 0.000734 * delta_T**2
    else:
        raise NotImplementedError("'source' must be one of  ['air', 'soil']")


def celsius_to_kelvin(t_celsius: Union[float, xr.DataArray, np.array]) -> Union[float, xr.DataArray, np.array]:
    if (np.asarray(t_celsius) > 200).any():
        raise ValueError("t_celsius > 200. Are you sure you are using the right units?")
    return t_celsius + 273.15


def logarithmic_mean(t_hot: Union[float, xr.DataArray, np.ndarray], t_cold: Union[float, xr.DataArray, np.ndarray]) -> Union[float, xr.DataArray, np.ndarray]:
    if (np.asarray(t_hot <= t_cold)).any():
        raise ValueError("t_hot must be greater than t_cold")
    return (t_hot - t_cold) / np.log(t_hot / t_cold)

class CopDistrictHeating:

    def __init__(
        self,
        forward_temperature_celsius: Union[xr.DataArray, np.array],
        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
        return_temperature_celsius: Union[xr.DataArray, np.array],
        source_outlet_temperature_celsius: Union[xr.DataArray, np.array],
        delta_t_pinch_point: float = 5,
        isentropic_compressor_efficiency: float = 0.8,
        heat_loss: float = 0.0,
    ) -> None:
        """
        Initialize the COPProfileBuilder object.

        Parameters:
        ----------
        forward_temperature_celsius : Union[xr.DataArray, np.array]
            The forward temperature in Celsius.
        return_temperature_celsius : Union[xr.DataArray, np.array]
            The return temperature in Celsius.
        source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
            The source inlet temperature in Celsius.
        source_outlet_temperature_celsius : Union[xr.DataArray, np.array]
            The source outlet temperature in Celsius.
        delta_t_pinch_point : float, optional
            The pinch point temperature difference, by default 5.
        isentropic_compressor_efficiency : float, optional
            The isentropic compressor efficiency, by default 0.8.
        heat_loss : float, optional
            The heat loss, by default 0.0.
        """
        self.t_source_in = celsius_to_kelvin(source_inlet_temperature_celsius)
        self.t_sink_out = celsius_to_kelvin(forward_temperature_celsius)

        self.t_sink_in = celsius_to_kelvin(return_temperature_celsius)
        self.t_source_out = celsius_to_kelvin(source_outlet_temperature_celsius)

        self.isentropic_efficiency_compressor = isentropic_compressor_efficiency
        self.heat_loss = heat_loss
        self.delta_t_pinch = delta_t_pinch_point

    def cop(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the coefficient of performance (COP) for the system.

        Returns:
            Union[xr.DataArray, np.array]: The calculated COP values.
        """
        return (
            self.ideal_lorenz_cop
            * (
                (
                    1
                    + (self.delta_t_refrigerant_sink + self.delta_t_pinch)
                    / self.t_sink_mean
                )
                / (
                    1
                    + (
                        self.delta_t_refrigerant_sink
                        + self.delta_t_refrigerant_source
                        + 2 * self.delta_t_pinch
                    )
                    / self.delta_t_lift
                )
            )
            * self.isentropic_efficiency_compressor
            * (1 - self.ratio_evaporation_compression_work)
            + 1
            - self.isentropic_efficiency_compressor
            - self.heat_loss
        )

    @property
    def t_sink_mean(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the logarithmic mean temperature difference between the cold and hot sinks.

        Returns
        -------
        Union[xr.DataArray, np.array]
            The mean temperature difference.
        """
        return logarithmic_mean(t_cold=self.t_sink_in, t_hot=self.t_sink_out)

    @property
    def t_source_mean(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the logarithmic mean temperature of the heat source.

        Returns
        -------
        Union[xr.DataArray, np.array]
            The mean temperature of the heat source.
        """
        return logarithmic_mean(t_hot=self.t_source_in, t_cold=self.t_source_out)

    @property
    def delta_t_lift(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the temperature lift as the difference between the logarithmic sink and source temperatures.

        Returns
        -------
        Union[xr.DataArray, np.array]
            The temperature difference between the sink and source.
        """
        return self.t_sink_mean - self.t_source_mean

    @property
    def ideal_lorenz_cop(self) -> Union[xr.DataArray, np.array]:
        """
        Ideal Lorenz coefficient of performance (COP).

        The ideal Lorenz COP is calculated as the ratio of the mean sink temperature
        to the lift temperature difference.

        Returns
        -------
        np.array
            The ideal Lorenz COP.

        """
        return self.t_sink_mean / self.delta_t_lift

    @property
    def delta_t_refrigerant_source(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the temperature difference between the refrigerant source inlet and outlet.

        Returns
        -------
        Union[xr.DataArray, np.array]
            The temperature difference between the refrigerant source inlet and outlet.
        """
        return self._approximate_delta_t_refrigerant_source(
            delta_t_source=self.t_source_in - self.t_source_out
        )

    @property
    def delta_t_refrigerant_sink(self) -> Union[xr.DataArray, np.array]:
        """
        Temperature difference between the refrigerant and the sink based on approximation.

        Returns
        -------
        Union[xr.DataArray, np.array]
            The temperature difference between the refrigerant and the sink.
        """
        return self._approximate_delta_t_refrigerant_sink()

    @property
    def ratio_evaporation_compression_work(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the ratio of evaporation to compression work based on approximation.

        Returns
        -------
        Union[xr.DataArray, np.array]
            The calculated ratio of evaporation to compression work.
        """
        return self._ratio_evaporation_compression_work_approximation()
    
    @property
    def delta_t_sink(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the temperature difference at the sink.

        Returns
        -------
        Union[xr.DataArray, np.array]
            The temperature difference at the sink.
        """
        return self.t_sink_out - self.t_sink_in

    def _approximate_delta_t_refrigerant_source(
        self, delta_t_source: Union[xr.DataArray, np.array]
    ) -> Union[xr.DataArray, np.array]:
        """
        Approximates the temperature difference between the refrigerant and the source.

        Parameters
        ----------
        delta_t_source : Union[xr.DataArray, np.array]
            The temperature difference for the refrigerant source.

        Returns
        -------
        Union[xr.DataArray, np.array]
            The approximate temperature difference for the refrigerant source.
        """
        return delta_t_source / 2

    def _approximate_delta_t_refrigerant_sink(
        self,
        refrigerant: str = "ammonia",
        a: float = {"ammonia": 0.2, "isobutane": -0.0011},
        b: float = {"ammonia": 0.2, "isobutane": 0.3},
        c: float = {"ammonia": 0.016, "isobutane": 2.4},
    ) -> Union[xr.DataArray, np.array]:
        """
        Approximates the temperature difference at the refrigerant sink.

        Parameters:
        ----------
        refrigerant : str, optional
            The refrigerant used in the system. Either 'isobutane' or 'ammonia. Default is 'ammonia'.
        a : float, optional
            Coefficient for the temperature difference between the sink and source, default is 0.2.
        b : float, optional
            Coefficient for the temperature difference at the sink, default is 0.2.
        c : float, optional
            Constant term, default is 0.016.

        Returns:
        -------
        Union[xr.DataArray, np.array]
            The approximate temperature difference at the refrigerant sink.

        Notes:
        ------
        This function assumes ammonia as the refrigerant.

        The approximate temperature difference at the refrigerant sink is calculated using the following formula:
        a * (t_sink_out - t_source_out + 2 * delta_t_pinch) + b * delta_t_sink + c

        """
        if refrigerant not in a.keys():
            raise ValueError(
                f"Invalid refrigerant '{refrigerant}'. Must be one of {a.keys()}"
            )
        return (
            a[refrigerant]
            * (self.t_sink_out - self.t_source_out + 2 * self.delta_t_pinch)
            + b[refrigerant] * self.delta_t_sink
            + c[refrigerant]
        )

    def _ratio_evaporation_compression_work_approximation(
        self,
        refrigerant: str = "ammonia",
        a: float = {"ammonia": 0.0014, "isobutane": 0.0035},
        b: float = {"ammonia": -0.0015, "isobutane": -0.0033},
        c: float = {"ammonia": 0.039, "isobutane": 0.053},
    ) -> Union[xr.DataArray, np.array]:
        """
        Calculate the ratio of evaporation to compression work approximation.

        Parameters:
        ----------
        refrigerant : str, optional
            The refrigerant used in the system. Either 'isobutane' or 'ammonia. Default is 'ammonia'.
        a : float, optional
            Coefficient 'a' in the approximation equation. Default is 0.0014.
        b : float, optional
            Coefficient 'b' in the approximation equation. Default is -0.0015.
        c : float, optional
            Coefficient 'c' in the approximation equation. Default is 0.039.

        Returns:
        -------
        Union[xr.DataArray, np.array]
            The calculated ratio of evaporation to compression work.

        Notes:
        ------
        This function assumes ammonia as the refrigerant.

        The approximation equation used is:
        ratio = a * (t_sink_out - t_source_out + 2 * delta_t_pinch) + b * delta_t_sink + c
        """
        if refrigerant not in a.keys():
            raise ValueError(
                f"Invalid refrigerant '{refrigerant}'. Must be one of {a.keys()}"
            )
        return (
            a[refrigerant]
            * (self.t_sink_out - self.t_source_out + 2 * self.delta_t_pinch)
            + b[refrigerant] * self.delta_t_sink
            + c[refrigerant]
        )


if __name__ == "__main__":
    if "snakemake" not in globals():
        from _helpers import mock_snakemake

        snakemake = mock_snakemake(
            "build_cop_profiles",
            simpl="",
            clusters=48,
        )

    set_scenario_config(snakemake)

    for source in ["air", "soil"]:
        source_T = xr.open_dataarray(snakemake.input[f"temp_{source}_total"])

        delta_T = snakemake.params.heat_pump_sink_T_individual_heating - source_T

        cop_individual_heating = coefficient_of_performance_individual_heating(delta_T, source)
        cop_individual_heating.to_netcdf(snakemake.output[f"cop_{source}_decentral_heating"])

        cop_district_heating = CopDistrictHeating(
            forward_temperature_celsius=snakemake.params.forward_temperature_district_heating,
            return_temperature_celsius=snakemake.params.return_temperature_district_heating,
            source_inlet_temperature_celsius=source_T,
            source_outlet_temperature_celsius=source_T - snakemake.params.heat_source_cooling_district_heating,
        ).cop()

        cop_district_heating.to_netcdf(snakemake.output[f"cop_{source}_central_heating"])