LSTM_SISO_MPC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Emulating a SISO MPC with an LSTM Network\n",
    "\n",
    "The purpose of this exercise is to automate a temperature control process with an LSTM network. The LSTM network is trained from a 1x1 SISO MPC (Single Input Single Output, Model Predictive Control). LSTM (Long Short Term Memory) networks are a special type of RNN (Recurrent Neural Network) that is structured to remember and predict based on long-term dependencies that are trained with time-series data. An LSTM repeating module has four interacting components.\n",
    "\n",
    "<img align=left width=400px src='https://apmonitor.com/pds/uploads/Main/lstm_mpc_1x1.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The LSTM is trained (parameters adjusted) with an input window of prior data and minimized difference between the predicted and next measured value. Sequential methods predict just one next value based on the window of prior data. In this case, the error between the set point and measured value is the feature and the heater value is the output label.\n",
    "\n",
    "See [Automation with LSTM Networks](https://apmonitor.com/pds/index.php/Main/LSTMAutomation) for additional details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:25:23.638976Z",
     "start_time": "2021-02-17T16:25:20.302787Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import random\n",
    "import time\n",
    "from tqdm import tqdm # Progress bar\n",
    "\n",
    "# For scaling, feature selection\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from sklearn.feature_selection import SelectKBest, f_regression\n",
    "from sklearn.model_selection import train_test_split \n",
    "\n",
    "# For LSTM model\n",
    "from keras.models import Sequential\n",
    "from keras.layers import LSTM, Dropout, Dense\n",
    "from keras.callbacks import EarlyStopping\n",
    "from tqdm.keras import TqdmCallback\n",
    "from keras.models import load_model\n",
    "\n",
    "# For Gekko and TCLab\n",
    "import tclab\n",
    "from gekko import GEKKO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use [TCLab](https://apmonitor.com/heat.htm) Microcontroller or Emulator (Digital Twin)\n",
    "\n",
    "Change `tclab_hardware` to `True` if TCLab device is connected.\n",
    "\n",
    "<img width=550px align=left src='https://apmonitor.com/pdc/uploads/Main/tclab_connect.png'>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tclab_hardware = False\n",
    "if tclab_hardware:\n",
    "    mlab = tclab.TCLab      # Physical hardware\n",
    "else:\n",
    "    speedup = 100\n",
    "    mlab = tclab.setup(connected=False, speedup = speedup) # Emulator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PID controller \n",
    "\n",
    "See [LSTM Network Replaces PID](https://github.com/APMonitor/pds/blob/main/LSTM_Automation.ipynb) for code with a PID controller instead of a Model Predictive Controller (MPC)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MPC (Model Predictive Control)\n",
    "\n",
    "See the last example for [MPC with TCLab](https://apmonitor.com/pds/notebooks/12_time_series.html) or the [Linear MPC with TCLab](https://apmonitor.com/do/index.php/Main/TCLabF). This controller is implemented instead of a PID controller. The steps to create an MPC application are:\n",
    "\n",
    "- generate step test data\n",
    "- identify model\n",
    "- initialize controller\n",
    "- create MPC function to return Q1\n",
    "\n",
    "**Generate Step Test Data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate step test data on TCLab\n",
    "filename = 'tclab_data.csv'\n",
    "\n",
    "n = 500\n",
    "tm = np.linspace(0,n*2,n+1)\n",
    "T1 = np.zeros(n+1)\n",
    "T2 = np.zeros(n+1)\n",
    "\n",
    "# heater steps\n",
    "Q1d = np.zeros(n+1)\n",
    "Q1d[10:150]=80; Q1d[150:300]=20; Q1d[300:450]=70; Q1d[450:]=50\n",
    "\n",
    "Q2d = np.zeros(n+1)\n",
    "Q2d[50:150]=35; Q2d[150:250]=95; Q2d[250:350]=25; Q2d[350:]=100\n",
    "\n",
    "p1 = 1 if tclab_hardware else 100\n",
    "# Connect to TCLab\n",
    "with mlab() as lab:\n",
    "    # run step test\n",
    "    i = 0\n",
    "    for t in tclab.clock(tm[-1]+1, 2):\n",
    "        # set heater values\n",
    "        lab.Q1(Q1d[i])\n",
    "        lab.Q2(Q2d[i])\n",
    "        T1[i] = lab.T1\n",
    "        T2[i] = lab.T2\n",
    "        if i%p1==0:\n",
    "            print('Time: ' + str(2*i) + \\\n",
    "              ' Q1: ' + str(Q1d[i]) + \\\n",
    "              ' Q2: ' + str(Q2d[i]) + \\\n",
    "              ' T1: ' + str(round(T1[i],2)) + \\\n",
    "              ' T2: ' + str(round(T2[i],2)))\n",
    "        i+=1\n",
    "\n",
    "# write data to file\n",
    "fid = open(filename,'w')\n",
    "fid.write('Time,Q1,Q2,T1,T2\\n')\n",
    "for i in range(n+1):\n",
    "    fid.write(str(tm[i])+','+str(Q1d[i])+','+str(Q2d[i])+',' \\\n",
    "              +str(T1[i])+','+str(T2[i])+'\\n')\n",
    "fid.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('tclab_data.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Identify Model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#########################################################\n",
    "# Initialize Model\n",
    "#########################################################\n",
    "# load data and parse into columns\n",
    "data = pd.read_csv('tclab_data.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; nb=2 # use 2nd order model\n",
    "print('Identify model')\n",
    "yp,p,K = m.sysid(t,u,y,na,nb,objf=10000,scale=False,diaglevel=0,pred='model')\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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initialize Controller**\n",
    "\n",
    "Create 2x2 (2 Heaters, 2 Temperature Sensors) MPC but only use 1x1 (`Q1`, `T1`) for this application. The LSTM learning can be extended to the 2x2 MPC (MVs: `Q1` and `Q2`, CVs: `T1` and `T2`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################################################\n",
    "# create control ARX model\n",
    "m = GEKKO(remote=False)\n",
    "m.y = m.Array(m.CV,2)\n",
    "m.u = m.Array(m.MV,2)\n",
    "m.arx(p,m.y,m.u)\n",
    "\n",
    "# rename CVs\n",
    "m.TC1 = m.y[0]\n",
    "m.TC2 = m.y[1]\n",
    "\n",
    "# rename MVs\n",
    "m.Q1 = m.u[0]\n",
    "m.Q2 = m.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 = 2 # Objective type\n",
    "m.options.NODES   = 2 # Collocation nodes\n",
    "m.options.SOLVER  = 1 # APOPT\n",
    "m.time=np.linspace(0,60,31)\n",
    "\n",
    "# Manipulated variables\n",
    "m.Q1.STATUS = 1  # manipulated\n",
    "m.Q1.FSTATUS = 0 # not measured\n",
    "m.Q1.DMAX = 100.0\n",
    "m.Q1.DCOST = 2.0\n",
    "m.Q1.UPPER = 100.0\n",
    "m.Q1.LOWER = 0.0\n",
    "\n",
    "m.Q2.STATUS = 0  # manipulated, turn off Q2\n",
    "m.Q2.FSTATUS = 1 # use measured value\n",
    "m.Q2.DMAX = 100.0\n",
    "m.Q2.DCOST = 2.0\n",
    "m.Q2.UPPER = 100.0\n",
    "m.Q2.LOWER = 0.0\n",
    "m.Q2.MEAS = 0    # set Q2=0\n",
    "\n",
    "# Controlled variables\n",
    "m.TC1.STATUS = 1     # drive to set point\n",
    "m.TC1.FSTATUS = 1    # receive measurement\n",
    "m.TC1.TAU = 8        # response speed (time constant)\n",
    "m.TC1.TR_INIT = 2    # reference trajectory\n",
    "m.TC1.TR_OPEN = 5\n",
    "\n",
    "m.TC2.STATUS = 0     # drive to set point\n",
    "m.TC2.FSTATUS = 1    # receive measurement\n",
    "m.TC2.TAU = 8        # response speed (time constant)\n",
    "m.TC2.TR_INIT = 2    # dead-band\n",
    "m.TC2.TR_OPEN = 1    # for CV_TYPE=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-18T15:06:22.503795Z",
     "start_time": "2021-02-18T15:06:22.495792Z"
    }
   },
   "outputs": [],
   "source": [
    "def mpc(m,T1,T1sp,T2,T2sp):    \n",
    "    # Insert measurements\n",
    "    m.TC1.MEAS = T1\n",
    "    m.TC2.MEAS = T2\n",
    "\n",
    "    # Adjust setpoints\n",
    "    db1 = 1.0 # dead-band\n",
    "    m.TC1.SP   = T1sp\n",
    "    m.TC1.SPHI = T1sp + db1\n",
    "    m.TC1.SPLO = T1sp - db1\n",
    "\n",
    "    db2 = 0.2\n",
    "    m.TC2.SP   = T2sp\n",
    "    m.TC2.SPHI = T2sp + db2\n",
    "    m.TC2.SPLO = T2sp - db2\n",
    "    \n",
    "    # Adjust heaters with MPC\n",
    "    m.solve(disp=False) \n",
    "\n",
    "    if m.options.APPSTATUS == 1:\n",
    "        # Retrieve new values\n",
    "        Q1  = m.Q1.NEWVAL\n",
    "        Q2  = m.Q2.NEWVAL\n",
    "    else:\n",
    "        # Solution failed\n",
    "        Q1  = 0.0\n",
    "        Q2  = 0.0    \n",
    "    return [Q1,Q2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate data for training LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:25:24.035744Z",
     "start_time": "2021-02-17T16:25:23.687213Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "##### Set up run parameters #####\n",
    "run_time = 90.0 # minutes\n",
    "\n",
    "loops = int(30.0*run_time+1) # cycles (2 sec each)\n",
    "\n",
    "# arrays for storing data\n",
    "T1 = np.zeros(loops) # measured T (degC)\n",
    "T2 = np.zeros(loops) # measured T (degC)\n",
    "Q1 = np.zeros(loops) # Heater values\n",
    "Q2 = np.zeros(loops) # Heater values\n",
    "tm = np.linspace(0,2*(loops-1),loops) # Time\n",
    "\n",
    "# Temperature set point (degC)\n",
    "with mlab() as lab:\n",
    "    Tsp1 = np.ones(loops) * lab.T1\n",
    "    Tsp2 = np.ones(loops) * lab.T2\n",
    "\n",
    "# vary temperature setpoint\n",
    "end = 2 # leave first couple cycles of temp set point as room temp\n",
    "while end <= loops:\n",
    "    start = end\n",
    "    # keep new temp set point value for anywhere from 3 to 5 min\n",
    "    end += random.randint(90,150) \n",
    "    Tsp1[start:end] = random.randint(30,70)\n",
    "Tsp1[-120:] = Tsp1[0] # last 4 minutes at room temperature\n",
    "    \n",
    "if tclab_hardware:\n",
    "    # print every cycle with hardware\n",
    "    p1 = 10; p2 = 1\n",
    "else:\n",
    "    # print 20x less with emulator\n",
    "    p1 = 200; p2 = 20   \n",
    "\n",
    "# Plot\n",
    "plt.plot(tm,Tsp1,'b.-')\n",
    "plt.xlabel('Time',size=14)\n",
    "plt.ylabel(r'Temp SP ($^oC$)',size=14)\n",
    "plt.xticks(size=12)\n",
    "plt.yticks(size=12)\n",
    "plt.savefig('SP_profile.png');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:26:54.869018Z",
     "start_time": "2021-02-17T16:25:24.035744Z"
    }
   },
   "outputs": [],
   "source": [
    "# Data collection\n",
    "with mlab() as lab:\n",
    "    # Find current T1, T2\n",
    "    print('Temperature 1: {0:0.2f} °C'.format(lab.T1))\n",
    "    print('Temperature 2: {0:0.2f} °C'.format(lab.T2))\n",
    "    \n",
    "    i = 0\n",
    "    for t in tclab.clock(tm[-1]+1, 2):\n",
    "        # Read temperatures in Celcius \n",
    "        T1[i] = lab.T1; T2[i] = lab.T2\n",
    "\n",
    "        # Calculate MPC output every 2 sec\n",
    "        try:\n",
    "            [Q1[i],Q2[i]] = mpc(m,T1[i],Tsp1[i],T2[i],Tsp2[i])\n",
    "        except:\n",
    "            Q1[i] = 0; Q2[i] = 0 # catch any failure to converge\n",
    "        # Write heater output (0-100)\n",
    "        lab.Q1(Q1[i])\n",
    "        lab.Q2(Q2[i])\n",
    "         \n",
    "        if i%p1==0:            \n",
    "            print('  Time_____Q1___Tsp1_____T1______Q2____Tsp2_____T2')\n",
    "        if i%p2==0:\n",
    "            print(('{:6.1f} {:6.2f} {:6.2f} {:6.2f}  {:6.2f}  {:6.2f} {:6.2f}').format( \\\n",
    "                      tm[i],Q1[i],Tsp1[i],T1[i],Q2[i],Tsp2[i],T2[i]))\n",
    "        i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:26:55.171468Z",
     "start_time": "2021-02-17T16:26:54.870445Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Save csv file\n",
    "df = pd.DataFrame()\n",
    "df['time'] = tm[:i]\n",
    "df['Q1'] = Q1[:i]\n",
    "df['T1'] = T1[:i]\n",
    "df['Tsp'] = Tsp1[:i]\n",
    "df.set_index('time',inplace=True)\n",
    "df.to_csv('MPC_train_data.csv',index=False)\n",
    "\n",
    "# Plot\n",
    "df[['Q1','T1','Tsp']].plot()\n",
    "plt.savefig('MPC_train.png')\n",
    "plt.show()\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature engineering: create SP error feature and analyze feature importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:26:55.327215Z",
     "start_time": "2021-02-17T16:26:55.177474Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Create new feature: setpoint error\n",
    "df['err'] = df['Tsp'] - df['T1']\n",
    "\n",
    "# Load possible features\n",
    "X = df[['T1','Tsp','err']]\n",
    "y = np.ravel(df[['Q1']])\n",
    "\n",
    "# SelectKBest feature selection\n",
    "bestfeatures = SelectKBest(score_func=f_regression, k='all')\n",
    "fit = bestfeatures.fit(X,y)\n",
    "plt.bar(x=X.columns,height=fit.scores_);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature selection, scaling, and formatting data to LSTM input format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:26:55.373213Z",
     "start_time": "2021-02-17T16:26:55.344216Z"
    }
   },
   "outputs": [],
   "source": [
    "X = df[['Tsp','err']].values\n",
    "y = df[['Q1']].values\n",
    "\n",
    "# Scale data\n",
    "s_x = MinMaxScaler()\n",
    "Xs = s_x.fit_transform(X)\n",
    "\n",
    "s_y = MinMaxScaler()\n",
    "ys = s_y.fit_transform(y)\n",
    "\n",
    "# Each input uses last 'window' number of Tsp and err to predict the next Q1\n",
    "window = 15\n",
    "X_lstm = []\n",
    "y_lstm = []\n",
    "for i in range(window,len(df)):\n",
    "    X_lstm.append(Xs[i-window:i])\n",
    "    y_lstm.append(ys[i])\n",
    "\n",
    "# Reshape data to format accepted by LSTM\n",
    "X_lstm, y_lstm = np.array(X_lstm), np.array(y_lstm)\n",
    "\n",
    "# Split into train and test \n",
    "Xtrain, Xtest, ytrain, ytest = train_test_split(X_lstm,y_lstm,test_size=0.2,shuffle=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Keras LSTM "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:40:26.987200Z",
     "start_time": "2021-02-17T16:33:55.613757Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Keras LSTM model\n",
    "model = Sequential()\n",
    "\n",
    "# First layer specifies input_shape and returns sequences\n",
    "model.add(LSTM(units=100, return_sequences=True, \n",
    "               input_shape=(Xtrain.shape[1],Xtrain.shape[2])))\n",
    "# Dropout layer to prevent overfitting\n",
    "model.add(Dropout(rate=0.1))\n",
    "\n",
    "# Last layer doesn't return sequences (middle layers should return sequences)\n",
    "model.add(LSTM(units=100))\n",
    "model.add(Dropout(rate=0.1))\n",
    "\n",
    "# Dense layer to return prediction\n",
    "model.add(Dense(1))\n",
    "\n",
    "# Compile model; adam optimizer, mse loss\n",
    "model.compile(optimizer='adam', loss='mean_squared_error')\n",
    "\n",
    "es = EarlyStopping(monitor='val_loss',mode='min',verbose=1,patience=25)\n",
    "\n",
    "result = model.fit(Xtrain, ytrain, verbose=0, validation_split=0.2,\n",
    "                   callbacks = [TqdmCallback(verbose=1)],#es\n",
    "                   batch_size=100,\n",
    "                   epochs=300)\n",
    "\n",
    "# Plot loss and save model\n",
    "epochs = es.stopped_epoch\n",
    "plt.semilogy(result.history['loss'],label='loss')\n",
    "plt.semilogy(result.history['val_loss'],label='val_loss')\n",
    "plt.legend()\n",
    "\n",
    "model.save('lstm_control.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preliminary model performance assessment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:40:28.617490Z",
     "start_time": "2021-02-17T16:40:26.990200Z"
    }
   },
   "outputs": [],
   "source": [
    "# Predict using LSTM\n",
    "yp_s = model.predict(Xtest)\n",
    "\n",
    "# Unscale data\n",
    "Xtest_us = s_x.inverse_transform(Xtest[:,-1,:])\n",
    "ytest_us = s_y.inverse_transform(ytest)\n",
    "yp = s_y.inverse_transform(yp_s)\n",
    "\n",
    "# Derive Tsp (sp) and T1 (pv) from X data\n",
    "sp = Xtest_us[:,0]\n",
    "pv = Xtest_us[:,0] + Xtest_us[:,1]\n",
    "\n",
    "# Plot SP, MPC response, and LSTM response\n",
    "tm = np.linspace(0,2*(len(sp)-1)/60,len(sp))\n",
    "plt.plot(tm,sp,'k-',label='$SP$ $(^oC)$')\n",
    "plt.plot(tm,pv,'r-',label='$T_1$ $(^oC)$')\n",
    "plt.plot(tm,ytest_us,'b-',label='$Q_{MPC}$ (%)')\n",
    "plt.plot(tm,yp,'g-',label='$Q_{LSTM}$ (%)')\n",
    "plt.legend()\n",
    "plt.xlabel('Time (min)')\n",
    "plt.ylabel('Value')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate SP data for test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:40:28.976494Z",
     "start_time": "2021-02-17T16:40:28.622493Z"
    }
   },
   "outputs": [],
   "source": [
    "# Run time in minutes\n",
    "run_time = 45.0\n",
    "\n",
    "# Number of cycles\n",
    "loops = int(30.0*run_time)\n",
    "\n",
    "# arrays for storing data\n",
    "T1 = np.zeros(loops) # measured T (degC)\n",
    "T2 = np.zeros(loops)\n",
    "Q1mpc = np.zeros(loops) # Heater values for MPC controller\n",
    "Q2mpc = np.zeros(loops)\n",
    "Qlstm = np.zeros(loops) # Heater values for LSTM controller\n",
    "tm = np.linspace(0,2*(loops-1),loops) # Time\n",
    "\n",
    "# Temperature set point (degC)\n",
    "with mlab() as lab:\n",
    "    Tsp1 = np.ones(loops) * lab.T1\n",
    "    Tsp2 = np.ones(loops) * lab.T2\n",
    "\n",
    "# vary temperature setpoint\n",
    "end = window + 5 # leave 1st window + 10 seconds of temp set point as room temp\n",
    "while end <= loops: \n",
    "    start = end\n",
    "    # keep new temp set point value for anywhere from 3 to 5 min\n",
    "    end += random.randint(90,150) \n",
    "    Tsp1[start:end] = random.randint(30,70)\n",
    "    \n",
    "# leave last 120 seconds as room temp\n",
    "Tsp1[-60:] = Tsp1[0]\n",
    "plt.plot(Tsp1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part A: Run the controller with the MPC, but also output the LSTM for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:40:29.006489Z",
     "start_time": "2021-02-17T16:40:28.983493Z"
    }
   },
   "outputs": [],
   "source": [
    "# LSTM Controller\n",
    "def lstm(T1_m, Tsp_m):\n",
    "    # Calculate error (necessary feature for LSTM input)\n",
    "    err = Tsp_m - T1_m\n",
    "    \n",
    "    # Format data for LSTM input\n",
    "    X = np.vstack((Tsp_m,err)).T\n",
    "    Xs = s_x.transform(X)\n",
    "    Xs = np.reshape(Xs, (1, Xs.shape[0], Xs.shape[1]))\n",
    "    \n",
    "    # Predict Q for controller and unscale\n",
    "    Q1c_s = model.predict(Xs)\n",
    "    Q1c = s_y.inverse_transform(Q1c_s)[0][0]\n",
    "    \n",
    "    # Ensure Q1c is between 0 and 100\n",
    "    Q1c = np.clip(Q1c,0.0,100.0)\n",
    "    return Q1c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:45:48.349921Z",
     "start_time": "2021-02-17T16:44:02.026886Z"
    }
   },
   "outputs": [],
   "source": [
    "# Run test\n",
    "with mlab() as lab:\n",
    "    # Find current T1, T2\n",
    "    print('Temperature 1: {0:0.2f} °C'.format(lab.T1))\n",
    "    print('Temperature 2: {0:0.2f} °C'.format(lab.T2))\n",
    "    \n",
    "    i = 0\n",
    "    for t in tclab.clock(tm[-1]+1, 2):\n",
    "               \n",
    "        # Read temperatures in Celcius \n",
    "        T1[i] = lab.T1; T2[i] = lab.T2\n",
    "\n",
    "        # Calculate MPC output every 2 sec\n",
    "        try:\n",
    "            [Q1mpc[i],Q2mpc[i]] = mpc(m,T1[i],Tsp1[i],T2[i],Tsp2[i])\n",
    "        except:\n",
    "            Q1mpc[i]=0; Q2mpc[i]=0\n",
    "        # Write heater output (0-100)\n",
    "        lab.Q1(Q1mpc[i])\n",
    "        lab.Q2(Q2mpc[i])\n",
    "         \n",
    "        if i%p1==0:            \n",
    "            print('  Time_____Q1___Tsp1_____T1')\n",
    "        if i%p2==0:\n",
    "            print(('{:6.1f} {:6.2f} {:6.2f} {:6.2f}').format( \\\n",
    "                      tm[i],Q1mpc[i],Tsp1[i],T1[i]))\n",
    "        \n",
    "        # Run LSTM model to get Q1 value for control\n",
    "        if i >= window:\n",
    "            # Load data for model\n",
    "            T1_m = T1[i-window:i]\n",
    "            Tsp_m = Tsp1[i-window:i]\n",
    "            # Predict and store LSTM value for comparison\n",
    "            Qlstm[i] = lstm(T1_m,Tsp_m)\n",
    "        i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:45:48.536872Z",
     "start_time": "2021-02-17T16:45:48.352873Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(tm[:i],Tsp1[:i],'k-',label='SP $(^oC)$')\n",
    "plt.plot(tm[:i],T1[:i],'r-',label='$T_1$ $(^oC)$')\n",
    "plt.legend(loc='upper right',fontsize=14)\n",
    "plt.ylim((0,100))\n",
    "plt.xlabel('Time (s)',size=14)\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(tm[:i],Q1mpc[:i],'b-',label='$Q_{MPC}$ (%)')\n",
    "plt.plot(tm[:i],Qlstm[:i],'g-',label='$Q_{LSTM}$ (%)')\n",
    "plt.legend(loc='upper right',fontsize=14)\n",
    "plt.ylim((0,100))\n",
    "plt.xlabel('Time (s)',size=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part B: Run the controller with just the LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-17T16:47:32.690137Z",
     "start_time": "2021-02-17T16:45:48.538872Z"
    }
   },
   "outputs": [],
   "source": [
    "# Run test\n",
    "with mlab() as lab:\n",
    "    # Find current T1, T2\n",
    "    print('Temperature 1: {0:0.2f} °C'.format(lab.T1))\n",
    "    print('Temperature 2: {0:0.2f} °C'.format(lab.T2))\n",
    "    \n",
    "    i = 0\n",
    "    for t in tclab.clock(tm[-1]+1, 2):\n",
    "        # Read temperatures in Celcius \n",
    "        T1[i] = lab.T1; T2[i] = lab.T2\n",
    "\n",
    "        # Run LSTM model to get Q1 value for control\n",
    "        if i >= window:\n",
    "            # Load data for model\n",
    "            T1_m = T1[i-window:i]\n",
    "            Tsp_m = Tsp1[i-window:i]\n",
    "            # Predict and store LSTM value for comparison\n",
    "            Qlstm[i] = lstm(T1_m,Tsp_m)\n",
    "            \n",
    "        if i%p1==0:            \n",
    "            print('  Time_____Q1___Tsp1_____T1')\n",
    "        if i%p2==0:\n",
    "            print(('{:6.1f} {:6.2f} {:6.2f} {:6.2f}').format( \\\n",
    "                      tm[i],Qlstm[i],Tsp1[i],T1[i]))\n",
    "\n",
    "        # Write heater output (0-100)\n",
    "        lab.Q1(Qlstm[i])\n",
    "        i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-18T15:13:49.761601Z",
     "start_time": "2021-02-18T15:13:49.565402Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "plt.plot(tm[:i],Tsp1[:i],'k-',label='SP $(^oC)$')\n",
    "plt.plot(tm[:i],T1[:i],'r-',label='$T_1$ $(^oC)$')\n",
    "plt.plot(tm[:i],Qlstm[:i],'g-',label='$Q_{LSTM}$ (%)')\n",
    "plt.legend()\n",
    "plt.ylim((0,100))\n",
    "plt.xlabel('Time (s)')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part C\n",
    "\n",
    "Use the LSTM controller to control the [TCLab microcontroller](https://apmonitor.com/heat.htm) or [Emulated TCLab](https://tclab.readthedocs.io/en/latest/notebooks/04_Emulation_of_TCLab_for_Offline_Use.html) from a separate program. Demonstrate that you can [package and deploy the LSTM model](https://apmonitor.com/pds/index.php/Main/DeployMachineLearning) so that it could be used on a different computer than the training computer. The application does not necessarily need to run on a separate computer, but show that it is a separate program and independent of this training notebook.\n",
    "\n",
    "<img align=left width=150px src='https://apmonitor.com/pds/uploads/Main/deploy_machine_learning.png'>\n",
    "\n",
    "<img align=left width=400px src='https://apmonitor.com/pdc/uploads/Main/tclab_connect.png'>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "# export model name and scaling\n",
    "z = ['lstm_control.h5',s_x,s_y,window]\n",
    "pickle.dump(z,open('lstm_control.pkl','wb'))"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}