Getting accustomed with MPC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Getting accustomed with MPC\n",
    "author: Vasken Dermardiros\n",
    "date: 2019-05-09\n",
    "Confidential. For internal use only.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mise en Situation\n",
    "An advanced alien civilization has landed on the blue pebble in the vastness of space we call Earth, which literaly means dirt, and we, scientists, have presented them a token of humble peace: a net-zero energy building!\n",
    "\n",
    "Alas, EnergyPlus got a new update that made the model of the building completely unuseable. Years of work, wasted. A moment of silence please. Let us grieve.\n",
    "\n",
    "But hoorah! A young engineer, full of hope and full of intelimaginences exclaims, \"why the sad faces my utterly depressed PhD friends?\" The PhD students who worked tiredlessly replied, \"we have spent an eternity measuring the thermal conductance of the tile floors in this building to have the most accurate model ever, my research compatriot there has even proclaimed his love of the building by marrying it in holy matrimony, but with the recent update to version 9.0 and the introduction of JDF formating, the IDF file no longer loads and we can no longer use the predictions to optimize or quasi-heuristico-polonialy optimize the building.\" #sad\n",
    "\n",
    "\"My friends, I have no idea what the hell any of that means,\" asserted the young engineer, \"but have you not been monitoring this institution?\" \"But indeed we have,\" replied the researcher befuddled, \"do you really think we are here doing Mickey Mouse stuff?,\" slightly in annoyance. \"I do not wish to strain your delicate pysche, I merely believe you can do away with the complex modeling and use a simplified low-order approximation, like an RC model, but I will need some help with fitting the parameters.\"\n",
    "\n",
    "This is where you come in! Help the engineer fit an RC model to then present it to the aliens. Did you forget the aliens? Is there anybody out there? Just nod if you can hear me, is there anyone home?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make jupyter notebook full-width\n",
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Default packages\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "import matplotlib as mpl\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from datetime import datetime, timedelta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save and display settings\n",
    "mpl.rc('figure', figsize=(4, 4))\n",
    "# mpl.rcParams.update({'figure.subplot.left': 0.20, 'figure.subplot.bottom': 0.15})\n",
    "mpl.rcParams.update({'figure.subplot.left': 0.15, 'figure.subplot.bottom': 0.11})\n",
    "# figext, figdpi = \".svg\", 600\n",
    "figext, figdpi = \".png\", 300\n",
    "pd.set_option('display.max_columns', 8)\n",
    "pd.set_option('display.width', 200)\n",
    "RED, BLUE, PINK = \"#f80d1c\", \"#00bdff\", \"#fb3983\"\n",
    "\n",
    "np.set_printoptions(edgeitems=30, linewidth=100000,)\n",
    "#     formatter=dict(float=lambda x: \"%.3g\" % x))\n",
    "np.set_printoptions(suppress=True) # suppress scientific notation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Exploratory Data Analysis\n",
    "After importing all the necessary packages, we will begin by reading the csv file with the data and doing some exploratory data analysis (pretty much making graphs and looking at them)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data from csv file using pandas\n",
    "df = pd.read_csv('intro_data.csv', parse_dates=True, infer_datetime_format=True, index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# So what's in the dataframe?\n",
    "df.head() # print the first few rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Check if there is missing or corrupt data\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some statistics perhaps?\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the data!\n",
    "# df.plot()\n",
    "# plt.show()\n",
    "\n",
    "# ... actually don't because there's too much of it, what about something a little more useful?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Like seeing what are correlations between the data using seaborn + the distribution of the data in the diagonal\n",
    "# Read through the documentations to learn all the options\n",
    "sns.pairplot(df, diag_kind=\"kde\", kind='reg', plot_kws={'line_kws':{'color':RED}, 'scatter_kws': {'color':BLUE,'s':1,'alpha': 0.1}})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the previous plot shows what exactly? How one variable is related to another for that time. Could you think of a way to do this with lagged values? What would you be trying to learn from that?\n",
    "\n",
    "Back to the previous plots, the red lines are best-fit regression lines, do you agree with them?\n",
    "\n",
    "On the diagonal, we see the density of the parameters. Look at GHI, it appears it's mostly 0 most of the time, why? Consumption has 2 bumps. Heating Total has a few bumps, what do you think those represent?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Do a scatter plot of production versus GHI and fit a trendline; what can you observe?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train an RC Model\n",
    "Draw a 1-R 1-C model of this building. Which variables go where? Heating Total and Cooling Total is in kW (average for that timestep), Consumption and Production are electrical and in kWh, GHI is global horizontal irradiance in W/m^2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Draw RC model on paper, write down the equation (1) here in the form:\n",
    "# T_int^{t+1} = b_0*T_int^{t} + b_1*T_ext^{t} + ...\n",
    "# where you replace b_0 with what they correspond to, like U/dt\n",
    "# use conductance (U=1/R) instead of resistance to avoid dividing by 0 (if it happens)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before finding the best parameters for the data, it is good practice to split the data into a training and a testing set. The training set is what the model learns on, the testing set is what we use to observe it's performance on an out-of-sample situation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['T_int_plusone'] = df['T_int'].shift(-1) # search on google \"pandas shift\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split data to train/test sets\n",
    "df_train = df['2018-01-01':'2018-02-14'].copy()\n",
    "df_test = df['2018-02-14':'2018-02-28'].copy()\n",
    "df_test.interpolate(method='linear', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Complete this to match your equation\n",
    "# Inputs used in model\n",
    "X_columns = ['T_int','T_ext',...]\n",
    "# Output\n",
    "y_column = ['T_int_plusone']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Transfer the equation you wrote from above here where X is a matrix containing the training data\n",
    "dt = 15*60. # s\n",
    "def model(X,a,b,c,d,e,f...): # returns T_int{t+1}\n",
    "    return a*X[0,] + b*X[1,] + ... # you will copy this equation in the MPC part"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Fit your model using the training set!\n",
    "X_data_curvefit = df_train[X_columns].values.T\n",
    "y_data_curvefit = df_train[y_column].values.flatten()\n",
    "\n",
    "# Read the documentation for \"curve_fit\" and use it to fit your model coefficients to the data\n",
    "# The function returns the fitted parameters and the covariance\n",
    "popt, pcov = curve_fit(model, X_data_curvefit, y_data_curvefit)\n",
    "a,b,c,d,e,f... = popt # <-- update!\n",
    "perr = np.sqrt(np.diag(pcov)) # covariance matrix diagonal is the variance\n",
    "print(popt, perr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make predictions!\n",
    "y_pred = model(df_test[X_columns].values.T, popt[0], popt[1], popt[2], popt[3], ...) # <-- update!\n",
    "y_test = df_test[y_column].values.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can change the start and end of the plot\n",
    "# I prefer using Bokeh since you can interact with the plot\n",
    "start, end = 0, 200\n",
    "y_test_plot = y_test[start:end]\n",
    "y_pred_plot = y_pred[start:end]\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.plot(y_test_plot, color='k', label='Truth')\n",
    "plt.plot(y_pred_plot, color=PINK, label='Predicted')\n",
    "plt.legend(frameon=False)\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(y_test_plot-y_pred_plot, label='Residual')\n",
    "plt.legend(frameon=False)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Augmented Dickey-Fuller test to see if residuals are stationary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "residuals = y_test - y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adf = sm.tsa.stattools.adfuller(residuals)\n",
    "print('if p-value <= 0.05: Reject the null hypothesis (H0), the data does not have a unit root and is stationary.')\n",
    "print('ADF Statistic: %f for y_%i, p-value: %f' %(adf[0], 1, adf[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Durbin-Watson test to capture autocorrelation in residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dw = sm.stats.stattools.durbin_watson(residuals)\n",
    "# The Durbin Watson statistic is a number that tests for autocorrelation in the \n",
    "# residuals from a statistical regression analysis. The Durbin-Watson statistic \n",
    "# is always between 0 and 4. A value of 2 means that there is no autocorrelation \n",
    "# in the sample. Values from 0 to less than 2 indicate positive autocorrelation \n",
    "# and values from more than 2 to 4 indicate negative autocorrelation.\n",
    "print(dw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] What's all that about?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Create a linear model for the electrical production using GHI; what would be a better approach?\n",
    "# Pretend you had also readings for diffuse horizontal irradiance (DHI) and direct normal irradiance (DNI)\n",
    "# and you know the number of panels at what orientation and tilt\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization\n",
    "Read Boyd's book on convex optimization. I won't be going anywhere. Take your time. https://web.stanford.edu/~boyd/cvxbook/\n",
    "\n",
    "Done? Good! I'm a big fan of cvxpy and that's what we will be using. Scipy isn't great here. Do your own bloody research.\n",
    "\n",
    "CVXPY examples ('Control' example is MPC): https://www.cvxpy.org/examples/index.html\n",
    "\n",
    "You probably didn't read the whole book or followed his online course *because you are lazy and want to be spoon-fed everything.*\n",
    "\n",
    "In a nutshell, convex optimization is a special type of optimization where the cost function and all constraints are convex. Everything being convex has a very nice property that when you find an optimal, it is *the* optimal. Period. What's convex? if the 2nd-derivative of the function or contraints is non-negative. Linear programming is a simpler case of convex optimization. Generally, convex optimization can be performed extremely fast for the number of variables we tend to have in our field. You can solve for a 24h prediction horizon with a 15 minute time step in <<1 second. (The actual optimization takes milli-seconds, most of the time is spent translating the optimization problem to the form needed by the actual solver.)\n",
    "\n",
    "Linear functions are convex. Quadratic functions are convex. Log is convex. Adding two convex functions is convex. Multiplying two convex functions is convex. Anyway. Do your research and really understand this. There are cases where your function is not convex but can be approximated convex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import cvxpy as cp\n",
    "print(cp.__version__) # make sure it's >1.0\n",
    "print(cp.installed_solvers()) # which solvers are available; you can also use commercial solvers if you have a license key"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPC Example from CVXPY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may wish to look at the other examples from the website, but let's focus here on the MPC example because we will do something similar with the model we trained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Replicate the 'control' example from the cvxpy website:\n",
    "# https://colab.research.google.com/github/cvxgrp/cvx_short_course/blob/master/intro/control.ipynb\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Change the sum of squares cost to a norm-1, how does the behaviour of the system change? What is norm-1 penalizing?\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Keep the sum of squares cost. Add a constraint to reduce variable \"u\" from changing more than a\n",
    "# certain amount from one timestep to the next, eg. |u^{t+1} - u^{t}| <= some value\n",
    "# What happens?\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPC Using Physics-Based Approach, 1st-order model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I hope it's safe to assume that you included heating_total in your model and this is what we will control. Cooling_total is practically 0 for the whole dataset.\n",
    "\n",
    "Traditionally, the variable 'x' represents the system variables (temperature, pressure) which are properties that are measured. The variable 'u' represents controllable variables (heater output, fan state, valve position) and is what MPC will change to attempt to minimize the cost function.\n",
    "\n",
    "Everything has to be expressed in the cost function and constraints. The MPC doesn't know a battery can't have a -24% charge level, you will need to constrain it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick the day from where to test it on\n",
    "df_opti = df_test['2018-02-20'::]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Modify my MPC code to include your calibrated model\n",
    "nx = len(X_columns)-1 # dynamics variables: 'T_ext', 'T_int', ...\n",
    "nu = 1    # controllable variables: 'heating_total'\n",
    "st = 4    # steps per hour\n",
    "ph = 24*st # 1 day\n",
    "ct = 1    # time step where to start optimization {>0}\n",
    "\n",
    "# For simplicity, I have these here like this; it's easier to have these within the df_test dataframe\n",
    "heating_min, heating_max = 0., 90. # u0_min, u0_max\n",
    "Tint_min, Tint_max = 18.*np.ones(ph), 26.*np.ones(ph)\n",
    "Tint_min[9*st:21*st] = 20.\n",
    "\n",
    "x = cp.Variable((nx, ph))\n",
    "u = cp.Variable((nu, ph))\n",
    "c = np.ones((nu, ph)) # unit cost\n",
    "# b = cp.Parameter(nx+nu)\n",
    "\n",
    "obj = cp.Minimize(c*u.T) # vectorized form\n",
    "\n",
    "constraints = cons = []\n",
    "cons.append(x[???,:] >= Tint_min,) # replace ??? with the index corresponding to T_int in your X_columns variable\n",
    "cons.append(x[???,:] <= Tint_max,) # replace ??? with the index corresponding to T_int in your X_columns variable\n",
    "cons.append(u >= heating_min,)\n",
    "cons.append(u <= heating_max,)\n",
    "\n",
    "cons.append(x[???,:] == df_opti['T_ext'][ct:ct+ph],) # replace ??? with the index corresponding to T_ext in your X_columns variable, if you used it\n",
    "cons.append(x[???,:] == df_opti['GHI'][ct:ct+ph],) # replace ??? with the index corresponding to GHI in your X_columns variable, if you used it\n",
    "cons.append(x[???,:] == df_opti[???][ct:ct+ph],) # other; you can make your own too; internal gains ?= cons - heating_total/COP, up to you\n",
    "T_int_0 = df_opti['T_int'][ct-1] # initial temperature\n",
    "\n",
    "# NOTE: x[1,t] is actually @ t=t+1 meaning at the next step, \n",
    "# otherwise we need to calculate Q_ph+1 which is useless\n",
    "# TODO: change the following section to match your calibrated equation format\n",
    "cons.append(x[1,0] == b[0]*u[0,0] + b[1]*x[0,0] + b[2]*T_int_0 + b[3]*x[2,0] + b[4]*x[3,0],)\n",
    "for t in range(1,ph):\n",
    "    cons.append(x[1,t] == b[0]*u[0,t] + b[1]*x[0,t] + b[2]*x[1,t-1] + b[3]*x[2,t] + b[4]*x[3,t],)\n",
    "constraints.extend(cons)\n",
    "\n",
    "prob = cp.Problem(obj, constraints)\n",
    "\n",
    "timer = time.time()\n",
    "# b.value = b_grey # grey-box corresponding parameter values\n",
    "prob.solve(verbose=True, solver=cp.ECOS)\n",
    "Tint_MPC = x.value[???,:] # replace ??? with the index corresponding to T_int in your X_columns variable\n",
    "Q_MPC = u.value[0,:]\n",
    "print(\"Elapsed time: %.5f seconds\" % (time.time()-timer))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tickhrs = 3\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "\n",
    "plt.plot(Tint_min,'--', c=RED, label='Heating SP')\n",
    "plt.plot(Tint_max,'--',c=BLUE, label='Cooling SP')\n",
    "plt.plot(df_opti['T_int'][ct:ct+ph].values,'k--', label='Interior, Current')\n",
    "plt.plot(Tint_MPC,'k', label='Interior, MPC')\n",
    "\n",
    "plt.xlim([0,ph])\n",
    "ax.set_xticks(np.arange(0,ph+1,4*tickhrs))\n",
    "ax.set_xticklabels(np.arange(0,24+1,tickhrs))\n",
    "ax.set_yticks(np.arange(18,27,2))\n",
    "# plt.legend(['Heating SP','Cooling SP'], loc='upper left')\n",
    "plt.legend(loc='upper left')\n",
    "plt.xlabel('Time, h')\n",
    "plt.ylabel('Temperature, degC')\n",
    "# plt.title('Setpoint Temperatures')\n",
    "plt.grid()\n",
    "plt.savefig('figures/MPC_better_than_normal'+figext, dpi=figdpi)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "\n",
    "plt.plot(df_opti['heating_total'][ct:ct+ph].values,'k--', label='Heating, Current')\n",
    "plt.plot(Q_MPC,'k', label='Heating, MPC')\n",
    "\n",
    "plt.xlim([0,ph])\n",
    "ax.set_xticks(np.arange(0,ph+1,4*tickhrs))\n",
    "ax.set_xticklabels(np.arange(0,24+1,tickhrs))\n",
    "ax.set_ylim([0,None])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Time, h')\n",
    "plt.ylabel('Heating, kW')\n",
    "# plt.title('Setpoint Temperatures')\n",
    "plt.grid()\n",
    "plt.savefig('figures/MPC_better_than_normal_Q'+figext, dpi=figdpi)\n",
    "plt.show()\n",
    "\n",
    "print(\"Energy consumed by current strategy: %.0f kWh\" % (15/60*df_opti['heating_total'][ct:ct+ph].sum()))\n",
    "print(\"Energy consumed by MPC strategy: %.0f kWh\" % (15/60*Q_MPC.sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Follow Up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Congrats! You just did 1 step of MPC for a 24h prediction horizon! What you would actually do from here is to apply the first action and 15 minutes later, run the optimization again and apply the first action then too. You could use a \"warm start\" approach in the optimization and use the previous optimal results as the initial guess to speed conversion. You could also use the parameter approach like where I have variable 'b' be weights that I can chose and run the optimization multiple times (same 24 hour period) and have the 'b' weights vary a little.\n",
    "\n",
    "Limitations\n",
    "+ T_ext, GHI, consumption, production are given for future timesteps with perfect certainty\n",
    "+ MPC may decide to apply actions that have never happened and the model's response might be far from reality\n",
    "+ I'm too sleepy to think of others"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Change the sum of energy cost to penalize the peak power, how does the behaviour of the system change?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Change the sum of energy cost to have varying cost at different times, how does the behaviour of the system change?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Change the sum of energy cost to replicate Hydro-Quebec Rate M, how does the behaviour of the system change?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Did you include the electrical consumption of the building? Did you make sure you're not double counting the HVAC system?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Change the sum of energy cost to replicate Hydro-Quebec Rate G, how does the behaviour of the system change?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Change the sum of energy cost to minimize grid interaction. Include the consumption and production \n",
    "# (was the production fault free?)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Change the start date and increase the prediction horizon to 2 days to have a sunny day followed by a cloudy day during the week\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [EXERCISE] Add in an ideal battery (fix a max kWh charge level and a max kW power output) and minimize grid interaction\n",
    "# How is the battery being used? [use the sunny day then cloudy day period]\n",
    "# Vary the initial state of charge.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you did all this and you feel confident about yourself then tap yourself on the shoulder because you are now initiated! As a future project, you may consider having a more detailed RC model, using an ARX model, using PCM as a storage solution, using Canmeteo predictions, adding noise to T_ext to make it more realistic, adapting the fitted RC model parameters at every timestep (fixed when running predictions). The world is now yours my young padawans. May the force be with you.\n",
    "\n",
    ".  \n",
    ".  \n",
    "...  \n",
    "DUN DUN DUN DUNNN DUNNNNNNNNNNN DUN DUN DUN DUNNNNNNNNNNN DUN  \n",
    "DUN DUN DUN DUNNNNNNNNNNNN DUN  \n",
    "DUN DUN DUN DUNNNNNNN  \n",
    "D.DUN DUNNNNNNNNN DUN DUN DUN DUNNNNNNNNNNNN DUN  \n",
    "DUN DUN DUN DUNNNNNNNNN DUN  \n",
    "DUN DUN DUN DUNNNNNNNN  \n",
    "Dididoooooooo do do do do dididoo didoodoo  \n",
    "dooooooo do do do do DOO DIDOODOO  \n",
    "do do dooooooo doo doo didoodoo didoodoo  \n",
    "DUN DUNDUN DUNDUN DUNDUN DUNDUN DUNDUNNNNNNNNNNNNNNNNNNNNNN  \n",
    "DUNDIDIDUN DUNDIDIDUNNNNNNN DUN DUN DUN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}