Py.Cafe

# run locally with
# pip install solara
# solara run app.py
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import typing

# data downloaded from https://www.knmi.nl/nederland-nu/klimatologie/daggegevens
# (this file is for Lauwersoog)
# Read the KNMI daily data
file_path = "etmgeg_277.txt"  # Replace this with your file path
data_all = pd.read_csv(file_path, skiprows=51, delimiter=",")  # Read the file, skipping the header rows
data_all.columns = data_all.columns.str.strip()

# Convert temperature to degrees Celsius
data_all["TG"] = data_all["TG"] / 10
data_all["TN"] = data_all["TN"] / 10
data_all["Date"] = pd.to_datetime(data_all["YYYYMMDD"], format="%Y%m%d")


class Device(typing.TypedDict):
    Prated: float
    T: typing.List[float]
    COP35: typing.List[float]
    Power35: typing.List[float]
    COP55: typing.List[float]
    Power55: typing.List[float]
    url: typing.Optional[str]


devices = {
    "Itho Amber 65": {
        "Prated": 6.5,  # Rated power in kW
        "T": [-7, 7],
        "COP35": [3.1, 4.4],
        "Power35": [4.3, 6.8],
        "COP55": [2.2, 2.9],
        "Power55": [3.9, 6.1],
        "url": "https://ithodaalderop.compano.com/Data/Environments/000001/Attachment/Bijlage/A02_Warmtepompen/A02_02_Lucht/A02_02_07_Amber/B01_06_MTE/01-05519_002%20handleiding%20Amber.pdf",
    },
    "WeHeat Flint P40": {
        "Prated": 4.0,  # Rated power in kW
        "T": [-10, -7, 2, 7, 12],  # Temperatures in degrees Celsius
        "COP35": [2.9, 3.2, 4.7, 6.4, 8.9],  # Corresponding COP values
        "Power35": [3.1, 3.6, 2.4, 2.1, 2.6],  # Power in kW
        "COP55": [2.0, 2.3, 3.6, 4.6, 6.9],  # Corresponding COP values
        "Power55": [3.7, 4.1, 2.5, 2.2, 2.7],  # Power in kW
        "url": "https://content.weheat.nl/weheat-flint-product-specification",
    },
    "WeHeat Sparrow P60": {
        "Prated": 6.0,  # Rated power in kW
        "T": [-10, -7, 2, 7, 12],  # Temperatures in degrees Celsius
        "COP35": [2.8, 2.9, 4.5, 6.8, 8.6],  # Corresponding COP values
        "Power35": [5.5, 4.7, 2.9, 2.2, 2.2],  # Power in kW
        "COP55": [2.1, 2.3, 3.6, 5.0, 6.9],  # Corresponding COP values
        "Power55": [5.6, 5.0, 3.0, 2.3, 2.2],  # Power in kW
        "url": "https://content.weheat.nl/weheat-sparrow-product-specification",
    },
    "WeHeat Blackbird P80": {
        "Prated": 8.0,  # Rated power in kW
        "T": [-10, -7, 2, 7, 12],  # Temperatures in degrees Celsius
        "COP35": [2.7, 3.2, 4.5, 6.4, 8.4],  # Corresponding COP values
        "Power35": [7.5, 6.7, 4.4, 2.9, 3.4],  # Power in kW
        "COP55": [1.9, 2.2, 3.4, 4.3, 5.7],  # Corresponding COP values
        "Power55": [7.5, 6.6, 4.0, 2.6, 1.15],  # Power in kW
        "url": "https://content.weheat.nl/weheat-blackbird-product-specification",
    },
    "Quatt Hybrid": {
        "Prated": 4.5,  # Rated power in kW
        "T": [-7, 2, 7],
        "COP35": [3.1, 4, 4.8],
        "Power35": [3.3, 4.2, 4.2],  # read from graph, not very accurate
        "COP55": [1.9, 2.3, 2.9],
        "Power55": [3.4, 4.0, 3.8],  # read from graph, not very accurate
        "url": "https://support.quatt.io/nl/articles/6121833-wat-zijn-de-specificaties-van-een-quatt-warmtepomp",
    },
    "Quatt Hybrid Duo": {
        "Prated": 9.0,  # Rated power in kW
        "T": [-7, 2, 7],
        "COP35": [3.1, 4, 4.8],
        "Power35": [3.3 * 2, 4.2 * 2, 4.2 * 2],  # read from graph, not very accurate
        "COP55": [1.9, 2.3, 2.9],
        "Power55": [3.4 * 2, 4.0 * 2, 3.8 * 2],  # read from graph, not very accurate
        "url": "https://support.quatt.io/nl/articles/6121833-wat-zijn-de-specificaties-van-een-quatt-warmtepomp",
    },
}


def compute_all(
    year: int,
    price_per_kwh: float,
    price_per_cubic_meter: float,
    electric_power_kwh: float,
    gas_usage_per_degree_day: float,
    indoor_temperature: float,
    caloric_value: float,
    device: Device,
    T_water: float,
    daily_hours: int,
    use_minumum: bool,
):
    # Filter for the required columns and process temperature
    # Create a column for the date

    # Select only the data for the year 2024
    data = data_all.copy()
    data = data[data["Date"].dt.year == year]
    T = data["TN"] if use_minumum else data["TG"]

    # Add a COP column using interpolation

    COP = [np.interp(T_water, [35, 55], [device["COP35"][i], device["COP55"][i]]) for i in range(len(device["T"]))]
    print("T", T)
    data["COP"] = np.interp(T, device["T"], COP)

    # Parameters
    daily_hours = daily_hours  # Operating hours per day

    # Calculate DeltaT and energy usage
    data["DeltaT"] = np.maximum(indoor_temperature - T, 0)  # Degree days
    data["Heat_Energy_Usage_m3"] = data["DeltaT"] * gas_usage_per_degree_day  # Gas usage in m³
    data["Heat_Energy_Usage_kWh"] = data["Heat_Energy_Usage_m3"] * caloric_value  # Convert to kWh

    data["Cumulative_Heat_Energy_Usage"] = data["Heat_Energy_Usage_m3"].cumsum()
    data["Cumulative_Heat_Energy_Usage_kWh"] = data["Cumulative_Heat_Energy_Usage"] * caloric_value  # Convert to kWh

    # Calculations
    data["Max_Electricity_Usage_kWh"] = electric_power_kwh * daily_hours  # Daily electricity usage
    data["Max_HP_Energy_kWh"] = data["Max_Electricity_Usage_kWh"]  # Max energy by heat pump
    # data["Cumulative_Max_HP_Energy_kWh"] = data["Max_HP_Energy_kWh"].cumsum()

    # Energy delivered by the heat pump, which is the minimum of Max_HP_Energy_kWh and Heat_Energy_Usage_kWh
    data["Heat_Energy_Delivered_HP_kWh"] = data[["Max_HP_Energy_kWh", "Heat_Energy_Usage_kWh"]].min(axis=1)
    # Energy not delivered by the heat pump
    data["Heat_Energy_Not_Delivered_HP_kWh"] = data["Heat_Energy_Usage_kWh"] - data["Heat_Energy_Delivered_HP_kWh"]
    data["Electrical_Energy_Usage_HP_kWh"] = data["Heat_Energy_Delivered_HP_kWh"] / data["COP"]  # Energy usage by heat pump
    data["Heat_Energy_Usage_Not_HP_kWh"] = data["Heat_Energy_Usage_kWh"] - data["Heat_Energy_Delivered_HP_kWh"]
    data["Energy_Usage_kWh"] = data["Electrical_Energy_Usage_HP_kWh"] + data["Heat_Energy_Usage_Not_HP_kWh"]

    data["Electricity_Cost"] = data["Electrical_Energy_Usage_HP_kWh"] * price_per_kwh
    data["Gas_Cost"] = data["Heat_Energy_Not_Delivered_HP_kWh"] / caloric_value * price_per_cubic_meter
    data["Total_Cost"] = data["Electricity_Cost"] + data["Gas_Cost"]
    data["Gas_Hypothetical_Cost"] = data["Heat_Energy_Usage_kWh"] / caloric_value * price_per_cubic_meter

    # For each temperature range, e.g., [6, 7), calculate the sum of Heat_Energy_Delivered_HP_kWh
    dft = data.groupby(pd.cut(T, bins=range(-10, 20)))["Heat_Energy_Delivered_HP_kWh"].sum()

    # View a sample of the results
    print(data[["YYYYMMDD", "TG", "COP"]].head())
    return data


def plot_stats(data, caloric_value: float, device: Device, T_water: float, use_minumum: bool, indoor_temperature: float, gas_usage_per_degree_day: float, electric_power_kwh: float):
    # filter out temperatures above 18 and below -10
    T = data["TN"] if use_minumum else data["TG"]
    data = data[T < 18]
    data = data[T > -10]
    T = data["TN"] if use_minumum else data["TG"]
    # Plotting
    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 12))

    # Define discrete bins and group data
    bins = range(-10, 20)  # Discrete temperature bins
    bin_labels = list(bins[:-1])  # Use discrete values as labels
    # dft = data.groupby(T.round().astype(int))["Heat_Energy_Delivered_HP_kWh"].sum()
    dft = data.groupby(T.round().astype(int))[["Heat_Energy_Delivered_HP_kWh", "Electrical_Energy_Usage_HP_kWh"]].sum()

    print("dft", dft)

    # Plot bars in the center of the bins
    ax1.bar(dft.index + 0.5, dft["Heat_Energy_Delivered_HP_kWh"].values, width=0.8, align="center", label="Heat Energy Delivered by heat pump")
    ax1.bar(dft.index + 0.5, dft["Electrical_Energy_Usage_HP_kWh"].values, width=0.8, align="center", label="Electrical Energy Used by heat pump")

    # Customize x-axis to show discrete values
    ax1.set_xticks(dft.index)
    ax1.set_xticklabels(dft.index, rotation=45, ha="right")

    # Add a new axis for the COP curve
    ax_cop = ax1.twinx()  # Create a twin Y-axis for COP
    temp = np.arange(-10, 20, 1.0) + 0.5  # Temperature range
    #
    COP_at_T_water = [np.interp(T_water, [35, 55], [device["COP35"][i], device["COP55"][i]]) for i in range(len(device["T"]))]
    Power_at_T_water = [np.interp(T_water, [35, 55], [device["Power35"][i], device["Power55"][i]]) for i in range(len(device["T"]))]
    data["COP"] = np.interp(T, device["T"], COP_at_T_water)
    cop = np.interp(temp, device["T"], COP_at_T_water)

    ax_cop.plot(temp, cop, label="COP", color="tab:green", linestyle="solid")
    ax_cop.scatter(temp, cop, label="COP", color="tab:green")
    ax_cop.set_ylabel("COP", color="tab:green")
    ax_cop.tick_params(axis="y", labelcolor="tab:green")
    ax_cop.spines["right"].set_position(("outward", 60))  # Offset the secondary Y-axis outward by 60 points

    # Add a legend for both
    handles1, labels1 = ax1.get_legend_handles_labels()
    handles2, labels2 = ax_cop.get_legend_handles_labels()
    ax1.legend(handles1 + handles2, labels1 + labels2, loc="upper left")

    # Add labels and title
    ax1.set_ylabel("Energy (kWh)")
    ax1.set_xlabel("Outdoor Temperature (°C)")
    plt.title("Energy Delivered and COP Across Temperature Bins")
    # if we change the parameters, the plot y max will continously change, making it difficult to compare different parameters
    # so we change the y limit max by a discretly rounded up value, by 30%
    max = ax1.get_ylim()[1]
    base = 10 ** int(np.log10(max))
    multiple = base // 2
    discete_max = int(max / multiple + 1) * multiple
    ax1.set_ylim(0, discete_max)
    print("max", max, discete_max)

    # On the right Y-axis, show the energy in m³ (1 m³ = caloric_value kWh)
    ax1b = ax1.twinx()
    ax1b.set_ylim(ax1.get_ylim())
    caloric_value = 8.8  # Example conversion factor
    formatter = plt.FuncFormatter(lambda x, _: f"{x / caloric_value:.0f}")
    ax1b.yaxis.set_major_formatter(formatter)
    ax1b.set_ylabel("Energy (m³)")

    ax2.scatter(temp, cop, label="Interpolated COP used for model", color="tab:green")
    ax2.scatter(device["T"], device["COP35"], label="COP at 35 deg from spec sheet", marker="x", color="tab:orange")
    ax2.scatter(device["T"], device["COP55"], label="COP at 55 deg from spec sheet", marker="x", color="tab:green")
    ax2.scatter(device["T"], COP_at_T_water, label="Interpolated COP from spec sheet", marker="x", color="tab:red")
    
    # DeltaT = np.maximum(indoor_temperature - temp, 0)  # Degree days
    # Thermal_load = DeltaT * gas_usage_per_degree_day * caloric_value/24 / electric_power_kwh

    # # https://www.nrel.gov/docs/fy23osti/85081.pdf
    # print("Thermal_load", Thermal_load)
    # print("DeltaT", DeltaT)
    # COP = -7.46*(Thermal_load - 0.0047 *temp  - 0.477)**2 + 0.0941*temp + 4.34
    # ax2.scatter(temp, COP, label="COP at load (theoretical modal)", marker="x", color="tab:red")

    # for Thermal_load in [0.2, 0.5, 0.8, 1.0]:
    #     # https://www.nrel.gov/docs/fy23osti/85081.pdf
    #     COP = -7.46*(Thermal_load - 0.0047 *temp  - 0.477)**2 + 0.0941*temp + 4.34
    #     ax2.scatter(temp, COP, label=f"COP at load={Thermal_load} (theoretical modal)", marker="x", color="tab:blue")


    ax2.set_ylabel("COP")
    ax2.set_xlabel("Temperature (°C)")
    ax2.legend()
    ax2.grid()

    Power = data["Heat_Energy_Delivered_HP_kWh"] / 24
    # ax3.scatter(T, Power/electric_power_kwh, label="Power (in form of heat) delivered by heat pump", color="tab:blue")
    # ax3.scatter(device["T"], device["Power35"]/electric_power_kwh, label="Power at 35 deg from spec sheet", marker="x", color="tab:orange")
    # ax3.scatter(device["T"], device["Power55"]/electric_power_kwh, label="Power at 55 deg from spec sheet", marker="x", color="tab:green")
    ax3.scatter(T, Power, label="Power (in form of heat) delivered by heat pump", color="tab:blue")
    ax3.scatter(device["T"], device["Power35"], label="Power at 35 deg from spec sheet", marker="x", color="tab:orange")
    ax3.scatter(device["T"], device["Power55"], label="Power at 55 deg from spec sheet", marker="x", color="tab:green")
    ax3.scatter(device["T"], Power_at_T_water, label="Power at water temperature (interpolated)", marker="x", color="tab:red")
    ax3.set_ylabel("Electric Power (kWh)")
    ax3.set_xlabel("Temperature (°C)")
    ax3.legend()

    plt.tight_layout()
    plt.show()


def plot_usage(data, caloric_value=8.8, gas_usage_per_degree_day=1.2, indoor_temperature: float = 20, use_minumum: bool = False):
    T = data["TN"] if use_minumum else data["TG"]
    # Create subplots: 2 rows, 1 column
    fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(10, 16))

    ax1, ax2, ax3, ax4 = ax3, ax1, ax4, ax2

    # Left Y-axis: Temperature (°C)
    ax1.step(data["Date"], T, label="Outdoor temperature (°C)", where="post")
    ax1.set_xlabel("Date")
    ax1.set_ylabel("Temperature (°C)")
    ax1.tick_params(axis="y")
    ax1.grid()

    # Second Y-axis (right): COP
    ax1c = ax1.twinx()  # Create a twin Y-axis
    ax1c.spines["right"].set_position(("axes", 1.0))  # Attach to the right edge of the axes
    ax1c.step(data["Date"], data["COP"], label="COP", color="tab:green", where="post")
    ax1c.set_ylabel("COP", color="tab:green")
    ax1c.tick_params(axis="y", labelcolor="tab:green")

    # Combine legends
    handles1, labels1 = ax1.get_legend_handles_labels()
    handles2, labels2 = ax1c.get_legend_handles_labels()
    ax1.legend(handles1 + handles2, labels1 + labels2, loc="upper left")

    # Plot stacked areas
    ax2.fill_between(data["Date"], 0, data["Heat_Energy_Delivered_HP_kWh"], label="Heat delivered by heat pump", color="tab:green", step="mid")
    ax2.fill_between(
        data["Date"], data["Heat_Energy_Delivered_HP_kWh"], data["Heat_Energy_Usage_kWh"], label="Heat Energy required", color="tab:blue", step="mid"
    )

    # Add a separate line for Max_HP_Energy_kWh
    ax2.plot(data["Date"], data["Max_HP_Energy_kWh"], label="Max delivered by heat pump", color="tab:red", linestyle="--", linewidth=2)

    # Add a line for Energy Usage not by heat pump
    ax2.step(data["Date"], data["Heat_Energy_Not_Delivered_HP_kWh"], label="Heat *NOT* delivered by heat pump", color="tab:red", linewidth=2, where="mid")

    ax2b = ax2.twinx()
    ax2b.set_ylim(ax2.get_ylim())
    formatter = plt.FuncFormatter(lambda x, _: f"{x / caloric_value:.0f}")
    ax2b.yaxis.set_major_formatter(formatter)
    ax2b.set_ylabel("Daily Heat Energy Usage (m³)")

    ax2.set_xlabel("Date")
    ax2.set_ylabel("Daily Heat Energy Usage (kWh)")
    ax2.set_title("Daily Heat Energy Usage")
    ax2.legend()
    ax2.grid()

    # Cumulative or daily energy data
    ax4.plot(data["Date"], data["Heat_Energy_Delivered_HP_kWh"].cumsum(), label="Heat delivered by heat pump", color="tab:green")
    ax4.plot(data["Date"], data["Heat_Energy_Not_Delivered_HP_kWh"].cumsum(), label="Heat/Energy(gas?) not delivered/used by heat pump", color="tab:red")
    ax4.plot(data["Date"], data["Electrical_Energy_Usage_HP_kWh"].cumsum(), label="Electricity used by heat pump", color="tab:green", linestyle="--")
    # same as Heat_Energy_Not_Delivered_HP_kWh since COP=1
    # ax4.plot(data["Date"], data["Heat_Energy_Usage_Not_HP_kWh"].cumsum(), label="Energy not used by heat pump", color="tab:orange", linestyle="--")
    ax4.plot(data["Date"], data["Energy_Usage_kWh"].cumsum(), label="Total (electrical+gas?) Energy Usage", color="tab:blue", linestyle="--")

    ax4c = ax4.twinx()  # Create a secondary Y-axis for cost
    ax4c.plot(data["Date"], data["Total_Cost"].cumsum(), label="Total Cost (€)", color="tab:purple", linestyle="-")
    ax4c.plot(data["Date"], data["Gas_Hypothetical_Cost"].cumsum(), label="(Pure gas based scenario) Gas Cost (€)", color="tab:red", linestyle=":")
    ax4c.set_ylabel("Total Cost (€)", color="tab:purple")
    ax4c.tick_params(axis="y", labelcolor="tab:purple")
    ax4c.spines["right"].set_position(("outward", 60))  # Offset the secondary Y-axis outward by 60 points

    # Adjust legend to include both axes
    handles1, labels1 = ax4.get_legend_handles_labels()
    handles2, labels2 = ax4c.get_legend_handles_labels()
    ax4.legend(handles1 + handles2, labels1 + labels2, loc="upper left")

    ax4.set_ylabel("Cumulative Heat Energy Usage (kWh)")
    ax4.set_xlabel("Date")
    ax4.set_title("Cumulative Heat Energy Usage")
    ax4.grid()

    ax4b = ax4.twinx()
    ax4b.set_ylim(ax4.get_ylim())
    formatter = plt.FuncFormatter(lambda x, _: f"{x / caloric_value:.0f}")
    ax4b.yaxis.set_major_formatter(formatter)
    ax4b.set_ylabel("CumulativeHeat Energy Usage (m³)")

    # predicted house temperature if only heat pump is used
    deltaT = data["Heat_Energy_Delivered_HP_kWh"] / caloric_value / gas_usage_per_degree_day
    data["Predicted_Temperature_HP"] = T + deltaT

    # Plot the predicted temperature
    ax3.plot(data["Date"], T, label="Outdoor Temperature (°C)")
    ax3.plot(data["Date"], data["Predicted_Temperature_HP"], label="Indoor Temperature (°C) - Heat Pump Only")
    # indoor temperature
    ax3.axhline(y=indoor_temperature, color="tab:red", linestyle="--", label="Desired Indoor Temperature")
    ax3.set_xlabel("Date")
    ax3.set_ylabel("Temperature (°C)")
    ax3.set_title("Predicted indoor Temperature if Only Heat Pump is Used")
    ax3.legend()

    # Title and layout adjustment
    fig.tight_layout()
    plt.show()


# plot_usage(data)

import solara


year_min = solara.reactive(data_all["Date"].dt.year.min().item())
year_max = solara.reactive(data_all["Date"].dt.year.max().item())
caloric_value = solara.reactive(8.8)  # Caloric value of gas in kWh/m³
year = solara.reactive(2024)
gas_usage_per_degree_day = solara.reactive(0.8)  # Gas consumption in m³ per degree day
device_name = solara.reactive("WeHeat Sparrow P60")

price_per_kwh = solara.reactive(0.28292)  # Price per kWh
price_per_cubic_meter = solara.reactive(1.35)
daily_hours = solara.reactive(21)
use_minumum = solara.reactive(False)


electric_power_kwh = solara.reactive(devices[device_name.value]["Prated"])  # Rated power in kW
indoor_temperature = solara.reactive(20.0)  # in °C
water_temperature = solara.reactive(45.0)  # in °C


@solara.component
def Parameters():
    with solara.AppBar():
        solara.lab.ThemeToggle()
    with solara.Card("Input Parameters"):
        with solara.Column(gap="28px"):
            solara.Select("Device", value=device_name.value, on_value=set_device, values=list(devices.keys()))
            with solara.Tooltip("Max number of hours the HP operates, by default 3 hours are reserved for defrost."):
                solara.SliderInt("Heating hours", value=daily_hours, min=10, max=24, thumb_label="always")
            solara.SliderInt("Year", min=year_min.value, max=year_max.value, value=year, thumb_label="always")
            power_comment = (
                "Using default power for the selected device" if electric_power_kwh.value == devices[device_name.value]["Prated"] else "Using custom power"
            )
            with solara.Tooltip("Water temperature in °C. This is the temperature of the water in the heating system."):
                solara.SliderFloat("Water Temperature (°C)", min=35.0, max=55.0, value=water_temperature, thumb_label="always")
            with solara.Tooltip(power_comment):
                solara.SliderFloat("Heatpump Power (kW)", min=0.0, max=12.0, value=electric_power_kwh, thumb_label="always")
            with solara.Tooltip(
                "Gas Usage per Degree Day (m³), this is how much gas is consumed per day, per degree Celsius difference from the reference temperature. mindergas.nl is a good source for this information."
            ):
                solara.SliderFloat("GUDD (m³/deg/day)", min=0.01, max=2.0, value=gas_usage_per_degree_day, thumb_label="always")
            with solara.Tooltip("Caloric value of gas in kWh/m³. This is the energy content of the gas in kWh per m³. The default value is 8.8 kWh/m³."):
                solara.SliderFloat("Caloric Value (kWh/m³)", min=5.0, max=12.0, value=caloric_value, thumb_label="always")
            with solara.Tooltip("Indoor temperature in °C. This is the desired indoor temperature."):
                solara.SliderFloat("Indoor Temperature (°C)", min=15.0, max=25.0, value=indoor_temperature, thumb_label="always")
            with solara.Tooltip("Use minimum temperature instead of average temperature."):
                solara.Switch(label="Use minimum temperature", value=use_minumum)
            solara.InputFloat("Price per kWh (€)", value=price_per_kwh, continuous_update=True)
            solara.InputFloat("Price per m³ (€)", value=price_per_cubic_meter, continuous_update=True)


@solara.component
def Page():
    dark_effective = solara.lab.use_dark_effective()
    plt.style.use("dark_background" if dark_effective else "default")
    device = devices[device_name.value]
    data = compute_all(
        year=year.value,
        electric_power_kwh=electric_power_kwh.value,
        gas_usage_per_degree_day=gas_usage_per_degree_day.value,
        indoor_temperature=indoor_temperature.value,
        caloric_value=caloric_value.value,
        device=device,
        T_water=water_temperature.value,
        price_per_kwh=price_per_kwh.value,
        price_per_cubic_meter=price_per_cubic_meter.value,
        daily_hours=daily_hours.value,
        use_minumum=use_minumum.value,
    )
    print("Heat_Energy_Delivered_HP_kWh", data["Heat_Energy_Delivered_HP_kWh"].min(), data["Heat_Energy_Delivered_HP_kWh"].max())

    solara.Title("Heat Pump Hybrid System")
    solara.AppBarTitle("Energy and cost analysis")
    COP_mean = (data["Heat_Energy_Delivered_HP_kWh"].sum() + data["Heat_Energy_Not_Delivered_HP_kWh"].sum()) / (
        data["Electrical_Energy_Usage_HP_kWh"].sum() + data["Heat_Energy_Not_Delivered_HP_kWh"].sum()
    )
    with solara.Sidebar():
        Parameters()
        with solara.Card():
            # Percentage of energy delivered by heat pump
            percent_by_heat_pump = data["Heat_Energy_Delivered_HP_kWh"].sum() / data["Heat_Energy_Usage_kWh"].sum() * 100
            print("Percentage of energy delivered by heat pump:", percent_by_heat_pump)
            solara.Markdown(
                f"""
            
            ## Energy
            - **Total Heat Energy Usage**: {data["Heat_Energy_Usage_kWh"].sum():,.0f} kWh / {data["Heat_Energy_Usage_kWh"].sum() / caloric_value.value:,.0f} m³
            - **Heat Energy Delivered by Heat Pump**: {data["Heat_Energy_Delivered_HP_kWh"].sum():,.0f} kWh / {data["Heat_Energy_Delivered_HP_kWh"].sum() / caloric_value.value:,.0f} m³
            - **Heat Energy Not Delivered by Heat Pump**: {data["Heat_Energy_Not_Delivered_HP_kWh"].sum():,.0f} kWh / {data["Heat_Energy_Not_Delivered_HP_kWh"].sum() / caloric_value.value:,.0f} m³
            - **Percentage of Heat Energy Delivered by Heat Pump**: {percent_by_heat_pump:.1f}%
            - **Electrical Energy Usage by Heat Pump**: {data["Electrical_Energy_Usage_HP_kWh"].sum():,.0f} kWh / {data["Electrical_Energy_Usage_HP_kWh"].sum() / caloric_value.value:,.0f} m³
            - **Average COP**: {COP_mean}

            ## Cost

             - **Total Cost**: €{data["Total_Cost"].sum():,.2f}
             - **Total Electricity Cost**: €{data["Electricity_Cost"].sum():,.2f}
             - **Total Gas Cost**: €{data["Gas_Cost"].sum():,.2f}
             - **Cost when only using Gas**: €{data["Gas_Hypothetical_Cost"].sum():,.2f}
             - **Total savings**: €{data["Gas_Hypothetical_Cost"].sum() - data["Total_Cost"].sum():,.2f}
             - **Saving percentage**: {(data["Gas_Hypothetical_Cost"].sum() - data["Total_Cost"].sum()) / data["Gas_Hypothetical_Cost"].sum() * 100:.1f}%

"""
            )

    print(
        year_min.value,
        year_max.value,
        year.value,
    )
    plot_usage(
        data,
        caloric_value=caloric_value.value,
        gas_usage_per_degree_day=gas_usage_per_degree_day.value,
        indoor_temperature=indoor_temperature.value,
        use_minumum=use_minumum.value,
    )


def set_device(name):
    device = devices[name]
    device_name.value = name
    electric_power_kwh.value = device["Prated"]


@solara.component
def COP():
    dark_effective = solara.lab.use_dark_effective()
    plt.style.use("dark_background" if dark_effective else "default")
    device = devices[device_name.value]
    data = compute_all(
        year=year.value,
        electric_power_kwh=electric_power_kwh.value,
        gas_usage_per_degree_day=gas_usage_per_degree_day.value,
        indoor_temperature=indoor_temperature.value,
        caloric_value=caloric_value.value,
        device=device,
        T_water=water_temperature.value,
        price_per_kwh=price_per_kwh.value,
        price_per_cubic_meter=price_per_cubic_meter.value,
        daily_hours=daily_hours.value,
        use_minumum=use_minumum.value,
    )
    solara.Title("Heat Pump Hybrid System")
    solara.AppBarTitle("Energy and cost analysis")
    with solara.Sidebar():

        Parameters()

    # plot_stats_plotly(data, caloric_value=caloric_value.value, device=device, T_water=water_temperature.value)
    plot_stats(
        data,
        caloric_value=caloric_value.value,
        device=device,
        T_water=water_temperature.value,
        use_minumum=use_minumum.value,
        indoor_temperature=indoor_temperature.value,
        gas_usage_per_degree_day=gas_usage_per_degree_day.value,
        electric_power_kwh=electric_power_kwh.value,
    )

    solara.Markdown(
        """

These figures compare values used for the COP in the modelling to the values obtained from the spec sheet.
                    

"""
    )


@solara.component
def Layout(children):
    dark_effective = solara.lab.use_dark_effective()
    return solara.AppLayout(children=children, toolbar_dark=dark_effective, color=None)


routes = [
    solara.Route("/", component=Page, layout=Layout),
    solara.Route("COP", component=COP),
]