least_squares_fit_exact.py

import casadi
import numpy as np
from casadi import Opti, vertcat, sum1
import datetime
import pytz

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from pathlib import Path
from tqdm import tqdm

from bems.helpers.data import load_simulation_results
from bems.config import EMSConfig
from bems.models.full import Model as BemsModel
from bems.helpers.data import prepare_fitting_data

from model_fitting.optimal_fit import solve_exact_constrained_least_squares_fit_model
from model_fitting.interfaces import Topology
from model_fitting.utility import build_state_space_model

from sklearn.model_selection import train_test_split

topology = Topology(
    rooms=list(range(0, 9)),
    outside=list(range(0, 9)),
    couplings=[(1, 8), (2, 3), (4, 5), (4, 7), (5, 7)]
)
# results_path = Path(r"C:\Users\re902174\workspaces\Hierarchical MPC\robustification\bems-cosimos\analysis\journal_publication_study\final_results\year_monolithic_dt_240412_163803")
# results_path = Path(r"C:\Users\re902174\workspaces\Hierarchical MPC\robustification\bems-cosimos\analysis\journal_publication_study\final_results\year_hierarchical_no_neg_dt_240117_095157")
# results_path = Path(r"C:\Users\re902174\workspaces\Hierarchical MPC\robustification\bems-cosimos\simulations\digital_twin\results\year_simple_rbc_dt_240830_130239")
# results_path = Path(r"C:\Users\re902174\workspaces\Hierarchical MPC\robustification\bems-cosimos\simulations\digital_twin\results\year_simple_rbc_dt_fixed_240909_111607")
results_path = Path(r"C:\Users\re902174\workspaces\Hierarchical MPC\robustification\bems-cosimos\simulations\model_fitting\results\year_2021_simple_rbc_dt_fixed_windows_240925_100119")
# results_path = Path(r"C:\Users\re902174\workspaces\Hierarchical MPC\robustification\bems-cosimos\simulations\model_fitting\results\year_2021_simple_rbc_dt_fixed_no_wall_inertia_240924_152709")
# results_path = Path(r"C:\Users\re902174\workspaces\Hierarchical MPC\robustification\bems-cosimos\simulations\digital_twin\results\year_rbc_dt_240411_164508")
data = load_simulation_results(results_path)

dump_directory = Path(r"C:\Users\re902174\workspaces\Hierarchical MPC\robustification\bems-cosimos\data\fitted_models")

features, labels = prepare_fitting_data(
    data,
    use_controller_predictions=False
)

experiment_name = "simple_rbc_fixed_windows_TEST"
# experiment_name = "p_rbc"

masks = [
    {"summer": False, "internal_heating": False},
    {"summer": False, "internal_heating": True},
    # {"summer": True, "internal_heating": False},
    # {"summer": True, "internal_heating": True},
]

split_train_test = True

dump = False

config = EMSConfig()
bems_model = BemsModel(config)

dump_directory = dump_directory.joinpath(experiment_name)

#%%
test_size = 0.25
train_test_random_state = 2342
# stratify_on_columns = ["time_of_day", "day_of_week"]
stratify_on_columns = None
shuffle_split = True

temperature_cols = [f"theta_{i}" for i in range(1, 10)]
power_cols = [f"Q_total_{i}" for i in range(1, 10)]
label_cols = np.array(temperature_cols)
ambient_col = "theta_air"
load_demand_col = "P_dem"

for mask in masks:
    if "summer" in mask:
        filter_only_summer = mask["summer"]
    else:
        filter_only_summer = False

    if "internal_heating" in mask:
        fit_internal_heating = mask["internal_heating"]
    else:
        fit_internal_heating = False

    if "weekends" in mask:
        filter_only_weekends = mask["weekends"]
    else:
        filter_only_weekends = False

    if "winter" in mask:
        filter_only_winter = mask["winter"]
    else:
        filter_only_winter = False

    if "Q_other" in mask:
        fit_Q_other = mask["Q_other"]
    else:
        fit_Q_other = True

    # try fitting only weekend days
    index_filter = np.repeat(True, len(features.index))
    if filter_only_weekends:
        index_filter &= features.index.dayofweek > 4
    if filter_only_summer:
        tz = pytz.timezone("Europe/Berlin")
        index_filter &= features["datetime"].between(
            datetime.datetime(2021, 6, 1, tzinfo=tz),
            datetime.datetime(2021, 8, 1, tzinfo=tz),
        )
    if filter_only_winter:
        tz = pytz.timezone("Europe/Berlin")
        index_filter &= features["datetime"].between(
            datetime.datetime(2021, 1, 1, tzinfo=tz),
            datetime.datetime(2021, 4, 1, tzinfo=tz),
        )

    if index_filter is not None:
        features = features[index_filter]
        labels = labels[index_filter]

    if split_train_test:
        X_train, X_test, y_train, y_test = train_test_split(
            features, labels, test_size=test_size, random_state=train_test_random_state, shuffle=shuffle_split, stratify=stratify_on_columns
        )
    else:
        X_train, X_test, y_train, y_test = features, features, labels, labels

    # weights = np.array([0.1, 0.25, 0.10, 0.05, 0.2, 0.2, 0.05, 0.001, 0.002])
    weights = None
    fitted_parameters = solve_exact_constrained_least_squares_fit_model(
        topology=topology,
        data=(X_train, y_train),
        temperature_columns=temperature_cols,
        power_columns=power_cols,
        ambient_temperature_column=ambient_col,
        load_demand_column=load_demand_col,
        fit_constant_offset=fit_Q_other,
        fit_internal_heating=fit_internal_heating,
        approximate_capacity_proportions=weights,
        allow_slack=False
    )

    fitted_state_space_model = build_state_space_model(topology, fitted_parameters)

    if dump:
        dump_directory.mkdir(parents=True, exist_ok=True)

        import pickle
        parametrization_str = f"{experiment_name}_"
        parametrization_str += f"{'other_' if fit_Q_other else ''}"
        parametrization_str += f"{'internal_' if fit_internal_heating else ''}"
        parametrization_str += f"{'weekends_' if filter_only_weekends else ''}"
        parametrization_str += f"{'summer_' if filter_only_summer else ''}"

        # file_name = f"fitted_state_space_{parametrization_str}_2021.pkl"
        # with open(dump_directory.joinpath(file_name), "wb") as fh:
        #     pickle.dump(fitted_state_space_model, fh)

        file_name = f"fitted_parameters_{parametrization_str}_2021.pkl"
        with open(dump_directory.joinpath(file_name), "wb") as fh:
            pickle.dump(fitted_parameters, fh)

    #%%
    C_th_fit = fitted_parameters.C_thermal(config.dt)
    C_th_ex = config.C_thermal

    C_comp = np.c_[C_th_ex, C_th_fit]

    #%%
    c_val = fitted_parameters.c
    beta_diag_fit = c_val
    beta_diag_ex = np.array([config.beta_coupling[x, y] for x, y in topology.couplings])

    beta_comp = np.c_[beta_diag_ex, beta_diag_fit]

    #%%
    coupling_ext_fit = fitted_parameters.b
    coupling_ext_ex = bems_model.S_dis[:, 0]

    coupling_ext_comp = np.c_[coupling_ext_ex, coupling_ext_fit]