{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 0,
     "link": "[](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html)",
     "section": ""
    }
   },
   "source": [
    "<!--NOTEBOOK_HEADER-->\n",
    "*This notebook contains material from [CBE32338](https://jckantor.github.io/CBE32338);\n",
    "content is available [on Github](https://github.com/jckantor/CBE32338.git).*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 0,
     "link": "[](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html)",
     "section": ""
    }
   },
   "source": [
    "<!--NAVIGATION-->\n",
    "< [2.3 First Order Model for a Single Heater](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html) | [Contents](toc.html) | [2.5 Two State Model for a Single Heater](https://jckantor.github.io/CBE32338/02.05-Two-State-Model-for-a-Single-Heater.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE32338/blob/master/docs/02.04-Two-Input-Two-Output-Model.ipynb\"> <img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.ipynb\"> <img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[2.4 Two-Input, Two-Output Model](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4-Two-Input,-Two-Output-Model)",
     "section": "2.4 Two-Input, Two-Output Model"
    }
   },
   "source": [
    "# 2.4 Two-Input, Two-Output Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.1 Interacting Heaters](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.1-Interacting-Heaters)",
     "section": "2.4.1 Interacting Heaters"
    }
   },
   "source": [
    "## 2.4.1 Interacting Heaters\n",
    "\n",
    "Our next model considers the interaction between the two heaters. From the step response it is obvious that the heaters interact, and we need to understand how that affects the dynamic performance.  Will we see qualitative differences in the step response once we've included the second heater?\n",
    "\n",
    "\\begin{align}\n",
    "C_p\\frac{dT_1}{dt} & = U_a(T_{amb} - T_1) + U_b(T_2 - T_1) + P_1u_1\\\\\n",
    "C_p\\frac{dT_2}{dt} & = U_a(T_{amb} - T_2) + U_b(T_1 - T_2) + P_2u_2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.1 Interacting Heaters](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.1-Interacting-Heaters)",
     "section": "2.4.1 Interacting Heaters"
    }
   },
   "source": [
    "As a first step we isolate the bare derivatives on the left-hand side of the equations. This is good practice in situations where you expect to be doing numerical simulation.\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dT_1}{dt} & = \\frac{U_a}{C_p}(T_{amb} - T_1) + \\frac{U_b}{C_p}(T_2 - T_1) + \\frac{P_1}{C_p}u_1\\\\\n",
    "\\frac{dT_2}{dt} & = \\frac{U_a}{C_p}(T_{amb} - T_2) + \\frac{U_b}{C_p}(T_1 - T_2) + \\frac{P_2}{C_p}u_2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.2 Deviation Variables](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.2-Deviation-Variables)",
     "section": "2.4.2 Deviation Variables"
    }
   },
   "source": [
    "## 2.4.2 Deviation Variables\n",
    "\n",
    "As before, the equations appear somewhat simpler if we introduce deviation variables from a nominal state. In this case we'll use the ambient temperature as the reference, so the deviation variable will express deviation from ambient.\n",
    "\n",
    "\\begin{align}\n",
    "T_1' & = T_1 - T_{amb} \\\\\n",
    "T_2' & = T_2 - T_{amb}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.2 Deviation Variables](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.2-Deviation-Variables)",
     "section": "2.4.2 Deviation Variables"
    }
   },
   "source": [
    "Then gathering terms\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dT_1'}{dt} & = -\\frac{U_a+U_b}{C_p}T_1' + \\frac{U_b}{C_p}T_2' + \\frac{P_1}{C_p}u_1\\\\\n",
    "\\frac{dT_2'}{dt} & = -\\frac{U_a+U_b}{C_p}T_2' + \\frac{U_b}{C_p}T_1' + \\frac{P_2}{C_p}u_2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.2 Deviation Variables](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.2-Deviation-Variables)",
     "section": "2.4.2 Deviation Variables"
    }
   },
   "source": [
    "An alternative way to express these equations is in matrix/vector form.  We won't use this today, but later will make extensive use of this 'state-space' model of a linear system.\n",
    "\n",
    "\\begin{align}\n",
    "\\left[\\begin{array}{c}\\frac{dT_1'}{dt} \\\\ \\frac{dT_2'}{dt}\\end{array}\\right] & = \n",
    "\\left[\\begin{array}{cc}-\\frac{U_a+U_b}{C_p} & \\frac{U_b}{C_p} \\\\ \\frac{U_b}{C_p} & -\\frac{U_a+U_b}{C_p} \\end{array}\\right]\n",
    "\\left[\\begin{array}{c}T_1' \\\\ T_2'\\end{array}\\right] +\n",
    "\\left[\\begin{array}{cc}\\frac{P_1}{C_p} & 0 \\\\ 0 & \\frac{P_2}{C_p} \\end{array}\\right]\n",
    "\\left[\\begin{array}{c}u_1 \\\\ u_2\\end{array}\\right]\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.3 Steady State](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.3-Steady-State)",
     "section": "2.4.3 Steady State"
    }
   },
   "source": [
    "## 2.4.3 Steady State\n",
    "\n",
    "The steady state is found by setting the derivatives to zero. We'll use a bar over the dependent variables to designate steady state values.\n",
    "\n",
    "\\begin{align}\n",
    "0 & = -(U_a+U_b)\\bar{T}_{1}' + U_b\\bar{T}_{2}' + P_1\\bar{u}_{1}\\\\\n",
    "0 & = -(U_a+U_b)\\bar{T}_{2}' + U_b\\bar{T}_{1}' + P_2\\bar{u}_{1}\n",
    "\\end{align}\n",
    "\n",
    "Our challenge is to solve the following 2x2 system of linear equations for the steady state temperatures.\n",
    "\n",
    "\\begin{align} \n",
    "\\left[\\begin{array}{cc}(U_a+U_b) & -U_b \\\\ -U_b & (U_a+U_b) \\end{array}\\right]\n",
    "\\left[\\begin{array}{c}\\bar{T}_1' \\\\ \\bar{T}_2'\\end{array}\\right] & = \n",
    "\\left[\\begin{array}{c} P_1\\bar{u}_1 \\\\ P_2\\bar{u}_2\\end{array}\\right]\n",
    "\\end{align}\n",
    "\n",
    "The solution to these equations is a bit of work, but you should be able to verify\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{T}_1' & = \\frac{(U_a+U_b)}{U_a(U_a+2U_b)}P_1\\bar{u}_1 + \\frac{U_b}{U_a(U_a+2U_b)}P_2\\bar{u}_2 \\\\\n",
    "\\bar{T}_2' & = \\frac{U_b}{U_a(U_a+2U_b)}P_1\\bar{u}_1 + \\frac{(U_a+U_b)}{U_a(U_a+2U_b)}P_2\\bar{u}_2\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.4 Estimating $U_a$ and $U_b$](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.4-Estimating-$U_a$-and-$U_b$)",
     "section": "2.4.4 Estimating $U_a$ and $U_b$"
    }
   },
   "source": [
    "## 2.4.4 Estimating $U_a$ and $U_b$\n",
    "\n",
    "Experimental step response data gives us the information needed to estimate parameters $U_a$ and $U_b$. Rearranging the steady state equations, we find a pair of linear equations for $U_a$ and $U_b$.\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{T}_{1}'U_a  + (\\bar{T}_{1}'-\\bar{T}_{2}')U_b & = P_1\\bar{u}_{1}\\\\\n",
    "\\bar{T}_{2}'U_a + (\\bar{T}_{2}'-\\bar{T}_{1}')U_b & = P_2\\bar{u}_{1}\n",
    "\\end{align}\n",
    "\n",
    "There is an analytical solution to these equations\n",
    "\n",
    "\\begin{align}\n",
    "U_a & = \\frac{P_1\\bar{u}_1 + P_2\\bar{u}_2}{\\bar{T}_1' + \\bar{T}_2'}\\\\\n",
    "U_b & = \\frac{\\bar{T}_2'P_1\\bar{u}_1 - \\bar{T}_1'P_2\\bar{u}_2}{(\\bar{T}_1' - \\bar{T}_2')(\\bar{T}_1' + \\bar{T}_2')}\n",
    "\\end{align}\n",
    "\n",
    "The step response data gives us information about steady-state that can be used to estimate values for $U_a$ and $U_b$. First we plot the step response data (see Part 1 if this cell cannot find the data set.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.4 Estimating $U_a$ and $U_b$](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.4-Estimating-$U_a$-and-$U_b$)",
     "section": "2.4.4 Estimating $U_a$ and $U_b$"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VOXZ//HPlR0SIKxh3wRZVBYT\nAUUtwX15FH1cq9atpdbaR2tbrdqqbe1TfWprrb+6VavWaqOl4oIbFokimwKyRJF9JywhhCyQda7f\nH+dknyQzISczmbner9e8mHPPWb5JhrnmbPctqooxxpjoFRPqAMYYY0LLCoExxkQ5KwTGGBPlrBAY\nY0yUs0JgjDFRzgqBMcZEOSsExhgT5awQGGNMlLNCYIwxUS4u1AEC0atXLx06dGirli0pKSE5Oblt\nA7UByxUcyxWccM0F4ZstEnMtX748T1V7tzijqob9Iz09XVtr/vz5rV7WS5YrOJYrOOGaSzV8s0Vi\nLmCZBvAZa4eGjDEmylkhMMaYKGeFwBhjopwVAmOMiXJWCIwxJspZITDGmChnhcAYY6Jch7ihzBhj\n2tv6vUXsKyxjT2EpZ41No1un+JrXlm/LZ9P+EnbmH663TN9unfj25MGN1jVn9W4E4YJx/TzP3RpW\nCIwxxo+zH/u05vmMCf3501UTa6b/+6nFNc9FnH+rh3+fOqInQ3rW3glcWlHFba9+CcAZY84lKT7W\nw9St42khEJGtQBFQBVSqaoaIPAh8D9jvznavqr7nZQ5jTNuq8ik/fm0lW/JK6NstqcX5/2t8f4b0\n6MyT2RvJP1DKv3avwOdT7jjzWOJjhUfnrqOiSmvmV1U27S/hr99JZ0SfLkFle+yj9YjAHWceG9Ry\nH22r4JW/L/P72psrd1NSXuVmq22/eEJ/HncLxNe7Czn/zwu447WV9EpJrJmnpKyy5vnNL31B54Ta\nj92teSX06ZrIry46jriYGG566QuG9UwmJkZq5pnarSqon6M12mOPIFNV8xq0Paaqj7bDto0xHliz\n6xBvr9oNQEWVDxFpct7dBUfYfuAwEwal8uFXe53GPbkADOrRmaS4GN7P2cPovl1rllmbWwjAD1/5\nkg9/fHrAuQ4dqeDxeRsAuPnUYXRJim9hCYfPp7yxoZzkpIP07uIUtl4pieQVlwEwum8Xdh48UjN/\nfKyQGBfL5emDatqOTUvhzDFp7Co4Um9egDH9nJ8tv6SC/JKKmm1u2FfMhn3F/HvFLj7fks/m/SVs\n3l9SMz9ARXB1sFXs0JAxBoDC0gqu+etSDh4ub3HeI+W131I/uKP5D+rHPlrP4/M2sPVACeMHdmPV\nzkOA8+H48uJtiMCotC68f/tpNcvM+MtCVu4oYN3eIk595ONG6xw/MJWfnjOK7/19GRv3FTMgtRMi\nTlGqduYfP6F75wRevHESd/97NZv2F3PT1GF8sTWfNbsO1VufKhyphF+dM5orThrUcHMBiYuN4bnr\nMwKe3+dTht/rHAx5ceFWjlQ4v9PhvZPr/S6ys7NblScYonX3c9p65SJbgIOAAs+o6rPuoaEbgEJg\nGfATVT3oZ9mZwEyAtLS09KysrFZlKC4uJiUlpVXLeslyBcdyNa28SomPod638uLiYhI7J1Ppa3q5\npFgo94HP/QjIyaviLyvLmNA7luT4pr/hV9tT4uPcYfGc1Lf575P7D/t4Z3MFVT44bWAca/cdoUun\nRPp0FpbmOh9+GX1jmdgnrt4yz+eU0SMphoZJ9pT42HTIxwXD4nl3i/PtOj4GJrk5qg/BF5cry/ZW\nceHweOZsdubrkSTklyojUmNI69zgoklfBVePTSYloeWfva0s2l3J9kIfReWKAgeO+Ljx+ET6Jtdm\nO5r3WGZm5nJVbbE6eV0I+qvqbhHpA3wE/AhYB+ThFIffAP1U9abm1pORkaHLlvk/dteS7Oxspk2b\n1qplvWS5gmO5Gpv71R5mvrwcgOumDOE3M46vee2dufO5d2E5RXWOTzc0flAqa3YW1BQCgNgYYc2D\nZ9c7jt3WjvZ3tmTzAa56dkm9thtOGcqDFx1Xr63Kp0z41dya38H00X34+Jt9ALz+/ZOZNKxHm+by\nytHkEpGACoGnh4ZUdbf77z4RmQ1MUtWaU/Ei8ldgjpcZjIkkPp+yOa+YiirlD3PX17S/snRbvcsW\nv9hTSVFZJd89dZjfk7mzlu9k1Y4CYmOE+84bXXPly7BeyZ4WgbYwaWgPHvnvEygqrWT8oFQOFJcx\nbVSfRvPFxghPX5fO2txCundOYNqo3ry1cjcpiXFkDOkeguThy7O/uIgkAzGqWuQ+Pxv4tYj0U9Vc\nd7ZLgByvMhgTaWYt38ld/17dqN2ncN7jC+q1JcTG8OOzjiU5sfF/88LSSr7ZU8TYfl353unDPcvr\nhZgY4cqTGl+r78/UEb2YOqJXzfRNpw7zKlaH5mXpTwNmu8ct44BXVfUDEXlZRCbgHBraCnzfwwzG\ntJsDxWUs3HSAKl/jA/NVPoj1cx9/jAgDu3cGlO0Nbk7y582Vu+iRnMD/XlJ9GEjo2y2J/UVl9bab\nk/MVZ5+a4bcIAMw8fTjH9+9a7+oUE708KwSquhkY76f9Oq+2aUwoXfb0YrbklXi+nQvG9ePc45u/\nQzUpbx3jBqY2+XpKYhxnH9e3raOZDiq8DwYaE+bKq5THPlrPkYqqmiJw1tg07jt/TM08X2zN52ez\nVjMgtROvfHdyveWveGYx+4qca9UvSx/IbZkjWtxm/9RObfgTGGOFwJijsnJ/FU+u3EBiXO1xnzvP\nOpahvWq7GOiZksCfP97A7y4ZV68dnJuefvf+N4Bzl2rD141pD1YIjAlClU/54Ssr+Dq3kCE9O7N+\nVzkJcTGsefAcEuL8d+bbJSmeBXdN9/va9791DN//1jFeRjamRVYIjAnC2txCPvhqDwA9khPokiBc\nNnlYk0XAmI7ACoExdZRWVJHx0H8oLqus13FYtbLK2q4VZt96Cp988gnTpo1uz4jGtDkrBCbq+XxK\nudtHzfJtByl270Q9a2wfv52p7Tp4hGsmD262ozVjOhIrBCbqXfLUIlbtKGjU/rtLx4UgjTHtzwqB\niSoVVT7W5hZS5Xawc7i8ilU7Cjh7bBoTB3d32yo5Y0xaKGMa066sEJio8uLCrfz2vbWN2meePpyM\noT38LGFM5LNCYCLa/HX7yC0orZmes3o3A7t3qtdTZ3JCHOnWCZmJYlYITMTaV1jKjS980aj92imD\nyfTTW6Ux0coKgYkoqsrzn22hpKyKnN3OKFTPX5/B8QO61czT289locZEMysEJqLk7CrkoXedcwCd\nE2IZ268rU0f0Iql62CpjTCNWCEyH98XWfJ7O3oRPlf3uYOMAf/n2iWSOtkNAxrTE7os3Hd4rS7ax\ncFMeB0rKiXFv8jp1RC9OGmZXARkTCNsjMB1Ozq5D/PDVFZS7I7MfKC4nc3RvnrmuxaFZjTF+WCEw\nYa+iykdxuVJwuByAd9fksj3/MJedOBAREIQrJw0KcUpjOi4rBCasqSrnPb6AjfsOw8cf1bSPSuvC\n7y9vNACeMaYVrBCYsFJWWcWGvcUUllYAUHikgo37ipnSL5ZzMkbVzHeS3QVsTJvxtBCIyFagCKgC\nKlU1Q0R6AK8BQ3EGr79CVQ96mcN0HHfPWs2bK3c3aj9vWDzXTx0WgkTGRL722CPIVNW8OtM/B+ap\n6sMi8nN3+u52yGHCWElZJXO/3lNTBAakduIPVziHflIS48jb8GUo4xkT0UJxaOhiYJr7/CUgGysE\nUe+fn2+vuREM4PYzRjJleM+a6ewNoUhlTHTwuhAoMFdEFHhGVZ8F0lQ1F0BVc0XE7viJYks2H+DD\nr/awcGMeA1I78dr3p9AzOZFOCXYnsDHtRVTVu5WL9FfV3e6H/UfAj4C3VTW1zjwHVbVR148iMhOY\nCZCWlpaelZXVqgzFxcWkpKS0alkvWS7HbxYfYWuhj4RYmD4onstHJYRFrkBZruCFa7ZIzJWZmblc\nVVu+wUZV2+UBPAj8FFgH9HPb+gHrWlo2PT1dW2v+/PmtXtZL0Z7rF7PX6FXPLNYR976rv3336xbn\nj/bfV7DCNZdq+GaLxFzAMg3g89mzLiZEJFlEulQ/B84GcoC3gevd2a4H3vIqgwlPew6V8vKSbewt\nKiVjSA8uGt8/1JGMiWpeniNIA2a7A3zHAa+q6gci8gXwuojcDGwHLvcwgwkTFVU+LnlyITsPHqGy\nyjkc+fvLxtuAMMaEAc8KgapuBhrd+qmqB4AzvNquCT8Hisv4anchObsKOXNMHwakdqJ7cgITBqW2\nvLAxxnN2Z7Hx1MKNeVzz3NKa6Qf+6zgG9egcwkTGmIasEJg2p6ps3FdMwZEK3lq5i9gY4aEZx5PW\nNdGKgDFhyAqBaXPPfLqZh9//pmZ63MBuXD1pcAgTGWOaY4XAtJnySh9vrtzF059sqmn7x82TGZkW\nftdmG2NqWSEwbeY/a/dy16zVNdO3nzGSU0f2CmEiY0wgrBCYVvP5lMfnbWBfkTNO8Fe7D5EQF8OC\nuzLpmZxAbIyEOKExJhBWCEyrrd1TyOPzNtCtUzwJcc69iTMm9Ceta1KIkxljgmGFwATtr59u5t01\nuRw64gwe8/ZtUxnSMznEqYwxrWWFwATt+c+2EBsjHNMnhakjejLYLgk1pkOzQmACUlZZxRVPL2ZX\nwRHyisv5xQVj+O5pw0MdyxjTBqwQmBYdOlzBqp0FrNp5iOmj+zC4R2cunjAg1LGMMW3ECoFp1qfr\n9/Odv31eM/3QjOPpn9ophImMMW3NCoFp0o78w7yxYicJsTH88sIx9O3WyYqAMRHICoFp0tV/XcLO\ng0c4aWh3rjt5aKjjGGM8YoXANLJyRwEL1u9n58EjfO+0Yfxg2ohQRzLGeMgKgWnkZ/9axYZ9xcTH\nCpelD6JHsv9xhI0xkcEKgamRX1LOCzllbNxfyY+mj+C26SNIjIsNdSxjjMc8G7PYdDwf5Ozhk52V\nDOuVzHnH97MiYEyUaHKPQER6A71V9esG7ccB+1R1v9fhTPtZtaOAe2evoXMczLvzW7hjTRtjokBz\newRPAL39tA8EHg90AyISKyJfisgcd/pFEdkiIivdx4TgIhsvzP5yFwAXDo+3ImBMlGnuHMEJqvpJ\nw0ZV/VBE/hDENm4H1gJd67T9TFVnBbEO45EV2w/yvZeWUXCkgknDenD+8LJQRzLGtLPm9gjiW/la\nDREZCFwAPBdMKNM+8orLeHvlbgqOVHD9yUO565xRoY5kjAmB5grBBhE5v2GjiJwHbA5w/X8C7gJ8\nDdp/KyKrReQxEUkMcF2mDa3aUUDGQ//hxUVbGdOvC/f/11gyhvYIdSxjTAiIqvp/QeRYYA6wCFju\nNmcAJwMXqur6ZlcsciFwvqreKiLTgJ+q6oUi0g/YAyQAzwKbVPXXfpafCcwESEtLS8/KymrFjwfF\nxcWkpITfmLmhzLX/sI93Nlfw6c5Krh+bwKgesfRPiQl5ruZYruCEay4I32yRmCszM3O5qma0OKOq\nNvkAEoEbgT+4j5uApOaWqbPs74CdwFacD/7DwD8azDMNmNPSutLT07W15s+f3+plvRTKXBf+eYEO\nuXuOnvmH7Eav2e8rOJYreOGaLRJzAcs0gM/rZm8oU9Uy4IXWVCJVvQe4B6DOHsG1ItJPVXPFuTRl\nBpDTmvWb4G3YW8SnG/L4OreQb08ezM/OtnMCxpjm7yM4DbjNnXxS/VxB1EqvuPcoCLASuKWN1mta\ncP9bX7F48wFE4MJx/ehuXUcYY2j+8tE/AKe7zxcAJ7V2I6qaDWS7z6e3dj2mdY6UV/HIB9/w5Y6D\nXD1pML+4YAzJida7iDHG0dynQTbwK5wrfha3SxrjiU837OfFRVsZkNqJi8b3tyJgjKmnyU8EVb1L\nREYBsdqgmwnTcTy3YDPPf7aFuBhh3k++RVK89R9kjKmvpZPF69oriPHGCwu3EiPC/5wx0oqAMcYv\n6300Qq3fW8QJD3zIroIj3HDKUP7njJGhjmSMCVNWCCJQeaWPD3P2UFRWyW2ZI7gsfWCoIxljwlhz\nl4/eA3ygql+2Yx7TBm568Qs+25hHv25J/NT6DzLGtKC5cwRbgNtFZDywCngfmKuqB9slmQlalU9Z\nsvkAn2/J54zRfbj9TDscZIxpWXNXDWUBWQAiMhE4F3hDRGKB/+DsLXzeLilNQD78ag+3vrICgKsm\nDWbcwNQQJzLGdAQBXVDuHh76EvidiHQFzgK+C1ghCBP7ikp58O2vSIyLYfatUxnTr0uoIxljOoig\n7yxS1ULg3+7DhIkXFm5lX1EZZ41NY2z/ri0vYIwxLrtqKAIs2pTHU9mbGN23C89elx7qOMaYDsYK\nQQR49lNnnKBrJg+28YaNMUFrVSEQkdFtHcS0zhPzNpC9bj9XZAzkupOHhjqOMaYDau0ewdw2TWFa\n7ZWl2wG4dsqQECcxxnRUzd1Q9uemXgLsusQwcOUzi9lTWMovLhhjl4oaY1qtuauGbgR+ApT5ee1q\nb+KYQG3NK2HplnwGpHbi0hOtCwljTOs1Vwi+AHJUdVHDF0TkQc8SmRaVlFUy7dFsAB69fDw9bKQx\nY8xRaK4QXAaU+ntBVYd5E8e0RFV55pNNAFw7ZTBThvcIcSJjTEfXXBcT+e0ZxARmza5D/PnjjQDc\nedYou1zUGHPU7D6CDqSotIKf/msVAG/9cKodEjLGtAnPC4GIxIrIlyIyx50eJiJLRWSDiLwmIvZp\nFqB3VuWyfm8xI/qkMH6QXSVkjGkbARcCEUlu5TZuB9bWmX4EeExVRwIHgZtbud6osmFvEffOXkNq\n53g++vHpoY5jjIkgLRYCETlFRL7G/TAXkfEi8mQgKxeRgcAFwHPutADTgVnuLC8BM1qRO+pU3zh2\nZcYgOy9gjGlToqrNzyCyFOcKordVdaLblqOqx7e4cpFZwO+ALsBPgRuAJao6wn19EPC+v3WJyExg\nJkBaWlp6VlZWED9WreLiYlJSUlq1rJeCyVVQ5uOO+UcY1T2GeyZ3Cptc7clyBSdcc0H4ZovEXJmZ\nmctVNaPFGVW12Qew1P33yzptqwJY7kLgSff5NGAO0BvYWGeeQcCaltaVnp6urTV//vxWL+ulYHL9\na9kOHXL3HP374q3eBXJFwu+rPVmu4IVrtkjMBSzTFj5fVTWg8Qh2iMgpgLondv+H+sf8mzIVuEhE\nzgeSgK7An4BUEYlT1UpgILA7gHVFrc825PHTf62iU3ws10waHOo4xpgIFMjJ4luAHwIDgJ3ABHe6\nWap6j6oOVNWhwFXAx6p6DTAf51ATwPXAW63IHTVeW7YDgMeunEBMjJ0bMMa0vWb3CNzxia9zP8Db\nyt1Alog8hDP85fNtuO6IsiP/MO+s2s24gd049/i+oY5jjIlQzRYCVa0SkYuBx45mI6qaDWS7zzcD\nk45mfdGgssrHw+9/A8AdZ44McRpjTCQL5BzBQhH5f8BrQEl1o6qu8CyV4T9r9/HumlwS4mI4fWTv\nUMcxxkSwQArBKe6/v67Tpjj3AxgP5BWX8cDbOQAs+vl04mKtJxBjjHdaLASqmtkeQUytfyzZxt7C\nMjJH9aZXSmKo4xhjIlyLhUBE7vfXrqq/9tdujs6SzQf40382MLJPCn+74aRQxzHGRIFADg2V1Hme\nhHOjWCD3EZhW+NtnWwC4ZvJg60rCGNMuAjk09Ie60yLyKPC2Z4mi2JPZG5n79V5mTOjPDVNt7B9j\nTPtozVnIzsDwtg4S7Xw+5eXF2wCsCBhj2lUg5wjW4FwlBBCL01/Qb7wMFY0uf2YxuYdKuee80Uyw\nsQaMMe0okHMEF9Z5XgnsdfsJMm1kR/5hlm87yPBeyVx1kvUnZIxpX4EcGnpIVbe5j12qWikiL3ue\nLEqUVlRx2v/NB+DBi46jW+f4ECcyxkSbQArBcXUnRCQOSPcmTvR5PycXgKtOGsSpI3qFOI0xJho1\nWQhE5B4RKQLGiUih+ygC9mI9hraZp7I3AXDn2cda76LGmJBoshCo6u+AbsDfVbWr++iiqj1V9Z72\nixi51uYWsn5vMd+ePJg+XZJCHccYE6WaPTSkqj5gfDtliTovLHRuHrt04oAQJzHGRLNAzhEsERHr\n66CNfb6nkteX7eRbx/YmY2iPUMcxxkSxQC4fzQRuEZGtON1NCKCqOs7LYJHu/S0VgHOS2BhjQimQ\nQnCe5ymizNLNB9hyyMf3Tx/OeSf0C3UcY0yUa/HQkKpuAwYB093nhwNZzjTt7VW7Afjv9IEhTmKM\nMQF8oIvIAzjjDFdfKRQP/MPLUJFsR/5hXlm6ndE9Yjg2rUuo4xhjTEDf7C8BLsLtjlpVdwMtfoKJ\nSJKIfC4iq0TkKxH5ldv+oohsEZGV7mPC0fwAHc1bK3cBMH2w3UFsjAkPgZwjKFdVFREFEJHkANdd\nhnM4qVhE4oHPROR997WfqeqsVuTt0CqrfDw6dz2JcTFkpMWGOo4xxgCB7RG8LiLPAKki8j3gP8Bf\nW1pIHcXuZLz70GYWiXjPfLoZgFu+dQwxNuiMMSZMBHKy+FFgFvBv4FjgflV9IpCVi0isiKwE9gEf\nqepS96XfishqEXlMRKJiUN7N+4v5/YfrALjCLhk1xoQRUW35S7qI9AUm4Xyj/0JV9wS1EZFUYDbw\nI+AAsAdIAJ4FNvkb/1hEZgIzAdLS0tKzsrKC2WSN4uJiUlJSWrVsW3pmdSmLd1dxz6QkRvWIDZtc\nDVmu4Fiu4IVrtkjMlZmZuVxVM1qcUVWbfQDfBbYDLwIvAVuBm1pazs96HgB+2qBtGjCnpWXT09O1\ntebPn9/qZdtKVZVPx/zyfR1y9xytrPKpanjk8sdyBcdyBS9cs0ViLmCZBvD5HMjJ4p8BE1X1AICI\n9AQWAX9rbiER6Q1UqGqBiHQCzgQeEZF+qporzsjsM4CcADJ0aP/8YjuHy6t49PLxxFoPo8aYMBNI\nIdgJFNWZLgJ2BLBcP+AlEYnFORfxuqrOEZGP3SIhwErgliAzdzjvuDeQnT7SxhswxoSfQArBLmCp\niLyFc47gYuBzEbkTQFX/6G8hVV0NTPTTPr31cTue5dvyWbI5n+tPHkKfrtbVtDEm/ARSCDa5j2rV\ng9LYbbEBeGOFcwPZlTYWsTEmTLVYCFT1V+0RJBI9t2AzryzdzomDUxnbv2uo4xhjjF8tFgIRyQDu\nA4bUnV+tG+pmqSr/WLINgIdmnBDiNMYY07RADg29gnPl0BrA522cyPHxN/vYeuAwV2YMsr0BY0xY\nC6QQ7FfVtz1PEmFeXLQVgLvPGx3aIMYY04JACsEDIvIcMA+nIzkAVPUNz1J1cLmHjrBgQx5Denam\nR3JCqOMYY0yzAikENwKjcTqNqz40pIAVgiY884nTudzDl9ppFGNM+AukEIxXVTvbGSBVrRmB7MQh\nqSFOY4wxLQukG+olIjLW8yQR4o0Vu8gvKec3M44nMc7GHDDGhL9A9ghOBa4XkS045wgEZ7gBO+7h\nxzurnb2BM0b3CXESY4wJTCCF4FzPU0SIlTsKyF63nysyBtI/tVOo4xhjTEACubN4m4icCoxU1Rfc\nDuPCr9PuMFDdudy1U4aEOImJWFWVUFUGvkqoLKv/WqfuENtgLOyKUqgsBfU5ywDEd4bEBv+FfT44\nnOdd7k7doeIw+Kq824ZptUDuLH4AyABGAS/gXD30D2Cqt9E6lrziMp7/bAuj0rowbqCdJDYe+cel\nsOUT/68N+xZcX+eWH58P/nQClOyrP198Z7gjB5J71ra9+2NY/mKbx/Wn13H34AxFYsJFIIeGLsHp\nRXQFgKruFhHrcK6BD79yBm27YerQ0AYxHZMqFOVCWTGUFTY9T90icO7DtXsA6+fCpo9hx+cgMXQp\nXAffHKpfBC74AxzaCZ89Bmv+BQPrDFy18WMYkAETrm77n+2De529GICYOHrvXwg7lwW/nrgkZ+8G\ngb7HQ5yfUW5VYW8OdO4JXfsfVewm5W+BsiLQKpAY55GUCsV7nde7D4Wqcqdt39cQEweJXeBIAXTu\n7sx/+ICTFaD3aGcPzVcFe9bU7rlV/9gVxXgtkEJQrqoqIgogIskeZ+qQHn7vGwBmTBgQ4iSmQ/ry\nZXj7R4HPn5QKU35QO52SBhs+hOfPAiAd3K9uruQ+cNJ3nUKz+En44O7G65w805mnrRXsgIV/cp4P\nOYW0LZ/Cc2cc3TpPvRPOfKBx+7r3IOvbkNgV7t4GMYFcGBmEnH/DrJuanyc2wSkEgRp3JVz6rPMe\neOf2Ri93PcHPz9nGAikEr4vIM0CqiHwPuAl4zttYHcsXW/MpKqvkuilD6JRgl4yaOnw+yJkFFUec\n4/RV5c4HRUOLnqh9fulfnWPq/sQlQlI351HXqPPh+ndqzhusXr2acePGOevpOqD223NiCsycD4W7\n6y8fEwuDT2nlD9mC6b+E0Rc4OWJiWT33FSdbMFa+Al/NhlEXQME25wM/1U/X7uved/4tK3SKT1O/\nRz/67V4Py7Y0P1Nzh89OvRNK9jsf6NUGnwzbF/uf/5pZsOjPsDkblr3g7KUl94YZT9WbrWjL4YDy\nH41AThY/KiJnAYU45wnuV9WPPE/WgTw5fyMA355sYw6YBlb9E966NfD5uw6AEy4HCXJI05hYGHZ6\nzWT+rngYOc3/vGnHOY/2EhsHgybVTOb3TG86W1M6dXcKwSk/cg6Bffp/MOcO//N2H+ocApsXXA/6\nowDWB7HAgAzYvw7K3QEcM26Ewtz6hWDsxc4hoB1LYOK1kDMbKkqcvYCRZ0HhLmcvoPpnqW6vo2JX\ndlA/R2sEcrL4EVW9G/jIT1vUU1XW7Crk3OP6Mqaf9TIaVb55jxNWPwol70DBdkj1c7XYjqW1z6f/\nAj5+yPn2foGfgf0Skp1j/sEWgWgwMAN+mef8fgZPgZNurj3G3lByLygvcfbCgrBo8SJOOTmAvaKk\nbs5enYibQZ0ropK6OXspd21xzgtUHIGUPjBppnPcPzYBzn3EOUcQ545WmH6D836ovpoqJTT3HwVy\naOgsoOGH/nl+2qLS7kOl5BWXccqIni3PbCJL1tX0BMhfXtvWqYf/ec/5Xzj+Mucb4Rn3Q9d+7ZEw\nslSfGBeBLn2bn7dTqvMIQnnardJuAAAVKUlEQVRiz9b/XWLrHKrr7L4HktwvhhLr7LFB48t2IWQf\n/nU1WQhE5AfArcBwEVld56UuwEKvg3UUj7zvnCQeb5eMRh6fD978gXNMuiF/30aHnAo3vtv8Om9d\n1DbZjGlDze0RvAq8D/wO+Hmd9iJVzW9pxSKSBHwKJLrbmaWqD4jIMCAL6IFzXcN1qhrEKfbwoaos\n3JhHYlyMDT4TifI3weos6HNc/Wvuq404i7yDhfS6+DfOycxT/qf9MxrTBposBKp6CDgEtPbC4jJg\nuqoWi0g88JmIvA/cCTymqlki8jRwM/BUcysKV68s3c6BknJ+c/FxxMe28WVqpn0c2ASvXuH/eHJl\nqfPvpc861637kZOdzbTBk2HwZA9DGuOtQM4RtIqqKlB9J0S8+1BgOvBtt/0l4EE6YCHw+ZS3Vu4C\n4LwT7Hhv2Ks44v/Dfu07cGAjjLvKubqloZS+0Mc63zWRTbSpM+9tsXKRWGA5MAL4C/B7YImqjnBf\nHwS8r6qNvm6JyExgJkBaWlp6VlZWqzIUFxeTktL2XSM9v6aMBbsqOW1AHDef4OcOxxDlOlqRmCuh\nLJ/JS79PrM//EciyhB4sPvlvrbpaJxJ/X14L12yRmCszM3O5qma0NJ9newQAqloFTBCRVGA2MMbf\nbE0s+yzwLEBGRoZOmzatVRmys7Np7bLNuWfxPLomwSPXnd6qnka9ynW0OmyuA5ugvNi59T8mvv6d\nnTs2ga/cueEnJa3Roon9JzqHd7zIFSLhmgvCN1s05/K0EFRT1QIRyQam4NyhHKeqlcBAYHezC4eh\nvYWl5B4q5f4Lx1p30+FgTw483UIfiAkpMO0eiLMxpI1pyLNC4HZXXeEWgU7AmcAjwHzgMpwrh64H\n3vIqg1f+scS5nHD8oG4tzGnaVFUFrHwVivfVvx47d1XjeS96AroPq53u2t+KgDFN8HKPoB/wknue\nIAZ4XVXniMjXQJaIPAR8CTzvYYY2p6r88/MdABzX3wpBu9rwEbzTwiWanXs6PTuOu9J/75TGmEa8\nvGpoNU731Q3bNwOTGi/RMcz9ei95xWXcd/4YkuKtgzlPrcqCLQsAGLUnF9bXuX3lu/Og5zG10/HJ\n9o3fmFZql3MEkeT1L5y9gbPGNj7paNqQKnxwj3M4KKkb3ctKITHJ6U65/0TnEWOF2Ji2YIUgCJv2\nFzPvm31cOnEAQ3vZsAxtpqwY3rzFGbijmq8KjuQ7nbOddDNLwvSKDmMigd0OG4R3V+cCcPFEG3ym\nTW1b5NzYVVboFIDqnhiPmQ7HnhvabMZEAdsjCNChIxX88aP1jOiTwreO7R3qOB1Xzhsw60any15x\nv4dUljnPb3jPf++MxhhPWSEI0IptBwGYMcGjcVAjla+qdixXcIoAOB/+J36ntr3PWCsCxoSIFYIA\nrdxRQIzAjVOHtTyzqTX7FljzeuP2jJvh3P9t/zzGmEasEATA51Nmf7mLkX26kJxov7Jm5W2E/M21\n05vnw5CpMO6K2raYeDjhsvbPZozxyz7VAvDphv1szz/MpXaSuHmq8OL59Q8FgTNEY/oNIYlkjGmZ\nFYIAPP/ZFgB+eaF1R1zPjs/rd+9QXuwUgdN+6ozDCk7Xzmn++/I3xoQHKwQt2F9UxoINefTrlkT3\nZLtztYYqZF0DJfvqt0uscxio96jQ5DLGBM0KQQuezN4IwB+uGB/iJGFgy6eQ82/neVWFUwTO+jWM\n/3btPPFJkNglNPmMMa1ihaAF761xbiI7cXD3ECcJA/P/F3atgE6pznT3oTDmvyDF7qswpiOzQtCM\nN7/cxd7CMu6/cGz0dTBXkgdv3eZ086A+p233Spj0PTj3d6HNZoxpU1YImjHH7VLinOP7hjhJCHz+\nLKx/33k+IN25E/iYTBh/VWhzGWPanBWCJqzaUcB/1u7l0okDGBANo5BVHIGnpnJawQ74LBaqympf\nuzoLUvqELpsxxlNWCJqQvW4/ANefMjS0QbxUWQ6lh5znuSshfxMHep9Kn5EnOm1FuTDkFCsCxkQ4\nKwRNWL2zgBF9Uhg/KDXUUbzz/JmNhnncOOIm+pz93yEKZIwJBSsEfpRWVDnjDpwYYXcSlxyAPe4H\nf0WpUwSOu8TpAgKg20DKc6PgMJgxph4rBH48t8DpKydjSI8QJ2ljb98G696r3zb5BzB4cu10bna7\nRjLGhJ5nhUBEBgF/B/oCPuBZVX1cRB4Evgfsd2e9V1Xf87+W9ldZ5eOZTzeTkhjHFRkDQx3n6K18\nFYr2QEwcbFsIoy+EU37kvBbfGfqNC20+Y0zIeblHUAn8RFVXiEgXYLmIfOS+9piqPurhtlttwYY8\nikorufTEAcTFdvAB3HYuhzd/UL/thMtg8JTQ5DHGhCXPCoGq5gK57vMiEVkLhP1B90c++AaAh2Z0\n0I7Sdi6HL/7q9AV0aEdt+z27nMHe4+0cgDGmvnb5yisiQ4GJwFK36TYRWS0ifxORsOm7YV9hKd/s\nKeLYtBQ6J3TQ0yeL/5/TH9D2xXBop9M27V5n9C8rAsYYP0RVvd2ASArwCfBbVX1DRNKAPECB3wD9\nVPUmP8vNBGYCpKWlpWdlZbVq+8XFxaSkBDYE4jOrSlmcW8V9k5MY2d3bLiWCydWUmKpSTljzW+Ir\nCmvaOh3ZzYGeGXx93N0hy+UFyxWccM0F4ZstEnNlZmYuV9WMFmdUVc8eQDzwIXBnE68PBXJaWk96\nerq21vz58wOed9JvP9Jj73tPyyurWr29QAWTq0kbP1Z9oKvqCxeo/vPbtY8tn4U2lwcsV3DCNZdq\n+GaLxFzAMg3gs9rLq4YEeB5Yq6p/rNPeT53zBwCXADleZQhG1ufb2VtYxi8vHEt8uJ4kVoUXL4A9\na5zpqnLn3ytfhk5hc4TNGNPBeHkgfCpwHbBGRFa6bfcCV4vIBJxDQ1uB73uYIWAff+MMsDJjQv8Q\nJ6mjvMTpBbRaSZ5zCegx06H3aKet10grAsaYo+LlVUOfAeLnpbC5Z6CuVTsLmDGhPz1TEkMdxaEK\nT58G+Zsavzb9F06PoMYY0wY66KUxbStn1yH2FpaFR79CRw7C9iVQWugUgQnXOh2/VeuUCv1PDF0+\nY0zEsUIA/N+H6wA4aWgYdCkx79ew7G+105O/b3f/GmM8FfWFIK+4jE/X72fqiJ4cP6Bb+2780E74\n4vn61/evnwuDT4ZzH3bG/u15TPtmMsZEnagvBK8s2Q7AlScNbv+Nv3A+FGxr3D71dug/of3zGGOi\nUlQXgtKKKv788QaO6Z3MRePb4WqhkgPw5i3QuRej9u6rLQLHTIdrZtXOFxNl4yMbY0IqqgvB4k0H\nqPIpJx/Ts302OPc+2DAXgO6JPSGxK5QVwtkP2Ye/MSZkoroQ3P+2cy/bPeeN8WYDZUXw2rVwpMCZ\nznVvp5AYlpz8N6ZNm+bNdo0xJghhegut9/YXlbEj/wgTBqWSnOhRPdy2CDZnQ0IypKTByHOg6wC4\n6UNvtmeMMa0QtXsEv3jT6abh3vPbYG9g79fw1MnuhNQe5vFVgcTANf9yikFdm7KPfrvGGNMGorIQ\nFJVW8PE3++iZnED6kCC6Z1CF0kNQXux0/1Dtsz/WPo9LgpNvrZ3uPaZxETDGmDASlYXgkicXUVGl\n/Gj6CGJj/PWC0YTlL8KcO5qfZ/JMOOP+o8pnjDHtKeoKwby1e9m4r5gpw3tw1SQ/9w7kroa89f4X\n/uSR2ucXPVH/m35sAsQmwjGZbRvYGGM8FlWFoLiskptfWgbAXeeOJim+wSWbZcXwzGktr6jrADjx\nOx4kNMaY9hfZhWDjPAZvmw0LVgDw9dZ8fhC7n7OPS2PitnXQ8Kbeoj21z2dmQ7yfY/uJKXbM3xgT\nUSK7EKx7j+FbXoYtzuQkYFI8sN59NCWpG/Sf6H0+Y4wJA5FdCM59mE+TzsHXdzRvrdzNezm53Dpt\nBLdPH9n0MjFxEBvZvxZjjKkrsj/xYuPxxSbw+Cc7WL+niGFpPTnrhCEQnxTqZMYYEzYiuhA8MW8D\n/1x8mL2H4aapQ7nvgrGhjmSMMWEnogtB7y6J9E+OIf2YNK7IGBTqOMYYE5Y8KwQiMgj4O9AX8AHP\nqurjItIDeA0YijN4/RWqetCLDFdNGkzfw5uZNs1O/BpjTFO87HSuEviJqo4BpgA/FJGxwM+Beao6\nEpjnThtjjAkRzwqBquaq6gr3eRGwFhgAXAy85M72EjDDqwzGGGNa1i7dUIvIUGAisBRIU9VccIoF\n0Kc9MhhjjPFPVNXbDYikAJ8Av1XVN0SkQFVT67x+UFUbdQEqIjOBmQBpaWnpWVlZrdp+cXExKSkp\nrQvvIcsVHMsVnHDNBeGbLRJzZWZmLlfVjBZnVFXPHkA88CFwZ522dUA/93k/YF1L60lPT9fWmj9/\nfquX9ZLlCo7lCk645lIN32yRmAtYpgF8Vnt2aEhEBHgeWKuqdTrs523gevf59cBbXmUwxhjTMi/v\nI5gKXAesERF3sF7uBR4GXheRm4HtwOUeZjDGGNMCzwqBqn4GNDXqyxlebdcYY0xwPD9Z3BZEZD+N\nO40OVC8grw3jtBXLFRzLFZxwzQXhmy0Scw1R1d4tzdQhCsHREJFlGshZ83ZmuYJjuYITrrkgfLNF\nc652uY/AGGNM+LJCYIwxUS4aCsGzoQ7QBMsVHMsVnHDNBeGbLWpzRfw5AmOMMc2Lhj0CY4wxzYjY\nQiAi54rIOhHZKCLt3tW1iPxNRPaJSE6dth4i8pGIbHD/7e62i4j82c26WkRO9CjTIBGZLyJrReQr\nEbk9HHK520oSkc9FZJWb7Vdu+zARWepme01EEtz2RHd6o/v6UA+zxYrIlyIyJ1wyudvbKiJrRGSl\niCxz28Lhb5kqIrNE5Bv3vXZyqHOJyCj391T9KBSRO0Kdy93Wj933fI6I/NP9v9C+77FA+qHoaA8g\nFtgEDAcSgFXA2HbOcDpwIpBTp+3/gJ+7z38OPOI+Px94H+cGvCnAUo8y9QNOdJ93AdYDY0Ody92W\nACnu83icnmqnAK8DV7ntTwM/cJ/fCjztPr8KeM3DbHcCrwJz3OmQZ3K3sRXo1aAtHP6WLwHfdZ8n\nAKnhkKtOvlhgDzAk1LlwuubfAnSq8966ob3fY57+wkP1AE4GPqwzfQ9wTwhyDKV+IfDb4R7wDHC1\nv/k8zvcWcFYY5uoMrAAm49xIE9fw74rTmeHJ7vM4dz7xIMtAnAGUpgNz3A+GkGaqk20rjQtBSP+W\nQFf3g03CKVeDLGcDC8MhF04h2AH0cN8zc4Bz2vs9FqmHhqp/udV2um2h1tRYDO2eVwIbI6Jdc7mH\nYFYC+4CPcPbqClS10s/2a7K5rx8CenoQ60/AXTjDreJuI9SZqikwV0SWi9NtO4T+bzkc2A+84B5O\ne05EksMgV11XAf90n4c0l6ruAh7F6XctF+c9s5x2fo9FaiHw18dROF8e1a55xRkj4t/AHapa2Nys\nfto8y6WqVao6Aedb+CRgTDPb9zybiFwI7FPV5XWbQ5mpgamqeiJwHs5QsKc3M297ZYvDOST6lKpO\nBEpofjja9n7vJwAXAf9qaVY/bW2eyz0ncTEwDOgPJOP8PZvatie5IrUQ7AQG1ZkeCOwOUZa69opI\nPwD3331ue7vlFZF4nCLwiqq+ES656lLVAiAb59hsqohUd45Yd/s12dzXuwH5bRxlKnCRiGwFsnAO\nD/0pxJlqqOpu9999wGyc4hnqv+VOYKeqLnWnZ+EUhlDnqnYesEJV97rToc51JrBFVferagXwBnAK\n7fwei9RC8AUw0j3znoCzK/h2iDNB02MxvA18x71SYQpwqHp3tS2JBD1GRLvkcrP1FpFU93knnP8g\na4H5wGVNZKvOfBnwsboHTtuKqt6jqgNVdSjOe+hjVb0mlJmqiUiyiHSpfo5z3DuHEP8tVXUPsENE\nRrlNZwBfhzpXHVdTe1ioevuhzLUdmCIind3/n9W/r/Z9j3l5UiaUD5yz/utxjjPfF4Lt/xPnmF8F\nThW/GedY3jxgg/tvD3deAf7iZl0DZHiU6VSc3cjVwEr3cX6oc7nbGgd86WbLAe5324cDnwMbcXbn\nE932JHd6o/v6cI//ntOovWoo5JncDKvcx1fV7/Ew+VtOAJa5f8s3ge5hkqszcADoVqctHHL9CvjG\nfd+/DCS293vM7iw2xpgoF6mHhowxxgTICoExxkQ5KwTGGBPlrBAYY0yUs0JgjDFRzgqB6VDcni1v\nrTPdX0RmebSteBFZ3vKc7UtEhkqdXm2NOVpWCExHk4rTAyPg3F2rqpc1M//ROBVY5NG6jQkbVghM\nR/MwcIzbp/zv6347FpEbRORNEXlHRLaIyG0icqfb+dkSEenhzneMiHzgdta2QERGN7Gtc3G6Iq7h\ndoz3ott3/BoR+XFz6xSRNBGZLc44C6tE5BS3/U53HTkicofbNlSc/vv/Kk7/9HPdu6wRkXR3+cXA\nD+vkOU6ccRxWitNv/si2/GWbKOHV3XL2sIcXDxp37V0zjdOP+0acsRZ64/TMeIv72mM4neyBcwfp\nSPf5ZJzb9P1t63Ogc4O2dOCjOtOpza0TeK3OdmNx+oZJx7lbNRlIwbkzeKL7s1QCE9z5XweudZ+v\nBr7lPv99nZ/5CeAa93kCbr/29rBHMI/qTo2MiRTzVbUIKBKRQ8A7bvsaYJzb8+opwL+crl0A55b+\nekSkP5CvqocbvLQZGC4iTwDv4nQD3dw6pwPfAad3VeCQiJwKzFbVEndbbwCn4fQjs0VVV7rLLgeG\nikg3nILzidv+MrU9VC4G7hORgcAbqroh0F+UMdWsEJhIU1bnua/OtA/n/R6D09f7hBbWcx7OICD1\nqOpBERmPM3jID4ErgDsCXGc1f10JV6ubvwro5M7vty8YVX1VRJYCFwAfish3VfXjAHMYA9g5AtPx\nFOEc+mkVdcZf2CIil0PN2LTj/cza6PyAO38vIEZV/w38Emfoz+bWOQ/4gdseKyJdgU+BGW6Pk8nA\nJcCCZjIXULsnAXBNnTzDgc2q+mecPYpxgf4ujKlmhcB0KKp6AFjonmT9fStXcw1ws4hU99x5cd0X\nRSQW53j/N36WHQBkizOS2os4w6A2t87bgUwRWYNzqOc4VV3hLvs5zghxz6nqly1kvhH4i3uy+Eid\n9iuBHDfPaODvLazHmEas91FjGnC/eV+rqreEOosx7cEKgTHGRDk7NGSMMVHOCoExxkQ5KwTGGBPl\nrBAYY0yUs0JgjDFRzgqBMcZEOSsExhgT5f4/U1b5M+iuv3EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10df40588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = pd.read_csv('Step_Test_Data.csv').set_index('Time')[1:]\n",
    "t = data.index\n",
    "T1 = data['T1'].values\n",
    "T2 = data['T2'].values\n",
    "\n",
    "plt.plot(t, T1, t, T2)\n",
    "plt.xlabel('time / seconds')\n",
    "plt.ylabel('temperture 1 / °C')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.4 Estimating $U_a$ and $U_b$](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.4-Estimating-$U_a$-and-$U_b$)",
     "section": "2.4.4 Estimating $U_a$ and $U_b$"
    }
   },
   "source": [
    "**Exercise:** Roughly speaking, the has an increase of 35 degrees and T2 and increase of 10 degrees due to a setting u1 to 50% of the maximum power of 4 watts.  Use that information and the above equations to estimate values for $U_a$ and $U_b$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.5 Simulation with odeint](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.5-Simulation-with-odeint)",
     "section": "2.4.5 Simulation with odeint"
    }
   },
   "source": [
    "## 2.4.5 Simulation with odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.5 Simulation with odeint](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.5-Simulation-with-odeint)",
     "section": "2.4.5 Simulation with odeint"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VPW9//HXZyYrWQhLiGEHQQRB\nkQTFrRK3KnqrvVWv2lq7ib29rW3totZfd72tt962tvXWaq1Vq1KrFS1tbakSF0CEsMkiyBJ2AiQh\n+zYzn98f5wQChmQmycyZzHyej8d5nHO+c+acN5NwPpmzfI+oKsYYY5KXz+sAxhhjvGWFwBhjkpwV\nAmOMSXJWCIwxJslZITDGmCRnhcAYY5KcFQJjjElyVgiMMSbJWSEwxpgkl+J1gHAMHTpUx44d26P3\nNjQ0kJWV1beB+oDliozliozliky85oLeZSsrKzukqvndLqiqcT8UFRVpTy1atKjH740myxUZyxUZ\nyxWZeM2l2rtswAoNYx9rh4aMMSbJWSEwxpgkZ4XAGGOSnBUCY4xJclYIjDEmyVkhMMaYJGeFwBhj\nkly/uKHMGGN6S1UJhJRgyB0HlUAoRDCkVDaF2FXVeOS1kCqBoDt23xNSdxxSgtqxDYIh53r89nZV\njryuCkF1pkMhJaQ404o7f7RN3WntMD+8LRT1z8YKgTGmzwRDSnNb0BkCIZrbgrS0hWgJBGkJhFh7\nMEDL+v20BEK0HhmCR+eD7hAI0RYM0RbQI21t7W1BdcdHpwMhJeDOB0IhAm57MKS0uTvyYKib57O/\nvig2H1IEROCrM9Kjvh0rBMYkoUAwRENLkPrWAPXNAepb3KE5QENrgMaWAA2tQZpag+68M25uC9LY\nGqSpzXmtfbq5NUhzIEhbsJudLUBZWZcvp6X4SPf7SE3xkeb3kZoipPrdab+PVL8zn5WeQopPSHHb\nUnw+UvxCavvY7yPFJ/jdNr9PSPULPt/R+RS/4PcJWzZvZvLkU53XxVmX3wd+d+wTZzm/OO/3++S4\nNo60tbf7xHmfz532i4A79vsEEbfdJwjOOoSj7xMBEaG0tLRvfuhdsEJgTD/VEghyuLGN6sZWqhva\nqGlqpaapjZqmNtZubmVh9bvUNLVR1xygttkZ17njxtZgWNvwCWSlp5CVlsKAND+ZaX4yU/1kp6eQ\nn53OgDQ/GanOkJnmJyPFT0aqz21zxukpPtLd8Ya1a5h1VrHTluInLcXn7PjdcYq7g4y10qbtzC4e\nFfPtxouoFgIRKQfqgCAQUNViEfkecCtw0F3sW6r6t2jmMKY/CIaUyoYWDtS2cKi+hcr6ViobnPGh\nDtPOjr+Vhi525j6BvIr9DMxMJTcjhZyMVE7KzSA3I5Ucdz47I4Wc9BRnR5/uJycj5chOPyvd2fGn\np/j6dMfcusvP1BED+2x9pm/E4htBiaoeOq7tZ6r6QAy2bUxcaGwNsK+mmX2Hm9lX0+RM1zRzoLaZ\nirpmDtS2UNnQ2ulx7LQUH0Oz0hiSnc7grDQmDMtm0IA0Bg1IZVBW2pHpgQNSGZiZSt6ANJYveZOS\nkhIP/qWmP7JDQ8b0gabWIDuqGthd1cTu6kZ2VTvj3dVN7K5uoqap7QPvGZKVRkFuBsNy05lSmMuw\nHGd6WE46Q7OdYUh2GtnpKRH/Ve7F4RXTf0W7ECjwTxFR4Deq+ojb/kUR+SSwAviaqlZHOYcxvdYa\nCLGjsoGtB+vZfqiR8kMNrN7axJ1L/kVFbcsxy2ak+hg1aAAjB2Vy5ug8hudlMnxgJicNzGD4wEyG\n5aaTker36F9izLHE6bI6SisXGa6qe0VkGLAQ+BKwCTiEUyR+CBSq6mc6ee9cYC5AQUFB0bx583qU\nob6+nuzs7B7+C6LHckUmlrkCIWVfg7KnLsSehhB7653hQKPS8aKY3DQYmq4U5qRSkCUMG+BjWKYw\nNNNHTpq3f5XbzzEy8ZoLepetpKSkTFWLu1suqoXgmA05J4nrO54bEJGxwAJVndrVe4uLi3XFihU9\n2m5paSmzZ8/u0XujyXJFJlq5ahrbWL+vho376tiwt5aN+2rZcqCe1qBzE4/fJ4wZMoAJ+dlMLMhm\n4rAcTs7PZszQAeRmpCbd59VblityvckmImEVgqgdGhKRLMCnqnXu9GXAD0SkUFX3uYt9FFgXrQzG\ndNTcFmTDvlrW7DrsDLtr2H6o4cjr+TnpTC7M5UOn5DO5MIdTT8pl7NABpKfYIRyT2KJ5jqAAeNH9\nepwCPKOqr4jIUyIyHefQUDlwWxQzmCR2qL6FFeVVLNtexYryajbuqyXgXpVzUm4GZ4wayLVFI5k2\nYiCTC3PJz4n+HZzGxKOoFQJV3Qac0Un7zdHapkluFbXNLNl6iHe2V/HO9iq2HnT+2s9I9XHmqEHc\nduF4zhiZxxmj8ijIzfA4rTHxwy4fNf1Wc1uQZdureHPzQd58/xCbKuoAyMlIYebYwVxXPIqzxg1m\n6vCBpKVYR7vGnIgVAtOv7KpqZOGGChZtOsCy7VW0BkKk+X3MHDeIj844lfMnDGVyYS5+n11Hb0y4\nrBCYuKaqbNxXx/wtrfzPmjfZsK8WgAnDsrl51hgumDiUs8cNITPNTuga01NWCEzcUVU27Ktl/qo9\n/H3dfnZXNyFA8dhs7pkzmUunFDB2aJbXMY1JGFYITNzYXd3IS6v3Mn/VHt4/UE+qXzh/wlC+WDKB\nzOqtXP3hc72OaExCskJgPNXQEuDlNXt5ceUe3imvAqB4zCDuvWYqV04rZFBWGgClpdu8jGlMQrNC\nYDyxaX8df3h7By+u2kN9S4CT87P4+mWncPX0EYwaPMDreMYkFSsEJmZaAkFeWbefP7y9g+Xl1aSl\n+LhyWiGfmDWaGaMHWY+ZxnjECoGJuprGNp5cWs7vl5RT2dDKmCED+NacU7m2aBSD3UM/xhjvWCEw\nUbO/ppnH3trGM8t20tAapGRSPp8+bxznTxiKz67zNyZuWCEwfW7bwXp+8/o2/rxqNyGFfzu9kNsu\nPJnJhbleRzPGdMIKgekzu6oa+dnCzby4eg9pfh83zBzNrReMZ/QQO/lrTDyzQmB6raqhlYcWbeGp\npTsQgbkXjOdzF4y33jyN6SesEJgea2wN8Pjich4u3UpDa4DrikbxlUsnUjgw0+toxpgIWCEwEVNV\nni/bzU/+sYkDdS1cMrmAb14+iVMKcryOZozpASsEJiLvV9Rxz/x1vLO9ijNH5/HQx2cwc+xgr2MZ\nY3rBCoEJS3NbkF++9j6PvLGNAWkp/Pjfp3F98Si7DNSYBGCFwHTr9c0H+fb8deysauTfZ4zgnjmT\nGZJtJ4KNSRRWCMwJ1TS18Z2X1vHS6r2Mz8/imVvP5tyTh3odyxjTx6JaCESkHKgDgkBAVYtFZDDw\nR2AszsPrr1fV6mjmMJFbtq2SO55bQ0VtM1+5ZCL/Oftk0lPs4S/GJKJYPMi1RFWnq2qxO38X8Kqq\nTgRededNnAiElAf+sYkbH32bVL/wwn+ey1cuOcWKgDEJzItDQ1cDs93pJ4BS4E4Pcpjj7Khs4L+X\nNbOtZgvXF4/ku/92GlnpdvTQmEQX7f/lCvxTRBT4jao+AhSo6j4AVd0nIsOinMGE4YWy3XznpXVo\nKMRDN83gytMLvY5kjIkRUdXorVxkuKrudXf2C4EvAS+ral6HZapVdVAn750LzAUoKCgomjdvXo8y\n1NfXk52d3aP3RlO85AqElGffa+XVnQEmDfLxiQlBRg3xPtfx4uXzOp7liozlilxvspWUlJR1OCx/\nYqoakwH4HvB1YBNQ6LYVApu6e29RUZH21KJFi3r83miKh1zVDS164yNLdcydC/TeBes1EAzFRa7O\nWK7IWK7IxGsu1d5lA1ZoGPvnqJ0sFpEsEclpnwYuA9YBLwO3uIvdArwUrQzmxN6vqOPqhxazorya\nB647g3uunILfbg4zJilF8xxBAfCi+/jBFOAZVX1FRJYDz4nIZ4GdwHVRzGA68dp7Fdz+7GoyUv08\nO3cWRWM+cGTOGJNEolYIVHUbcEYn7ZXAxdHarjkxVeWRN7bx41feY0phLo9+spjhedZTqDHJzq4N\nTBKhkPKDBRv4/ZJyrpxWyE+uO50BafbjN8ZYIUgKwZBy1wtr+VPZbj57/jj+35WTcQ/ZGWOMFYJE\n1xYM8dU/rmbB2n3cfvFEvnrJRCsCxphjWCFIYM1tQb74zEr+tfEAd19xKrddeLLXkYwxccgKQYJq\naAkw96kVLN5SyQ+vPo2bzxnrdSRjTJyyQpCAapvb+PTjy1m107lH4NqikV5HMsbEMSsECaa5Lchn\nf7+cNbsO86ubZjBnmvUZZIzpmhWCBBIIhvjiM6tYsaOaX954phUBY0xYYvE8AhMDqsr/m7+Of22s\n4LtXTeGq04d7HckY009YIUgQ//vPzcxbvosvlkzgU+eN8zqOMaYfsUKQAH6/eDu/WrSFG2aO4muX\nneJ1HGNMP2OFoJ/7y5q9fH/BBi6dUsC910y1m8WMMRGzQtCPLd5yiDueW83MMYP55Y1nkuK3H6cx\nJnK25+indlU18oWnVzJuaBaP3lJMRqo9XN4Y0zNWCPqh5rYgn/9DGarKo58sZmBmqteRjDH9mN1H\n0M+0Xya6fm8tv/tUMWOGZHkdyRjTz9k3gn7mmXd28nzZbm6/eCIXnVrgdRxjTAKwQtCPrNpZzfde\nXs+Fp+Tz5Ysneh3HGJMgTlgIRORHInJbJ+1fFZH7ohvLHK+yvoUvPL2SgtwMHrxhuj1o3hjTZ7r6\nRnA18Ggn7b9wXwuLiPhFZJWILHDnfy8i20VktTtMjyxy8gkEQ3zp2VVUNrTy8CeKyBuQ5nUkY0wC\n6epkcUhVQ8c3qmpQIrtr6cvARiC3Q9s3VPX5CNaR1P534WaWbK3kf649nakjBnodxxiTYLr6RtAs\nIh94pJXb1hzOykVkJHAl8NuexTPLtlXy8OtbuWHmKK4vHuV1HGNMAuqqEHwX+JuIfEJEJrvDzcBf\n3dfC8XPgm8Dx3yzuE5G1IvIzEUmPPHZyaGgJ8I3n1zJq0AC+fdUUr+MYYxKUqOqJXxQ5A2dHPtVt\nWgf8RFVXd7tikauAOar6BRGZDXxdVa8SkUJgP5AGPAJsVdUfdPL+ucBcgIKCgqJ58+ZF9A9rV19f\nT3Z2do/eG03h5HpqQwuv7Qxw11kZTBocmzuH+/Pn5QXLFRnLFbneZCspKSlT1eJuF1TVqAzAj4Dd\nQDnOjr8R+MNxy8wGFnS3rqKiIu2pRYsW9fi90dRdrsXvH9Qxdy7Q77+8PjaBXP318/KK5YqM5Ypc\nb7IBKzSM/XVXl48OE5EfusNJkVYiVb1bVUeq6ljgBuA1Vf2E+40A94TzNTjfMkwHdc1tfOP5tYwb\nmsU3PjzJ6zjGmATX1VVDfwC+Aag7fUkfbfNpEckHBFgNfL6P1psw/vtv77Gvpok/ff5cMtOsMzlj\nTHR1VQhagGE4hSDQm42oailQ6k5f1Jt1Jbo3Nh/k2Xd2ctuHxlM0ZpDXcYwxSaCrq4ZuAk4BJuMc\n2jFRVtPUxp0vrOXk/Cy+eqk9acwYExsn/EagqnXAQzHMkvTuXbCBitpm/vyF8+z5AsaYmLFO5+LE\n0q2V/KlsN5+/8GSmj8rzOo4xJolYIYgDwZDy/b+sZ0ReJrdbr6LGmBjr6vLRmbEMkszmLd/Je/vr\n+NacyXZIyBgTc119I7jN7TX0D243E/kxS5VEahrbeOAfmzhr3GDmTIv4dg1jjOm1rk4Wfw5ARKYC\nVwDPikgG8BrwCvC2dtI7qYnMg6++z+GmNr77b1OIrFNXY4zpG92eI1DVdar6E1W9BPgwsAK42R2b\nXthyoJ4nl5Zzw8zRnDbcupc2xngjoofXq2oD8LI7mF66968byEzz8/XL7J4BY4x37Kohj6w5GKB0\n00G+fPFEhmRbT9zGGO9YIfBAayDEsxtbGZ+fxSfPGet1HGNMkgurEIjISBEpcafTRSQrurES25NL\ny9nfqHz7qimkpVgtNsZ4q9u9kIh8BuecQPvjJscAL0UzVCI7VN/Cg/96n9Pz/ZRMGuZ1HGOMCesb\nwe3ALKAWQFU34/RKanrgN69vpaE1wA2T0ryOYowxQHiFoFlVW9tnRMSP8ywBE6FD9S089fYOrpk+\nguHZdkjIGBMfwtkbLRaRbwIZ7nmCPwILohsrMT365jZaAyG+eNEEr6MYY8wR4RSCbwJ1wHvAl4FX\ngXuiGSoRVTW08tTSHXzkjOGMz4/Ph2QbY5JTlzeUuYeBfqeqtwC/jk2kxPTom9toagvatwFjTNzp\n8huBqgaBQhFJjVGehFTd0MqTS8q56vThTBiW43UcY4w5RjhdTGwD3hSRl4CG9kZV/UU4G3C/VawA\n9qjqVSIyDpgHDAZWAjd3PBmdiB57azuNbUFut28Dxpg4FM45goPAQmAAkN9hCNeXgY0d5u8Hfqaq\nE4Fq4LMRrKvfOdzYyu+XlDNnWiETC+zbgDEm/nT7jUBVv93TlYvISOBK4D7gDnH6Wb4IuMld5Ang\neyTw+YffvbWd+pYAX7JvA8aYONVtIRCRhYAe366ql4Wx/p/jXHXU/qfwEOCwqgbc+d3AiPCi9j81\njW08vricK6aexKkn5XodxxhjOiWqH9jHH7uAyNkdZjOAjwEtqvqNbt53FTBHVb8gIrOBrwOfBpaq\n6gR3mVHA31R1WifvnwvMBSgoKCiaN29e2P+ojurr68nO9uZyzflbWpm/pY0fnJvB6NxjH0HpZa6u\nWK7IWK7IWK7I9SZbSUlJmaoWd7ugqkY8AK+HscyPcP7iLwf2A43A08AhIMVd5hzgH92tq6ioSHtq\n0aJFPX5vb9Q0teq0776itz6xvNPXvcrVHcsVGcsVGcsVud5kA1ZoGPv0cDqdy+0w5InIxUBhGAXm\nblUdqapjgRuA11T148Ai4Fp3sVtI0A7snlq6g9rmALdfPNHrKMYY06VwLh9dj3OOQIAAsB24tRfb\nvBOYJyL3AquAx3qxrrjUFgzx1NIdXDBxKFNH2CMojTHxLZxCMF5V2zo2iEikj7gsBUrd6W3AWZG8\nv79ZuKGC/bXN3HvNVK+jGGNMt8K5j2BZJ23v9HWQRPLEknJGDsqk5FTrrdsYE/9O+Je9iAzDOReQ\nKSLTONr1dC7OzWWmExv31bJsexXfmnMqfp/11m2MiX9dHeK5EvgMMBL4vw7tdUCPbzJLdE8u3UF6\nio/ri0d5HcUYY8JywkKgqo8Dj4vI9ar6XAwz9Vs1jW3MX7WHa6aPIG+APYHMGNM/hNPFxHMi8mHg\nNJwbytrb/zuawfqjP5XtoqktyCfPHeN1FGOMCVs4XUz8H5AHfAh4HOfO4rejnKvfCYWUp97ewcyx\ngzhtuF0yaozpP8K5auh8Vb0JqFSnA7qzcc4bmA5e33yQHZWNfPKcsV5HMcaYiIT18Pr2sYic5M6P\njVqifuqJpeUMy0nn8qkneR3FGGMiEk4h+JuI5AEPAKtx+g56Ppqh+pvthxoo3XSQj589hlR/OB+p\nMcbEj+6eWewD/q6qh4E/icgCIFNVq2KSrp94aukOUv3CjWfbJaPGmP6nu2cWh4AHO8w3WRE4VkNL\ngD+V7eKKqYUMy8no/g3GGBNnwjmOsVBEro56kn5q/uo91DUHuMUuGTXG9FPhdB73RWCgiLQATThd\nTaiqDo5qsn7imWU7OW14LjNGD/I6ijHG9Eg43wiGAqlANs5D64cS2cPrE9bmijrW763l2qKROI9j\nNsaY/qfbQqCqQeA64E53uhCYHu1g/cGLq/bg9wn/dsZwr6MYY0yPhfOEsl8BJcDNblMj8HA0Q/UH\noZDy0qo9fGjiUIZmp3sdxxhjeiycQ0PnquptuDeWuVcNJX2Pasu2V7G3pplrzhzhdRRjjOmVcApB\nm3s/gQKIyBAgFNVU/cD8VXvISvNz2RS7k9gY07+FUwgeAl4A8kXk+8BbwP1RTRXnmtuC/O3dfVw+\ntZDMNL/XcYwxplfC6Yb6SREpAy5xm65T1XXRjRXfXt14gLqWAB+1w0LGmAQQbsc4fqANaA33PSKS\nISLviMgaEVnvfptARH4vIttFZLU79LsrkF5ctYeC3HTOOXmI11GMMabXwrlq6B7gWWA4TvfTz4jI\n3WGsuwW4SFXPwLnc9HIRmeW+9g1Vne4Oq3uY3RNVDa2UbjrA1dNH2DOJjTEJIZw7iz8BFKlqI4CI\n3AeUAT/q6k2qqkC9O5vqDtrzqPHhr2v3Eggp10y3w0LGmMQgzv66iwVEXgGuV9Vadz4XeFZVr+x2\n5SJ+nKIxAXhIVe8Ukd8D5+B8Y3gVuEtVWzp571xgLkBBQUHRvHnzIvl3HVFfX092dnaP3tuZe99u\nojmg/PC8zF7dTdzXufqK5YqM5YqM5Ypcb7KVlJSUqWpxtwuqapcD8GdgF/Bb4FFgB/A08FPgp929\n311HHrAImIpzZ7IA6cATwHe6e39RUZH21KJFi3r83uNtP1ivY+5coL8u3dLrdfVlrr5kuSJjuSJj\nuSLXm2zACg1jHx3OoaG/ukO7iJ9XrKqHRaQUuFxVH3CbW0TkceDrka7PK/NX70EEPmJdShhjEkg4\nl48+1pMVi0g+0OYWgUycy0/vF5FCVd0nznGVa4B+cSmqqjJ/1R5mjRvC8LxMr+MYY0yf6bYQiMjl\nwA+BMe7y4XZDXQg84Z4n8AHPqeoCEXnNLRKC8+jLz/fmHxArq3YdpryykS+UTPA6ijHG9KlwDg39\nCrgeeJcIupZQ1bXAmZ20XxR2ujgyf9Ue0lN89nB6Y0zCCacQ7AZWq/PYyqQUCIZYsHYfl0wuIDcj\n1es4xhjTp8IpBN8E/uKe7D1ymaeq/iJaoeJN2Y5qqhpamTOt0OsoxhjT58IpBN/H6V4ijyTtdXTh\nhgrS/D4unGQPZjPGJJ5wCsEwVS2KepI4paos3FjBrJOHkJ0ezsdljDH9SzgdyL0qIv3yBG9f2HKg\nnh2VjVw6pcDrKMYYExXhFIJbgX+JSL2IVIlItYhURTtYvFi4sQKASyYP8ziJMcZERzjHOoZGPUUc\nW7ihgmkjBlI40G4iM8Ykpm6/EahqELgOuNOdLsTpVjrhHahrZvWuw3ZYyBiT0MJ5HsGvgBLgZrep\nEXg4mqHixWsbD6AKl0y2QmCMSVzhHBo6V1VniMgqAFWtEpG0KOeKC//aWMGIvEwmF+Z4HcUYY6Im\nnJPFbSLiw32ojIgMIQnuJ2hsDfDm+4e4dEpBr547YIwx8e6EhUBE2r8tPAS8AOS7zx1+C7g/Btk8\n9db7h2gJhOz8gDEm4XV1aOgdYIaqPikiZTjdSAtwnar2i66je2PhhgpyMlI4a1x3nawaY0z/1lUh\nOHI8RFXXA+ujHyc+BEPKa+8doGTSMFL94Rw9M8aY/qurQpAvInec6EVV/WkU8sSFVTurqWxotcNC\nxpik0FUh8APZdPhmkCwWbqwg1S/WyZwxsRYKQbDVHdo+OB1qc6c7zgeOtocCR8ehNggFj2tz59sH\nDUEowMRdO6DuRff1IGjw2OU11KHdHase1xZyp0NHpzXU4b0h5z3Hv3bMoMeNQww67W5gdlQ/9q4K\nwT5V/UFUtx6nFm6oYNb4IfbsAZO8VJ0dbWuDM7Q1ukNTh3GHIdAEbc3OONDitrVAoPmY8Yzqg7Ax\n3Z1v38m3uNMtzk43Fnyp4EtxBz/5wRDUZHRo84H4j7zuTLe3ufPic6dTO2nzfXBa/CBydNzxNdx5\nxG2TI681t0X/D9KwzhEkk60H69l2sIFPnTvW6yjGREbV2cE21zCgYTfsWg4tNdBS18lQCy31R3f0\nrXUdpt1Bg5Fn8KdBSiakZkBKujOdkgYpGZCSQVtqDuQVOm3+9A7jdOe9/jS37fgh1R3SnJ24P+XY\naZ/7+pH5lKM7++PnfR8877ektJTZs2f3/mcQBU2lpVHfRleF4OKobz0O/WuD08ncxXY3sfGKKrTW\nQ2MVNFZCU5U7XeVMN1VD02Fn3Hz46HxzjXM4BDgLYPkJ1p+aBenZkJ4DaVmQlgPZJznT6dnO62kD\nIHWA09ZxnJrZyTjzyI6+s51sR+/G8Q43mZ2wEKhqr3oYFZEM4A0g3d3O86r6XREZB8wDBgMrgZtV\ntbU32+pLCzdUcNrwXEbkWSdzpg+pOjvs+gpnqHPHDQegoRIaDrrDIWccbDnBigQyciFzkDNk5MHA\nUZCZBxkDIT0XMgayYdsepsw4x53PhTR3x5+e4x6CMOaoaD5ppQW4SFXrRSQVeEtE/g7cAfxMVeeJ\nyMPAZ4FfRzFH2CrrWyjbWc3tF030OorpT4IBqNsLNbuhdq8z1O2D2j1Qu8+Zrtt/5K/1Y6RkQFY+\nZA2F7GEwbIozPWCIM84c7EwPcMcZA8PakR9oKGXKxNl9/281CSlqhUBVFah3Z1PdQYGLgJvc9ieA\n7xEnhWDJ1kpUoeRUe/aA6SAYgNrdUF3uDjugZhcc3sWsA1vg9SrnCo+O0rIhpxByC2HMeZBzEmQX\nQE6BM84+ydnxp+c4JwaN8ZA4++sorVzED5QBE3C6qvgJ8LaqTnBfHwX8XVWndvLeucBcgIKCgqJ5\n8+b1KEN9fT3Z2dlhLfv4uhbe2R/goYsH4Ivyf85IcsVSsuaSUICM5goym/YxoHEPmU17yWzaR2ZT\nBRnNB5AO3WuFxE9L+lBa0vOpSxlEMLuQlvR8mjPyaUkfQkv6UIIpA6KWNRzJ+nPsqXjNBb3LVlJS\nUqaqxd0tF9WH8LrPL5guInnAi8DkzhY7wXsfAR4BKC4u1p6eYCqN4OTUd95ZxPmnDOGikm4/t16L\nJFcsJXyulno4tBkOboKDG53xofedv/Q7XiWTOQgGj4cR58OgsccMvtzhZPr8ZAKrE/3z6mOWK3Kx\nyBaTp7Gr6mERKQVmAXkikqKqAWAksDcWGbqzq6qRnVWNfPq8sV5HMX0hGICqrbD/XahYBxXr4cB7\nULPz6DL+NBgyAU6aBqd91JkeMgGGnOwckzcmSUStEIhIPtDmFoFMnE7r7gcWAdfiXDl0C/BStDJE\nYunWSgDOm5DUT+bsnwItzs5+7yrYu9qZPrDRuaYenOvH8yfB6LMh/5OQf6ozDBrnXGNuTJKL5v+C\nQuAJ9zyBD3hOVReIyAZgnohjX7KCAAASNUlEQVTcC6wCHotihrAt3nqIodnpTBwWn8cJjSsUdHby\ne1a4O/5VULHh6BU5mYOdv/Bnfg4KpsJJU2HoJOcmJWNMp6J51dBa4MxO2rfh3u8SL1SVJVsrOffk\nIfYQmnjTVA27VzB2+3Ow82ewu8y5Cxaca+iHnwnnfskZDz8TBo60q3CMiZB9Lwa2HKjnYF0L500Y\n4nUUU38QdiyG8rec4eBGAMbgc/66P+M/YNTZMLLYObRjO31jes0KAbB4yyEAzj3Zzg/EXGMVbH8d\ntr/pFICD7zntqQNg9CyY+jEYdRZvbWvkgkuu8DarMQnKCgGweGslowZnMmqwt9d+J4VgG+xeDltf\nc4Y9KwF1bsAaPQtO/w8YewEMn+50Itb+tp2lnkU2JtElfSEIhpS3t1Vy5bRCr6Mkrrr9sPkV2PxP\n2P6Gc4xffDByJsy+C8aXwIgZx+z4jTGxk/SFYN2eGuqaA5xrl432HVXYvxY2vQKb/+5c2QMwcDRM\nuxYmXOz81Z+Z521OYwxghYDFW53zA+eMtxPFvRIKwq5lsH4+vLfA6XANcU7qXvRtmHSF06Gandw1\nJu4kfSFYsqWSSQU55Oekex2l/wkGYOcS2PASbPyL061ySgZMuARKvgUTPwzZ9rhPY+JdUheClkCQ\n5eVV3HT2aK+j9B+qsOsdePc556//xkPOU6hOuQymXO3s/NPtpjxj+pOkLgQrdxymJRCyy0bDceh9\nWPucUwCqy52//CddAVOugYmXOk+wMsb0S0ldCJZsPYRP4Ozx1sFYpxqr4N0/wZp5sHelc6XPuAvh\nwjvh1KucJ18ZY/q9JC8ElZw+Mo/cDLts8QgNwbbXYeWTznH/YIvTd89l9zk3d+XaZbbGJJqkLQT1\nLQHW7DrM3A+N9zpKfKjdC6uf5uxlj8LrFc4jEYtugTNvhsLTvU5njImipC0E72yvJBDS5O52WhXK\n34R3HoX3/goapDlvGplz7oPJV0FqptcJjTExkLSFYPGWStJSfBSNGeR1lNhrqXOO+y//rdO3T+Yg\nOOe/oPjTrFm7k9mnz/Y6oTEmhpK2ECzZWknR6EFkpPq9jhI7lVth2W9g9TNONw+F0+Hq/4Op/97h\nr/+dXa7CGJN4krIQVNa3sHFfLV+/7BSvo0SfKux8G5b+yjn84091Hst41lwYUWR3+hpjkrMQvL2t\nCoBzEvn+gWAA3vsLLPkl7ClzDv9c8DWnAOQUeJ3OGBNHkrIQlO2oJj3Fx7QRA72O0vfammDVH2DJ\nL+DwThg8HuY8ANNvspu+jDGdSspCsGpXNaePHEhais/rKH2nuRZWPAZLH4KGgzDyLPjwf8OkOeBL\novMgxpiIRa0QiMgo4EngJCAEPKKqD4rI94BbgYPuot9S1b9FK8fxWgJB1u+p5VPnjY3VJqOroRKW\nPQzv/Aaaa+Dki5xDQGPOs+P/xpiwRPMbQQD4mqquFJEcoExEFrqv/UxVH4jitk9o/d5aWoMhZozu\n533h1x+EJQ/C8segrdHp8uGCrzkPeDHGmAhErRCo6j5gnztdJyIbgRHR2l64Vu6oBuDM0f30/oGO\nBSDQDFOvhQvugGGTvU5mjOmnRFWjvxGRscAbwFTgDuBTQC2wAudbQ3Un75kLzAUoKCgomjdvXo+2\nXV9fT3b20W6R/291M1sPh/jf2d4+n/j4XN1Jba1h1K75jNjzV3yhNioKPsSOMdfRNGCkp7lixXJF\nxnJFJl5zQe+ylZSUlKlqcbcLqmpUByAbKAP+3Z0vAPyAD7gP+F136ygqKtKeWrRo0THz5/7oVf3C\n02U9Xl9fOT7XCTVUqi78ruq9harfHaj6/GdVD2zyPleMWa7IWK7IxGsu1d5lA1ZoGPvpqF41JCKp\nwAvA06r6Z7fwVHR4/VFgQTQzdFRR28yew0185vxxsdpkz7XUwdsPO5eBttQ5d/9eeCfkT/I6mTEm\nwUTzqiEBHgM2qupPO7QXqnP+AOCjwLpoZTjeqp3t5wfi+ERxW7NzGeibP3We/jXpSrjoHig4zetk\nxpgEFc1vBOcBNwPvishqt+1bwI0iMh1QoBy4LYoZjrFy52HS/D5OGx6HD1QJtsHqp+H1/3Ee/D5+\ntvPQ95HdH94zxpjeiOZVQ28BnV3IHrN7Bo63amc1p43IJT0ljm6wCoVgw3x47V6o2gojZ8I1v4bx\nF3qdzBiTJJLmzuLWQIi1u2v4xKwxXkdxqMKWV+HV78O+NZA/GW541nkOsN0IZoyJoaQpBO/tr6Ul\nEIqP8wO7yzhjzbfh9Xdh4Gi45mE4/XrrCsIY44mkKQTtN5LN8PJGsoOb4LUfwsa/kJU6EC6/H4o/\nDSnp3mUyxiS9pCkEq3YdpiA3ncKBGbHf+OFdUPpjWPMMpGbB7LtZFjidC2bNiX0WY4w5TtIUgpU7\nq5kxehASy+PvDYecy0CXP+rMn/2fTncQWUMJlpbGLocxxnQhKQrBwboWdlU1cXOsThQ3HYa3f+10\nCd3W4DwL4MK7IG9UbLZvjDERSIpCsHrXYSAG5wda6pwuoZf80ukSevJH4KL/Z3cDG2PiWlIUgpU7\nq0nxCVOj9USy1kbn8M/iB6GxEk65AkruhsIzorM9Y4zpQ0lRCFbtrOa04blkpPbx5Zkt9bDyCXjr\n59BwAE6+GErugZFFfbsdY4yJooQvBMGQsmZXDf8xsw+PzzdUOk8EW/YbaD4MYy+A65+EMef03TaM\nMSZGEr4Q7K4P0dQW7JsbyQ7vdE4Alz0BgSbnqWDnfQVGzez9uo0xxiMJXwi2Hg4BvThRrAo7l8KK\nx2H9n5220/8Dzr0dhp3aRymNMcY7SVEIhmanMXJQZmRvbKiENc865wAObYa0HJh5K5zzX3YZqDEm\noSR8IdhyOMiZY4eEdyNZsA3K34SVT8F7CyDY6vQG+pFfOQ+GScuKfmBjjImxhC4E1Q2tVDRq1+cH\nqsudXkC3vArb34DWOsjIg+LPwIxboGBKzPIaY4wXEroQHLmRbNRAqKuAur1Q6w6HNsPW16Byi7Pw\nwNEw7VqYcDFMuBRSPeiTyBhjPJDQhSBt8QMsTn+B4U8fhlDg2BdTMmHs+TDzczDhEhgywZ4DYIxJ\nSgldCFLzhrP9wBRGzCyC3BGQOxxyCp3prHzw+byOaIwxnkvoQnDWx75Cael0mD3b6yjGGBO3ovYn\nsYiMEpFFIrJRRNaLyJfd9sEislBE3nfHHj4pxhhjTDSPjQSAr6nqZGAW8F8iMgW4C3hVVScCr7rz\nxhhjPBK1QqCq+1R1pTtdB2wERgBXA0+4iz0BXBOtDMYYY7oXk7OlIjIWOBNYBhSo6j5wigUwLBYZ\njDHGdE5UNbobEMkGXgfuU9U/i8hhVc3r8Hq1qn7gPIGIzAXmAhQUFBTNmzevR9uvr68nOzu7Z+Gj\nyHJFxnJFxnJFJl5zQe+ylZSUlKlqcbcLqmrUBiAV+AdwR4e2TUChO10IbOpuPUVFRdpTixYt6vF7\no8lyRcZyRcZyRSZec6n2LhuwQsPYV0fzqiEBHgM2qupPO7z0MnCLO30L8FK0MhhjjOleNO8jOA+4\nGXhXRFa7bd8Cfgw8JyKfBXYC10UxgzHGmG5E/RxBXxCRg8COHr59KHCoD+P0FcsVGcsVGcsVmXjN\nBb3LNkZV87tbqF8Ugt4QkRUazsmSGLNckbFckbFckYnXXBCbbNbZjjHGJDkrBMYYk+SSoRA84nWA\nE7BckbFckbFckYnXXBCDbAl/jsAYY0zXkuEbgTHGmC4kdCEQkctFZJOIbBGRmPZyKiK/E5EDIrKu\nQ1unXXCL4xduzrUiMiOKuSLqHjxW2UQkQ0TeEZE1bq7vu+3jRGSZm+uPIpLmtqe781vc18dGI5e7\nLb+IrBKRBfGSyd1euYi8KyKrRWSF2xYPv2N5IvK8iLzn/p6d43UuEZnkfk7tQ62IfMXrXO62vur+\nzq8TkWfd/wux/R0L5/bj/jgAfmArMB5IA9YAU2K4/Q8BM4B1Hdr+B7jLnb4LuN+dngP8HRCcLruX\nRTFXITDDnc4BNgNTvM7mrj/bnU7F6aBwFvAccIPb/jDwn+70F4CH3ekbgD9G8TO7A3gGWODOe57J\n3UY5MPS4tnj4HXsC+Jw7nQbkxUOuDvn8wH5gjNe5cHpk3g5kdvjd+lSsf8ei+oF7OQDnAP/oMH83\ncHeMM4zl2ELQaT9LwG+AGztbLgYZXwIujadswABgJXA2zo00Kcf/THH6sDrHnU5xl5MoZBmJ89yM\ni4AF7o7B00wdspXzwULg6c8RyHV3bBJPuY7LchmwOB5y4RSCXcBg93dmAfDhWP+OJfKhofYPuN1u\nt81LJ+qC25OsEl734DHL5h6CWQ0cABbifKM7rKqBTrZ9JJf7eg0wJAqxfg58Ewi580PiIFM7Bf4p\nImXi9NYL3v8cxwMHgcfdw2m/FZGsOMjV0Q3As+60p7lUdQ/wAE53O/twfmfKiPHvWCIXAumkLV4v\nkYp5VnG6B38B+Iqq1na1aCdtUcmmqkFVnY7zV/hZwOQuth31XCJyFXBAVcs6NnuZ6TjnqeoM4Aqc\nJwB+qItlY5UtBeeQ6K9V9Uygga6fQhjTz8w91v4R4E/dLdpJW5/ncs9JXA2MA4YDWTg/zxNtOyq5\nErkQ7AZGdZgfCez1KEu7ChEpBHDHB9z2mGYVkVScIvC0qv45nrIBqOphoBTn2GyeiLR3jthx20dy\nua8PBKr6OMp5wEdEpByYh3N46OceZzpCVfe64wPAizjF0+uf425gt6ouc+efxykMXudqdwWwUlUr\n3Hmvc10CbFfVg6raBvwZOJcY/44lciFYDkx0z76n4XwdfNnjTCfqgvtl4JPulQqzgJr2r6t9TSTi\n7sFjkk1E8kUkz53OxPkPshFYBFx7glztea8FXlP3wGlfUdW7VXWkqo7F+f15TVU/7mWmdiKSJSI5\n7dM4x73X4fHPUVX3A7tEZJLbdDGwwetcHdzI0cNC7dv3MtdOYJaIDHD/b7Z/XrH9HYvmSRmvB5wz\n/5txjjXfE+NtP4tzzK8Np4p/FudY3qvA++54sLusAA+5Od8FiqOY63ycr5JrgdXuMMfrbMDpwCo3\n1zrgO277eOAdYAvO1/l0tz3Dnd/ivj4+yj/P2Ry9asjzTG6GNe6wvv332+ufo7ut6cAK92c5HxgU\nJ7kGAJXAwA5t8ZDr+8B77u/9U0B6rH/H7M5iY4xJcol8aMgYY0wYrBAYY0ySs0JgjDFJzgqBMcYk\nOSsExhiT5FK6X8SY+CAi7Zf6AZwEBHG6MwBoVNVzo7TdMpz+XVqjsf6eEpF6Vc32Oofp/6wQmH5D\nVStxrlFHRL4H1KvqA9Hcptsf0554KwLG9CU7NGQSgojUu+PZIvK6iDwnIptF5Mci8nFxnnXwroic\n7C6XLyIviMhydzjvBKu+Anilk+39WEQ2uH3VP9DVOkUkW0Qed7e/VkQ+5rbf6LatE5H7O/5bROQ+\ncZ7N8LaIFLjt40RkqbvuH3ZYvlBE3hCnn/11InJB33yqJllYITCJ6Azgy8A04GbgFFU9C/gt8CV3\nmQeBn6nqTOBj7muduZzjCoGIDAY+CpymqqcD93azzm/jdFEwzV3+NREZDtyP03/RdGCmiFzjLp8F\nvK2qZwBvALd2WP+v3fXv7xDpJpxuiqe7//bVYXxGxhxhh4ZMIlqubr8wIrIV+Kfb/i5Q4k5fAkxx\nuncBIFdEclS1rr3B7aNqpKpuO279tUAz8FsR+StOH/InXKfbfkN7o6pWuz2FlqrqQXdbT+M8zGg+\n0NphnWU4z4sApxO8j7nTT+EUEnD61fqd25ngfFW1QmAiYoXAJKKWDtOhDvMhjv7O+3BOADd1sZ4L\ngLeOb1TVgIichdNB2A3AF3H+su90nW5nYsf35dJZd8Lt2vRo3y9Bjv1/+oE+YVT1DbewXAk8JSI/\nUdUnu1i/McewQ0MmWf0TZwcOgIhM72SZy3EeV3gMcZ7lMFBV/wZ8BfcEdhfrPL59EM7DgC4UkaEi\n4sfpFfP1bjIv5ug3i493WN8YnOcmPIrTs2zUnq9rEpMVApOsbgeK3ZO3G4DPd7LMbDrfOecAC0Rk\nrfv6V7tZ573AIPdE7hqgxD10dTdOd8NrcPrIf4mufRnnATTLcfqh75hztYiswjl09GA36zHmGNb7\nqDGdEJGRwKOq2tnTooxJKFYIjDEmydmhIWOMSXJWCIwxJslZITDGmCRnhcAYY5KcFQJjjElyVgiM\nMSbJWSEwxpgk9/8B9PYcT6vG9+gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b98a780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np                  # basic package for numerical calculations\n",
    "import matplotlib.pyplot as plt     # plotting package\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# parameter values and units\n",
    "P1 = 4                 # watts\n",
    "P2 = 2                 # watts\n",
    "Ua = 0.044             # watts/deg C\n",
    "Ub = 0.018             # watts/deg C\n",
    "Cp = 6                 # joules/deg C\n",
    "u1 = 0.5               # steady state value of u1 (fraction of total power)\n",
    "u2 = 0.0               # steady state value of u2 (fraction of total power)\n",
    "T_ambient = 21         # ambient temperature of the surroundings\n",
    "\n",
    "def deriv(T,t):\n",
    "    T1, T2 = T\n",
    "    dT1 = -(Ua + Ub)*T1/Cp + Ub*T2/Cp + P1*u1/Cp\n",
    "    dT2 = Ub*T1/Cp - (Ua + Ub)*T2/Cp + P2*u2/Cp\n",
    "    return [dT1, dT2]\n",
    "\n",
    "# compute the transient solution\n",
    "t = np.linspace(0,800)\n",
    "T = odeint(deriv, [0,0], t)\n",
    "\n",
    "# plot\n",
    "plt.plot(t, T + T_ambient)\n",
    "plt.xlabel('Time / seconds')\n",
    "plt.ylabel('Temperature / °C')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.6 Fitting the Model to Measured Step Response Data](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.6-Fitting-the-Model-to-Measured-Step-Response-Data)",
     "section": "2.4.6 Fitting the Model to Measured Step Response Data"
    }
   },
   "source": [
    "## 2.4.6 Fitting the Model to Measured Step Response Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.6 Fitting the Model to Measured Step Response Data](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.6-Fitting-the-Model-to-Measured-Step-Response-Data)",
     "section": "2.4.6 Fitting the Model to Measured Step Response Data"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dcdd10cdbb834532ba0216d18536017f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.compare>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "from ipywidgets import interact\n",
    "\n",
    "data = pd.read_csv('Step_Test_Data.csv').set_index('Time')[1:]\n",
    "t = data.index\n",
    "T1 = data['T1'].values\n",
    "T2 = data['T2'].values\n",
    "\n",
    "# known parameter values\n",
    "P1 = 4\n",
    "u1 = 0.5   # steady state value of u1 (fraction of total power)\n",
    "P2 = 2\n",
    "u2 = 0.0\n",
    "T_ambient = 21\n",
    "\n",
    "def compare(Ua, Ub, Cp):\n",
    "\n",
    "    def deriv(T,t):\n",
    "        T1, T2 = T\n",
    "        dT1 = -(Ua + Ub)*T1/Cp + Ub*T2/Cp + P1*u1/Cp\n",
    "        dT2 = Ub*T1/Cp - (Ua + Ub)*T2/Cp + P2*u2/Cp\n",
    "        return [dT1, dT2]\n",
    "\n",
    "    T = odeint(deriv, [0,0], t)\n",
    "\n",
    "    # plot\n",
    "    plt.plot(t, T + T_ambient, t, T1, t, T2)\n",
    "    plt.xlabel('Time / seconds')\n",
    "    plt.ylabel('Temperature / °C')\n",
    "    plt.grid()\n",
    "    plt.text(400,28,'Ua = ' + str(Ua))\n",
    "    plt.text(400,24,'Ub = ' + str(Ub))\n",
    "    plt.text(400,20,'Cp = ' + str(Cp))\n",
    "\n",
    "# parameter values and units\n",
    "P1 = 4                 # watts\n",
    "P2 = 2                 # watts\n",
    "Ua = 0.044             # watts/deg C\n",
    "Ub = 0.018             # watts/deg C\n",
    "Cp = 7                 # joules/deg C\n",
    "\n",
    "interact(compare, Ua=(0,0.06,0.001), Ub=(0,0.06,0.001), Cp=(2,10,0.002))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.6 Fitting the Model to Measured Step Response Data](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.6-Fitting-the-Model-to-Measured-Step-Response-Data)",
     "section": "2.4.6 Fitting the Model to Measured Step Response Data"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.4.6 Fitting the Model to Measured Step Response Data](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html#2.4.6-Fitting-the-Model-to-Measured-Step-Response-Data)",
     "section": "2.4.6 Fitting the Model to Measured Step Response Data"
    }
   },
   "source": [
    "<!--NAVIGATION-->\n",
    "< [2.3 First Order Model for a Single Heater](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html) | [Contents](toc.html) | [2.5 Two State Model for a Single Heater](https://jckantor.github.io/CBE32338/02.05-Two-State-Model-for-a-Single-Heater.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE32338/blob/master/docs/02.04-Two-Input-Two-Output-Model.ipynb\"> <img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.ipynb\"> <img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}