Adding Building Energy Statistical Models

building_energy_statistical_models/12-The-pySIP-library.ipynb

@@ -0,0 +1,111 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "d8558e9f",
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "ModuleNotFoundError",
+     "evalue": "No module named 'pysip'",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[1], line 4\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mpandas\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpd\u001b[39;00m\n\u001b[1;32m      2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[0;32m----> 4\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mpysip\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mstatespace\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m TwTi_RoRi\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mpysip\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mregressors\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m FreqRegressor \u001b[38;5;28;01mas\u001b[39;00m Regressor\n\u001b[1;32m      7\u001b[0m \u001b[38;5;66;03m# Reading the data\u001b[39;00m\n",
+      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'pysip'"
+     ]
+    }
+   ],
+   "source": [
+    "import pandas as pd\n",
+    "import numpy as np\n",
+    "\n",
+    "from pysip.statespace import TwTi_RoRi\n",
+    "from pysip.regressors import FreqRegressor as Regressor\n",
+    "\n",
+    "# Reading the data\n",
+    "df = pd.read_csv('data/statespace.csv').set_index('Time')\n",
+    "df.drop(df.index[-1], axis=0, inplace=True)\n",
+    "\n",
+    "# pySIP's fit() method takes a dataframe df as argument\n",
+    "# It also needs to be passed the labels of inputs and outputs in this dataframe\n",
+    "inputs = ['T_ext', 'P_hea']\n",
+    "outputs = 'T_int'\n",
+    "\n",
+    "# The time can be scaled on daily units instead of seconds.\n",
+    "# This brings heat capacities to a range closer to other variables\n",
+    "sT = 3600.0 * 24.0\n",
+    "df.index /= sT\n",
+    "\n",
+    "# Specification of parameters\n",
+    "parameters = [\n",
+    "    dict(name='Ro', value=1.0, scale=0.1, transform='log'),\n",
+    "    dict(name='Ri', value=1.0, scale=0.01, transform='log'),\n",
+    "    dict(name='Cw', value=1.0, scale=1e7 / sT, transform='log'),\n",
+    "    dict(name='Ci', value=1.0, scale=1e6 / sT, transform='log'),\n",
+    "    dict(name='sigw_w', value=1.0, scale=0.01 * sT ** 0.5, transform='log'),\n",
+    "    dict(name='sigw_i', value=1.0, scale=0.01 * sT ** 0.5, transform='log'),\n",
+    "    dict(name='sigv', value=1.0, scale=0.01, transform='log'),\n",
+    "    dict(name='x0_w', value=1.0, scale=25.0, transform='log'),\n",
+    "    dict(name='x0_i', value=26.7),\n",
+    "    dict(name='sigx0_w', value=1.0, transform='fixed'),\n",
+    "    dict(name='sigx0_i', value=1.0, transform='fixed'),\n",
+    "]\n",
+    "\n",
+    "# Definition of the regressor, and model fitting\n",
+    "reg = Regressor(TwTi_RoRi(parameters, hold_order=1))\n",
+    "out = reg.fit(df=df, inputs=inputs, outputs=outputs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2ccc1e37",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Specification of parameters\n",
+    "parameters = [\n",
+    "    dict(name='Ro', value=1.0, scale=0.1, transform='log'),\n",
+    "    dict(name='Ri', value=1.0, scale=0.01, transform='log'),\n",
+    "    dict(name='Cw', value=1.0, scale=1e7 / sT, transform='log'),\n",
+    "    dict(name='Ci', value=1.0, scale=1e6 / sT, transform='log'),\n",
+    "    dict(name='sigw_w', value=1.0, scale=0.01 * sT ** 0.5, transform='log'),\n",
+    "    dict(name='sigw_i', value=1.0, scale=0.01 * sT ** 0.5, transform='log'),\n",
+    "    dict(name='sigv', value=1.0, scale=0.01, transform='log'),\n",
+    "    dict(name='x0_w', value=1.0, scale=25.0, transform='log'),\n",
+    "    dict(name='x0_i', value=26.7),\n",
+    "    dict(name='sigx0_w', value=1.0, transform='fixed'),\n",
+    "    dict(name='sigx0_i', value=1.0, transform='fixed'),\n",
+    "]\n",
+    "\n",
+    "# Definition of the regressor, and model fitting\n",
+    "reg = Regressor(TwTi_RoRi(parameters, hold_order=1))\n",
+    "out = reg.fit(df=df, inputs=inputs, outputs=outputs)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

building_energy_statistical_models/4-Ordinary linear regression.qmd

@@ -0,0 +1,92 @@
+---
+title: "4-Ordinary linear regression"
+format: html
+editor: visual
+---
+
+```{r}
+library(tidyverse)
+library(lubridate)
+
+df <- read_csv("data/linearregression.csv") %>% transform(TIMESTAMP = ymd(TIMESTAMP))
+head(df)
+```
+
+```{r}
+df <- df %>% mutate(tite = ti - te)
+lm1.fit = lm(e_hp ~ tite + 0, data=df)
+summary(lm1.fit)
+```
+
+```{r}
+conf.int = as.data.frame( predict(lm1.fit, interval="confidence") )
+pred.int = as.data.frame( predict(lm1.fit, interval="prediction") )
+
+ggplot() +
+  geom_point(data = df, aes(tite, e_hp)) +
+  geom_line(aes(x=df$tite, y=pred.int$fit)) +
+  geom_ribbon(aes(x=df$tite, ymin=pred.int$lwr, ymax=pred.int$upr), alpha=0.2, fill='blue') +
+  geom_ribbon(aes(x=df$tite, ymin=conf.int$lwr, ymax=conf.int$upr), alpha=0.2, fill='red')
+```
+
+```{r}
+df <- df %>% mutate(titg = ti - tg,
+                    tits = ti - ts,
+                    vtite = wind_speed * (ti-te))
+lm2.fit = lm(e_hp ~ tite + titg + i_sol + tits + vtite + 0, data=df)
+summary(lm2.fit)
+```
+
+```{r}
+lm3.fit = lm(e_hp ~ tite + i_sol + tits + 0, data=df)
+summary(lm3.fit)
+```
+
+```{r}
+par(mfrow=c(2,2))
+plot(lm3.fit)
+```
+
+
+```{r}
+library(rstan)
+rstan_options(auto_write = TRUE)
+rstan_options(threads_per_chain = 1)
+options(mc.cores = parallel::detectCores())
+
+lr_model= "
+data {
+  int<lower=0> N;   // number of data items
+  int<lower=0> K;   // number of predictors
+  matrix[N, K] x;   // predictor matrix
+  vector[N] y;      // outcome vector
+}
+parameters {
+  vector[K] theta;       // coefficients for predictors
+  real<lower=0> sigma;  // error scale
+}
+model {
+  y ~ normal(x * theta, sigma);  // likelihood
+}
+"
+
+model_data <- list(
+  N = nrow(df),
+  K = 3,
+  x = df %>% select(tite, i_sol, tits),
+  y = df$e_hp
+)
+
+fit1 <- stan(
+  model_code = lr_model,    # Stan program
+  data = model_data,        # named list of data
+  chains = 4,               # number of Markov chains
+  warmup = 1000,            # number of warmup iterations per chain
+  iter = 4000,              # total number of iterations per chain
+  cores = 2,                # number of cores (could use one per chain)
+)
+
+print(fit1)
+traceplot(fit1)
+pairs(fit1)
+```