main_kpis.py

import pandas as pd
import numpy as np
from pathlib import Path
from typing import Optional, Callable

from sklearn.metrics import mean_squared_error, mean_absolute_error

from bems.config import EMSConfig
from cosimos.simulator.history import History, HistoryData


def main_results_kpis(history: History or HistoryData, results_dir: Optional[Path], config: EMSConfig, store: bool = True):
    if isinstance(history, HistoryData):
        data = history
    else:
        # don't do anything if simulation stopped immediately
        if len(history.time) == 1:
            return pd.DataFrame()

        data = history.to_dataframe()


    # left out: monetary_costs_gas_consumption, monetary_costs_peak, max_peak_kw_used,
    list_kpis = [
        monetary_costs,
        max_peak_kw,
        comfort_building_rmse,
        comfort_weighted_rmse,
        tracking_error_rmse,
        tracking_error_std,
        residual_error_mae,
        residual_error_mean
    ]

    df_kpis = calculate_kpis(data, config, list_kpis)

    if store and results_dir is not None:
        df_kpis.to_csv(results_dir.joinpath('main_kpis.csv'), index=False)

    return df_kpis


def calculate_kpis(history_data: HistoryData, config: EMSConfig, kpis: list[Callable]) -> Optional[pd.DataFrame]:
    dict_kpis = {}
    for kpi_function in kpis:
        dict_kpis[kpi_function.__name__] = kpi_function(history_data, config)

    df_kpis = pd.DataFrame(dict_kpis.items(), columns=['kpi_name', 'kpi_value'])

    return df_kpis


def monetary_costs_grid_gas(data: HistoryData, config: EMSConfig) -> float:
    dt = (data.state.index[1] - data.state.index[0]).total_seconds() / 60 / 60
    E_gas_chp = (1 + 1 / config.c_chp) / config.eta_chp * sum(data.applied_controller_output['P_chp']) * dt
    E_gas_rad = dt * sum(data.applied_controller_output['Q_rad']) / config.eps_rad
    E_gas = E_gas_rad + E_gas_chp

    E_grid_in = sum(np.maximum(data.applied_controller_output['P_grid'], 0)) * dt
    E_grid_out = -1 * sum(np.minimum(data.applied_controller_output['P_grid'], 0)) * dt

    costs_grid = E_grid_in * config.c_grid_buy - E_grid_out * config.c_grid_sell
    costs_gas = E_gas * config.c_gas

    return costs_grid + costs_gas


def monetary_costs(data: HistoryData, config: EMSConfig) -> float:
    costs_grid_gas = monetary_costs_grid_gas(data, config)

    P_peak = data.applied_controller_output['P_grid'].max()
    costs_peak = P_peak * config.c_grid_peak

    return costs_grid_gas + costs_peak


def max_peak_kw(data: HistoryData, config: EMSConfig) -> float:
    return data.state['P_grid_peak_current'].values[-1]


def comfort_weighted_rmse(data: HistoryData, config: EMSConfig) -> float:
    """
    Weighted sum of the RMSE of individual zones from the reference temperature

    :param data:
    :param config:
    :return:
    """
    errors = []
    for i in range(1, 8):
        errors.append(mean_squared_error(
            data.state[f'theta_{i}'][1:],
            np.ones(len(data.state.index) - 1) * config.theta_ref,
            squared=False
        ))

    rmse = config.building_zone_weights @ errors
    return rmse

def comfort_weighted_mse(data: HistoryData, config: EMSConfig) -> float:
    """
    Weighted sum of the MSE of individual zones from the reference temperature,
    equivalent to J_{comf,dis} objective of monolithic or naive distributor

    :param data:
    :param config:
    :return:
    """
    errors = []
    for i in range(1, 8):
        errors.append(mean_squared_error(
            data.state[f'theta_{i}'][1:],
            np.ones(len(data.state.index) - 1) * config.theta_ref,
            squared=True
        ))

    rmse = config.building_zone_weights @ errors
    return rmse

def comfort_building_rmse(data: HistoryData, config: EMSConfig) -> float:
    """
    RMSE of the average building temperature from the reference temperature

    :param data:
    :param config:
    :return:
    """

    theta_b = config.building_zone_weights @ data.state[[f'theta_{i}' for i in range(1, 8)]].T
    rmse = mean_squared_error(
        theta_b,
        np.ones(len(theta_b)) * config.theta_ref,
        squared=False
    )

    return rmse


def comfort_weighted_mae(data: HistoryData, config: EMSConfig) -> float:
    """
    Weighted sum of the MAE of individual zones from the reference temperature

    :param data:
    :param config:
    :return:
    """
    errors = []
    for i in range(1, 8):
        errors.append(mean_absolute_error(
            data.state[f'theta_{i}'][1:],
            np.ones(len(data.state.index) - 1) * config.theta_ref
        ))

    mae = config.building_zone_weights @ errors
    return mae


def tracking_error_std(data: HistoryData, config: EMSConfig) -> float:
    """ calculate "standard deviation" of individual zone temperatures from average of all zones """
    errors = []

    theta_b = config.building_zone_weights @ data.state[[f'theta_{i}' for i in range(1, 8)]].T
    for i in range(1, 8):
        error = mean_squared_error(
            data.state[f'theta_{i}'],
            theta_b,
            squared=False
        )
        errors.append(error)

    rmse = config.building_zone_weights @ errors

    return rmse


def tracking_error_rmse(data: HistoryData, config: EMSConfig) -> float:
    """
    calculate total (weighted) rmse of tracking error, i.e. deviation of individual zone temperatures from aggregator reference
    returns nan if no aggregator was used
    """

    rmse_zones = tracking_error_zones_rmse(data, config)
    if not isinstance(rmse_zones, list) and np.isnan(rmse_zones):
        return rmse_zones

    rmse = config.building_zone_weights @ rmse_zones

    return rmse


def tracking_error_zones_rmse(data: HistoryData, config: EMSConfig) -> list[float] | float:
    """
    calculate rmse of tracking error, i.e. deviation of individual zone temperatures from aggregator reference
    returns nan if no aggregator was used
    """
    errors = []

    if "theta_b_ref" not in data.logged_data.columns:
        return np.nan

    for i in range(1, 8):
        error = mean_squared_error(
            data.state[f'theta_{i}'][1:],
            data.logged_data["theta_b_ref"],
            squared=False
        )
        errors.append(error)

    return errors


def tracking_error_zones_rmse_config_ref(data: HistoryData, config: EMSConfig) -> list[float] | float:
    """
    calculate rmse of tracking error, i.e. deviation of individual zone temperatures from config.theta_ref
    """
    errors = []

    for i in range(1, 8):
        error = mean_squared_error(
            data.state[f'theta_{i}'][1:],
            np.ones(len(data.state.index)-1) * config.theta_ref,
            squared=False
        )
        errors.append(error)

    return errors


def tracking_error_zones_std(data: HistoryData, config: EMSConfig) -> list[float] | float:
    """
    calculate rmse of tracking error, i.e. deviation of individual zone temperatures from aggregator reference
    returns nan if no aggregator was used
    """
    errors = []

    if "theta_b" not in data.state.columns:
        return np.nan

    for i in range(1, 8):
        error = mean_squared_error(
            data.state[f'theta_{i}'][1:],
            data.state["theta_b"][1:],
            squared=False
        )
        errors.append(error)

    return errors


def residual_error_mae(data: HistoryData, config: EMSConfig) -> list[float]:
    """
    Calculates the thermal residual model error for each zone.

    :param data:
    :param config:
    :return:
    """

    theta = [f"theta_{i}" for i in range(1, 10)]

    mae_zones = []
    for zone in theta:
        mae_zones.append(
            mean_absolute_error(
                data.state[zone].shift(-1)[:-1],
                data.expected_state[zone]
            )
        )

    return mae_zones

def residual_error_std(data: HistoryData, config: EMSConfig) -> list[float]:
    """
    Calculates the thermal residual model error for each zone.

    :param data:
    :param config:
    :return:
    """

    theta = [f"theta_{i}" for i in range(1, 10)]

    mae_zones = []
    for zone in theta:
        mae_zones.append(
            np.std(np.abs(
                data.state[zone].shift(-1)[:-1] - data.expected_state[zone]
            ))
        )

    return mae_zones


def residual_error_mean(data: HistoryData, config: EMSConfig) -> list[float]:
    """
    Calculates the thermal residual model error for each zone.

    :param data:
    :param config:
    :return:
    """

    theta = [f"theta_{i}" for i in range(1, 10)]

    mae_zones = []
    for zone in theta:
        mae_zones.append(
            (data.state[zone].shift(-1)[:-1] - data.expected_state[zone]).mean()
        )

    return mae_zones