From 633c4d12d35e9019be09b8f9aa9265cd6bcb2524 Mon Sep 17 00:00:00 2001 From: "Radu C. Martin" Date: Sat, 31 Jul 2021 12:57:47 +0200 Subject: [PATCH] Updated Notebook comments --- .../38_gp_hyperparameter_estimation.ipynb | 83 +- .../39_svgp_hyperparameter_estimation.ipynb | 65 +- Notebooks/41_casadi_gp_test.ipynb | 42 + Notebooks/50_mpc_formulation.ipynb | 2375 ----------------- Notebooks/70_Server_result_analysis.ipynb | 64 + 5 files changed, 216 insertions(+), 2413 deletions(-) delete mode 100644 Notebooks/50_mpc_formulation.ipynb diff --git a/Notebooks/38_gp_hyperparameter_estimation.ipynb b/Notebooks/38_gp_hyperparameter_estimation.ipynb index b88418f..c3f23c4 100644 --- a/Notebooks/38_gp_hyperparameter_estimation.ipynb +++ b/Notebooks/38_gp_hyperparameter_estimation.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Bayesian Optimisation of starting Gaussian Process hyperparameters" + "# Gaussian Process Model Training and Performance Evaluation" ] }, { @@ -213,7 +213,7 @@ "id": "0aba0df5-b0e3-4738-bb61-1dad869d1ea3" }, "source": [ - "## Load previously exported data" + "## Load previously exported CARNOT 'experimental' data" ] }, { @@ -226,6 +226,13 @@ "dfs_test = []" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Separate into training and testing data sets:" + ] + }, { "cell_type": "code", "execution_count": 16, @@ -249,6 +256,13 @@ " dfs_test.append(pd.read_csv(f\"../Data/Good_CARNOT/{exp}_table.csv\").rename(columns = {'Power': 'SimulatedHeat'}))" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Separate columns into exogenous inputs, controlled inputs and outputs:" + ] + }, { "cell_type": "code", "execution_count": 18, @@ -262,6 +276,13 @@ "y_cols = ['SimulatedTemp']" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Impose the autoregressive lags for each input group:" + ] + }, { "cell_type": "code", "execution_count": 19, @@ -511,6 +532,13 @@ " return df_gpr" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Merge all the training dataframes:" + ] + }, { "cell_type": "code", "execution_count": 28, @@ -667,13 +695,19 @@ "df_gpr_train.head()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Select all points in the training dataset:" + ] + }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ - "train_dataset_size = 15 * 96\n", "train_dataset_size = -1" ] }, @@ -1401,6 +1435,13 @@ "y_range = np.arange(1,6)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Iterate over all combination of lags and compute for each the RMSE, SMSE, LPD and MSLL errors:" + ] + }, { "cell_type": "code", "execution_count": null, @@ -1545,6 +1586,13 @@ "## Multistep prediction" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Select the dataset which will be used for multistep prediction:" + ] + }, { "cell_type": "code", "execution_count": 54, @@ -1556,6 +1604,13 @@ "df_output = dfs_gpr_test[test_dataset_idx][dict_cols['y'][1]]" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Select the starting index in the test dataset and the number of consecutive points to simulate:" + ] + }, { "cell_type": "code", "execution_count": 55, @@ -1617,26 +1672,6 @@ "plt.title(f\"Multi step prediction over {N_pred} steps for Test dataset {test_dataset_idx}\")\n", "plt.savefig(f\"../Thesis/Plots/GP_{w_lags}{u_lags}{y_lags}_{train_dataset_size}pts_test_prediction_{N_pred}_steps.pdf\")" ] - }, - { - "cell_type": "code", - "execution_count": 141, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "TensorShape([2612, 7])" - ] - }, - "execution_count": 141, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "m.data[0].shape" - ] } ], "metadata": { diff --git a/Notebooks/39_svgp_hyperparameter_estimation.ipynb b/Notebooks/39_svgp_hyperparameter_estimation.ipynb index a546e46..7b7112d 100644 --- a/Notebooks/39_svgp_hyperparameter_estimation.ipynb +++ b/Notebooks/39_svgp_hyperparameter_estimation.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Bayesian Optimisation of starting Gaussian Process hyperparameters" + "# Sparse and Variational Gaussian Process Model Training and Performance Evaluation" ] }, { @@ -78,9 +78,6 @@ "cell_type": "markdown", "metadata": { "id": "90fdac33-eed4-4ab4-b2b1-de0f1f27727b", - "jupyter": { - "source_hidden": true - }, "tags": [] }, "source": [ @@ -199,7 +196,7 @@ "id": "0aba0df5-b0e3-4738-bb61-1dad869d1ea3" }, "source": [ - "## Load previously exported data" + "## Load previously exported CARNOT 'experimental' data" ] }, { @@ -212,6 +209,13 @@ "dfs_test = []" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Separate into training and testing data sets:" + ] + }, { "cell_type": "code", "execution_count": 13, @@ -235,6 +239,13 @@ " dfs_test.append(pd.read_csv(f\"../Data/Good_CARNOT/{exp}_table.csv\").rename(columns = {'Power': 'SimulatedHeat'}))" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Separate columns into exogenous inputs, controlled inputs and outputs:" + ] + }, { "cell_type": "code", "execution_count": 15, @@ -248,6 +259,13 @@ "y_cols = ['SimulatedTemp']" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Impose the autoregressive lags for each input group:" + ] + }, { "cell_type": "code", "execution_count": 16, @@ -497,6 +515,13 @@ " return df_gpr" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Merge all the training dataframes:" + ] + }, { "cell_type": "code", "execution_count": 25, @@ -647,15 +672,6 @@ "df_gpr_train.head()" ] }, - { - "cell_type": "code", - "execution_count": 26, - "metadata": {}, - "outputs": [], - "source": [ - "#df_gpr_train = df_gpr_train.sample(n = 500)" - ] - }, { "cell_type": "code", "execution_count": 27, @@ -1245,6 +1261,13 @@ "y_lags = 5" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Iterate over all combination of lags and compute for each the RMSE, SMSE, LPD and MSLL errors:" + ] + }, { "cell_type": "code", "execution_count": 54, @@ -1470,6 +1493,13 @@ "## Multistep prediction" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Select the dataset which will be used for multistep prediction:" + ] + }, { "cell_type": "code", "execution_count": 47, @@ -1481,6 +1511,13 @@ "df_output = dfs_gpr_test[test_dataset_idx][dict_cols['y'][1]]" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Select the starting index in the test dataset and the number of consecutive points to simulate:" + ] + }, { "cell_type": "code", "execution_count": 48, diff --git a/Notebooks/41_casadi_gp_test.ipynb b/Notebooks/41_casadi_gp_test.ipynb index ccb71de..99da151 100644 --- a/Notebooks/41_casadi_gp_test.ipynb +++ b/Notebooks/41_casadi_gp_test.ipynb @@ -75,6 +75,13 @@ "## GP model" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the underlying function that will be modeled:" + ] + }, { "cell_type": "code", "execution_count": 8, @@ -108,6 +115,13 @@ "})" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Sample `n` points and add measurement noise:" + ] + }, { "cell_type": "code", "execution_count": 11, @@ -142,6 +156,13 @@ "Y_sampled = Y_sampled + noise" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the function and the sampled points:" + ] + }, { "cell_type": "code", "execution_count": 14, @@ -167,6 +188,13 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Train the GP" + ] + }, { "cell_type": "code", "execution_count": 15, @@ -362,6 +390,13 @@ "## CasADi part" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the CasADi callback:" + ] + }, { "cell_type": "code", "execution_count": 23, @@ -1074,6 +1109,13 @@ "grads" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define another casadi callback, which directly uses the GPflow gradients:" + ] + }, { "cell_type": "code", "execution_count": 29, diff --git a/Notebooks/50_mpc_formulation.ipynb b/Notebooks/50_mpc_formulation.ipynb deleted file mode 100644 index 03f0077..0000000 --- a/Notebooks/50_mpc_formulation.ipynb +++ /dev/null @@ -1,2375 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "from pathlib import Path\n", - "from shutil import copyfile\n", - "import pickle" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Load the general math/data manipulation packages" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import pandas as pd" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Load the packages related to the Gaussian Process Regressor:" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "import gpflow\n", - "import tensorflow as tf" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "from gpflow.utilities import print_summary" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "gpflow.config.set_default_summary_fmt(\"notebook\")" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "tf.config.set_visible_devices([], 'GPU')" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Load the CasADi package used for optimization:" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "import casadi" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Import MATLAB engine and start it in the background since this takes a while:" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [], - "source": [ - "import matlab.engine" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "eng = matlab.engine.start_matlab()" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "eng.load_system(\"../Simulink/polydome\", background = True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Copy the experimental data set to the CARNOT input location: " - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "exp_id = 'Exp5'" - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "'../Data/input_WDB.mat'" - ] - }, - "execution_count": 13, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "copyfile(f\"../Data/Experimental_data_WDB/{exp_id}_WDB.mat\", \"../Data/input_WDB.mat\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Load the existing GP model" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "model_path = Path(Path.cwd(), 'model')" - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "metadata": {}, - "outputs": [], - "source": [ - "def load_gpr(model_path):\n", - " x_scaler = pickle.load(open(Path(model_path, 'x_scaler.pkl'), 'rb'))\n", - " m_params = pickle.load(open(Path(model_path, 'gp_params.gpf'), 'rb'))\n", - " m_data = pickle.load(open(Path(model_path, 'gp_data.gpf'), 'rb'))\n", - "\n", - " k = gpflow.kernels.SquaredExponential(lengthscales=([1] * m_data[0].shape[1])) + gpflow.kernels.Constant()\n", - "\n", - " m = gpflow.models.GPR(\n", - " data = m_data, \n", - " kernel = k, \n", - " mean_function = None\n", - " )\n", - " \n", - " gpflow.utilities.multiple_assign(m, m_params)\n", - " \n", - " return x_scaler, m" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "
name class transform prior trainable shape dtype value
GPR.kernel.kernels[0].variance ParameterSoftplus True () float64284790.175758825
GPR.kernel.kernels[0].lengthscalesParameterSoftplus True (7,) float64[8.23155895e+04, 1.32876975e+05, 2.21167354e+02...
GPR.kernel.kernels[1].variance ParameterSoftplus True () float64220815.62255492577
GPR.likelihood.variance ParameterSoftplus + Shift True () float640.1263230455247145
" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "x_scaler, model = load_gpr(model_path)\n", - "print_summary(model)" - ] - }, - { - "cell_type": "code", - "execution_count": 85, - "metadata": {}, - "outputs": [], - "source": [ - "def get_gpr_horizon_array(W, x0, u):\n", - " \n", - " n_rows = W.shape[0]\n", - " \n", - " n_w = W.shape[1]\n", - " n_inputs = u.shape[1]\n", - " n_cols = n_w + n_inputs + x0.shape[1]\n", - " \n", - " X = np.zeros([n_rows, n_cols])\n", - " X[:, :2] = W\n", - " X[0, (n_w + n_inputs):] = x0\n", - " X[:, 2] = u.reshape((-1,))\n", - " \n", - " # autoregressive inputs\n", - " X[1: , 3] = X[:-1, 2]\n", - " \n", - " for idx in range(1, n_rows):\n", - " x_sc = x_scaler.transform(X[idx-1, :].reshape((1, -1)))\n", - " X[idx, 4], _ = model.predict_y(x_sc)\n", - " X[idx, 5:] = X[idx -1, 4:-1]\n", - " \n", - " return X" - ] - }, - { - "cell_type": "code", - "execution_count": 86, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[ 1. , 1. , 3. , 2. , 2. ,\n", - " 2. , 2. ],\n", - " [ 1. , 1. , 3. , 3. , 3.16887448,\n", - " 2. , 2. ],\n", - " [ 1. , 1. , 3. , 3. , 6.02687796,\n", - " 3.16887448, 2. ],\n", - " [ 1. , 1. , 3. , 3. , 11.22601695,\n", - " 6.02687796, 3.16887448],\n", - " [ 1. , 1. , 3. , 3. , 19.6304083 ,\n", - " 11.22601695, 6.02687796]])" - ] - }, - "execution_count": 86, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "test_W = np.ones([5, 2])\n", - "test_x0 = 2 * np.ones([1, 4])\n", - "test_u0 = 3 * np.ones([5, 1])\n", - "\n", - "test_X = get_gpr_horizon_array(test_W, test_x0, test_u0)\n", - "test_X" - ] - }, - { - "cell_type": "code", - "execution_count": 87, - "metadata": {}, - "outputs": [], - "source": [ - "n_states = model.data[0].shape[1]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Set up the CasADi optimization problem" - ] - }, - { - "cell_type": "code", - "execution_count": 88, - "metadata": {}, - "outputs": [], - "source": [ - "# Package the regression model in a CasADi callback\n", - "class GPR(casadi.Callback):\n", - " def __init__(self, name, opts={}):\n", - " casadi.Callback.__init__(self)\n", - " self.construct(name, opts)\n", - " \n", - " # Number of inputs and outputs\n", - " def get_n_in(self): return 1\n", - " def get_n_out(self): return 1\n", - "\n", - " def get_sparsity_in(self,i):\n", - " return casadi.Sparsity.dense(n_states,1)\n", - "\n", - " def eval(self, arg):\n", - " x_input = x_scaler.transform(np.array(arg[0]).reshape(1, -1))\n", - " #x_input = np.array(arg[0]).reshape(1, -1)\n", - " [mean, _] = model.predict_y(x_input)\n", - " return [mean.numpy()]" - ] - }, - { - "cell_type": "code", - "execution_count": 89, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "GPR:(i0[7])->(o0) CallbackInternal\n" - ] - } - ], - "source": [ - "# Instantiate the Callback (make sure to keep a reference to it!)\n", - "gpr = GPR('GPR', {\"enable_fd\":True})\n", - "print(gpr)" - ] - }, - { - "cell_type": "code", - "execution_count": 90, - "metadata": {}, - "outputs": [], - "source": [ - "T_set = 20\n", - "N_horizon = 5" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Define optimization variables" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\"Gaussian" - ] - }, - { - "cell_type": "code", - "execution_count": 91, - "metadata": {}, - "outputs": [], - "source": [ - "X = casadi.MX.sym(\"X\", N_horizon, n_states)\n", - "W = casadi.MX.sym(\"W\", N_horizon, 2)\n", - "x0_lags = casadi.MX.sym(\"lags\", 1, n_states - 3)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Impose initial lags:" - ] - }, - { - "cell_type": "code", - "execution_count": 92, - "metadata": {}, - "outputs": [], - "source": [ - "g = casadi.vec(X[0,3:] - x0_lags)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Impose disturbances:" - ] - }, - { - "cell_type": "code", - "execution_count": 93, - "metadata": {}, - "outputs": [], - "source": [ - "g = casadi.vertcat(\n", - " g,\n", - " casadi.vec(X[:, :2] - W)\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Compute objective:" - ] - }, - { - "cell_type": "code", - "execution_count": 94, - "metadata": {}, - "outputs": [], - "source": [ - "Y = gpr(X.T)\n", - "J = casadi.dot(Y - T_set, Y - T_set)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Fix input/output lags between time steps (equality constraints):" - ] - }, - { - "cell_type": "code", - "execution_count": 95, - "metadata": {}, - "outputs": [], - "source": [ - "for idx in range(1, N_horizon):\n", - " g = casadi.vertcat(\n", - " g,\n", - " X[idx, 3] - X[idx-1, 2],\n", - " X[idx, 4] - gpr(X[idx-1,:]),\n", - " X[idx, 5] - X[idx-1, 4],\n", - " X[idx, 6] - X[idx-1, 5]\n", - " )" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Impose input inequality constraints:" - ] - }, - { - "cell_type": "code", - "execution_count": 96, - "metadata": {}, - "outputs": [], - "source": [ - "g = casadi.vertcat(\n", - " g,\n", - " X[:, 2]\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Compile the optimization problem" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Compile the parameter vector:" - ] - }, - { - "cell_type": "code", - "execution_count": 97, - "metadata": {}, - "outputs": [], - "source": [ - "p = casadi.vertcat(\n", - " casadi.vec(W),\n", - " casadi.vec(x0_lags)\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 98, - "metadata": {}, - "outputs": [], - "source": [ - "prob = {\"x\": casadi.vec(X), \"f\": J, \"p\": p, \"g\": g}\n", - "options = {\"ipopt\": {\"hessian_approximation\": \"limited-memory\", \"max_iter\": 10,\n", - " #\"acceptable_tol\": 1e-6, \"tol\": 1e-6,\n", - " #\"linear_solver\": \"SPRAL\",\n", - " #\"acceptable_obj_change_tol\": 1e-5, \n", - " #\"mu_strategy\": \"adaptive\",\n", - " \"expect_infeasible_problem\": \"yes\"\n", - " }}" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Get the solver object:" - ] - }, - { - "cell_type": "code", - "execution_count": 99, - "metadata": {}, - "outputs": [], - "source": [ - "solver = casadi.nlpsol(\"solver\",\"ipopt\",prob, options)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Compute the `lbg` `ubg` since they are always the same:" - ] - }, - { - "cell_type": "code", - "execution_count": 100, - "metadata": {}, - "outputs": [], - "source": [ - "Pel_max = 6300\n", - "COP_heat = 5\n", - "COP_cool = 5\n", - "u_min = - COP_cool * Pel_max\n", - "u_max = COP_heat * Pel_max" - ] - }, - { - "cell_type": "code", - "execution_count": 101, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "#u_min_sc = x_scaler.transform(np.array([1, 1, u_min, u_min, 1, 1, 1]).reshape((1, -1)))[0, 2]\n", - "#u_max_sc = x_scaler.transform(np.array([1, 1, u_max, u_max, 1, 1, 1]).reshape((1, -1)))[0, 2]" - ] - }, - { - "cell_type": "code", - "execution_count": 102, - "metadata": {}, - "outputs": [], - "source": [ - "real_lbg = [0] * (4 + 2 * N_horizon + 4 * (N_horizon - 1)) + [u_min] * (N_horizon)\n", - "real_ubg = [0] * (4 + 2 * N_horizon + 4 * (N_horizon - 1)) + [u_max] * (N_horizon)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Load CARNOT building with MATLAB backend" - ] - }, - { - "cell_type": "code", - "execution_count": 103, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
PowerSetpointOutsideTempSupplyTempInsideTempSolRadHeatSimulatedTemp
timestamp
2017-06-30 20:00:00+02:008.62069022.519.022.122.283333178.6182678.62069022.283285
2017-06-30 20:05:00+02:008.62069022.519.022.122.233333135.1221008.62069022.266734
2017-06-30 20:10:00+02:007.56666722.519.022.122.366667126.2199677.56666722.258667
2017-06-30 20:15:00+02:008.17241422.519.022.122.366667123.7724678.17241422.255146
2017-06-30 20:20:00+02:008.37931022.517.522.422.566667106.163600-25.13793122.255204
\n", - "
" - ], - "text/plain": [ - " Power Setpoint OutsideTemp SupplyTemp \\\n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 8.620690 22.5 19.0 22.1 \n", - "2017-06-30 20:05:00+02:00 8.620690 22.5 19.0 22.1 \n", - "2017-06-30 20:10:00+02:00 7.566667 22.5 19.0 22.1 \n", - "2017-06-30 20:15:00+02:00 8.172414 22.5 19.0 22.1 \n", - "2017-06-30 20:20:00+02:00 8.379310 22.5 17.5 22.4 \n", - "\n", - " InsideTemp SolRad Heat SimulatedTemp \n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 22.283333 178.618267 8.620690 22.283285 \n", - "2017-06-30 20:05:00+02:00 22.233333 135.122100 8.620690 22.266734 \n", - "2017-06-30 20:10:00+02:00 22.366667 126.219967 7.566667 22.258667 \n", - "2017-06-30 20:15:00+02:00 22.366667 123.772467 8.172414 22.255146 \n", - "2017-06-30 20:20:00+02:00 22.566667 106.163600 -25.137931 22.255204 " - ] - }, - "execution_count": 103, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df = pd.read_pickle(f\"../Data/CARNOT_output/{exp_id}_full.pkl\")\n", - "df.head()" - ] - }, - { - "cell_type": "code", - "execution_count": 104, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
timetimestampzenithazimuthdnidhiOutsideTempTsky_radrelative_humidityprecipitationcloud_indexpressurewind_speedwind_directionaoiincidence_mainincidence_secondpoa_directpoa_diffuse
0020170630200077.058396290.138593461.58739475.75999919.013.050-99990.5963000-999977.058396-9999-9999103.37612275.759999
130020170630200577.859576290.978550206.01476192.04208719.013.050-99990.5963000-999977.859576-9999-999943.32662892.042087
260020170630201078.655742291.820500205.63205386.03480619.013.050-99990.5963000-999978.655742-9999-999940.44855686.034806
390020170630201579.446620292.664742254.54169277.50275619.013.050-99990.5963000-999979.446620-9999-999946.61969577.502756
4120020170630202080.231464293.511567186.92863174.72604117.511.550-99990.5963000-999980.231464-9999-999931.71587174.726041
\n", - "
" - ], - "text/plain": [ - " time timestamp zenith azimuth dni dhi \\\n", - "0 0 201706302000 77.058396 290.138593 461.587394 75.759999 \n", - "1 300 201706302005 77.859576 290.978550 206.014761 92.042087 \n", - "2 600 201706302010 78.655742 291.820500 205.632053 86.034806 \n", - "3 900 201706302015 79.446620 292.664742 254.541692 77.502756 \n", - "4 1200 201706302020 80.231464 293.511567 186.928631 74.726041 \n", - "\n", - " OutsideTemp Tsky_rad relative_humidity precipitation cloud_index \\\n", - "0 19.0 13.0 50 -9999 0.5 \n", - "1 19.0 13.0 50 -9999 0.5 \n", - "2 19.0 13.0 50 -9999 0.5 \n", - "3 19.0 13.0 50 -9999 0.5 \n", - "4 17.5 11.5 50 -9999 0.5 \n", - "\n", - " pressure wind_speed wind_direction aoi incidence_main \\\n", - "0 96300 0 -9999 77.058396 -9999 \n", - "1 96300 0 -9999 77.859576 -9999 \n", - "2 96300 0 -9999 78.655742 -9999 \n", - "3 96300 0 -9999 79.446620 -9999 \n", - "4 96300 0 -9999 80.231464 -9999 \n", - "\n", - " incidence_second poa_direct poa_diffuse \n", - "0 -9999 103.376122 75.759999 \n", - "1 -9999 43.326628 92.042087 \n", - "2 -9999 40.448556 86.034806 \n", - "3 -9999 46.619695 77.502756 \n", - "4 -9999 31.715871 74.726041 " - ] - }, - "execution_count": 104, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df_wdb = pd.read_pickle(f\"../Data/Experimental_python/{exp_id}_WDB.pkl\")\n", - "df_wdb.head()" - ] - }, - { - "cell_type": "code", - "execution_count": 105, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Experiment runtime: 208500\n" - ] - } - ], - "source": [ - "runtime = df_wdb['time'].iloc[-1]\n", - "print(f\"Experiment runtime: {runtime}\")" - ] - }, - { - "cell_type": "code", - "execution_count": 106, - "metadata": {}, - "outputs": [], - "source": [ - "eng.workspace['t0'] = float(df['InsideTemp'][0])" - ] - }, - { - "cell_type": "code", - "execution_count": 107, - "metadata": {}, - "outputs": [], - "source": [ - "day_air_exchange_rate = 2.75\n", - "night_air_exchange_rate = 2.75" - ] - }, - { - "cell_type": "code", - "execution_count": 108, - "metadata": {}, - "outputs": [], - "source": [ - "air_exchange_rate = np.zeros((df_wdb.shape[0], 2))\n", - "air_exchange_rate[:, 0] = df_wdb['time']\n", - "air_exchange_rate[:, 1] = np.where(df['Power'] < 100, day_air_exchange_rate, night_air_exchange_rate)\n", - "eng.workspace['air_exchange_rate'] = matlab.double(air_exchange_rate.tolist())" - ] - }, - { - "cell_type": "code", - "execution_count": 109, - "metadata": {}, - "outputs": [], - "source": [ - "power = np.array([df_wdb['time'], df['Heat']]).T\n", - "eng.workspace['power'] = matlab.double(power.tolist())" - ] - }, - { - "cell_type": "code", - "execution_count": 110, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
PowerSetpointOutsideTempSupplyTempInsideTempSolRadHeatSimulatedTemptime
timestamp
2017-06-30 20:00:00+02:008.62069022.519.022.122.283333178.6182678.62069022.2832850
2017-06-30 20:05:00+02:008.62069022.519.022.122.233333135.1221008.62069022.266734300
2017-06-30 20:10:00+02:007.56666722.519.022.122.366667126.2199677.56666722.258667600
2017-06-30 20:15:00+02:008.17241422.519.022.122.366667123.7724678.17241422.255146900
2017-06-30 20:20:00+02:008.37931022.517.522.422.566667106.163600-25.13793122.2552041200
\n", - "
" - ], - "text/plain": [ - " Power Setpoint OutsideTemp SupplyTemp \\\n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 8.620690 22.5 19.0 22.1 \n", - "2017-06-30 20:05:00+02:00 8.620690 22.5 19.0 22.1 \n", - "2017-06-30 20:10:00+02:00 7.566667 22.5 19.0 22.1 \n", - "2017-06-30 20:15:00+02:00 8.172414 22.5 19.0 22.1 \n", - "2017-06-30 20:20:00+02:00 8.379310 22.5 17.5 22.4 \n", - "\n", - " InsideTemp SolRad Heat SimulatedTemp \\\n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 22.283333 178.618267 8.620690 22.283285 \n", - "2017-06-30 20:05:00+02:00 22.233333 135.122100 8.620690 22.266734 \n", - "2017-06-30 20:10:00+02:00 22.366667 126.219967 7.566667 22.258667 \n", - "2017-06-30 20:15:00+02:00 22.366667 123.772467 8.172414 22.255146 \n", - "2017-06-30 20:20:00+02:00 22.566667 106.163600 -25.137931 22.255204 \n", - "\n", - " time \n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 0 \n", - "2017-06-30 20:05:00+02:00 300 \n", - "2017-06-30 20:10:00+02:00 600 \n", - "2017-06-30 20:15:00+02:00 900 \n", - "2017-06-30 20:20:00+02:00 1200 " - ] - }, - "execution_count": 110, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df = pd.DataFrame(df).assign(time = df_wdb['time'].values)\n", - "df.head()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Control loop" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Initiation setup" - ] - }, - { - "cell_type": "code", - "execution_count": 111, - "metadata": {}, - "outputs": [], - "source": [ - "current_timestamp = 1500" - ] - }, - { - "cell_type": "code", - "execution_count": 112, - "metadata": {}, - "outputs": [], - "source": [ - "df_power = df['Heat']\n", - "df_power = pd.DataFrame(df_power).assign(time = df_wdb['time'].values)\n", - "df_power.loc[df_power['time'] >= current_timestamp, 'Heat'] = np.NaN" - ] - }, - { - "cell_type": "code", - "execution_count": 113, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Heattime
timestamp
2017-06-30 20:00:00+02:008.6206900
2017-06-30 20:05:00+02:008.620690300
2017-06-30 20:10:00+02:007.566667600
2017-06-30 20:15:00+02:008.172414900
2017-06-30 20:20:00+02:00-25.1379311200
.........
2017-07-03 05:35:00+02:00NaN207300
2017-07-03 05:40:00+02:00NaN207600
2017-07-03 05:45:00+02:00NaN207900
2017-07-03 05:50:00+02:00NaN208200
2017-07-03 05:55:00+02:00NaN208500
\n", - "

696 rows × 2 columns

\n", - "
" - ], - "text/plain": [ - " Heat time\n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 8.620690 0\n", - "2017-06-30 20:05:00+02:00 8.620690 300\n", - "2017-06-30 20:10:00+02:00 7.566667 600\n", - "2017-06-30 20:15:00+02:00 8.172414 900\n", - "2017-06-30 20:20:00+02:00 -25.137931 1200\n", - "... ... ...\n", - "2017-07-03 05:35:00+02:00 NaN 207300\n", - "2017-07-03 05:40:00+02:00 NaN 207600\n", - "2017-07-03 05:45:00+02:00 NaN 207900\n", - "2017-07-03 05:50:00+02:00 NaN 208200\n", - "2017-07-03 05:55:00+02:00 NaN 208500\n", - "\n", - "[696 rows x 2 columns]" - ] - }, - "execution_count": 113, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df_power" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Compute input to apply at current time step" - ] - }, - { - "cell_type": "code", - "execution_count": 114, - "metadata": {}, - "outputs": [], - "source": [ - "u_1 = float(df_power.loc[df['time'] == (current_timestamp - 300 * 1), 'Heat'])\n", - "\n", - "y_0 = float(df.loc[df['time'] == (current_timestamp - 300 * 0), 'SimulatedTemp'])\n", - "y_1 = float(df.loc[df['time'] == (current_timestamp - 300 * 1), 'SimulatedTemp'])\n", - "y_2 = float(df.loc[df['time'] == (current_timestamp - 300 * 2), 'SimulatedTemp'])" - ] - }, - { - "cell_type": "code", - "execution_count": 115, - "metadata": {}, - "outputs": [], - "source": [ - "real_x0 = np.array([u_1, y_0, y_1, y_2])" - ] - }, - { - "cell_type": "code", - "execution_count": 116, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([-25.13793103, 22.25542067, 22.25520408, 22.25514583])" - ] - }, - "execution_count": 116, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "real_x0" - ] - }, - { - "cell_type": "code", - "execution_count": 117, - "metadata": {}, - "outputs": [], - "source": [ - "#real_x0 = x_scaler.transform(np.hstack([np.zeros((3)), real_x0]).reshape((1, -1)))[0, 3:]" - ] - }, - { - "cell_type": "code", - "execution_count": 118, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([-25.13793103, 22.25542067, 22.25520408, 22.25514583])" - ] - }, - "execution_count": 118, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "real_x0" - ] - }, - { - "cell_type": "code", - "execution_count": 119, - "metadata": {}, - "outputs": [], - "source": [ - "iter_idx = (df['time'] >= current_timestamp)\n", - "real_W = df[iter_idx].iloc[:N_horizon, [5, 2]].to_numpy()" - ] - }, - { - "cell_type": "code", - "execution_count": 120, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[92.377 , 17.5 ],\n", - " [54.69793333, 16. ],\n", - " [51.83236667, 16. ],\n", - " [40.2653 , 16. ],\n", - " [29.8469 , 16. ]])" - ] - }, - "execution_count": 120, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "real_W = casadi.DM(real_W)" - ] - }, - { - "cell_type": "code", - "execution_count": 123, - "metadata": {}, - "outputs": [], - "source": [ - "real_p = casadi.vertcat(\n", - " casadi.vec(real_W),\n", - " casadi.vec(real_x0)\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 124, - "metadata": {}, - "outputs": [], - "source": [ - "real_X0 = get_gpr_horizon_array(real_W, real_x0.reshape((1, -1)), np.zeros((N_horizon, 1)))" - ] - }, - { - "cell_type": "code", - "execution_count": 125, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[ 92.377 , 17.5 , 0. , -25.13793103,\n", - " 22.25542067, 22.25520408, 22.25514583],\n", - " [ 54.69793333, 16. , 0. , 0. ,\n", - " 22.27224362, 22.25542067, 22.25520408],\n", - " [ 51.83236667, 16. , 0. , 0. ,\n", - " 22.31347714, 22.27224362, 22.25542067],\n", - " [ 40.2653 , 16. , 0. , 0. ,\n", - " 22.39125122, 22.31347714, 22.27224362],\n", - " [ 29.8469 , 16. , 0. , 0. ,\n", - " 22.52311219, 22.39125122, 22.31347714]])" - ] - }, - "execution_count": 125, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "real_X0 = casadi.DM(real_X0)" - ] - }, - { - "cell_type": "code", - "execution_count": 126, - "metadata": { - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "This is Ipopt version 3.13.4, running with linear solver mumps.\n", - "NOTE: Other linear solvers might be more efficient (see Ipopt documentation).\n", - "\n", - "Number of nonzeros in equality constraint Jacobian...: 70\n", - "Number of nonzeros in inequality constraint Jacobian.: 5\n", - "Number of nonzeros in Lagrangian Hessian.............: 0\n", - "\n", - "Total number of variables............................: 35\n", - " variables with only lower bounds: 0\n", - " variables with lower and upper bounds: 0\n", - " variables with only upper bounds: 0\n", - "Total number of equality constraints.................: 30\n", - "Total number of inequality constraints...............: 5\n", - " inequality constraints with only lower bounds: 0\n", - " inequality constraints with lower and upper bounds: 5\n", - " inequality constraints with only upper bounds: 0\n", - "\n", - "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 0 8.5409843e+03 6.54e+01 2.57e+01 0.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", - " 1 5.4574789e+02 7.59e+00 1.87e+02 -1.5 6.54e+01 - 9.89e-01 1.00e+00f 1\n", - " 2 3.4479029e+01 4.23e-01 6.34e+01 -3.5 7.50e+00 - 9.90e-01 1.00e+00f 1\n", - " 3 1.7089734e+01 2.18e-01 7.77e+00 -1.4 1.43e+00 - 1.00e+00 1.00e+00f 1\n", - " 4 1.3407625e+01 1.02e-01 5.52e+00 -3.2 1.47e+00 - 1.00e+00 1.00e+00f 1\n", - " 5 2.9545703e+01 1.38e-01 1.14e+01 -5.1 1.45e+00 - 1.00e+00 1.00e+00h 1\n", - " 6 2.6625453e+01 2.30e-02 1.46e+01 -7.0 3.08e-01 - 1.00e+00 1.00e+00f 1\n", - " 7 1.1836963e+01 4.04e-01 1.68e+01 -8.9 2.76e+00 - 1.00e+00 1.00e+00f 1\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "CasADi - 2021-04-15 12:50:16 WARNING(\"solver:nlp_jac_g failed:Error in Function::operator() for 'nlp_jac_g' [MXFunction] at .../casadi/core/function.cpp:1368:\n", - "Error in Function::operator() for 'fwd7_GPR' [CentralDiff] at .../casadi/core/function.cpp:1368:\n", - "Error in Function::operator() for 'GPR' [CallbackInternal] at .../casadi/core/function.cpp:1368:\n", - ".../casadi/core/function_internal.cpp:3366: Failed to evaluate 'eval_dm' for GPR:\n", - ".../casadi/core/callback_internal.cpp:122: Error calling \"eval\" for object GPR:\n", - "KeyboardInterrupt\") [.../casadi/core/oracle_function.cpp:223]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Number of Iterations....: 8\n", - "\n", - "Number of objective function evaluations = 9\n", - "Number of objective gradient evaluations = 9\n", - "Number of equality constraint evaluations = 9\n", - "Number of inequality constraint evaluations = 9\n", - "Number of equality constraint Jacobian evaluations = 8\n", - "Number of inequality constraint Jacobian evaluations = 9\n", - "Number of Lagrangian Hessian evaluations = 0\n", - "Total CPU secs in IPOPT (w/o function evaluations) = 8.005\n", - "Total CPU secs in NLP function evaluations = 71.463\n", - "\n", - "EXIT: Invalid number in NLP function or derivative detected.\n", - " solver : t_proc (avg) t_wall (avg) n_eval\n", - " nlp_f | 1.36 s (150.85ms) 778.94ms ( 86.55ms) 9\n", - " nlp_g | 1.07 s (118.45ms) 616.40ms ( 68.49ms) 9\n", - " nlp_grad | 8.54 s ( 8.54 s) 5.05 s ( 5.05 s) 1\n", - " nlp_grad_f | 45.21 s ( 4.52 s) 26.23 s ( 2.62 s) 10\n", - " nlp_jac_g | 33.52 s ( 3.35 s) 19.31 s ( 1.93 s) 10\n", - " total | 89.70 s ( 89.70 s) 52.01 s ( 52.01 s) 1\n" - ] - } - ], - "source": [ - "res = solver(x0 = real_X0.reshape((-1, )), p = real_p, lbg = real_lbg, ubg = real_ubg)" - ] - }, - { - "cell_type": "code", - "execution_count": 59, - "metadata": {}, - "outputs": [], - "source": [ - "x_res = np.array(res['x'].reshape((N_horizon, -1)))\n", - "#x_res_sc = x_scaler.inverse_transform(x_res)" - ] - }, - { - "cell_type": "code", - "execution_count": 60, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "DM(\n", - "[[92.377, 17.5, 0.415039, -25.1379, 22.2554, 22.2552, 22.2551], \n", - " [54.6979, 16, 144.993, 0.415039, 22.2437, 22.2554, 22.2552], \n", - " [51.8324, 16, -103.495, 144.993, 22.272, 22.2437, 22.2554], \n", - " [40.2653, 16, -48.8636, -103.495, 23.3874, 22.272, 22.2437], \n", - " [29.8469, 16, 6.33578, -48.8636, 24.9003, 23.3874, 22.272]])" - ] - }, - "execution_count": 60, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "casadi.DM(x_res)" - ] - }, - { - "cell_type": "code", - "execution_count": 61, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([ 0.41503917, 144.99327115, -103.49457825, -48.8635556 ,\n", - " 6.33578304])" - ] - }, - "execution_count": 61, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "x_res[:, 2]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Apply the first computed input as the next input:" - ] - }, - { - "cell_type": "code", - "execution_count": 62, - "metadata": {}, - "outputs": [], - "source": [ - "df_power.loc[df_power['time'] == current_timestamp, 'Heat'] = res['x'].reshape((N_horizon, -1))[0, 2]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Simulate the building including the current input:" - ] - }, - { - "cell_type": "code", - "execution_count": 63, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
Heattime
timestamp
2017-06-30 20:00:00+02:008.6206900
2017-06-30 20:05:00+02:008.620690300
2017-06-30 20:10:00+02:007.566667600
2017-06-30 20:15:00+02:008.172414900
2017-06-30 20:20:00+02:00-25.1379311200
2017-06-30 20:25:00+02:000.4150391500
\n", - "
" - ], - "text/plain": [ - " Heat time\n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 8.620690 0\n", - "2017-06-30 20:05:00+02:00 8.620690 300\n", - "2017-06-30 20:10:00+02:00 7.566667 600\n", - "2017-06-30 20:15:00+02:00 8.172414 900\n", - "2017-06-30 20:20:00+02:00 -25.137931 1200\n", - "2017-06-30 20:25:00+02:00 0.415039 1500" - ] - }, - "execution_count": 63, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df_power.dropna()" - ] - }, - { - "cell_type": "code", - "execution_count": 64, - "metadata": {}, - "outputs": [], - "source": [ - "power = np.array(df_power[['time', 'Heat']].dropna())" - ] - }, - { - "cell_type": "code", - "execution_count": 65, - "metadata": {}, - "outputs": [], - "source": [ - "eng.workspace['power'] = matlab.double(power.tolist())" - ] - }, - { - "cell_type": "code", - "execution_count": 66, - "metadata": {}, - "outputs": [], - "source": [ - "eng.set_param('polydome', 'StopTime', str(current_timestamp + 300), nargout = 0)" - ] - }, - { - "cell_type": "code", - "execution_count": 67, - "metadata": {}, - "outputs": [], - "source": [ - "eng.workspace['result'] = eng.sim('polydome')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Interpret the MATLAB results as python:" - ] - }, - { - "cell_type": "code", - "execution_count": 68, - "metadata": {}, - "outputs": [], - "source": [ - "dict_simulation = {}\n", - "dict_simulation['values'] = np.asarray(eng.eval('result.SimulatedTemp.Data')).reshape(-1)\n", - "dict_simulation['time'] = np.asarray(eng.eval('result.SimulatedTemp.Time')).reshape(-1)" - ] - }, - { - "cell_type": "code", - "execution_count": 69, - "metadata": {}, - "outputs": [], - "source": [ - "df_simulation = pd.DataFrame(dict_simulation)\n", - "#df_simulation['time'] = df_simulation['time'].astype(int)\n", - "df_simulation.set_index('time', inplace = True, drop = True)" - ] - }, - { - "cell_type": "code", - "execution_count": 70, - "metadata": {}, - "outputs": [], - "source": [ - "df_simulation['timestamp'] = df.index[0] + df_simulation.index.map(lambda x: pd.Timedelta(seconds = x))" - ] - }, - { - "cell_type": "code", - "execution_count": 71, - "metadata": {}, - "outputs": [], - "source": [ - "df_simulation = df_simulation.reset_index().set_index('timestamp')" - ] - }, - { - "cell_type": "code", - "execution_count": 72, - "metadata": {}, - "outputs": [], - "source": [ - "df_resampled_5 = df_simulation['values'].resample('5min').mean().pad()" - ] - }, - { - "cell_type": "code", - "execution_count": 73, - "metadata": {}, - "outputs": [], - "source": [ - "df_simulation = pd.concat([df['time'], df_resampled_5], axis = 1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Update the simulated temperature (output data) with the new info from this step:" - ] - }, - { - "cell_type": "code", - "execution_count": 74, - "metadata": {}, - "outputs": [], - "source": [ - "df.loc[:, 'SimulatedTemp'] = df_simulation['values']" - ] - }, - { - "cell_type": "code", - "execution_count": 75, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
PowerSetpointOutsideTempSupplyTempInsideTempSolRadHeatSimulatedTemptime
timestamp
2017-06-30 20:00:00+02:008.62069022.519.022.122.283333178.6182678.62069022.2587760
2017-06-30 20:05:00+02:008.62069022.519.022.122.233333135.1221008.62069021.924169300
2017-06-30 20:10:00+02:007.56666722.519.022.122.366667126.2199677.56666721.781290600
2017-06-30 20:15:00+02:008.17241422.519.022.122.366667123.7724678.17241421.651448900
2017-06-30 20:20:00+02:008.37931022.517.522.422.566667106.163600-25.13793121.4714031200
..............................
2017-07-03 05:35:00+02:00-22.26666722.516.021.521.0166673.911200-22.266667NaN207300
2017-07-03 05:40:00+02:00-21.31034522.516.021.520.8500004.535500-21.310345NaN207600
2017-07-03 05:45:00+02:00-22.00000022.516.021.520.8500005.259500-22.000000NaN207900
2017-07-03 05:50:00+02:00-23.03448322.516.021.520.8666676.644067-23.034483NaN208200
2017-07-03 05:55:00+02:00-20.66666722.516.021.520.8500009.509900-20.666667NaN208500
\n", - "

696 rows × 9 columns

\n", - "
" - ], - "text/plain": [ - " Power Setpoint OutsideTemp SupplyTemp \\\n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 8.620690 22.5 19.0 22.1 \n", - "2017-06-30 20:05:00+02:00 8.620690 22.5 19.0 22.1 \n", - "2017-06-30 20:10:00+02:00 7.566667 22.5 19.0 22.1 \n", - "2017-06-30 20:15:00+02:00 8.172414 22.5 19.0 22.1 \n", - "2017-06-30 20:20:00+02:00 8.379310 22.5 17.5 22.4 \n", - "... ... ... ... ... \n", - "2017-07-03 05:35:00+02:00 -22.266667 22.5 16.0 21.5 \n", - "2017-07-03 05:40:00+02:00 -21.310345 22.5 16.0 21.5 \n", - "2017-07-03 05:45:00+02:00 -22.000000 22.5 16.0 21.5 \n", - "2017-07-03 05:50:00+02:00 -23.034483 22.5 16.0 21.5 \n", - "2017-07-03 05:55:00+02:00 -20.666667 22.5 16.0 21.5 \n", - "\n", - " InsideTemp SolRad Heat SimulatedTemp \\\n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 22.283333 178.618267 8.620690 22.258776 \n", - "2017-06-30 20:05:00+02:00 22.233333 135.122100 8.620690 21.924169 \n", - "2017-06-30 20:10:00+02:00 22.366667 126.219967 7.566667 21.781290 \n", - "2017-06-30 20:15:00+02:00 22.366667 123.772467 8.172414 21.651448 \n", - "2017-06-30 20:20:00+02:00 22.566667 106.163600 -25.137931 21.471403 \n", - "... ... ... ... ... \n", - "2017-07-03 05:35:00+02:00 21.016667 3.911200 -22.266667 NaN \n", - "2017-07-03 05:40:00+02:00 20.850000 4.535500 -21.310345 NaN \n", - "2017-07-03 05:45:00+02:00 20.850000 5.259500 -22.000000 NaN \n", - "2017-07-03 05:50:00+02:00 20.866667 6.644067 -23.034483 NaN \n", - "2017-07-03 05:55:00+02:00 20.850000 9.509900 -20.666667 NaN \n", - "\n", - " time \n", - "timestamp \n", - "2017-06-30 20:00:00+02:00 0 \n", - "2017-06-30 20:05:00+02:00 300 \n", - "2017-06-30 20:10:00+02:00 600 \n", - "2017-06-30 20:15:00+02:00 900 \n", - "2017-06-30 20:20:00+02:00 1200 \n", - "... ... \n", - "2017-07-03 05:35:00+02:00 207300 \n", - "2017-07-03 05:40:00+02:00 207600 \n", - "2017-07-03 05:45:00+02:00 207900 \n", - "2017-07-03 05:50:00+02:00 208200 \n", - "2017-07-03 05:55:00+02:00 208500 \n", - "\n", - "[696 rows x 9 columns]" - ] - }, - "execution_count": 75, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Keep track of the all the prediction horizon, to add to the graph:" - ] - }, - { - "cell_type": "code", - "execution_count": 76, - "metadata": {}, - "outputs": [], - "source": [ - "gpr_horizon = np.array(gpr(res['x'].reshape((N_horizon, -1)).T)).flatten()" - ] - }, - { - "cell_type": "code", - "execution_count": 77, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([22.25542067, 22.2722139 , 22.23139534, 22.34675811, 25.0869712 ,\n", - " 27.19437496])" - ] - }, - "execution_count": 77, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "T_sim_horizon = np.hstack([np.array(y_0), gpr_horizon])\n", - "T_sim_horizon" - ] - }, - { - "cell_type": "code", - "execution_count": 78, - "metadata": {}, - "outputs": [], - "source": [ - "simul_idx = (df_simulation['time'] >= current_timestamp) & (df_simulation['time'] <= (current_timestamp + N_horizon * 300))" - ] - }, - { - "cell_type": "code", - "execution_count": 79, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
timevalues
timestamp
2017-06-30 20:25:00+02:00150022.255421
2017-06-30 20:30:00+02:00180022.272214
2017-06-30 20:35:00+02:00210022.231395
2017-06-30 20:40:00+02:00240022.346758
2017-06-30 20:45:00+02:00270025.086971
2017-06-30 20:50:00+02:00300027.194375
\n", - "
" - ], - "text/plain": [ - " time values\n", - "timestamp \n", - "2017-06-30 20:25:00+02:00 1500 22.255421\n", - "2017-06-30 20:30:00+02:00 1800 22.272214\n", - "2017-06-30 20:35:00+02:00 2100 22.231395\n", - "2017-06-30 20:40:00+02:00 2400 22.346758\n", - "2017-06-30 20:45:00+02:00 2700 25.086971\n", - "2017-06-30 20:50:00+02:00 3000 27.194375" - ] - }, - "execution_count": 79, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df_T_sim_horizon = df_simulation[simul_idx].copy()\n", - "df_T_sim_horizon.loc[:, 'values'] = T_sim_horizon.reshape((-1, ))\n", - "df_T_sim_horizon" - ] - }, - { - "cell_type": "code", - "execution_count": 80, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAE/CAYAAADCCbvWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5CUlEQVR4nO3dd5wdVf3/8ddnWzqQQIKEjoBCCiEFKQIJUQGRL1gAFSF0aVJsoH75EvwK0pSmEkJRVL4qSBULNaEIpEEAIRSRAIH8IIQaUrad3x9zd3N3s5st2d17d/N6Ph7z2Jkzc2bOvXsM+/bMnImUEpIkSZKkwiopdAMkSZIkSYYzSZIkSSoKhjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSeohImJJ3lIbEcvytg8tdPvaIyLmR8RnCt2O1YmI6RFxTCeef2pEPJ/7nR7RaN8REVHT6Hc/Pm//oIi4NSI+iohXIuLrjepPjIjnImJpREyLiM0763NIklpmOJOkHiKl1L9uAV4F9s8ru6HQ7WssIsp6wjW6wJPAicDjzex/NP93n1Kanrfvl0AlsCFwKHBlRAwDiIgNgFuAs4BBwGzgT53zESRJrWE4k6QeLiJKIuLMiHgpIhZHxI0RMSi3b4uISBFxZES8FhHvRsTxETEuIp6KiPci4hd55zoiIv4ZEVdExPu5UZeJefvXjYhrI2JhRLweET+JiNJGdS+JiHeAyRHx8Yi4P9eutyPihohYL3f874DNgL/kRoS+HxHjI2JBo89XP7oWEZMj4s8R8fuI+AA4ooU2bR0RD+Q+y9sR0WQ4iYjeuXMuzn0nsyJiw4g4F9gd+EWujb/IHf/JiLgnIt7JjXodnHeu30TElNz+D3PXb3bEKqX0y5TSfcDy1v/WISL6AV8GzkopLUkpPQzcARyWO+RLwDMppZtSSsuBycAOEfHJtlxHktRxDGeS1POdAhwI7AkMBd4lG1HJ9ylgG+AQ4FLgR8BngGHAwRGxZ6Nj/wNsAJwN3FIX9oDrgWpga2BH4HPAMU3UHQKcCwTw01y7tgM2JQsJpJQOo+EI4IWt/LwHAH8G1gNuaKFN/wvcDQwENgGuaOack4B1c+1bHzgeWJZS+hHwEHByro0n50LRPcD/5T7n14Bf1Y1Y5Ryau/YGwNxcO9trx1ywfCEizsobLdwWqEkpvZB37JNkv1NyP5+s25FS+gh4KW+/JKmLGc4kqef7JvCjlNKClNIKsvDzlUa3/P1vSml5Sulu4CPgDymlt1JKr5OFjx3zjn0LuDSlVJVS+hPwPLBfRGwI7AucllL6KKX0FnAJ8NW8um+klK5IKVWnlJallP6dUronpbQipbQI+DlZiFwTj6aUbksp1QLrtNCmKmBzYGju8z/czDmryELZ1imlmpTSnJTSB80c+wVgfkrp17nP+ThwM/CVvGP+mlJ6MPf7+BGwS0Rs2o7P+iAwnCwEfpksCH4vt68/8H6j498HBrRyvySpi/WEe/ElSau3OXBrRNTmldWQPYdU58289WVNbPfP2349pZTytl8hG/naHCgHFkZE3b4S4LW8Y/PXiYghwOVktwYOyB3/bqs+VfPyr9FSm75PNoI1MyLeBX6WUrquiXP+jmzU7I+52y5/TxZ4q5o4dnPgUxHxXl5ZWe4cq7QxpbQkd5vn0EZtb1FK6T95m09HxI/JwtlPgSVk4TTfOsCHufWW9kuSupgjZ5LU870G7JtSWi9v6Z0bFWuPjSMv6ZA9F/ZG7jorgA3yrrNOSin/Nrn8UAdZiEjAyJTSOsA3yG51bO74j4C+dRu5Z8cGNzomv85q25RS+n8ppWNTSkPJRhh/FRFbN/7AuVHCc1JK2wO7ko2OHd5MG18DHmj0ffdPKZ2Qd0z9KFlE9CebkOONxtdth8TK7+8FoCwitsnbvwPwTG79mdx2XTv6AR/P2y9J6mKGM0nq+aYA59ZNOhERgyPigDU43xDglIgoj4iDyJ4V+1tKaSHZ81s/i4h1chORfLzR82qNDSAbwXkvIjZm5S15dd4EtsrbfgHoHRH7RUQ58N9Ar+ZO3lKbIuKgiNgkd/i7ZOGmpvF5ImJCRIzIhcEPyG5zrDuucRvvBLaNiMNy31F5boKV7fKO+XxEfDoiKshG7maklJocNYuIiojoTRa6ynOTk5Tk9u2bu52U3EQeZwG35z77R2SzMf44IvpFxG5kz+PVjeDdCgyPiC/nzv8/wFMppeea+z4lSZ3LcCZJPd9lZLP03R0RHwKPkU3M0V4zyCYPeZtsUo+vpJQW5/YdDlQAz5KFnT8DG63mXOcAo8medforWZjI91Pgv3MzJH43pfQ+2bTy1wCvk42kLWD1VtemccCMiFhC9h2dmlJ6uYlzfCxX7wNgHvAA2a2NkH2/X4lspsvLU0ofkk068lWy0bD/B1xAwxD5f2STqbwDjCGbIKQ5d5PdWrorMDW3vkdu30TgqYj4CPgb2fd3Xl7dE4E+ZM8J/gE4IaX0DEDuGb8vk/0O3yXrE/nPB0qSulg0fGxAkqTmRfYS5GNSSp8udFu6q4j4DbAgpfTfhW6LJKm4OHImSZIkSUWgxXCWu7d9ZkQ8GRHPRMQ5ufJBuRdovpj7ObDzmytJkiRJPVOLtzXmZuTql5vqtxx4GDgV+BLwTkrp/Ig4ExiYUjqj01ssSZIkST1QiyNnKbMkt1meWxLZjE/X58qvBw7sjAZKkiRJ0tqgVc+cRURpRMwlm+3pnpTSDGDD3BTFdVMVD+m0VkqSJElSD1fWmoNSSjXAqIhYD7g1Ioa39gIRcRxwHEC/fv3GfPKTn2xPOyVJkiSp25szZ87bKaXBTe1rVTirk1J6LyKmA/sAb0bERimlhRGxEdmoWlN1ppK9l4WxY8em2bNnt6nxkiRJktRTRMQrze1rzWyNg3MjZkREH+AzwHNkL+uclDtsEnD7GrdUkiRJktZSrRk52wi4PiJKycLcjSmlOyPiUeDGiDgaeBU4qBPbKUmSJEk9WovhLKX0FLBjE+WLgYmd0ShJkiRJWtu06ZkzSZIkqaeoqqpiwYIFLF++vNBNUQ/Uu3dvNtlkE8rLy1tdx3AmSZKktdKCBQsYMGAAW2yxBRFR6OaoB0kpsXjxYhYsWMCWW27Z6nqtes+ZJEmS1NMsX76c9ddf32CmDhcRrL/++m0elTWcSZIkaa1lMFNnaU/fMpxJkiRJBXLuuecybNgwRo4cyahRo5gxYwYAxxxzDM8++2yHXGOLLbbg7bffXu0x5513XpvP+5vf/IaTTz65Qdmvf/1rRo0axahRo6ioqGDEiBGMGjWKM888s83n7wqXXnopS5cuLXQz6vnMmSRJklQAjz76KHfeeSePP/44vXr14u2336ayshKAa665pkvbct555/HDH/5wjc9z5JFHcuSRRwJZKJw2bRobbLDBGp+3vVJKpJQoKWl6TOrSSy/lG9/4Bn379m31Oaurqykr65wY5ciZJEmSVAALFy5kgw02oFevXgBssMEGDB06FIDx48cze/ZsAPr3788ZZ5zBmDFj+MxnPsPMmTMZP348W221FXfccQew6ijWF77wBaZPn77KNQ888EDGjBnDsGHDmDp1KgBnnnkmy5YtY9SoURx66KEA/P73v2ennXZi1KhRfPOb36SmpgbIRsa23XZb9txzT/75z3+2+rNedNFFjBs3jpEjR3L22WcDMH/+fD75yU9yzDHHMHz4cA499FDuvfdedtttN7bZZhtmzpwJwOTJkznssMPYa6+92Gabbbj66qtbPO92223HiSeeyOjRo3nttdc44YQTGDt2LMOGDas/7vLLL+eNN95gwoQJTJgwof67rvPnP/+ZI444AoAjjjiCb3/720yYMIEzzjiDl156iX322YcxY8aw++6789xzz7X6u1itujTZFcuYMWOSJEmSVAyeffbZgl7/ww8/TDvssEPaZptt0gknnJCmT59ev2/PPfdMs2bNSimlBKS//e1vKaWUDjzwwPTZz342VVZWprlz56YddtghpZTSr3/963TSSSfV199vv/3StGnTUkopbb755mnRokUppZQWL16cUkpp6dKladiwYentt99OKaXUr1+/+rrPPvts+sIXvpAqKytTSimdcMIJ6frrr09vvPFG2nTTTdNbb72VVqxYkXbdddcG12ys7rp33XVXOvbYY1NtbW2qqalJ++23X3rggQfSyy+/nEpLS9NTTz2Vampq0ujRo9ORRx6Zamtr02233ZYOOOCAlFJKZ599dho5cmRaunRpWrRoUdpkk03S66+/vtrzRkR69NFH69tS97mrq6vTnnvumZ588slVvpvG38NNN92UJk2alFJKadKkSWm//fZL1dXVKaWU9tprr/TCCy+klFJ67LHH0oQJE5r8DprqY8Ds1Exe8rZGSZIkrfXO+cszPPvGBx16zu2HrsPZ+w9rdn///v2ZM2cODz30ENOmTeOQQw7h/PPPrx+tqVNRUcE+++wDwIgRI+jVqxfl5eWMGDGC+fPnt6lNl19+ObfeeisAr732Gi+++CLrr79+g2Puu+8+5syZw7hx4wBYtmwZQ4YMYcaMGYwfP57BgwcDcMghh/DCCy+0eM27776bu+++mx133BGAJUuW8OKLL7LZZpux5ZZbMmLECACGDRvGxIkTiYhVPtsBBxxAnz596NOnDxMmTGDmzJk8/PDDzZ538803Z+edd66vf+ONNzJ16lSqq6tZuHAhzz77LCNHjmzTd3fQQQdRWlrKkiVLeOSRRzjooIPq961YsaJN52qO4UySJEkqkNLSUsaPH8/48eMZMWIE119//SrhrLy8vH7mv5KSkvrbIEtKSqiurgagrKyM2tra+jpNTeE+ffp07r33Xh599FH69u3L+PHjmzwupcSkSZP46U9/2qD8tttua9cMhCklfvCDH/DNb36zQfn8+fPrP8vqPhusOvNhRKz2vP369avffvnll7n44ouZNWsWAwcO5Igjjmh2ivv86zQ+pu6ctbW1rLfeesydO7elj95mhjNJkiSt9VY3wtVZnn/+eUpKSthmm20AmDt3Lptvvnm7zrXFFlvwq1/9itraWl5//fX657Xyvf/++wwcOJC+ffvy3HPP8dhjj9XvKy8vp6qqivLyciZOnMgBBxzA6aefzpAhQ3jnnXf48MMP+dSnPsWpp57K4sWLWWeddbjpppvYYYcdWmzb3nvvzVlnncWhhx5K//79ef311ykvL2/T57v99tv5wQ9+wEcffcT06dM5//zz6dOnT6vO+8EHH9CvXz/WXXdd3nzzTf7+978zfvx4AAYMGMCHH35YP2nJhhtuyLx58/jEJz7BrbfeyoABA1Y53zrrrMOWW27JTTfdxEEHHURKiaeeeqpV30VLDGeSJElSASxZsoRvfetbvPfee5SVlbH11lvXT9LRVrvttlv9LYLDhw9n9OjRqxyzzz77MGXKFEaOHMknPvGJBrf9HXfccYwcOZLRo0dzww038JOf/ITPfe5z1NbWUl5ezi9/+Ut23nlnJk+ezC677MJGG23E6NGj6ycKWZ3Pfe5zzJs3j1122QXIbuf8/e9/T2lpaas/30477cR+++3Hq6++yllnncXQoUMZOnRoq867ww47sOOOOzJs2DC22mordttttwafe99992WjjTZi2rRpnH/++XzhC19g0003Zfjw4SxZsqTJ9txwww2ccMIJ/OQnP6GqqoqvfvWrHRLOInsmrWuMHTs21c06I0mSJBXSvHnz2G677QrdDLVg8uTJ9O/fn+9+97uFbkqbNdXHImJOSmlsU8c7lb4kSZIkFQFva5QkSZJUtCZPnlzoJnQZR84kSZIkqQgYziRJkiSpCBjOJEmSJKkIGM4kSZIkqQgYziRJkqQCKS0tZdSoUQwfPpyDDjqIpUuXtvtcRxxxBH/+858BOOaYY3j22WebPXb69Ok88sgj9dtTpkzht7/9bbuvXWf+/PkMHz68QdnkyZO5+OKL23SejmpPd+NsjZIkSVKB9OnTh7lz5wJw6KGHMmXKFL797W/X76+pqWnTy5rrXHPNNavdP336dPr378+uu+4KwPHHH9/ma3SW6urqompPV3LkTJIkSWrJhRfCtGkNy6ZNy8o7yO67786///1vpk+fzoQJE/j617/OiBEjqKmp4Xvf+x7jxo1j5MiRXHXVVQCklDj55JPZfvvt2W+//XjrrbfqzzV+/Hhmz54NwD/+8Q9Gjx7NDjvswMSJE5k/fz5TpkzhkksuYdSoUTz00EMNRrfmzp3LzjvvzMiRI/niF7/Iu+++W3/OM844g5122oltt92Whx56qM2fcXXn/uEPf8iee+7JZZddVt+eN954g1GjRtUvpaWlvPLKK7zyyitMnDiRkSNHMnHiRF599VUgGz085ZRT2HXXXdlqq63qRxK7C8OZJEmS1JJx4+Dgg1cGtGnTsu1x4zrk9NXV1fz9739nxIgRAMycOZNzzz2XZ599lmuvvZZ1112XWbNmMWvWLK6++mpefvllbr31Vp5//nmefvpprr766ga3KdZZtGgRxx57LDfffDNPPvkkN910E1tssQXHH388p59+OnPnzmX33XdvUOfwww/nggsu4KmnnmLEiBGcc845Ddo5c+ZMLr300gbl+V566aUGgWrKlCmtOvd7773HAw88wHe+8536sqFDhzJ37lzmzp3Lsccey5e//GU233xzTj75ZA4//HCeeuopDj30UE455ZT6OgsXLuThhx/mzjvv5Mwzz2zjb6KwvK1RkiRJOu00yN1e2KyhQ2HvvWGjjWDhQthuOzjnnGxpyqhRcOmlqz3lsmXLGDVqFJCNnB199NE88sgj7LTTTmy55ZYA3H333Tz11FP1o0Dvv/8+L774Ig8++CBf+9rXKC0tZejQoey1116rnP+xxx5jjz32qD/XoEGDVtue999/n/fee48999wTgEmTJnHQQQfV7//Sl74EwJgxY5g/f36T5/j4xz9ef6smrHyJdEvnPuSQQ5pt1z//+U+uueaa+tG6Rx99lFtuuQWAww47jO9///v1xx544IGUlJSw/fbb8+abb6728xYbw5kkSZLUGgMHZsHs1Vdhs82y7TWU/8xZvn79+tWvp5S44oor2HvvvRsc87e//Y2IWO35U0otHtMWvXr1ArKJTKqrqzvsvNDwM+dbuHAhRx99NHfccQf9+/dv8pj8z1jXRsg+f3diOJMkSZJaGOECVt7KeNZZcOWVcPbZMGFCpzdt77335sorr2SvvfaivLycF154gY033pg99tiDq666isMPP5y33nqLadOm8fWvf71B3V122YWTTjqJl19+mS233JJ33nmHQYMGMWDAAD744INVrrXuuusycOBAHnroIXbffXd+97vf1Y90ran2nLuqqoqDDz6YCy64gG233ba+fNddd+WPf/wjhx12GDfccAOf/vSnO6SNhWY4kyRJklpSF8xuvDELZBMmNNzuRMcccwzz589n9OjRpJQYPHgwt912G1/84he5//77GTFiBNtuu22TQWfw4MFMnTqVL33pS9TW1jJkyBDuuece9t9/f77yla9w++23c8UVVzSoc/3113P88cezdOlSttpqK37961932Gdp67kfeeQRZs2axdlnn83ZZ58NZCOGl19+OUcddRQXXXQRgwcP7tA2FlK0NNQXEZsCvwU+BtQCU1NKl0XEKGAK0BuoBk5MKc1c3bnGjh2b6maNkSRJkgpp3rx5bLfddq07+MILs8k/8oPYtGkwaxbkPe8k5Wuqj0XEnJTS2KaOb83IWTXwnZTS4xExAJgTEfcAFwLnpJT+HhGfz22PX6PWS5IkScWoqQBWN4ImdZAWw1lKaSGwMLf+YUTMAzYGErBO7rB1gTc6q5GSJEmS1NO16ZmziNgC2BGYAZwG3BURF5O9L23Xjm6cJEmSJK0tWv0S6ojoD9wMnJZS+gA4ATg9pbQpcDpwbTP1jouI2RExe9GiRR3RZkmSJKlDdLep1tV9tKdvtSqcRUQ5WTC7IaV0S654ElC3fhOwUzONmppSGptSGjt48OA2N1CSJEnqDL1792bx4sUGNHW4lBKLFy+md+/ebarX4m2Nkb3R7VpgXkrp53m73gD2BKYDewEvtunKkiRJUgFtsskmLFiwAO/uUmfo3bs3m2yySZvqtOaZs92Aw4CnI2JuruyHwLHAZRFRBiwHjmvTlSVJkqQCKi8vZ8sttyx0M6R6rZmt8WEgmtk9pmObI0mSJElrp1ZPCCJJkiRJ6jyGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEk9w4UXwrRpDcumTcvKuwHDmSRJkqSeYdw4OPhguP9+WLYsC2YHH5yVdwNlhW6AJEmSJHWI8ePhtNNg771hxx3h5ZfhxhthwoRCt6xVHDmTJEmS1L2lBHfemY2Q/fd/Q//+MGsWnHBCtwlmYDiTJEmS1F2lBH/9K+y0E+y/P7zzDnzve1BWBmedBVdeueozaEXM2xolSZIkdS8pwd//DpMnZyNkW24J114Lm24KX//6ylsZJ0zInjnrJrc2OnImSZIkqXuoC2U77wz77QeLFsE118Dzz8NRR8ETTzQMYhMmZNuzZhW23a0UKaUuu9jYsWPT7Nmzu+x6kiRJknqAlOCuu7KRshkzYPPNs2fLDj8cKioK3bo2iYg5KaWxTe1z5EySJElScaoLZbvuCvvuCwsXwtSp8MILcMwx3S6YtcRwJkmSJKm4pAR33w277Qb77ANvvAFXXQUvvgjHHtvjQlkdw5kkSZKk4pAS3HMPfPrT2bvKFiyAKVOyUHbccT02lNUxnEmSJEkqrJTg3nth993hc5+D117LpsF/8UX45jd7fCirYziTJEmSVBgpwX33wR57wGc/C6+8Ar/6VRbKjj8eevUqdAu7lOFMkiRJUtdKCe6/H/bcEz7zGXj5ZfjlL+Hf/4YTTljrQlkdw5kkSZKkrpESTJsG48fDxInwn//AL34BL70EJ5641oayOoYzSZIkSZ1v+vQslO21VzZCdsUV2c+TTlrrQ1kdw5kkSZKkzvPAA1komzAhe5bs8suzkbKTT4bevQvduqJiOJMkSZLU8R58MAtk48dnL42+7LLsNsZvfctQ1gzDmSRJkqSO89BD2a2Le+4Jzz0Hl16ajZSdcoqhrAWGM0mSJElr7uGHs0k+9tgD5s3LQtl//gOnngp9+hS6dd2C4UySJElS+/3zn9l0+LvvDs88A5dcYihrpxbDWURsGhHTImJeRDwTEafm7ftWRDyfK7+wc5sqSZIkqWg88kj24uhPfxr+9S/4+c+zUHbaaYaydiprxTHVwHdSSo9HxABgTkTcA2wIHACMTCmtiIghndlQSZIkSUXg0Ufh7LPhnntgyBD42c/g+OOhb99Ct6zbazGcpZQWAgtz6x9GxDxgY+BY4PyU0orcvrc6s6GSJEmSCujRR2HyZLj7bhg8GC6+OAtl/foVumU9RpueOYuILYAdgRnAtsDuETEjIh6IiHHN1DkuImZHxOxFixatcYMlSZIkdaHHHoN99oFdd4UnnoCLLoKXX4bvfMdg1sFaHc4ioj9wM3BaSukDslG3gcDOwPeAGyMiGtdLKU1NKY1NKY0dPHhwBzVbkiRJUqeaMQP23Rd22QXmzIELL8xC2Xe/ayjrJK0KZxFRThbMbkgp3ZIrXgDckjIzgVpgg85ppiRJkqQuMXMmfP7zsPPOMHs2XHBBFsq+9z1DWSdrzWyNAVwLzEsp/Txv123AXrljtgUqgLc7oY2SJEmSOtvMmbDffvCpT2Xr55+fhbLvfx/69y9069YKrZmtcTfgMODpiJibK/shcB1wXUT8C6gEJqWUUqe0UpIkSVLnmDULzjkH/vpXWH99+OlP4aSTYMCAQrdsrdOa2RofBlZ5liznGx3bHEmSJEldYvbsLJTdeScMGgTnnQcnn2woK6DWjJxJkiRJ6inmzMlC2V/+koWyc8+Fb33LUFYEDGeSJEnS2uDxx7NQdscdMHAg/OQnWShbZ51Ct0w5hjNJkiSpJ3viiSyU3X47rLce/O//wimnGMqKkOFMkiRJ6onmzs1C2W23ZaHsxz/OQtm66xa4YWqO4UySJEnqSfJD2brrZuunnmoo6wYMZ5IkSVJP8OSTWRC79dYsiE2enIWy9dYrdMvUSoYzSZIkqTt76qkslN1yS/Yc2dlnw2mnGcq6IcOZJEmS1B09/XQWym6+OQtl//M/WSgbOLDQLVM7lRS6AZIkSZKaceGFMG1aw7LrroORI7PlnnuyUDZ/fhbUDGbdmiNnkiRJUrEaNw4OPhhuvBGGDIGTToIHHoA+feCss7KRskGDCt1KdRDDmSRJklSsJkyAKVNg331hxQqIgG98Ay67zFDWA3lboyRJklSMUoLf/haOPRaqq7Oyb38bfvc7g1kPZTiTJEmSis2rr8LnPw+TJsHQodmEH2edBddfv+ozaOoxDGeSJElSsaithSuvhGHD4MEH4eST4c03sxkZf/zj7Nmzgw82oPVQhjNJkiSpGLz4YvaM2Yknws47w7/+BZtumgWyCROyYyZMyLZnzSpsW9UpIqXUZRcbO3Zsmj17dpddT5IkSSp6NTVwySXZbYu9esHPfgZHHZVN/qEeJyLmpJTGNrXP2RolSZKkQvnXv7IgNmsW/Nd/Zbc0Dh1a6FapQLytUZIkSepqlZXZS6NHj4aXX4Y//AFuu81gtpZz5EySJEnqSrNnZ6NlTz8NX/ta9s6ywYML3SoVAUfOJEmSpK6wbBmccQZ86lOweDHccQf83/8ZzFTPkTNJkiSpsz30EBx9dDYj4zHHwEUXwXrrFbpVKjKOnEmSJEmd5cMPs3eV7bEHVFfDvffC1VcbzNQkw5kkSZLUGe6+G4YPh1/9Ck49NXvGbOLEQrdKRcxwJkmSJHWkd9+FI4+EvfeGPn3g4Yfh0kuhX79Ct0xFznAmSZIkdZRbb4Xtt4ff/Q5++EOYOxd23bXQrVI34YQgkiRJ0pp680341rfgpptg1Cj4299gxx0L3Sp1M46cSZIkSe2VEvz+99lo2e23w09+AjNnGszULi2Gs4jYNCKmRcS8iHgmIk5ttP+7EZEiYoPOa6YkSZJUZF57DfbfHw47DLbdFp54An70IygvL3TL1E21ZuSsGvhOSmk7YGfgpIjYHrLgBnwWeLXzmihJkiQVkdpauOoqGDYMpk2DSy7JJv3YfvtCt0zdXIvhLKW0MKX0eG79Q2AesHFu9yXA94HUaS2UJEmSisVLL2XT4R9/PIwbl02Pf9ppUFpa6JapB2jTM2cRsQWwIzAjIv4LeD2l9GRnNEySJEkqGjU18POfw4gR8PjjMHVq9kLprbYqdMvUg7R6tsaI6A/cDJxGdqvjj4DPtaLeccBxAJtttlm7GilJkiQVzDPPwNFHw4wZ8IUvwJVXwiabFLpV6oFaNXIWEeVkweyGlNItwMeBLYEnI2I+sAnweER8rHHdlNLUlNLYlNLYwYMHd1zLJUmSpM5UVZXNvjh6NPz733DDDXDHHQYzdZoWR84iIoBrgXkppZ8DpJSeBobkHTMfGJtSeruT2ilJkiR1nccfh6OOgiefhEMOgcsvhyFDWq4nrYHWjJztBhwG7BURc3PL5zu5XZIkSVLXW74cfvAD2GkneOstuO02+OMfDWbqEi2OnKWUHgaihWO26KgGSZIkSQXxz39mz5Y9/3w2anbxxTBwYKFbpbVIm2ZrlCRJknqcJUvglFNg992zkbO77oJrrzWYqcu1erZGSZIkqce591449lh45RU4+WQ47zzo37/QrdJaypEzSZIkrX3eey+7hfGzn4WKCnjwwWzSD4OZCshwJkmSpLXL7bfD9tvD9dfDGWfA3Lnw6U8XulWStzVKkiRpLbFoUfZs2R//CCNHwl/+AmPGFLpVUj1HziRJktSzpQR/+EM2WnbzzfDjH8OsWQYzFR1HziRJktRzvf46HH883Hln9u6y666DYcMK3SqpSY6cSZIkqedJCa65Jhstu+8++NnP4JFHDGYqao6cSZIkqWd5+eVsevz77oPx4+Hqq2HrrQvdKqlFjpxJkiSpZ6ipgcsug+HDYeZMmDIlC2gGM3UTjpxJkiSp+5s3L3tv2aOPwuc/nwWzTTctdKukNnHkTJIkSd1XVRWcdx6MGgXPPw+/+102+YfBTN2QI2eSJEnqnp54Ihste+IJOOgguOIK2HDDQrdKajdHziRJktS9LF8OP/oRjBsHCxdm7y678UaDmbo9R84kSZLUfTz6KBx1FDz3HBxxBPz85zBwYKFbJXUIR84kSZJU/D76CE4/HXbbDZYuhX/8A379a4OZehRHziRJklTc7r8/e2/Zf/4DJ54I558PAwYUulVSh3PkTJIkScXp/ffhuONg4kQoLYUHHoBf/tJgph7LcCZJkqTic+edMGwYXHstfO978OSTsMcehW6V1KkMZ5IkSSoeb78Nhx4K+++fPU/22GNw4YXQp0+hWyZ1OsOZJEmSCi8l+NOfYPvt4aabYPJkmDMnmy5fWks4IYgkSZIK6403sok+br89C2PXXQfDhxe6VVKXc+RMkiRJXefCC2HatGw9pSyIbbst/PWvcNFF8MgjBjOttQxnkqSeLf8PwTrTpmXlaprfmTrT6NFw0EHw29/C3nvD0UdDZWUW0r77XSjzxi6tvez9kqSebdw4OPhguPFGmDAhCxl122qa31nbXHhh9p1NmLCybNo0mDULvv/9wrULspGpykpYsWLNluXL1/wcdUttbda2SZOgvBz694fbbsumy5fWcpFS6rKLjR07Ns2ePbvLrtcaN89ZwLtLKxnUr4KB/SpYv18FA/tWsH7/CvqUlxIRhW6iJGlNTZsGX/kK7LQTPPQQHHkkbL119odrRy+1tZ1z3q5uz+LF2dTlm20Gr70GY8bAkCHZu6ZKSlb92VRZV/0s5LXr3r311a9m4XWPPeDuu+Gww+Cqq2Ds2K4PP/lLZWXH/e+orAx69eq45f774a674Kyz4Mc/7rh2SkUuIuaklMY2uW9tD2dfvvIR5rzybpP7epWVMKhfRYNlYN+KVcrqygf2Laes1DtFJakoHXJI5478RGR/rEcUz7Km7VmwAF5/HTbaCD72MaipycJe459NlbX0s270RM3ryCDU2qV376bLKyqyINpR6kZjTzgBrrxy5SittBZYXThr8bbGiNgU+C3wMaAWmJpSuiwiLgL2ByqBl4AjU0rvdViru8ifj9+FD5ZV887SSt75KFve/aiywXbd8uo7S3lnSSUfrqhu9nzr9inPRt/qRuDyR+T6VTCoXzmD+vViUN8KBvWvoF+Fo3OS1OmmTcv+X/pvfxuuvx6uuSYb4eioENQT1f3xfNZZ2R/PP/tZx/7xXDdC19ZQ1xHBsKN+Ni67/354+GHYa6/sHV1rEozKy3t+36oLZBMmNNyW1mItjpxFxEbARimlxyNiADAHOBDYBLg/pVQdERcApJTOWN25inHkrD0qq2t5d+nKILf4o8r67cZLXXlVTdPfc0VpCQPrAlt9cCtvFOhyS99su9zROUlqvcZ/CDbe1qr8ztrOkaDWK+Zn9KQu0KG3NUbE7cAvUkr35JV9EfhKSunQ1dXtKeGsrVJKLFlR3SCwLV6S+1k3UvdRFe98tIJ3l1bxzkeVvL+sqtnzDehdtjK4NbrNsr6s/8qfA3qVOTonae3lH4Jt53fWNoZZSW3QYeEsIrYAHgSGp5Q+yCv/C/CnlNLvV1d/bQ1n7VFVU8t7uaBWH+jqg1zTIa+yuun798tKosFkJ3XBLX90bv285+kG9iunV1kH3lcuSVJPZpiV1AYdEs4ioj/wAHBuSumWvPIfAWOBL6UmThYRxwHHAWy22WZjXnnllbZ/ArUopcTSypqVt1QubRjkmgp57y2rorlff/9eZXkjcStvuxzYr4J1epfTu7yU3uUl9CrLfvYuL6V3WSm9ykvoXbZyX6/yEnqVlThyJ0mSJNEB4SwiyoE7gbtSSj/PK58EHA9MTCktbek8jpwVl5raxHtLV3+b5TtLqxqEvGVVNW2+TkQ282V+kOtV1jDQNbmvmfJeZc3vy18vLTEQSpIkqbis6WyNAVwLzGsUzPYBzgD2bE0wU/EpLQnW79+L9fv3YushrauzrLKGD5dXsbyqlhXVNSyvqmV5dQ0rqmpZXlXD8lxZ/b6qGlZU17Kiqibb30S9JSuqeXtJJSvyz5OrV13b/lc9lJfGKuGvuUCXP9LXcARwdfXy6uSCZnlpOEooSZKkdmkxnAG7AYcBT0fE3FzZD4HLgV7APbk/Rh9LKR3fGY1U8ehTUUqfiq57Hq26ppbldeGuumFwW16Vhbn8ILhyX15orK5p8tgPl1c3DJS5azT37F5rRNDgts78QJd/m2dF3XZZtt2rPG8977j6Y8pLqChtujz/PL5nT5IkqftqMZyllB4GmhoK+FvHN0dqqKy0hP6lJfTv1Zr/H6Fj1NYmKmtq60f68sNgg/W80Lii8b5GQbCyppYVVbV8tKKadz6qzUYTq+sCY269urbZZwBbq7QkGoS8ivrw1ygMNhMA69YrVhMAGwTJvPWKshIqSn2+UJIkqb267i9eqZsoKQl6l2S3NHallBLVtan+NtAstDUR4qpqs7DXRHl+0KsbKawLhnX7lqyort9Xf40OuI20TutHA0upKC1Zpbx3eSl9ykvoW1FGn4pS+uZGa/tWlGXr5VlZ34oyepcbBiVJUs9hOJOKRERQXhqUd/FIYb6a2kRlEwGv2QDYKNytLihmIbGG95dVNXtcW28pjaA+rPWpKKVv+cpAl5WV0be8dNWyuvXy0lwILKFPeV55Lvw5qYwkSepKhjNJ9UpLosufK8xXd0vp0soallZWs6yyJrdew7Kq6pXr9T9zZVV1ZdX1+99bWsWyqoZlbR0ZrCgryQJb+crA1iDorRLoVgbCVUb86urlAmRFmc8HSpKkhgxnkopG/i2lg/pVdPj5K6trsxBXVd0w5FXlBb288qVV1XlBcGXQe+ejSha827BsRRtH/cpyQbjuFs0GI4C5Ub0Go3y5kLi62z1750JkbyeHkSSpWzKcSVprVOQmLlmX8g4/d01tqh+pazDiVzcKWJVf1igI5oXDD5dX89YHKxoGw6qaNk8WU1FaQu/ykmwkNC+49SmvC34Nfzbe37vBsXXPAq68fbR3hRPASJLU0QxnktQBSkuC/r3KOuV5wZQSy6tqV96iWR/0sgC3vKo2G/3Lhbxlldn28tztnnX76rbfX5bd8rm8wb62v0KiJPfMX5+K0gbhrX59lX0luYlcVm7XB8NmgqIvlJckrU0MZ5JU5CJWPgu4fiddo7Y2sbx6ZZjLglx+6MuVVTU+pmFZ3fp7SytZWL+dvVJiaWU17ZkQtKKsZOVtn02Et2y9pIkRv2x/37yyullE6179UFFWQnlpSf1rJypKSygxDEqSCsRwJkmipCRyz6913n8WUkpU1aRVgt3SyuZH+eq369drGwTFxR9VrnLsmrxIHrLnAfNDW0Vp3rv86spLV25XlJXQqzTv+Ebhr8GxbVnPKysrCW8hlaS1gOFMktQlIoKKsiz4rNun45/7q1NTmxqM8jUe8auszl7tUJl7fUP9en5Zbruqprb+NQ/5x1TV1LJ0aXW2L69eVU1+3TV/b2CdCCgvzUJg4wCYHwp75YW68kaBr1dZySqhM79eU6GzorSEzdfvy4Denff7kiStZDiTJPUopSVBv15l9CvQ+wLr1L0aorKmlqpGIa4u1DUur8yFwaomgmKzoTJvfcmK6obhMu+cdWVtnVzmuiPGstcnN+ycL0mS1IDhTJKkTpD/aohikVKiujY1CH/1QbFRIFyR2x6x8XqFbrYkrTUMZ5IkrSUigvLSoLy0hL4d/ypBSdIa8i2lkiRJklQEDGeSJEmSVAQMZ5IkSZJUBAxnkiRJklQEDGeSJEmSVAQMZ5IkSZJUBAxnkiRJklQEDGeSJEmSVAQMZ5IkSZJUBAxnkiRJklQEDGeSJEmSVAQMZ5IkSZJUBAxnkiRJklQEDGeSJEmSVAQMZ5IkSZJUBFoMZxGxaURMi4h5EfFMRJyaKx8UEfdExIu5nwM7v7mSJEmS1DO1ZuSsGvhOSmk7YGfgpIjYHjgTuC+ltA1wX25bkiRJktQOLYazlNLClNLjufUPgXnAxsABwPW5w64HDuykNkqSJElSj9emZ84iYgtgR2AGsGFKaSFkAQ4Y0uGtkyRJkqS1RKvDWUT0B24GTkspfdCGesdFxOyImL1o0aL2tFGSJEmSerxWhbOIKCcLZjeklG7JFb8ZERvl9m8EvNVU3ZTS1JTS2JTS2MGDB3dEmyVJkiSpx2nNbI0BXAvMSyn9PG/XHcCk3Pok4PaOb54kSZIkrR3KWnHMbsBhwNMRMTdX9kPgfODGiDgaeBU4qFNaKEmSJElrgRbDWUrpYSCa2T2xY5sjSZIkSWunNs3WKEmSJEnqHIYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAoYzSZIkSSoChjNJkiRJKgKGM0mSJEkqAi2Gs4i4LiLeioh/5ZWNiojHImJuRMyOiJ06t5mSJEmS1LO1ZuTsN8A+jcouBM5JKY0C/ie3LUmSJElqpxbDWUrpQeCdxsXAOrn1dYE3OrhdkiRJkrRWKWtnvdOAuyLiYrKAt2tzB0bEccBxAJtttlk7LydJkiRJPVt7JwQ5ATg9pbQpcDpwbXMHppSmppTGppTGDh48uJ2XkyRJkqSerb3hbBJwS279JsAJQSRJkiRpDbQ3nL0B7Jlb3wt4sWOaI0mSJElrpxafOYuIPwDjgQ0iYgFwNnAscFlElAHLyT1TJkmSJElqnxbDWUrpa83sGtPBbZEkSZKktVZ7b2uUJEmSJHUgw5kkSZIkFQHDmSRJkiQVAcOZJEmSJBUBw5kkSZIkFQHDmSRJkiQVAcOZJEmSJBUBw5kkSZIkFQHDmSRJkiQVAcOZJEmSJBUBw5kkSZIkFQHDmSRJkiQVAcOZJEmSJBUBw5kkSZIkFQHDmSRJkiQVAcOZJEmSJBUBw5kkSZIkFQHDmSRJkiQVAcOZJEmSJBUBw5kkSZIkFQHDmSRJkiQVAcOZJEmSJBUBw5kkSZIkFQHDmSRJkiQVAcOZJEmSJBUBw5kkSZIkFQHDmSRJkiQVgRbDWURcFxFvRcS/GpV/KyKej4hnIuLCzmuiJEmSJPV8rRk5+w2wT35BREwADgBGppSGARd3fNMkSZIkae3RYjhLKT0IvNOo+ATg/JTSitwxb3VC2yRJkiRprdHeZ862BXaPiBkR8UBEjOvIRkmSJEnS2qZsDeoNBHYGxgE3RsRWKaXU+MCIOA44DmCzzTZrbzslSZIkqUdr78jZAuCWlJkJ1AIbNHVgSmlqSmlsSmns4MGD29tOSZIkSerR2hvObgP2AoiIbYEK4O0OapMkSZIkrXVavK0xIv4AjAc2iIgFwNnAdcB1uen1K4FJTd3SKEmSJElqnRbDWUrpa83s+kYHt0WSJEmS1lrtva1RkiRJktSBDGeSJEmSVAQMZ5IkSZJUBAxnkiRJklQEDGeSJEmSVAQMZ5IkSZJUBAxnkiRJklQEDGeSJEmSVAQMZ5IkSZJUBAxnkiRJklQEDGeSJEmSVAQMZ5IkSZJUBAxnkiRJklQEIqXUdReLWAS80mUXbL0NgLcL3Qj1WPYvdSb7lzqbfUydyf6lzlSs/WvzlNLgpnZ0aTgrVhExO6U0ttDtUM9k/1Jnsn+ps9nH1JnsX+pM3bF/eVujJEmSJBUBw5kkSZIkFQHDWWZqoRugHs3+pc5k/1Jns4+pM9m/1Jm6Xf/ymTNJkiRJKgKOnEmSJElSEehW4SwiekfEzIh4MiKeiYhz8vYNioh7IuLF3M+BTdQfFRGP5uo+FRGH5O3bMiJm5Or/KSIqmmnDpNwxL0bEpLbWV/Eqkv71j4h4LyLubFRu/+rmCt2/OqJ/qngVQf/aPCLmRMTc3DmOb0t9FbdC96+8Y9eJiNcj4hftqa/iVQx9LCJqcv+GzY2IO9pav8OklLrNAgTQP7deDswAds5tXwicmVs/E7igifrbAtvk1ocCC4H1cts3Al/NrU8BTmii/iDgP7mfA3PrA1tb36W4l0L3r9y+icD+wJ2Nyu1f3XwpdP/qiP7pUrxLEfSvCqBXbr0/MB8Yav/qGUuh+1feeS4D/g/4RV6Z/asHLMXQx4AlzZR3aR8r+C9jDX6JfYHHgU/ltp8HNsqtbwQ834pzPAlsk+sQbwNlufJdgLuaOP5rwFV521flylpV36X7LIXoX3n1xpMXzuxfPW8pZP/qqPouxbsUun8B6wOv5v5Asn/1sKVQ/QsYA/wROIJcOLN/9cylgH1slXBWiD7WrW5rBIiI0oiYC7wF3JNSmpHbtWFKaSFA7ueQFs6zE9n/0/cS2X9I3kspVed2LwA2zh03NiKuyZVvDLyWd5q645qtr+6lwP2rOfavHqJY+ldr66t7KXT/iohNI+Ipsv9OXpBSemN19dW9FLJ/RUQJ8DPge41OZ//qQQr9bxjQOyJmR8RjEXFgrqzL+1i3C2cppZqU0ihgE2CniBje1nNExEbA74AjU0q1ZKl4lUvlrjc7pXRMXdVmjmu2vrqXAvevZk/ZXH11L8XQv9pSX91LoftXSum1lNJIYGtgUkRsuLr66l4K3L9OBP6WUnqt0bH2rx6k0P+GAZullMYCXwcujYiPr65+Z+l24axOSuk9YDqwT67ozdwvpO4X81ZT9SJiHeCvwH+nlB7LFb8NrBcRZbntTYA3mqi+ANg0b7vuuNbWVzdRoP7VHPtXD1Oo/tVJ/VNFptD/fuVGzJ4Bdm9PfRW3AvWvXYCTI2I+cDFweESc34b66kYK9W9Y7t8uUkr/yV1/x7bU7yjdKpxFxOCIWC+33gf4DPBcbvcdwKTc+iTg9ibqVwC3Ar9NKd1UV56ym0inAV9ZXX3gLuBzETEwN1PM58juO21tfRWxIuhfTbJ/9QyF7l+d1T9VHIqgf22Suy65/z7uRvZciP2rByh0/0opHZpS2iyltAXw3dx5zrR/9RyF7mO5v+175dY3IPs37NmC9LG2PKBW6AUYCTwBPAX8C/ifvH3rA/cBL+Z+Dmqi/jeAKmBu3jIqt28rYCbwb+AmVs46NRa4Ju8cR+WO+TfZkCmrq+/SfZYi6V8PAYuAZWQjtXvbv3rGUuj+1Z76Lt1nKYL+9dnctZ/M/Twu79z2r26+FLp/NTrXETScrdH+1QOWQvcxYFfg6dy/YU8DRxeqj0XuopIkSZKkAupWtzVKkiRJUk9lOJMkSZKkImA4kyRJkqQiYDiTJEmSpCJgOJMkSZKkImA4kyRJkqQiYDiTJEmSpCJgOJMkSZKkIvD/ATvMEK6NgptEAAAAAElFTkSuQmCC\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "plt.figure(figsize = (15, 5))\n", - "plt.plot(df_simulation.index, df_simulation['values'], label = 'Simulated Temperature')\n", - "plt.plot(df_T_sim_horizon.index, df_T_sim_horizon['values'], label = 'Prediction Horizon', color = 'red', marker = 'x')\n", - "#plt.plot(df.index, df['InsideTemp'], label = 'Inside Temperature')\n", - "#plt.plot(df.index, df['OutsideTemp'], label = 'Outside Temperature')\n", - "plt.title(f'Temperatures step {current_timestamp}')\n", - "plt.legend()\n", - "plt.ylim((15, 30))\n", - "plt.savefig(f\"sim_{current_timestamp}.png\")\n", - "plt.show()\n" - ] - }, - { - "cell_type": "code", - "execution_count": 81, - "metadata": { - "scrolled": true, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Timestamp 1500\n" - ] - }, - { - "ename": "RuntimeError", - "evalue": ".../casadi/core/function_internal.hpp:1257: Input 1 (p) has mismatching shape. Got 15-by-1. Allowed dimensions, in general, are:\n - The input dimension N-by-M (here 14-by-1)\n - A scalar, i.e. 1-by-1\n - M-by-N if N=1 or M=1 (i.e. a transposed vector)\n - N-by-M1 if K*M1=M for some K (argument repeated horizontally)\n - N-by-P*M, indicating evaluation with multiple arguments (P must be a multiple of 1 for consistency with previous inputs)", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 19\u001b[0m )\n\u001b[1;32m 20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreal_p\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlbg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreal_lbg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mubg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreal_ubg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 22\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0mdf_power\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdf_power\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'time'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mcurrent_timestamp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Heat'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'x'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN_horizon\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/usr/lib/python3.9/site-packages/casadi/casadi.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 8504\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8505\u001b[0m \u001b[0;31m# Named inputs -> return dictionary\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 8506\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 8507\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8508\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mbuffer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/usr/lib/python3.9/site-packages/casadi/casadi.py\u001b[0m in \u001b[0;36mcall\u001b[0;34m(self, *args)\u001b[0m\n\u001b[1;32m 7677\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7678\u001b[0m \"\"\"\n\u001b[0;32m-> 7679\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0m_casadi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFunction_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 7680\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7681\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mmapsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m->\u001b[0m \u001b[0;34m\"std::vector< casadi::MX,std::allocator< casadi::MX > >\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mRuntimeError\u001b[0m: .../casadi/core/function_internal.hpp:1257: Input 1 (p) has mismatching shape. Got 15-by-1. Allowed dimensions, in general, are:\n - The input dimension N-by-M (here 14-by-1)\n - A scalar, i.e. 1-by-1\n - M-by-N if N=1 or M=1 (i.e. a transposed vector)\n - N-by-M1 if K*M1=M for some K (argument repeated horizontally)\n - N-by-P*M, indicating evaluation with multiple arguments (P must be a multiple of 1 for consistency with previous inputs)" - ] - } - ], - "source": [ - "for loop_idx in range(100):\n", - " current_timestamp = 1500 + 300*loop_idx\n", - " print(f\"Timestamp {current_timestamp}\")\n", - " \n", - " u_1 = float(df_power.loc[df['time'] == (current_timestamp - 300 * 1), 'Heat'])\n", - " u_2 = float(df_power.loc[df['time'] == (current_timestamp - 300 * 2), 'Heat'])\n", - "\n", - " y_1 = float(df.loc[df['time'] == (current_timestamp - 300 * 1), 'SimulatedTemp'])\n", - " y_2 = float(df.loc[df['time'] == (current_timestamp - 300 * 2), 'SimulatedTemp'])\n", - " y_3 = float(df.loc[df['time'] == (current_timestamp - 300 * 3), 'SimulatedTemp'])\n", - " \n", - " real_x0 = np.array([u_1, u_2, y_1, y_2, y_3])\n", - " iter_idx = (df['time'] >= current_timestamp)\n", - " real_W = df[iter_idx].iloc[:N_horizon, [5, 2]].to_numpy()\n", - "\n", - " real_p = casadi.vertcat(\n", - " casadi.vec(real_W),\n", - " casadi.vec(real_x0)\n", - " )\n", - "\n", - " res = solver(p = real_p, lbg = real_lbg, ubg = real_ubg)\n", - " \n", - " df_power.loc[df_power['time'] == current_timestamp, 'Heat'] = res['x'].reshape((N_horizon, -1))[1, 2]\n", - " \n", - " power = np.array(df_power[['time', 'Heat']].dropna())\n", - " eng.workspace['power'] = matlab.double(power.tolist())\n", - " eng.set_param('polydome', 'StopTime', str(current_timestamp + 300), nargout = 0)\n", - " eng.workspace['result'] = eng.sim('polydome')\n", - " \n", - " \n", - " dict_simulation = {}\n", - " dict_simulation['values'] = np.asarray(eng.eval('result.SimulatedTemp.Data')).reshape(-1)\n", - " dict_simulation['time'] = np.asarray(eng.eval('result.SimulatedTemp.Time')).reshape(-1)\n", - " \n", - " df_simulation = pd.DataFrame(dict_simulation)\n", - " #df_simulation['time'] = df_simulation['time'].astype(int)\n", - " df_simulation.set_index('time', inplace = True, drop = True)\n", - " \n", - " df_simulation['timestamp'] = df.index[0] + df_simulation.index.map(lambda x: pd.Timedelta(seconds = x))\n", - " df_simulation = df_simulation.reset_index().set_index('timestamp')\n", - " df_resampled_5 = df_simulation['values'].resample('5min').mean().pad()\n", - " df_simulation = pd.concat([df['time'], df_resampled_5], axis = 1)\n", - " \n", - " df.loc[:, 'SimulatedTemp'] = df_simulation['values']\n", - " T_sim_horizon = np.array(gpr(res['x'].reshape((N_horizon, -1)).T))\n", - " simul_idx = (df_simulation['time'] >= current_timestamp) & (df_simulation['time'] < (current_timestamp + N_horizon * 300))\n", - " \n", - " \n", - " df_T_sim_horizon = df_simulation[simul_idx].copy()\n", - " df_T_sim_horizon.loc[:, 'values'] = T_sim_horizon.reshape((-1, ))\n", - " \n", - " plt.figure(figsize = (15, 5))\n", - " plt.plot(df_simulation.index, df_simulation['values'], label = 'Simulated Temperature')\n", - " plt.plot(df_T_sim_horizon.index, df_T_sim_horizon['values'], label = 'Prediction Horizon', color = 'red', marker = 'x')\n", - " #plt.plot(df.index, df['InsideTemp'], label = 'Inside Temperature')\n", - " #plt.plot(df.index, df['OutsideTemp'], label = 'Outside Temperature')\n", - " plt.title(f'Temperatures step {current_timestamp}')\n", - " plt.ylim((15, 30))\n", - " plt.legend()\n", - " plt.savefig(f\"sim_{current_timestamp}.png\")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.6" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/Notebooks/70_Server_result_analysis.ipynb b/Notebooks/70_Server_result_analysis.ipynb index 53dce19..3c3aa57 100644 --- a/Notebooks/70_Server_result_analysis.ipynb +++ b/Notebooks/70_Server_result_analysis.ipynb @@ -169,6 +169,14 @@ "### Compute the time index" ] }, + { + "cell_type": "markdown", + "id": "0a39bb25-f673-488b-99a4-57ee0d789f54", + "metadata": {}, + "source": [ + "The time index is computed by adding the elapsed time (`sample nr.` * `Tsample`) to the dataset start time. Since the CARNOT weather set represents the year 2010, the starting time is taken as 2010-01-01 at midnight." + ] + }, { "cell_type": "code", "execution_count": 14, @@ -278,6 +286,14 @@ "### Get reference temperature" ] }, + { + "cell_type": "markdown", + "id": "8b53a7c0-966f-4c75-8857-1a96ec3135b0", + "metadata": {}, + "source": [ + "The reference temperature is computed according to the SIA norm, as a function of the last 48h of outside temperature (at a `Tsample` of 15min this comes out to 2 * 96 sample points)" + ] + }, { "cell_type": "code", "execution_count": 16, @@ -292,6 +308,14 @@ "df_tref = df_tref.shift(1) # The reference at time t is computed using info up to t-1" ] }, + { + "cell_type": "markdown", + "id": "5d2e03c9-044a-454b-8a6b-c1a2a6adf6cc", + "metadata": {}, + "source": [ + "Compute mean value and standard deviation of tracking error:" + ] + }, { "cell_type": "code", "execution_count": 17, @@ -358,6 +382,14 @@ "}" ] }, + { + "cell_type": "markdown", + "id": "b0fc6fbc-ed75-4724-a5cb-8896a89f5039", + "metadata": {}, + "source": [ + "Scale the dataset and compute the autoregressive inputs that are passed to the GP:" + ] + }, { "cell_type": "code", "execution_count": 20, @@ -548,6 +580,14 @@ "df_output = df_gpr[dict_cols['y'][1]]" ] }, + { + "cell_type": "markdown", + "id": "ddbfd08f-7d5d-4f2e-b43e-a282a0503931", + "metadata": {}, + "source": [ + "Load the only trained model in the GP case, and the first trained model in the SVGP case:" + ] + }, { "cell_type": "code", "execution_count": 23, @@ -571,6 +611,14 @@ " m = model" ] }, + { + "cell_type": "markdown", + "id": "c18135f0-06d8-44e3-b86c-4a84b2508f79", + "metadata": {}, + "source": [ + "Plot the multistep prediction performance for 25 consecutive points:" + ] + }, { "cell_type": "code", "execution_count": 25, @@ -1091,6 +1139,14 @@ "output_notebook()" ] }, + { + "cell_type": "markdown", + "id": "3ec059bc-9ce6-4271-b5cd-5003d1330a0b", + "metadata": {}, + "source": [ + "### Plot evolution of reference/measured temperature" + ] + }, { "cell_type": "code", "execution_count": 22, @@ -1393,6 +1449,14 @@ "plt.savefig(f\"../Thesis/Plots/{sim_id}_abserr.pdf\", bbox_inches='tight')" ] }, + { + "cell_type": "markdown", + "id": "4b0595d2-59d7-4536-b8a3-8c8aced19dc4", + "metadata": {}, + "source": [ + "### Plot evolution of hyperparameters" + ] + }, { "cell_type": "code", "execution_count": 29,