12. Time_Series.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. Time Series\n",
    "\n",
    "[![Python Data Science](https://apmonitor.com/che263/uploads/Begin_Python/DataScience12.png)](https://www.youtube.com/watch?v=yfgE0GheCWY&list=PLLBUgWXdTBDg1Qgmwt4jKtVn9BWh5-zgy \"Python Data Science\")\n",
    "\n",
    "**Time series** data is produced sequentially as new measurements are recorded. Models derived from the data give insight into what happens next. They also show how the system can be changed to achieved a different future outcome. Time series models are a representation of a dynamic system in discrete time. Putting a model into time series form is the basis for many methods in dynamics and control. A **digital twin** is a virtual representation of a process that runs in parallel to the physical system. A time series model can be considered a digital twin in the narrow definition of just specific inputs and outputs included in the model. Below is the time series model with a single input `u` and single output `y` with `k` as an index that refers to the time step.\n",
    "\n",
    "$y_{k+1} = \\sum_{i=1}^{n_a} a_i y_{k-i+1} + \\sum_{i=1}^{n_b} b_i u_{k-i+1}$\n",
    "\n",
    "Time series models are used for identification and control. There is additional information on specific types of time-series and dynamic models such as [ARX (Auto-Regressive eXogenous inputs)](https://apmonitor.com/wiki/index.php/Apps/ARXTimeSeries), [Discrete State Space](https://apmonitor.com/wiki/index.php/Apps/DiscreteStateSpace), and [Continuous State Space](https://apmonitor.com/wiki/index.php/Apps/LinearStateSpace) models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![idea](https://apmonitor.com/che263/uploads/Begin_Python/idea.png)\n",
    "\n",
    "### Predict Time Series Response `y` with Input Step `u`\n",
    "\n",
    "Similar to differential equation models, a time series model may have an input (feature) that changes from an external source such as actively changed by a measurement sensor, a person (manually), or selected by a computer.\n",
    "\n",
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "Calculate the response `y` when the input `u` changes from `0` to `100` at `k=5`. Use model size with $n_a$=3, $n_b$=1, $n_u$=1, and $n_y$=1. The time series model is: \n",
    "\n",
    "$y_{k+1} = a_1 \\, y_k + a_2 \\, y_{k-1} + a_3 \\, y_{k-2} + b_1 \\, u_k$\n",
    "\n",
    "| **Parameter** | **Value** |\n",
    "| ----------- | ----------- |\n",
    "| $a_1$ | 0.6 |\n",
    "| $a_2$ | -0.15 |\n",
    "| $a_3$ | 0.46 |\n",
    "| $b_1$ | 0.08 |\n",
    "\n",
    "The initial condition is $y_0, y_1, y_2 = 0$ and the solution should be calculated until $k=100$. Fill in the equation for the time series in the loop.\n",
    "\n",
    "```python\n",
    "    y[k+1] = y[k]  # fill in time-series equation here\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "n = 101\n",
    "t = np.linspace(0,100,101)\n",
    "u = np.zeros(n); u[5:]=100\n",
    "y = np.zeros(n)\n",
    "\n",
    "a = [0.6,-0.15,0.46]\n",
    "b = [0.08]\n",
    "\n",
    "for i in range(2,n-1):\n",
    "    k = int(t[i])\n",
    "    y[k+1] = y[k]  # fill in time-series equation here\n",
    "    \n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(t,y,t,u)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Time-Series Regression\n",
    "\n",
    "Now that you've simulated a time series model, the next step is to determine the coefficients from data. There is a Gekko function `sysid` that automates the process of system identification of a time series model.\n",
    "\n",
    "![gekko](https://apmonitor.com/che263/uploads/Begin_Python/gekko.png)\n",
    "\n",
    "You can access the help on this function with `help(m.sysid)`. Some of the help is shown below.\n",
    "\n",
    "    y,p,K = sysid(t,u,y,na=1,nb=1,shift='calc',pred='model')\n",
    "\n",
    "    Input:     t = time data\n",
    "               u = input data for the regression\n",
    "               y = output data for the regression   \n",
    "               na   = number of output coefficients (default=1)\n",
    "               nb   = number of input coefficients (default=1)\n",
    "               nk   = input delay steps (default=0)\n",
    "               shift (optional) = \n",
    "                  'none' (no shift)\n",
    "                  'init' (initial pt),\n",
    "                  'mean' (mean center)\n",
    "                  'calc' (calculate c)\n",
    "               pred (option) = \n",
    "                  'model' for output error regression form, implicit solution\n",
    "                  'meas' for ARX regression form, explicit solution\n",
    "                  Using 'model' favors an unbiased model prediction but\n",
    "                    can require more time to compute, especially for large\n",
    "                    data sets\n",
    "                  Using 'meas' computes the coefficients of the time series\n",
    "                    model with an explicit solution\n",
    "               \n",
    "    Output:    returns\n",
    "               ypred (predicted outputs)\n",
    "               p as coefficient dictionary with keys 'a','b','c'\n",
    "               K gain matrix\n",
    "               \n",
    "There are many options such as `pred` as `meas` or `model`. With `meas`, the time-series next step is predicted from prior measurements as in `ARX` form. With `model`, it uses the prior model predictions to predict the next time step. This is also called an Output Error (`OE`) model. You should use `pred=meas` if it is a big data set because it will be much faster to solve.\n",
    "\n",
    "The most important thing to decide is how many coefficients to include in the model by setting `na` and `nb`. It is important to start with small numbers and only add additional coefficients if more accuracy is needed to predict higher-order dynamics. Another factor is `shift` where `init` is preferred if you are starting from steady state conditions. Otherwise `mean` or `calc` are good options to create an unbiased model that does not have offset in the predictions.\n",
    "\n",
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "Change the number of coefficients `na` and `nb` and observe the accuracy of the predictions. Also, set `na=2` and `nb=2` and change `pred=model`.\n",
    "\n",
    "```python\n",
    "na = 2 # output coefficients\n",
    "nb = 2 # input coefficients\n",
    "yp,p,K = m.sysid(t,u,y,na,nb,pred='meas')\n",
    "```\n",
    "\n",
    "How much longer does it take to solve as you increase the number of coefficient or change `pred`? You can time a function with:\n",
    "\n",
    "```python\n",
    "import time\n",
    "start = time.time()\n",
    "### the function\n",
    "print('Elapsed time: ' + str(time.time()-start))\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gekko import GEKKO\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# load data and parse into columns\n",
    "url = 'http://apmonitor.com/pdc/uploads/Main/tclab_data4.txt'\n",
    "data = pd.read_csv(url)\n",
    "t = data['Time']\n",
    "u = data['Q1']\n",
    "y = data['T1']\n",
    "\n",
    "# generate time-series model\n",
    "m = GEKKO(remote=False)\n",
    "\n",
    "# system identification\n",
    "na = 2 # output coefficients\n",
    "nb = 2 # input coefficients\n",
    "yp,p,K = m.sysid(t,u,y,na,nb,pred='meas')\n",
    "\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(t,u)\n",
    "plt.legend([r'$Q_1$ (%)'])\n",
    "plt.ylabel('MV Heater (%)')\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(t,y,'b-',label=r'$T_{1,meas}$')\n",
    "plt.plot(t,yp,'r--',label=r'$T_{1,pred}$')\n",
    "plt.legend(); plt.ylabel('CV Temp (°C)')\n",
    "plt.xlabel('Time (sec)'); plt.savefig('12-sysid.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![idea](https://apmonitor.com/che263/uploads/Begin_Python/idea.png)\n",
    "\n",
    "### Simulate Time-Series\n",
    "\n",
    "There may also be multiple inputs and multiple outputs such as when $n_a$=2, $n_b$=1, $n_u$=2, and $n_y$=2.\n",
    "\n",
    "$y_{1,k+1} = a_{1,1} \\, y_{1,k} + a_{2,1} \\, y_{1,k-1} + b_{1,1} \\, u_{1,k} + b_{1,2} \\, u_{2,k}$\n",
    "\n",
    "$y_{2,k+1} = a_{1,2} \\, y_{2,k} + a_{2,2} \\, y_{2,k-1} + b_{2,1} \\, u_{1,k} + b_{2,2} \\, u_{2,k}$\n",
    "\n",
    "Gekko has the `arx` model that solves time series models once they are identified. It requires a Python dictionary with $A\\in\\mathbb{R}^{n_a \\, \\mathrm{x} \\, n_y}$, $B\\in\\mathbb{R}^{n_y \\, \\mathrm{x} \\, \\left(n_b \\mathrm{x} n_u\\right)}$, and $C\\in\\mathbb{R}^{n_y}$ coefficient matrices. This dictionary is automatically created with the Gekko `sysid` (system identification) function. Below is an example of creating the dictionary manually.\n",
    "\n",
    "```python\n",
    "    # python dictionary\n",
    "    p = {'a':A,'b':B,'c':C}\n",
    "```\n",
    "\n",
    "$A = \\begin{bmatrix}0.36788 & 0.36788 \\\\ 0.223 & -0.136\\end{bmatrix}$\n",
    "$B = \\begin{bmatrix}0.63212 & 0.18964 \\\\ 0.31606 & 1.2642\\end{bmatrix}$\n",
    "$C = \\begin{bmatrix}0 & 0\\end{bmatrix}$\n",
    "\n",
    "[Additional tutorials on Gekko](https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization) show how to solve other types of equations and optimization problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from gekko import GEKKO\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "na = 2 # Number of A coefficients\n",
    "nb = 1 # Number of B coefficients\n",
    "ny = 2 # Number of outputs\n",
    "nu = 2 # Number of inputs\n",
    "\n",
    "# A (na x ny)\n",
    "A = np.array([[0.36788,0.36788],\\\n",
    "              [0.223,-0.136]]) \n",
    "# B (ny x (nb x nu))\n",
    "B1 = np.array([0.63212,0.18964]).T\n",
    "B2 = np.array([0.31606,1.26420]).T\n",
    "B = np.array([[B1],[B2]])\n",
    "\n",
    "C = np.array([0,0])\n",
    "\n",
    "# create parameter dictionary\n",
    "p = {'a':A,'b':B,'c':C}\n",
    "\n",
    "# Create GEKKO model\n",
    "m = GEKKO(remote=False)\n",
    "\n",
    "# Build GEKKO ARX model\n",
    "y,u = m.arx(p)\n",
    "\n",
    "# load inputs\n",
    "tf = 20 # final time\n",
    "u1 = np.zeros(tf+1)\n",
    "u2 = u1.copy()\n",
    "u1[5:] = 3.0\n",
    "u2[10:] = 5.0\n",
    "u[0].value = u1\n",
    "u[1].value = u2\n",
    "\n",
    "# customize names\n",
    "mv1 = u[0]; mv2 = u[1]\n",
    "cv1 = y[0]; cv2 = y[1]\n",
    "\n",
    "# options\n",
    "m.time = np.linspace(0,tf,tf+1)\n",
    "m.options.imode = 4; m.options.nodes = 2\n",
    "\n",
    "# simulate\n",
    "m.solve(disp=False)\n",
    "\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(m.time,mv1.value,'r-',label=r'$MV_1$')\n",
    "plt.plot(m.time,mv2.value,'b--',label=r'$MV_2$')\n",
    "plt.ylabel('MV')\n",
    "plt.legend(loc='best')\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(m.time,cv1.value,'r:',label=r'$CV_1$')\n",
    "plt.plot(m.time,cv2.value,'b.-',label=r'$CV_2$')\n",
    "plt.ylabel('CV'); plt.xlabel('Time (sec)')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### TCLab Activity\n",
    "\n",
    "Run the script to generate data for the ARX system identification and the Model Predictive Controller.\n",
    "\n",
    "![connections](https://apmonitor.com/che263/uploads/Begin_Python/connections.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import tclab\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# generate step test data on Arduino\n",
    "filename = '12-tclab.csv'\n",
    "\n",
    "# heater steps\n",
    "Q1d = np.zeros(601)\n",
    "Q1d[10:100] = 80\n",
    "Q1d[100:200] = 20\n",
    "Q1d[200:300] = 70\n",
    "Q1d[300:400] = 50\n",
    "Q1d[400:500] = 100\n",
    "Q1d[500:] = 0\n",
    "\n",
    "Q2d = np.zeros(601)\n",
    "Q2d[50:150] = 35\n",
    "Q2d[150:250] = 95\n",
    "Q2d[250:350] = 25\n",
    "Q2d[350:450] = 100\n",
    "Q2d[450:550] = 45\n",
    "Q2d[550:] = 0\n",
    "\n",
    "# Connect to Arduino\n",
    "a = tclab.TCLab()\n",
    "fid = open(filename,'w')\n",
    "fid.write('Time,Q1,Q2,T1,T2\\n')\n",
    "fid.close()\n",
    "\n",
    "# run step test (20 min)\n",
    "for i in range(601):\n",
    "    # set heater values\n",
    "    a.Q1(Q1d[i])\n",
    "    a.Q2(Q2d[i])\n",
    "    print('Time: ' + str(2*i) + \\\n",
    "          ' Q1: ' + str(Q1d[i]) + \\\n",
    "          ' Q2: ' + str(Q2d[i]) + \\\n",
    "          ' T1: ' + str(a.T1)   + \\\n",
    "          ' T2: ' + str(a.T2))\n",
    "    # wait 2 seconds\n",
    "    time.sleep(2)\n",
    "    fid = open(filename,'a')\n",
    "    fid.write(str(2*i)+','+str(Q1d[i])+','+str(Q2d[i])+',' \\\n",
    "              +str(a.T1)+','+str(a.T2)+'\\n')\n",
    "    fid.close()\n",
    "# close connection to Arduino\n",
    "a.close()\n",
    "\n",
    "# read data file\n",
    "data = pd.read_csv(filename)\n",
    "\n",
    "# plot measurements\n",
    "plt.figure()\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(data['Time'],data['Q1'],'r-',label='Heater 1')\n",
    "plt.plot(data['Time'],data['Q2'],'b--',label='Heater 2')\n",
    "plt.ylabel('Heater (%)')\n",
    "plt.legend(loc='best')\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(data['Time'],data['T1'],'r.',label='Temperature 1')\n",
    "plt.plot(data['Time'],data['T2'],'b.',label='Temperature 2')\n",
    "plt.ylabel('Temperature (degC)')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('Time (sec)')\n",
    "plt.savefig('12-tclab.png')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Model Predictive Control with Time Series Model\n",
    "\n",
    "Run the following application with the TCLab connected. It uses the step test data from above to create a time series model. The application then uses this time series to create a Model Predictive Control application that optimizes the heaters to target temperature values. As the MPC runs, blow on the heaters to cause a disturbance. Observe how the predicted heater profile changes as the disturbance is applied.\n",
    "\n",
    "<video width=\"500\" height=\"350\" controls src=\"https://apmonitor.com/do/uploads/Main/tclab_arx_mpc.mp4\" />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import json\n",
    "# get gekko package with:\n",
    "#   pip install gekko\n",
    "from gekko import GEKKO\n",
    "# get tclab package with:\n",
    "#   pip install tclab\n",
    "from tclab import TCLab\n",
    "\n",
    "# Detect session is IPython\n",
    "try:\n",
    "    from IPython import get_ipython\n",
    "    from IPython.display import display,clear_output\n",
    "    get_ipython().run_line_magic('matplotlib', 'inline')\n",
    "    ipython = True\n",
    "    print('IPython Notebook')\n",
    "except:\n",
    "    ipython = False\n",
    "    print('Not IPython Notebook')\n",
    "\n",
    "# Connect to Arduino\n",
    "a = TCLab()\n",
    "\n",
    "# Final time\n",
    "tf = 10 # min\n",
    "# number of data points (every 2 seconds)\n",
    "n = tf * 30 + 1\n",
    "\n",
    "# Percent Heater (0-100%)\n",
    "Q1s = np.zeros(n)\n",
    "Q2s = np.zeros(n)\n",
    "\n",
    "# Temperatures (degC)\n",
    "T1m = a.T1 * np.ones(n)\n",
    "T2m = a.T2 * np.ones(n)\n",
    "# Temperature setpoints\n",
    "T1sp = T1m[0] * np.ones(n)\n",
    "T2sp = T2m[0] * np.ones(n)\n",
    "\n",
    "# Heater set point steps about every 150 sec\n",
    "T1sp[3:] = 50.0\n",
    "T2sp[40:] = 35.0\n",
    "T1sp[80:] = 30.0\n",
    "T2sp[120:] = 50.0\n",
    "T1sp[160:] = 45.0\n",
    "T2sp[200:] = 35.0\n",
    "T1sp[240:] = 60.0\n",
    "\n",
    "#########################################################\n",
    "# Initialize Model\n",
    "#########################################################\n",
    "# load data (20 min, dt=2 sec) and parse into columns\n",
    "data = pd.read_csv('12-tclab.csv')\n",
    "t = data['Time']\n",
    "u = data[['Q1','Q2']]\n",
    "y = data[['T1','T2']]\n",
    "\n",
    "# generate time-series model\n",
    "m = GEKKO(remote=False)\n",
    "\n",
    "##################################################################\n",
    "# system identification\n",
    "na = 2 # output coefficients\n",
    "nb = 2 # input coefficients\n",
    "print('Identify model')\n",
    "yp,p,K = m.sysid(t,u,y,na,nb,objf=10000,scale=False,diaglevel=1)\n",
    "\n",
    "##################################################################\n",
    "# plot sysid results\n",
    "plt.figure()\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(t,u)\n",
    "plt.legend([r'$Q_1$',r'$Q_2$'])\n",
    "plt.ylabel('MVs')\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(t,y)\n",
    "plt.plot(t,yp)\n",
    "plt.legend([r'$T_{1meas}$',r'$T_{2meas}$',\\\n",
    "            r'$T_{1pred}$',r'$T_{2pred}$'])\n",
    "plt.ylabel('CVs')\n",
    "plt.xlabel('Time')\n",
    "plt.savefig('sysid.png')\n",
    "plt.show()\n",
    "\n",
    "##################################################################\n",
    "# create control ARX model\n",
    "y = m.Array(m.CV,2)\n",
    "u = m.Array(m.MV,2)\n",
    "m.arx(p,y,u)\n",
    "\n",
    "# rename CVs\n",
    "TC1 = y[0]\n",
    "TC2 = y[1]\n",
    "\n",
    "# rename MVs\n",
    "Q1 = u[0]\n",
    "Q2 = u[1]\n",
    "\n",
    "# steady state initialization\n",
    "m.options.IMODE = 1\n",
    "m.solve(disp=False)\n",
    "\n",
    "# set up MPC\n",
    "m.options.IMODE   = 6 # MPC\n",
    "m.options.CV_TYPE = 1 # Objective type\n",
    "m.options.NODES   = 2 # Collocation nodes\n",
    "m.options.SOLVER  = 3 # IPOPT\n",
    "m.time=np.linspace(0,120,61)\n",
    "\n",
    "# Manipulated variables\n",
    "Q1.STATUS = 1  # manipulated\n",
    "Q1.FSTATUS = 0 # not measured\n",
    "Q1.DMAX = 50.0\n",
    "Q1.DCOST = 0.1\n",
    "Q1.UPPER = 100.0\n",
    "Q1.LOWER = 0.0\n",
    "\n",
    "Q2.STATUS = 1  # manipulated\n",
    "Q2.FSTATUS = 0 # not measured\n",
    "Q2.DMAX = 50.0\n",
    "Q2.DCOST = 0.1\n",
    "Q2.UPPER = 100.0\n",
    "Q2.LOWER = 0.0\n",
    "\n",
    "# Controlled variables\n",
    "TC1.STATUS = 1     # drive to set point\n",
    "TC1.FSTATUS = 1    # receive measurement\n",
    "TC1.TAU = 20       # response speed (time constant)\n",
    "TC1.TR_INIT = 2    # reference trajectory\n",
    "TC1.TR_OPEN = 0\n",
    "\n",
    "TC2.STATUS = 1     # drive to set point\n",
    "TC2.FSTATUS = 1    # receive measurement\n",
    "TC2.TAU = 20        # response speed (time constant)\n",
    "TC2.TR_INIT = 2    # dead-band\n",
    "TC2.TR_OPEN = 1\n",
    "\n",
    "# Main Loop\n",
    "start_time = time.time()\n",
    "prev_time = start_time\n",
    "tm = np.zeros(n)\n",
    "\n",
    "# Create plot\n",
    "if not ipython:\n",
    "    plt.figure(figsize=(10,7))\n",
    "    plt.ion()\n",
    "    plt.show()\n",
    "\n",
    "try:\n",
    "    for i in range(1,n-1):\n",
    "        # Sleep time\n",
    "        sleep_max = 2.0\n",
    "        sleep = sleep_max - (time.time() - prev_time)\n",
    "        if sleep>=0.01:\n",
    "            time.sleep(sleep-0.01)\n",
    "        else:\n",
    "            time.sleep(0.01)\n",
    "\n",
    "        # Record time and change in time\n",
    "        t = time.time()\n",
    "        dt = t - prev_time\n",
    "        prev_time = t\n",
    "        tm[i] = t - start_time\n",
    "\n",
    "        # Read temperatures in Celsius \n",
    "        T1m[i] = a.T1\n",
    "        T2m[i] = a.T2\n",
    "\n",
    "        # Insert measurements\n",
    "        TC1.MEAS = T1m[i]\n",
    "        TC2.MEAS = T2m[i]\n",
    "\n",
    "        # Adjust setpoints\n",
    "        db1 = 1.0 # dead-band\n",
    "        TC1.SPHI = T1sp[i] + db1\n",
    "        TC1.SPLO = T1sp[i] - db1\n",
    "\n",
    "        db2 = 0.2\n",
    "        TC2.SPHI = T2sp[i] + db2\n",
    "        TC2.SPLO = T2sp[i] - db2\n",
    "\n",
    "        # Adjust heaters with MPC\n",
    "        m.solve() \n",
    "\n",
    "        if m.options.APPSTATUS == 1:\n",
    "            # Retrieve new values\n",
    "            Q1s[i+1]  = Q1.NEWVAL\n",
    "            Q2s[i+1]  = Q2.NEWVAL\n",
    "            # get additional solution information\n",
    "            with open(m.path+'//results.json') as f:\n",
    "                results = json.load(f)\n",
    "        else:\n",
    "            # Solution failed\n",
    "            Q1s[i+1]  = 0.0\n",
    "            Q2s[i+1]  = 0.0\n",
    "\n",
    "        # Write new heater values (0-100)\n",
    "        a.Q1(Q1s[i])\n",
    "        a.Q2(Q2s[i])\n",
    "\n",
    "        # Plot\n",
    "        # Plot\n",
    "        if ipython:\n",
    "            plt.figure(figsize=(8,5))\n",
    "        else:\n",
    "            plt.clf()\n",
    "        ax=plt.subplot(3,1,1)\n",
    "        ax.grid()\n",
    "        plt.plot(tm[0:i+1],T1sp[0:i+1]+db1,'k-',\\\n",
    "                 label=r'$T_1$ target',linewidth=3)\n",
    "        plt.plot(tm[0:i+1],T1sp[0:i+1]-db1,'k-',\\\n",
    "                 label=None,linewidth=3)\n",
    "        plt.plot(tm[0:i+1],T1m[0:i+1],'r.',label=r'$T_1$ measured')\n",
    "        plt.plot(tm[i]+m.time,results['v1.bcv'],'r-',\\\n",
    "                 label=r'$T_1$ predicted',linewidth=3)\n",
    "        plt.plot(tm[i]+m.time,results['v1.tr_hi'],'k--',\\\n",
    "                 label=r'$T_1$ trajectory')\n",
    "        plt.plot(tm[i]+m.time,results['v1.tr_lo'],'k--')\n",
    "        plt.ylabel('Temperature (degC)')\n",
    "        plt.legend(loc=2)\n",
    "        ax=plt.subplot(3,1,2)\n",
    "        ax.grid()        \n",
    "        plt.plot(tm[0:i+1],T2sp[0:i+1]+db2,'k-',\\\n",
    "                 label=r'$T_2$ target',linewidth=3)\n",
    "        plt.plot(tm[0:i+1],T2sp[0:i+1]-db2,'k-',\\\n",
    "                 label=None,linewidth=3)\n",
    "        plt.plot(tm[0:i+1],T2m[0:i+1],'b.',label=r'$T_2$ measured')\n",
    "        plt.plot(tm[i]+m.time,results['v2.bcv'],'b-',\\\n",
    "                 label=r'$T_2$ predict',linewidth=3)\n",
    "        plt.plot(tm[i]+m.time,results['v2.tr_hi'],'k--',\\\n",
    "                 label=r'$T_2$ range')\n",
    "        plt.plot(tm[i]+m.time,results['v2.tr_lo'],'k--')\n",
    "        plt.ylabel('Temperature (degC)')\n",
    "        plt.legend(loc=2)\n",
    "        ax=plt.subplot(3,1,3)\n",
    "        ax.grid()\n",
    "        plt.plot([tm[i],tm[i]],[0,100],'k-',\\\n",
    "                 label='Current Time',linewidth=1)\n",
    "        plt.plot(tm[0:i+1],Q1s[0:i+1],'r.-',\\\n",
    "                 label=r'$Q_1$ history',linewidth=2)\n",
    "        plt.plot(tm[i]+m.time,Q1.value,'r-',\\\n",
    "                 label=r'$Q_1$ plan',linewidth=3)\n",
    "        plt.plot(tm[0:i+1],Q2s[0:i+1],'b.-',\\\n",
    "                 label=r'$Q_2$ history',linewidth=2)\n",
    "        plt.plot(tm[i]+m.time,Q2.value,'b-',\n",
    "                 label=r'$Q_2$ plan',linewidth=3)\n",
    "        plt.plot(tm[i]+m.time[1],Q1.value[1],color='red',\\\n",
    "                 marker='.',markersize=15)\n",
    "        plt.plot(tm[i]+m.time[1],Q2.value[1],color='blue',\\\n",
    "                 marker='X',markersize=8)\n",
    "        plt.ylabel('Heaters')\n",
    "        plt.xlabel('Time (sec)')\n",
    "        plt.legend(loc=2)\n",
    "        if ipython:\n",
    "            clear_output(wait=True)\n",
    "            display(plt.gcf())\n",
    "        else:\n",
    "            plt.draw()\n",
    "            plt.pause(0.05)\n",
    "\n",
    "    # Turn off heaters and close connection\n",
    "    a.Q1(0)\n",
    "    a.Q2(0)\n",
    "    a.close()\n",
    "    # Save figure\n",
    "    plt.savefig('12-tclab_mpc.png')\n",
    "\n",
    "# Allow user to end loop with Ctrl-C           \n",
    "except KeyboardInterrupt:\n",
    "    # Turn off heaters and close connection\n",
    "    a.Q1(0)\n",
    "    a.Q2(0)\n",
    "    a.close()\n",
    "    print('Shutting down')\n",
    "    plt.savefig('12-tclab_mpc.png')\n",
    "\n",
    "# Make sure serial connection still closes when there's an error\n",
    "except:           \n",
    "    # Disconnect from Arduino\n",
    "    a.Q1(0)\n",
    "    a.Q2(0)\n",
    "    a.close()\n",
    "    print('Error: Shutting down')\n",
    "    plt.savefig('12-tclab_mpc.png')\n",
    "    raise"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}