{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 0,
     "link": "[](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.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.03-First-Order-Model-for-a-Single-Heater.html)",
     "section": ""
    }
   },
   "source": [
    "<!--NAVIGATION-->\n",
    "< [2.2 Fitting Step Test Data to Empirical Models](https://jckantor.github.io/CBE32338/02.02-Fitting-Step-Test-Data-to-Empirical-Models.html) | [Contents](toc.html) | [2.4 Two-Input, Two-Output Model](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE32338/blob/master/docs/02.03-First-Order-Model-for-a-Single-Heater.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.03-First-Order-Model-for-a-Single-Heater.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.3 First Order Model for a Single Heater](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3-First-Order-Model-for-a-Single-Heater)",
     "section": "2.3 First Order Model for a Single Heater"
    }
   },
   "source": [
    "# 2.3 First Order Model for a Single Heater"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[2.3 First Order Model for a Single Heater](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3-First-Order-Model-for-a-Single-Heater)",
     "section": "2.3 First Order Model for a Single Heater"
    }
   },
   "source": [
    "![](images/arduino_comsol.png)\n",
    "\n",
    "Mathematical modeling is an integral part of process control. The models come in many forms, ranging from the barest of information about a process to sophisticated simulation involving millions of computational nodes.  \n",
    "\n",
    "Our purpose here is to demonstrate certain basic approaches that will provide with insight about the qualitative nature of process dynamics, how to construct simple models from first-principles understanding of the processes, and how to develop models in those circumstances when its not possible to start from first principles.\n",
    "\n",
    "Our first example will be the Temperature Control Laboratory.\n",
    "\n",
    "\n",
    "* First-order lumped model for a single heater/sensor device.\n",
    "* Interacting first-order models for the dual heater/sensor.\n",
    "* Second-order model for a single heater/sensor device.\n",
    "* Interacting second-order model for the dual heater/sensor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.1 First-order lumped model for heater/sensor device.](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.1-First-order-lumped-model-for-heater/sensor-device.)",
     "section": "2.3.1 First-order lumped model for heater/sensor device."
    }
   },
   "source": [
    "## 2.3.1 First-order lumped model for heater/sensor device.\n",
    "\n",
    "We'll be begin by writing a model for one of the heater/sensor pairs. We will assume the whole heater/sensor pair is at single uniform temperature $T_1$. Then we write a dynamic energy balance\n",
    "\n",
    "\\begin{align}\n",
    "C_p\\frac{dT_1}{dt} & = U_a(T_{amb} - T_1) + P_1u_1 \\\\\n",
    "\\end{align}\n",
    "\n",
    "where $T_1$ is the combined temperature of heater/sensor one, $T_{amb}$ is the ambient temperature of the surroundings, and $u_1$ is a fraction of the maximum heater power $P_1$ being applied. The key parameters are the total heat capacity $C_p$ and the overall heat transfer coefficient with surroundings $U_a$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.3.1.1 Steady State](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.1.1-Steady-State)",
     "section": "2.3.1.1 Steady State"
    }
   },
   "source": [
    "### 2.3.1.1 Steady State\n",
    "\n",
    "We'll begin our analysis by investigating the steady-state response of this system to a steady-state input $\\bar{u}_{1}$.  At steady-state all variables are constant so $\\frac{dT_1}{dt} = 0$, which leaves \n",
    "\n",
    "\\begin{align}\n",
    "0 = U_a(T_{amb} - \\bar{T}_1) + P_1\\bar{u}_{1}\n",
    "\\end{align}\n",
    "\n",
    "Solving for $\\bar{T}_{1}$\n",
    "\n",
    "$$\\bar{T}_{1} = T_{amb} + \\frac{P_1}{U_a}\\bar{u}_{1}$$\n",
    "\n",
    "Next we'll load some experimental data that you can use to estimate the value of $U_a$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.2 Loading previously saved experimental data](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2-Loading-previously-saved-experimental-data)",
     "section": "2.3.2 Loading previously saved experimental data"
    }
   },
   "source": [
    "## 2.3.2 Loading previously saved experimental data\n",
    "\n",
    "Previously a step test was performed in which the temperature control laboratory was initially at steady state at an ambient temperature of 21°C. The power to heater 1 was set to 50% for $P_1$ = 4.0 watts. Temperature $T_1$ and $T_2$ were recorded every second for a period of 800 seconds.\n",
    "\n",
    "The data is in a file named \"Step_Test_Data.csv\". The following cell will attempt to read the data from the course github repository and store in a local file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.2 Loading previously saved experimental data](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2-Loading-previously-saved-experimental-data)",
     "section": "2.3.2 Loading previously saved experimental data"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "url = \"https://raw.githubusercontent.com/jckantor/CBE30338/master/notebooks/TCLab/Step_Test_Data.csv\"\n",
    "data = pd.read_csv(url)\n",
    "data.to_csv(\"Step_Test_Data.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.2 Loading previously saved experimental data](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2-Loading-previously-saved-experimental-data)",
     "section": "2.3.2 Loading previously saved experimental data"
    }
   },
   "source": [
    "If that fails to work, then manually download the dat from this link https://github.com/jckantor/CBE30338/blob/master/notebooks/TCLab/Step_Test_Data.csv and save to the same directory as this notebook before proceeding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.2 Loading previously saved experimental data](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2-Loading-previously-saved-experimental-data)",
     "section": "2.3.2 Loading previously saved experimental data"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VfX5wPHPkw0JEAgQ9hJkKSsIKoqAe1RxD9y2aFtbrbW1amu12lZ/tVU7bN1ateKoExdWCCACKspS9h5hJBAyIPv5/fE9SW6Sm+TekJt7kzzv1+u8cs/3rOeee3Oee77nnO9XVBVjjDGtV1S4AzDGGBNelgiMMaaVs0RgjDGtnCUCY4xp5SwRGGNMK2eJwBhjWjlLBCEmIioiA8Mdh6mdiPQUkc9FJFdEfh/ueIxpaq0iEYhIns9QJiKHfManefMcKSKvi0imiBwQkeUicpuIRItIP++AHhPiOJNF5FkR2eUdlNaKyB0+00OSVETkLp/9USAipT7j3x7Ges8QkfX1zDNDRAq995vr7ff7RSQpiO3sEpETGhon8CNgs6q2U9W7D2M95fHcJCIl3v7LEZGvReQMb1qiiPxXRLZ4n+exh7u9emIZLiJFIvJ0HfNEicgjIrLf+/4/EMqYvG1W+cxEZIiIlIR6u962ponIQu848FEIt3ONiGz1vgdviEiHWuZL8L4L+T7/d38PVVz+tIpEoKpJ5QOwFfieT9nLInIEsBjYBhytqh2Ai4GxQLsmDPURIAkYCnQAzgU2hHqjqvoHn/1zE7DQZ/8MD/X2gftVtR3QBfgBMBmYLyIJTbBtgL7Adw1ZsI4fB+ne/uwIvAK87iU3BeYCVwD7G7LNID0OfFHPPD8BTgWGAWOAS0Xk2hDHFU5ZwJ+Bv4RqAyIyGvgrcCnQHRDgsXoWG+zzf3dzqGLzS1Vb1QBsBk6pVvYS8H4dy/TD/QPH+Jk2DlgIZAMZwN+BOJ/pCvwU2AhkAn8ComrZzkpgai3T5nnrygfygEu98nOApd72PwdGVHuvd+IOcvuB54CEevbPtcBnfsqPAmZ761nlGydwHrAayMUl058CKcAhoMyLNw9I8bPeGcCvq5V1BPYC3/fGhwDpwD6v/AWgnTftdW8bB71t/BSIAf4L7Pb2yxzcP5m/9/sKUAwUesufCLQB/uF9ntu9zyzWm/8MYD3wG2/9T/lZ503A/3zGU7zP7qhq82UCx4bwu34t8G/gQeDpOub7GrjaZ/zHuETmb97FwNne61O89zXF57u4qIGf2R5vXeXfldHeOj4DDnjr+HeA77vK+/XWU+JnvpuBj0K07/8CPOszPtz7f6jx/wckeO+9V6i+C/UNreKMIACnAG80cNlS4GdAZ+A44GRcVYOv83FnF2NwB83ra1nXIuD3InKdiAzynaCqE72XI9X9YnhVRMYAzwI34g42TwDviki8z6LTgNOBI4AjgV8H+wZFpD3wCfCM9z6vBp71qaZ6FncgaQeMAuarapb3vjdq5a+crEC2p6r7cQfvE32Kfwd0A44GBgN3e/NejDuInOZt46/e/O9677kbLkm9UMu2Lscljfu95ecD9wEjvG2lAZOAX/os1g+IBXrjDmK18s4YbsAdzDbV++ZrLv+MiGTXMtT6S19EOuL20S9rm8fHMGCZz/gy3IHLn7m4/QEwEfcD5ySf8bk+8wbzmU0ESn2+K98AfwTeBpKBPrjvd5MSkUF17P9sEbmglkWH47NPVfVbIBr3nazNFyKSISKviUjvRnwb9bJE4KTgfv0FTVWXqOoiVS1R1c24L+tJ1WZ7SFX3qepW4FHg8lpW9xPgZdwvle9EZL2InFnH5n8APKGqi1W1VFVfwP2y9a13/ruqblPVfcDv69h2Xc4HVqrqy952vgTeAy70ppcAw0Wknapmef/Eh2sn0AlAVVer6mxVLVLVXbh9WH0fV/A+ixdUNU9VC3AH9nFBVDVNA36rqpmquht4ALjKZ3ohLnEUqeqhWtZxkohkA7twyX+qquYHuH3f93KDqibXMoyrY9EHcZ/9rrrWLyKxQDwuUZU7QO1VonOpeuD/o8/4Sd70oD+zWhTjkm43VT2kqguCXP6wqeq6OvZ/sqq+WcuiSVTdp+DOmP3t12Lcj56+uKScDbwjIk12fLZE4GTh6vGC5l1knuld/MoB/oD71exrm8/rLUAPf+vyvux/UNU0XHJ6DVe33KmWzfcFfu77CwX3K9V3/QFtux59gYnVtnMhlftsqje+VURmi8jYBmyjup64agVEpIe4C/k7vH38NDX3cQURiRGRh0Vkozf/alwdbUp9GxURwf2K3eJTvMWLp9wuVS2uZ1VzvQNFZ1WdoKrp9W27sYjIeGA8rnqrTt77KATa+xS3xx20/PkMGCkinXG/8l8ABnvjI73pQX9mtfgZ0Bb4xruJ4Moglw+nPKruU3BJoMZ+9X5cfaaqxd7Z8M24qti6zh4alSUC539U/roN1j9xB5pBqtoeuAt30PHle5rXB/drt06qWp5UEoH+tcy2Dfh9tV8obVX1lcPZdi3bmVVtO0mqeqsX60JVPQdIBWbh6t3B1XsGTUSScdUP872iP+GujRzl7ePvU3UfV9/OdcBpuIvOHXB1xFDzc6lBVRX3K76vT3EfYEcd2wsZEXleqt715jssqWWxycBAYLuI7MIdWK4QkYW1zP8d7iBebiTg924xVT2Au5Z1G7DESyRfeeMrve8tBP+Z1dinqrpDVa/H/eD4Ka46sk8t78FXPi6BlOsWwDJ+eT/0atv/eSJS23HjW3z2qYgMw1UjB3rzhxLA97WxWCJwfgscLyJ/EpFuACIyUERe8g5K5eK9W73Khyhcls8B8kRkCPBDP+v/hYh09Or9bgFe9ReEiPxGRI4RkTivGuMW3GniGm+W3cAAn0WeAm4SkfHiJIrI2SLie/r5YxHp5Z1V3FXbtuvxNjBaRC4VkVgvvmO9f5JEEbnMu45QjPvFU+oTb1cJ8FZQb5+OA97BJayXvEntcL+wcrwDwW3VFq2+X9oBBbgzvURc1U4wXgF+KyIpItIVV7f9Uj3LBExE4n2qqeKkjiorVb3Wp968+pBWy2J/wyWCUd7wHPAW8L1a5v837jvazfuO3go8X8dbmItLLuXXA9KrjUPwn9keINr3QO9933p4yTnbKy7xpi0SkV/VEt9SYLK450M6Anf4ThR3S3gC7qaCKO975/fuL1VdW8f+T1LV/9YSw0vAhd7/SRKuevJVr6qyChEZ4Q3R3v/RY8A63E0JTUPDdJU6XAN+7hryygfj7mbIwtXtLcP9Q0RTeddQ9eEUXD3patyXfj7uAtlnPuv1vWuo/La16Fpi+zXu11YOrlokHTjeZ/pNuGsZ2cAlXtkZwJdU3rX0OpV3Z2ym8q6hbNxpfNt69s+1+L9raDjwEe5Ol0zcWdRw3IF2Fu5uohzcXSXjvWUE9w+R5W2/k5/1zsBVTeR6+3Al7kyovc88o3D/3HnAEtw/9nqf6RfjzlqycQekDsD73vybvPdU610ZVLtzCfdr8p+4M4OduDtA4nz29/p69mGVu4b8TN/l57vULYTf+ep30ZwCZPqMR+FuXd7vfVa/r2d953kxl3/OY73x8xr6mXllD+HuDsr2ln/U+07n4Q6M1/osvx2YWEt8UbgfSQdwP6JuxOeuIe/zqb7//xWC/X6tF2c+7oaEDj7TZgO3ea9PB9Z68+3x5u0fqu+Dv0G8QEwLJCKbcbdg/i/csRjTWMTdrfaUqk4OdywtRUiflDXGmMamqutx10FMI7FrBMYY08pZ1ZAxxrRydkZgjDGtXLO4RtC5c2ft169fg5bNz88nMTGxcQNqBBZXcCyu4ERqXBC5sbXEuJYsWZKpql3qnbEpb1Fq6JCWlqYNNWfOnAYvG0oWV3AsruBEalyqkRtbS4wL+EoDOMZa1ZAxxrRylgiMMaaVs0RgjDGtXLO4WOxPcXEx27dvp6CgRtMdVXTo0IFVq1Y1UVSBSUhIwDVyaYwx4ddsE8H27dtp164d/fr1q/OgmpubS7t2TdnbZN1UlaysrIi8O8EY0zo126qhgoICUlJSmt0vaxEhJSWF6OjocIdijDFAM04EQLNLAuWaa9zGmJap2VYNGWNMKK3dncuenEJ25RRw6rBUOrSJrZi2ZMs+NuzNZ/u+g1WW6dahDVeMr9l3zszlOxGEs0c0qCPEkLNE0EBZWVmcfPLJAOzatYvo6Gi6dHEP8I0ZM4aZM2fStWtXVq5cGc4wjTENdNoj8ypeTx3Vg0cvG10xfuE/Kzt7Kz/BL2+2bcLAFPqmVF4DLCgu5eb/uG68Tx56BgmxkVctHNJE4LWHX95jVYmqjhWRe3Gdru/1ZrtLVT8IZRyhkJKSwtKlSwG49957SUpK4vbbbwdg3rx53HzzzVx99dXhDNGYkCktU3726lI2ZebTrUOtHaxV+N7IHvTt1JbH09ezL6uA13d+TVmZcuspRxIbLTw8aw3FpZUNYKoqG/bm89TVaQzsGtzNHo98shYRuPWUI4Na7pMtxbz876/8Tnt76U7yi0q92CrLzxvVg8e8BPHdzhzO+ut8bn11KZ2T4ivmyS8sqXh9wwtf0jau8rC7OTOfru3jue/c4cRERXH9C1/SPyWRqKjK6uMJHUoJtaY4I5isqpnVyh5R1YebYNthMXHiRDZv3hzuMIwJmRU7DvDuMtf9dXFpWZ3XvXZmH2Jr1kFG9U7m4293u8JdGQD07tSWhJgoPly5iyHdKvt6X5Xhuj7+8cvf8PHPJgYc14FDxTz26ToAbjihP+0SYutZwikrU95cV0Riwn66tHOJrXNSPJl5hQAM6daO7fsPVcwfGy3Ex0RzcVpll+BHpiZxytBUdmQfqjIvwNDu7r3tyy9mX35xxTbX7clj3Z48/vv1Dr7YtI+Ne/PZuDe/Yn6A4ia46bFFVA3d9963fLczx++00tLSBt2hM6xHe377veGHG5oxzUZOQTHTnlrM/oNF9c57qKjyV+pHt9Z9oH7kk7U89uk6NmflM7JXB5ZtPwC4g+OLC7cgAoNT2/HhLSdWLDP1HwtYui2bNbtzOeGh2TXWObJXMrefPpgf/Psr1u/Jo2dyG0RcUip3yl/m0rFtHM9fN447/rucDXvzuH5Cf77cvI8VOw5UWZ8qHCqB+04fwiXH9K6+uYDEREfx9DVjA56/rEwZcJerDHl+wWYOFbt9OqBLYpV9kZ6e3qB4ghHqRKDALBFR4AlVfdIrv1lErga+An6uqvurLygi04HpAKmpqTV2RocOHcjNzQWguKiY0lL/p0+qWuu0uhQXFVesvz6FhYXExsZWmT8vL4+ysrJa16GqTfIBBysvL8/iCkIkxFVUqsRGVb0bLS8vj09mz6GkrPblEqKhqAzKvKqOlZmlrNhRyKgu0STG1nNnWxvYFRPFGf1j633/vYrLmNgrhtIyOLFHIf1ilHZt4unatpDFXsxjuxVVWc+0/mUU5EfRKSEKobjK+nbll/H+igzIy2T9HjdtT84hxnVzh7OhHdzfvKISvtpVyAOvzmPuRjff3/+3in0FysDkKFLbVr1psl+C0nb/etLTN9T93hvR9BHxbM0pI7dIUWLIOlTGdUPKquyLpviOhToRTFDVnSLSFfhERFbjOgW/H5ck7sd15n599QW9pPEkwNixY3XSpElVpq9atariQbEHLhxVawBN8UBZfHw88fHxVbaTlJREVFRUrdsWEaq/p0iQnp5ucQUhnHHN+nYX019cAsBVx/bl/qlHVUx7b9YcbptbRK5P/XR1I3sns2J7dkUiAIiOEv7zk1Oq1GM3hot9Xge6zy4+y3/5oo1ZXPbkIt7fVJkgph3bj3vPrXoGX1qmjLpvFjO9JDBlSFdmr94DwB8uHc+4/p2qzB+OzzKQrTVFXCFNBKq60/u7R0TeAsapasWleBF5CpgZyhiMaUnKypSNmXkUlyp/nrW2ovzlxVuq3Lb45a4ScgtL+P4J/f1ezH1jyXaWbcsmOkq4+8whFXe+9O+c2OhJoLGN69eJhy48mtyCEkb2TiYrr5BJg7vWmC86SvjXVWmsysihY9s4Jg3uwjtLd5IUH8PYvh3DEHnkCtknLiKJQJSq5nqvTwN+JyLdVTXDm+18oMXdX3n55ZeTnp5OZmYmvXr14r777uOGG24Id1imBXhjyXZ++d/lNcrLFM58bH6VsrjoKH526pEkxtf8N88pKGH1rlyGdW/PDyYOCFm8oRAVJVx6TM179f2ZMLAzEwZ2rhi//oT+oQqrWQtl6k8F3vLqLWOA/6jqRyLyooiMwlUNbQZuDGEMTeLee++tMv7KK6+EJxATVll5hSzYkEVpWc2K+dIyiPbzHH+UCL06tgWUrdUeTvLn7aU76JQYxx/OL68GErp1SGBvbmGV7a5c+S2nnTDWbxIAmD5xAEf1aF/l7hTTeoUsEajqRmCkn/KrQrVNY8Lpon8tZFNmfsi3c/aI7pxxVN1PqCZkrmFEr+RapyfFx3Da8G6NHZpppiK7MtCYCFdUqjzyyVoOFZdWJIFTh6Vy91lDK+b5cvM+fvHGcnomt+Hl74+vsvwlTyxkT667V/2itF7cPHlgvdvskdymEd+BMZYIjDksS/eW8vjSdcTHVNb73HbqkfTrXNnEQEpSHH+dvY4/nj+iSjm4h57++OFqwD2lWn26MU3BEoExQSgtU3788td8l5FD35S2rN1RRFxMFCvuPZ24GP+N+bZLiGX+L6f4nXbjSUdw40lHhDJkY+plicCYIKzKyOGjb3cB0CkxjnZxwkXj+9eaBIxpDiwRGOOjoLiUsQ/8j7zCkioNh5UrLKl8Sv2tHx3P3LlzmTRpSFOGaEyjs0TQQLU1Q52bm0ufPn3YtWsXUVFRTJ8+nVtuuSXM0Zq6lJUpRV4bNUu27CfPexr31GFd/TamtmP/IaaN72MdDJkWwxJBA9XWDHVGRgYZGRmMGTOG3Nxc0tLSOPXUUxk2bFiYIza1Of+fn7NsW3aN8j9eMCIM0RjT9CwRNLLu3bvTvbu7x7tdu3YMHTqUHTt2WCKIEMWlZazKyKHUa2DnYFEpy7Zlc9qwVEb36eiVlXDy0NRwhmlMk2oZieDDX8GuFX4ntSktgegGvM1uR8OZDx5WWJs3b+abb75h/Pjx9c9smsTzCzbz+w9W1SifPnEAY/t18rOEMS1fy0gEESgvL48LL7yQRx99lPbt7TH+cJmzZg8Z2QUV4zOX76RXxzZVWupMjIshzRohM61Yy0gEdfxyP9QEzVBXV1xczIUXXsi0adO44IILmnTbptKenAKue+7LGuVXHtuHyX5aqzSmtWoZiSCCqCo33HADQ4cO5bbbbgt3OK2OqvLMZ5vILyxl5U7XC9Uz14zlqJ4dKubp4ue2UGNaM0sEjWzBggW8+OKLHH300Ywa5TrM+cMf/sBZZ9XSy4ZpVCt35PDA++4aQNu4aIZ1b8+EgZ1JiA2+u1JjWgtLBI3AtxnqE044AVWtfWbT6L7cvI9/pW+gTJW9XmfjAP+4YgyTh1gVkDH1sefiTbP38qItLNiQSVZ+EVHeQ14nDOzMMf3tLiBjAmFnBKbZWbnjAD/+z9cUeT2zZ+UVMXlIF564amyYIzOmeWrWiUBVm+Vj/lZ1FJzi0jLyipTsg0UAvL8ig637DnLRmF6IgCBcOq53mKM0pvlqtokgISGBrKwsUlJSmlUyUFWysrIoLS2tf2aDqnLmY/NZv+cgzP6konxwajv+dHGNDvCMMQ3QbBNBr1692L59O3v37q1zvoKCAhISEpooqsAkJCSQnx/6Lg2bo8KSUtbtziOnoBiAnEPFrN+Tx7Hdozl97OCK+Y6xp4CNaTQhTQQishnIBUqBElUdKyKdgFeBfrjO6y9R1f3Brjs2Npb+/fvXO196ejqjR48OdvUht2XLlnCHEJHueGM5by/dWaP8zP6xXDOh/s/bGBO8pjgjmKyqmT7jvwI+VdUHReRX3vgdTRCHiWD5hSXM+m5XRRLomdyGP1/iqn6S4mPIXPdNOMMzpkULR9XQecAk7/ULQDqWCFq9V77YWvEgGMAtJw/i2AEpFePp68IRlTGtQ6gTgQKzRESBJ1T1SSBVVTMAVDVDROyJn1Zs0cYsPv52FwvWZ9IzuQ2v3ngsKYnxtImzJ4GNaSoSylsZRaSHqu70DvafAD8B3lXVZJ959qtqjaYfRWQ6MB0gNTU1bcaMGQ2KIS8vj6SkpAYtG0oWl3P/wkNszikjLhqm9I7l4sFxERFXoCyu4EVqbC0xrsmTJy9R1fofsFHVJhmAe4HbgTVAd6+sO7CmvmXT0tK0oebMmdPgZUOptcf167dW6GVPLNSBd72vv3//u3rnb+37K1iRGpdq5MbWEuMCvtIAjs8ha2JCRBJFpF35a+A0YCXwLnCNN9s1wDuhisFEpl0HCnhx0RZ25xYwtm8nzh3ZI9whGdOqhfIaQSrwlvewVwzwH1X9SES+BF4TkRuArcDFIYzBRIji0jLOf3wB2/cfoqTUVUf+6aKR1iGMMREgZIlAVTcCNR79VNUs4ORQbddEnqy8Qr7dmcPKHTmcMrQrPZPb0DExjlG9k+tf2BgTcs32yWLTPCxYn8m0pxdXjP/2e8Pp3altGCMyxlRnicA0OlVl/Z48sg8V887SHURHCQ9MPYrU9vGWBIyJQJYITKN7Yt5GHvxwdcX4iF4duHxcnzBGZIypiyUC02iKSsp4e+kO/jV3Q0XZSzeMZ1Bq5N2bbYypZInANJr/rdrNL99YXjF+y8mDOGFQ5zBGZIwJhCUC02BlZcpjn65jT67rJ/jbnQeIi4li/i8nk5IYR3RU8+knwpjWzBKBabBVu3J47NN1dGgTS1yMezZx6qgepLaPrP4fjDF1s0RggvbUvI28vyKDA4dc5zHv3jyBvimJYY7KGNNQlghM0J75bBPRUcIRXZOYMDCFPnZLqDHNmiUCE5DCklIu+ddCdmQfIjOviF+fPZTvnzgg3GEZYxqBJQJTrwMHi1m2PZtl2w8wZUhX+nRqy3mjeoY7LGNMI7FEYOo0b+1ern72i4rxB6YeRY/kNmGMyBjT2CwRmFpt23eQN7/eTlx0FL85ZyjdOrSxJGBMC2SJwNTq8qcWsX3/IY7p15GrjusX7nCMMSFiicDUsHRbNvPX7mX7/kP84MT+/HDSwHCHZIwJIUsEpoZfvL6MdXvyiI0WLkrrTadE//0IG2NaBksEpsK+/CKeW1nI+r0l/GTKQG6eMpD4mOhwh2WMCbGQ9Vlsmp+PVu5i7vYS+ndO5MyjulsSMKaVsDMCA8Cybdnc9dYK2sbAp7edhNfXtDGmFaj1jEBEuojIMD/lw0WkS2jDMk3trW92AHDOgFhLAsa0MnVVDf0N8HfA7wU8FugGRCRaRL4RkZne+PMisklElnrDqOBCNo3p6637Sbv/E15ctIVx/Ttx1gC7MGxMa1NXIjhaVedWL1TVj4ERQWzjFmBVtbJfqOoob1gaxLpMI8rMK+TdpTvJPlTMNcf145enDw53SMaYMKjrGkFsA6dVEJFewNnA74HbgojLhNiybdmc948FABzVsz33fM/VAqZvDmNQxpiwEFX1P0HkfeAfqvpBtfIzgZ+q6pn1rlzkDeCPQDvgdlU9R0SeB44DCoFPgV+paqGfZacD0wFSU1PTZsyYEcz7qpCXl0dSUuT1mRvOuPYeLOO9jcXM217CNcPiGNwpmh5JUWGPqy4WV3AiNS6I3NhaYlyTJ09eoqpj651RVf0OwJHAWuB54Cfe8IJXdmRty/ksfw7wuPd6EjDTe90dECDeW9899a0rLS1NG2rOnDkNXjaUwhnXOX+dr33vmKmn/Dm9xjTbX8GxuIIXqbG1xLiAr7Se46uq1l41pKprReRo4ArgKK94LnCjqhYEkIwmAOeKyFlAAtBeRF5S1Su96YUi8hxwewDrMo1g3e5c5q3L5LuMHK4Y34dfnGbXBIwx9TxHoK7K5rmGrFhV7wTuBBCRSbiqoStFpLuqZoi7R3EqsLIh6zfBu+edb1m4MQsROGdEdzpa0xHGGOpIBCJyInCzN/q4+rmDqIFe9p5DEGApcFMjrdfU4lBRKQ99tJpvtu3n8nF9+PXZQ0mMt2cJjTFOXUeDPwMTvdfzgWMauhFVTQfSvddTGroe0zDz1u3l+c830zO5DeeO7GFJwBhTRV1HhHTgPqAMWNgk0ZhG9/T8jTzz2SZiooRPf34SCbHWfpAxpqq6Lhb/UkQGA9Gq+l0TxmQa0XMLNhMlwk9PHmRJwBjjV30Xi9c0VSCmca3dncuFj39ObmEJd581lB9MHBDukIwxEcqaoW6BikrK+HjlLnILS7h58kAuSusV7pCMMRHMrhq2QNc//yWfrc+ke4cEbrf2g4wx9airGeo7RWR0UwZjDk9pmbJgfSZfbNrHyUO68sRVaeEOyRjTDNR1RrAJuEVERgLLgA+BWaq6v0kiM0H7+Ntd/OjlrwG4bFwfRvRKDnNExpjmoK67hmYAMwC8M4MzgDdFJBr4H/CRqn7RJFGaeu3JLeDed78lPiaKt340gaHd24U7JGNMMxHQxWJV/UZV/6iqk3GNyX0LfD+kkZmgPLdgM3tyC5l4ZBeG9WhvvYwZYwIW9MViVc0B/usNJgJ8viGTf6ZvYEi3djxp1wWMMUGy20dbgCfnbQRg2vg+diZgjAmaJYJm7m+friN9zV4uGduLq47rF+5wjDHNUIMSgYgMaexATMO8vHgrAFce2zfMkRhjmquGnhHMatQoTINc+sRCduUU8Ouzh9qtosaYBqurP4K/1jYJsKNOmG3OzGfxpn30TG7DBWOsCQljTMPVddfQdcDPcZ3MV3d5aMIxgcgvLGHSw+kAPHzxSDpZT2PGmMNQVyL4Elipqp9XnyAi94YsIlMnVeWJuRsAuPLYPhw7oFOYIzLGNHd1JYKLAL+d1Ktq/9CEY+qzYscB/jp7PQC3nTrYbhc1xhy2upqY2NeUgZj65RYUc/vrywB458cTrErIGNMoQv4cgYhEi8g3IjLTG+8vIotFZJ2IvCoidjQL0HvLMli7O4+BXZMY2duu1xtjGkdTPFB2C7DKZ/wh4BFVHQTsB25oghiavXW7c7nrrRV0h0/tAAAgAElEQVQkt43lk59NDHc4xpgWJOBEICKJwa5cRHoBZwNPe+MCTAHe8GZ5AZga7Hpbo/IHxy4d29uuCxhjGpWoat0ziByPO5AnqWofr3+CG1X1R/WuXOQN4I9AO+B24FpgkaoO9Kb3Bj5U1aP8LDsdmA6QmpqaNmPGjGDeV4W8vDySkpIatGwoBRNXdmEZt845xOCOUdw5vk3ExNWULK7gRGpcELmxtcS4Jk+evERVx9Y7o6rWOQCLgd7ANz5lKwNY7hzgce/1JGAm0AVY7zNPb2BFfetKS0vThpozZ06Dlw2lYOJ6/att2veOmfrvhZtDF5CnJeyvpmRxBS9SY2uJcQFfaT3HV1UNrBlqVd1WrTqiNIDFJgDnishZQALQHngUSBaRGFUtAXoBOwOJobX6bF0mt7++jDax0Uwb1yfc4RhjWqBArhFs86qHVETiROR2ql789UtV71TVXqraD7gMmK2q04A5uGcUAK4B3mlY6K3Dq19tA+CRS0cRFWXXBowxjS+QRHAT8GOgJ7AdGOWNN9QdwG0ish5IAZ45jHW1aNv2HeS9ZTsZ0asDZxzVLdzhGGNaqDqrhrz+ia/yfsk3mKqmA+ne643AuMNZX2tQUlrGgx+uBuDWUwaFORpjTEtW5xmBqpYC5zVRLMbH/1bt4f0VGcTFRDFxUJdwh2OMacECuVi8QET+DrwK5JcXqurXIYuqlcvMK+S3764E4PNfTSEm2jqSM8aETiCJ4Hjv7+98yhT3YJgJgZcWbWF3TiGTB3ehc1J8uMMxxrRw9SYCVZ3cFIEYZ9HGLB793zoGdU3i2WuPCXc4xphWoN5EICL3+CtX1d/5KzeH59nPNgEwbXwfa0rCGNMkAqkayvd5nYB7Yrje5whM8B5PX8+s73YzdVQPrp1gXT4YY5pGIFVDf/YdF5GHgXdDFlErVVamvLhwC4AlAWNMk2rI7ShtgQGNHUhrd/ETC8k4UMCdZw5hlPU1YIxpQoFcI1iBu0sIIBrXcNz9oQyqtdm27yBLtuxnQOdELjvG2hMyxjStQK4RnOPzugTY7TUYZxpBQXEpJ/7fHADuPXc4HdrGhjkiY0xrE0jV0AOqusUbdqhqiYi8GPLIWokPV2YAcNkxvTlhYOcwR2OMaY0CSQTDfUdEJAZIC004rc8/0zcAcNtpR1rrosaYsKg1EYjInSKSC4wQkRxvyAV2Y01HN4pVGTms3Z3HFeP70LVdQrjDMca0UrUmAlX9I9AB+LeqtveGdqqaoqp3Nl2ILddzC9zDYxeM7hnmSIwxrVmdF4tVtczro9g0si92lfDa0u2cdGQXxvbrFO5wjDGtWCDXCBaJiDV608g+3FQMuIvExhgTToHcPjoZuElENuOamxBAVXVEKANryRZvzGLTgTJunDiAM4/uHu5wjDGtXCCJ4MyQR9HKvLtsJwAXpvUKcyTGGBNA1ZCqbgF6A1O81wcDWc74t23fQV5evJUhnaI4MrVduMMxxpj6D+gi8ltch/PldwrFAi+FMqiW7J2lOwCY0seeIDbGRIZAftmfD5yL1xy1qu4E6v0pKyIJIvKFiCwTkW9F5D6v/HkR2SQiS71h1OG8geakpLSMh2etJT4mirGp0eEOxxhjgMCuERSpqoqIAohIYoDrLsRVJ+WJSCzwmYh86E37haq+0YB4m7Un5m0E4KaTjiBKdoY5GmOMcQI5I3hNRJ4AkkXkB8D/gKfqW0idPG801hu0jkVatI178/jTx2sAuMRuGTXGRBBRrf/YLCKnAqd5o7NU9ZOAVi4SDSwBBgL/UNU7ROR54DjcGcOnwK9UtdDPstOB6QCpqalpM2bMCGSTNeTl5ZGUlNSgZRvTE8sLWLizlDvHJTC4U3TExFWdxRUciyt4kRpbS4xr8uTJS1R1bL0zqmq9A9ANd53ge0C3QJaptnwyMAc4CuiOexYhHngBuKe+5dPS0rSh5syZ0+BlG0tpaZkO/c2H2veOmVpSWqaqkRGXPxZXcCyu4EVqbC0xLuArDeAYHchdQ98HvgAuAC7CPWl8fTBZSVWzgXTgDFXN8GIsBJ4DxgWzrubolS+3crColIcvHkm0tTBqjIkwgVws/gUwWlWzAEQkBfgceLauhUSkC1Csqtki0gY4BXhIRLqraoaICDAVWHlY76AZeM97gGziIOtvwBgTeQJJBNuBXJ/xXGBbAMt1B17wrhNEAa+p6kwRme0lCQGWAjcFGXOzsmTLPhZt3Mc1x/Wla3tratoYE3kCSQQ7gMUi8g7urp/zgC9E5DYAVf2Lv4VUdTkw2k/5lIaH2/y8+bV7gOxS64vYGBOhAkkEG7yhXHmnNNY+Qj2enr+RlxdvZUyfZIb1aB/ucIwxxq96E4Gq3tcUgbQ0qspLi7YA8MDUo8McjTHG1K7eRCAiY4G7gb6+86s1Q12n2av3sDnrIJeO7W1nA8aYiBZI1dDLuDuHVgBloQ2n5Xj+880A3HHmkPAGYowx9QgkEexV1XdDHkkLknHgEPPXZdI3pS2dEuPCHY4xxtQpkETwWxF5GtccREVTEKr6ZsiiauaemOsal3vwAqs9M8ZEvkASwXXAEFyjceVVQwpYIvBDVSt6IBvTNznM0RhjTP0CSQQjVdVuewnQm1/vYF9+EfdPPYr4GOtzwBgT+QJphnqRiAwLeSQtxHvL3dnAyUO6hjkSY4wJTCBnBCcA14jIJtw1AsF1N2AV4NUs3ZZN+pq9XDK2Fz2S24Q7HGOMCUggieCMkEfRQpQ3LnflsX3DHIlpsUpLoLQQykqgpFo3Hm06QnS1vrCLC6CkALTMLQMQ2xbiq7VvX1YGBzNDF3ebjlB8EMpKQ7cN02CBPFm8RUROAAap6nNeg3GR13tDmGXmFfLMZ5sYnNqOEb3sIrEJkZcugE1z/U/rfxJc43Ond1kZPHo05O+pOl9sW7h1JSSmVJa9/zNY8nyjh+tP5+F3ApOaZFsmMIE8WfxbYCwwGNd/QCzwEjAhtKE1Lx9/uwuAayf0C28gpnlShdwMKMyDwpza5/FNAmc8WHkGsHYWbJgN274AiaJdzhpYfaBqEjj7z3BgO3z2CKx4HXr5dFy1fjb0HAujLm/89/bRXe4sBiAqhi57F8D2r4JfT0yCO7tBoNtREBNfcx5V2L0S2qZA+x6HFXat9m2CwlzQUpAoNyQkQ95uN71jPygtcmV7voOoGIhvB4eyoW1HN//BLBcrQJch7gytrBR2rag8cyt/28V5hFogVUPn41oR/RpAVXeKiDU4V82DH6wGYOqonmGOxDRL37wI7/4k8PkTkuHYH1aOJ6XCuo/hmVMBSAPvP9aT2BWO+b5LNAsfh4/uqLnO8dPdPI0texsseNS97ns8qZvmwdMnH946T7gNTvltzfI1H8CMKyC+PdyxBaICuR8mCCv/C2/U0y9XdJxLBIEacSlc8KT7Drx3S43J7Y/28z4bWSCJoEhVVUQUQEQSQxxTs/Pl5n3kFpZw1bF9aRNnt4waH2VlsPINKD7k6ulLi9yBorrP/1b5+oKnXJ26PzHxkNDBDb4GnwXXvFdx3WD58uWMGDHCrad9z8pfz/FJMH0O5OysunxUNPQ5voFvsh5TfgNDznZxREWzfNbLLrZgLH0Zvn0LBp8N2VvcAT/ZT9Puaz50fwtzXPKpbT/60X3nWvhqU90z1VV9dsJtkL/XHdDL9TkOti70P/+0N+Dzv8LGdPjqOXeWltgFpv6zymy5mw4GFP/hCCQRvCYiTwDJIvID4Hrg6dCG1bw8Pmc9AFeMtz4HTDXLXoF3fhT4/O17wtEXgwTZpWlUNPSfWDG6b0csDJrkf97U4W5oKtEx0LuyR9p9KWm1x1abNh1dIjj+J64KbN7/wcxb/c/bsZ+rAvs0uIaTBwOsDWKBnmNh7xoo8vrtGnsd5GRUTQTDznNVQNsWwegrYeVbUJzvzgIGnQo5O9xZQPl7KS/3UbwjPaj30RCBXCx+WEROBXJw++oeVf0k5JE1E6rKih05nDG8G0O7WyujrcrqDzh6+cOQ/x5kb4VkP3eLbVtc+XrKr2H2A+7X+9l++nOKS3R1/sEmgdag11j4TabbP32OhWNuqKxjry6xMxTlu7OwIHy+8HOOPy6As6KEDu6sTsSLQd0dUQkd3FnKLze56wLFhyCpK4yb7ur9o+PgjIfcNYIYr7fCtGvd96H8bqqk8Dx/FMjF4odU9Q7gEz9lrd7OAwVk5hVy/MCU+mc2LcuMy0kB2LeksqxNJ//znv4HOOoi94vw5HugffemiLBlKb8wLgLtutU9b5tkNwShKD6l4Z9LtE9VXVvvO5Dg/TCUaHfGBjVv24WwHfx9BVI1dCpQ/aB/pp+yVumhD91F4pF2y2jLU1YGb//Q1UlX5+/XaN8T4Lr3617njz5vnNiMaUS1XlIXkR+KyApgsIgs9xk2AcvrW7GIJIjIFyKyTES+FZH7vPL+IrJYRNaJyKsi0mzbaVZVFqzPJD4myjqfaYn2bYDlM6Agx/0a9R1i4mDgqWSmjIfrZ8GYa+B7j4U7YmMapK4zgv8AHwJ/BH7lU56rqvsCWHchMEVV80QkFvhMRD4EbgMeUdUZIvIv4Abgn3WtKFK9vHgrWflF3H/ecGKjG/k2NdM0sjbAfy7xX59cUuD+XvCku2/dj5Xp6UzqMx76jA9hkMaEVq2JQFUPAAeABj1hoqoKlD8JEesNCkwBrvDKXwDupRkmgrIy5Z2lOwA482ir7414xYf8H+xXvQdZ62HEZe7uluqSukFXa3PRtGyitV15b4yVi0QDS4CBwD+APwGLVHWgN7038KGq1vi5JSLTgekAqampaTNmzGhQDHl5eSQlNX6LGM+sKGT+jhJO7BnDDUf7ecIxTHEdrpYYV1zhPsYvvpHoMv8P+RTGdWLhcc826G6dlri/Qi1SY2uJcU2ePHmJqo6tb75ALhY3mKqWAqNEJBl4Cxjqb7Zaln0SeBJg7NixOmnSpAbFkJ6eTkOXrcudCz+lfQI8dNXEBrU0Gqq4DlezjStrAxTluUf/o2KrPtm5bQOUFbkHfpJSaywa32O0q94JRVxhEqlxQeTG1prjCmkiKKeq2SKSDhyLezAtRlVLgF7AzjoXjkC7cwrIOFDAPecMs+amI8GulfCvepq+ikuCSXe6i7zGmCpClgi8VkqLvSTQBjgFeAiYA1wEzACuAd4JVQyh8tIidzvhyN4d6pnTNKrSYlj6H8jbU/V+7IxlNec992/QsX/lePselgSMqUUozwi6Ay941wmigNdUdaaIfAfMEJEHgG+AZ0IYQ6NTVV75YhsAw3tYImhS6z6B935a9zxtU1zLjiMu9d86pTGmhpAlAlVdjmu1tHr5RmBczSWah1nf7SYzr5C7zxpKQqw1MBdSy2bApvkADN6VAWt97lr+/qeQckTleGyi/eI3poGa5BpBS/Lal+5s4NRhNS86mkakCh/d6aqDEjrQsbAA4hNcc8o9RrshyhKxMY3BEkEQNuzN49PVe7hgdE/6dbbWuBtNYR68fZPruKNcWSkc2ucaZzvmBhZF6B0dxrQE9jhsEN5fngHAeaOt85lGteVz92BXYY5LAOUtMR4xBY60LrONCTU7IwjQgUPF/OWTtQzsmsRJR3YJdzjN18o34Y3rXJO94v0OKSl0r6/9wH/rjMaYkLJEEKCvt+wHYOqoEPWD2lKVlVb25QouCYA7+I+5urK86zBLAsaEiSWCAC3dlk2UwHUT+tc/s6n01k2w4rWa5WNvgDP+0PTxGGNqsEQQgLIy5a1vdjCoazsS422X1SlzPezbWDm+cQ70nQAjLqksi4qFoy9q+tiMMX7ZUS0A89btZeu+g1xgF4nrpgrPn1W1KghcF41p14YlJGNM/SwRBOCZzzYB8JtzrDniKrZ9UbV5h6I8lwROvN31wwquaedU/235G2MigyWCeuzNLWT+uky6d0igY6I9uVpBFWZMg/w9Vcsl2lUDdRkcnriMMUGzRFCPx9PXA/DnS0aGOZIIsGkerPyve11a7JLAqb+DkVdUzhObAPHtwhOfMaZBLBHU44MV7iGyMX06hjmSCDDnD7Dja2iT7MY79oOh34Mke67CmObMEkEd3v5mB7tzCrnnnGGtr4G5/Ex452bXzIOWubKdS2HcD+CMP4Y3NmNMo7JEUIeZXpMSpx/VLcyRhMEXT8LaD93rnmnuSeAjJsPIy8IblzGm0VkiqMWybdn8b9VuLhjdk56toRey4kPwzwmcmL0NPouG0sLKaZfPgKSu4YvNGBNSlghqkb5mLwDXHN8vvIGEUkkRFBxwrzOWwr4NZHU5ga6Dxriy3Azoe7wlAWNaOEsEtVi+PZuBXZMY2Ts53KGEzjOn1Ojmcf3A6+l62oVhCsgYEw6WCPwoKC51/Q6MaWFPEudnwS7vwF9c4JLA8PNdExAAHXpRlNEKqsGMMVVYIvDj6fmurZyxfTuFOZJG9u7NsOaDqmXjfwh9xleOZ6Q3aUjGmPALWSIQkd7Av4FuQBnwpKo+JiL3Aj8A9nqz3qWqH/hfS9MrKS3jiXkbSYqP4ZKxvcIdzuFb+h/I3QVRMbBlAQw5B47/iZsW2xa6jwhvfMaYsAvlGUEJ8HNV/VpE2gFLROQTb9ojqvpwCLfdYPPXZZJbUMIFY3oSE93MO3DbvgTe/mHVsqMvgj7HhiceY0xEClkiUNUMIMN7nSsiq4CIr3R/6KPVADwwtZk2lLZ9CXz5lGsL6MC2yvI7d7jO3mPtGoAxpqom+ckrIv2A0cBir+hmEVkuIs+KSMS03bAnp4DVu3I5MjWJtnHN9PLJwr+79oC2LoQD213ZpLtc71+WBIwxfoiqhnYDIknAXOD3qvqmiKQCmYAC9wPdVfV6P8tNB6YDpKamps2YMaNB28/LyyMpKbAuEJ9YVsDCjFLuHp/AoI6hbVIimLhqE1VawNErfk9scU5FWZtDO8lKGct3w+8IW1yhYHEFJ1LjgsiNrSXGNXny5CWqOrbeGVU1ZAMQC3wM3FbL9H7AyvrWk5aWpg01Z86cgOcd9/tP9Mi7P9CiktIGby9QwcRVq/WzVX/bXvW5s1VfuaJy2PRZeOMKAYsrOJEal2rkxtYS4wK+0gCO1aG8a0iAZ4BVqvoXn/Lu6q4fAJwPrAxVDMGY8cVWducU8ptzhhEbqReJVeH5s2HXCjdeWuT+XvoitImYGjZjTDMTyorwCcBVwAoRWeqV3QVcLiKjcFVDm4EbQxhDwGavdh2sTB3VI8yR+CjKd62AlsvPdLeAHjEFugxxZZ0HWRIwxhyWUN419BkgfiZFzDMDvpZtz2bqqB6kJMWHOxRHFf51IuzbUHPalF+7FkGNMaYRNNNbYxrXyh0H2J1TGBntCh3aD1sXQUGOSwKjrnQNv5Vrkww9xoQvPmNMi2OJAPi/j9cAcEy/CGhS4tPfwVfPVo6Pv9Ge/jXGhFSrTwSZeYXMW7uXCQNTOKpnh6bd+IHt8OUzVe/vXzsL+hwHZzzo+v5NOaJpYzLGtDqtPhG8vGgrAJce06fpN/7cWZC9pWb5hFugx6imj8cY0yq16kRQUFzKX2ev44guiZw7sgnuFsrPgrdvgradGbx7T2USOGIKTHujcr6oVtY/sjEmrFp1Ili4IYvSMuW4I1KaZoOz7oZ1swDoGJ8C8e2hMAdOe8AO/saYsGnVieCed92zbHeeOTQ0GyjMhVevhEPZbjzDe5xColh03LNMmjQpNNs1xpggROgjtKG3N7eQbfsOMap3MonxIcqHWz6HjekQlwhJqTDodGjfE67/ODTbM8aYBmi1ZwS/fts103DXWY1wNrD7O/jncd6IVFbzlJWCRMG0110y8LUh/fC3a4wxjaBVJoLcgmJmr95DSmIcaX2DaJ5BFQoOQFGea/6h3Gd/qXwdkwDH/ahyvMvQmknAGGMiSKtMBOc//jnFpcpPpgwkOspfKxi1WPI8zLy17nnGT4eT7zms+Iwxpim1ukTw6ardrN+Tx7EDOnHZOD/PDmQsh8y1/hee+1Dl63P/VvWXfnQcRMfDEZMbN2BjjAmxVpUI8gpLuOGFrwD45RlDSIitdstmYR48cWL9K2rfE8ZcHYIIjTGm6bXsRLD+U/pseQvmfw3Ad5v38cPovZw2PJXRW9ZA9Yd6c3dVvp6eDrF+6vbjk6zO3xjTorTsRLDmAwZsehE2udFxwLhYYK031CahA/QYHfr4jDEmArTsRHDGg8xLOJ2ybkN4Z+lOPliZwY8mDeSWKYNqXyYqBqJb9m4xxhhfLfuIFx1LWXQcj83dxtpdufRPTeHUo/tCbEK4IzPGmIjRohPB3z5dxysLD7L7IFw/oR93nz0s3CEZY0zEadGJoEu7eHokRpF2RCqXjO0d7nCMMSYitehEcNm4PnQ7uJFJk+zCrzHG1CZkjc6JSG8RmSMiq0TkWxG5xSvvJCKfiMg6728QbTwYY4xpbKFsfbQE+LmqDgWOBX4sIsOAXwGfquog4FNv3BhjTJiELBGoaoaqfu29zgVWAT2B84AXvNleAKaGKgZjjDH1E1UN/UZE+gHzgKOAraqa7DNtv6rWqB4SkenAdIDU1NS0GTNmNGjbeXl5JCUlNWjZULK4gmNxBSdS44LIja0lxjV58uQlqjq23hlVNaQDkAQsAS7wxrOrTd9f3zrS0tK0oebMmdPgZUPJ4gqOxRWcSI1LNXJja4lxAV9pAMfpkPZQJiKxwH+Bl1X1Ta94t4h096Z3B/aEMgZjjDF1C+VdQwI8A6xSVZ+eW3gXuMZ7fQ3wTqhiMMYYU79QPkcwAbgKWCEiXq/t3AU8CLwmIjcAW4GLQxiDMcaYejTJxeLDJSJ7qdlodKA6A5mNGE5jsbiCY3EFJ1LjgsiNrSXG1VdVu9Q3U7NIBIdDRL7SQK6aNzGLKzgWV3AiNS6I3Nhac1whvVhsjDEm8lkiMMaYVq41JIInwx1ALSyu4FhcwYnUuCByY2u1cbX4awTGGGPq1hrOCIwxxtTBEoExxrRyLTYRiMgZIrJGRNaLSJM3dS0iz4rIHhFZ6VPmty8Gcf7qxbpcRMaEKKag+ohoqri8bSWIyBcissyL7T6vvL+ILPZie1VE4rzyeG98vTe9XwhjixaRb0RkZqTE5G1vs4isEJGlIvKVVxYJn2WyiLwhIqu979px4Y5LRAZ7+6l8yBGRW8Mdl7etn3nf+ZUi8or3v9C037FAGiRqbgMQDWwABgBxwDJgWBPHMBEYA6z0Kfs/4Ffe618BD3mvzwI+BATXd8PiEMXUHRjjvW4HrAWGhTsub1sCJHmvY4HF3jZfAy7zyv8F/NB7/SPgX97ry4BXQxjbbcB/gJneeNhj8raxGehcrSwSPssXgO97r+OA5EiIyye+aGAX0DfcceGa5t8EtPH5bl3b1N+xkO7wcA3AccDHPuN3AneGIY5+VE0Ea4Du3uvuwBrv9RPA5f7mC3F87wCnRmBcbYGvgfG4Jypjqn+uwMfAcd7rGG8+CUEsvXAdKE0BZnoHhrDG5BPbZmomgrB+lkB778AmkRRXtVhOAxZEQly4RLAN6OR9Z2YCpzf1d6ylVg2V79xy272ycEtV1QxwHfcAXb3yJo/XO6UcjfvlHRFxeVUwS3Et0n6CO6vLVtUSP9uviM2bfgBICUFYjwK/BMq88ZQIiKmcArNEZIm4/jsg/J/lAGAv8JxXnfa0iCRGQFy+LgNe8V6HNS5V3QE8jGt3LQP3nVlCE3/HWmoiED9lkXyfbJPGKyJJuObBb1XVnLpm9VMWsrhUtVRVR+F+hY8Dhtax/ZDHJiLnAHtUdYlvcThjqmaCqo4BzsR1BTuxjnmbKrYYXJXoP1V1NJBP3d3RNvV3Pw44F3i9vln9lDV6XN41ifOA/kAPIBH3eda27ZDE1VITwXagt894L2BnmGLxVVtfDE0WrwTXR0RY9qOqZgPpuLrZZBEpbyXXd/sVsXnTOwD7GjmUCcC5IrIZmIGrHno0zDFVUNWd3t89wFu45Bnuz3I7sF1VF3vjb+ASQ7jjKncm8LWq7vbGwx3XKcAmVd2rqsXAm8DxNPF3rKUmgi+BQd6V9zjcqeC7YY4Jau+L4V3gau9OhWOBA+Wnq41JJOg+IpokLi+2LiKS7L1ug/sHWQXMAS6qJbbymC8CZqtXcdpYVPVOVe2lqv1w36HZqjotnDGVE5FEEWlX/hpX772SMH+WqroL2CYig72ik4Hvwh2Xj8uprBYq334449oKHCsibb3/z/L91bTfsVBelAnngLvqvxZXz3x3GLb/Cq7OrxiXxW/A1eV9Cqzz/nby5hXgH16sK4CxIYrpBNxp5HJgqTecFe64vG2NAL7xYlsJ3OOVDwC+ANbjTufjvfIEb3y9N31AiD/PSVTeNRT2mLwYlnnDt+Xf8Qj5LEcBX3mf5dtAxwiJqy2QBXTwKYuEuO4DVnvf+xeB+Kb+jlkTE8YY08q11KohY4wxAbJEYIwxrZwlAmOMaeUsERhjTCtnicAYY1o5SwSmWfFatvyRz3gPEXkjRNuKFZEl9c/ZtESkn/i0amvM4bJEYJqbZFwLjIB7ulZVL6pj/sNxAvB5iNZtTMSwRGCamweBI7w25f/k++tYRK4VkbdF5D0R2SQiN4vIbV7jZ4tEpJM33xEi8pHXWNt8ERlSy7bOwDVFXMFrGO95r+34FSLys7rWKSKpIvKWuH4WlonI8V75bd46VorIrV5ZP3Ht9z8lrn36Wd5T1ohImrf8QuDHPvEMF9ePw1Jx7eYPasydbVqJUD0tZ4MNoRio2bR3xTiuHff1uL4WuuBaZrzJm/YIrpE9cE+QDvJej8c9pu9vW18AbauVpQGf+Iwn17VO4FWf7Ubj2oZJwz2tmggk4Z4MHu29lxJglDf/a8CV3uvlwEne6z/5vOe/AdO813F47drbYEMwQ3mjRsa0FHNUNRfIFZEDwHte+QpghNfy6vHA665pF8A90l+FiPQA9qnqwbUSFGsAAAGdSURBVGqTNgIDRORvwPu4ZqDrWucU4GpwrasCB0TkBOAtVc33tvUmcCKuHZlNqrrUW3YJ0E9EOuASzlyv/EUqW6hcCNwtIr2AN1V1XaA7yphylghMS1Po87rMZ7wM932PwrX1Pqqe9ZyJ6wSkClXdLyIjcZ2H/Bi4BLg1wHWW89eUcDnf+EuBNt78ftuCUdX/iMhi4GzgYxH5vqrODjAOYwC7RmCan1xc1U+DqOt/YZOIXAwVfdOO9DNrjesD3vydgShV/S/wG1zXn3Wt81Pgh155tIi0B+YBU70WJxOB84H5dcScTeWZBMA0n3gGABtV9a+4M4oRge4LY8pZIjDNiqpmAQu8i6x/auBqpgE3iEh5y53n+U4UkWhcff9qP8v2BNLF9aT2PK4b1LrWeQswWURW4Kp6hqvq196yX+B6iHtaVb+pJ+brgH94F4sP+ZRfCqz04hkC/Lue9RhTg7U+akw13i/vK1X1pnDHYkxTsERgjDGtnFUNGWNMK2eJwBhjWjlLBMYY08pZIjDGmFbOEoExxrRylgiMMaaV+3/qeTzN1gDSYAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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",
    "plt.plot(t,T1,t,T2)\n",
    "plt.xlabel('time / seconds')\n",
    "plt.ylabel('temperture 1 / °C')\n",
    "plt.title('TCLab Step Test Data for P1 = 4.0 watts, u1 = 0.5')\n",
    "plt.legend(['T1','T2'])\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.2 Loading previously saved experimental data](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2-Loading-previously-saved-experimental-data)",
     "section": "2.3.2 Loading previously saved experimental data"
    }
   },
   "source": [
    "**Exercise:**  From this step test data and the steady-state analysis described above, estimate the value of $U_a$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.3.2.1 Deviation Variables](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2.1-Deviation-Variables)",
     "section": "2.3.2.1 Deviation Variables"
    }
   },
   "source": [
    "### 2.3.2.1 Deviation Variables\n",
    "\n",
    "In examining the response of the temperature control laboratory, we see the temperature is a deviation from ambient temperature, i.e.,\n",
    "\n",
    "\\begin{align}\n",
    "T_1' = T_1 - T_{amb}\n",
    "\\end{align}\n",
    "\n",
    "For process control purposes, we are often interested in the deviation of a process variable from a nominal value. In this case the choice of deviation variable is clearly obvious which is designated $T_1'$. From the steady state equation we see\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{T}_1' = \\bar{T}_1 - T_{amb} = \\frac{P_1}{U_a}\\bar{u}_1\n",
    "\\end{align}\n",
    "\n",
    "which is a somewhat simpler expression.\n",
    "\n",
    "Let's see what happens to the transient model. Substituting $T_1 = T_{amb} + T_1'$ into the differential equation gives\n",
    "\n",
    "\\begin{align}\n",
    "C_p\\frac{d(T_{amb}+T_1')}{dt} & = U_a(T_{amb} - (T_{amb} + T_1')) + P_1u_1\n",
    "\\end{align}\n",
    "\n",
    "Expanding these terms\n",
    "\n",
    "\\begin{align}\n",
    "C_p\\underbrace{\\frac{dT_{amb}}{dt}}_{0} + C_p\\frac{dT_1'}{dt} & = U_a(\\underbrace{T_{amb} - T_{amb}}_{0} - T_1') + P_1u_1\n",
    "\\end{align}\n",
    "\n",
    "we see several terms drop out. The derivative of any constant is zero, and see a cancelation on the right hand side, leaving\n",
    "\n",
    "\\begin{align}\n",
    "C_p\\frac{dT_1'}{dt} & = - U_aT_1' + P_1u_1\n",
    "\\end{align}\n",
    "\n",
    "One last manipulation will bring this model into a commonly used standard form\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dT_1'}{dt} & = - \\frac{U_a}{C_p}T_1' + \\frac{P_1}{C_p}u_1\n",
    "\\end{align}\n",
    "\n",
    "Now let's recall some facts about first order differential equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.3.2.2 First Order Linear Differential Equations](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2.2-First-Order-Linear-Differential-Equations)",
     "section": "2.3.2.2 First Order Linear Differential Equations"
    }
   },
   "source": [
    "### 2.3.2.2 First Order Linear Differential Equations\n",
    "\n",
    "A standard form for a single differential equation is\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dx}{dt} & = ax + bu\n",
    "\\end{align}\n",
    "\n",
    "where $a$ and $b$ are model constants, $x$ is the dependent variable, and $u$ is an exogeneous input.  \n",
    "\n",
    "#### Steady State Response\n",
    "\n",
    "For a constant value $\\bar{u}$, the steady state response $\\bar{x}$ is given by solution to the equation\n",
    "\n",
    "\\begin{align}\n",
    "0 & = a\\bar{x} + b\\bar{u}\n",
    "\\end{align}\n",
    "\n",
    "which is\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{x} & = -\\frac{b}{a} \\bar{u}\n",
    "\\end{align}\n",
    "\n",
    "#### Transient Response\n",
    "\n",
    "The transient response is given by\n",
    "\n",
    "\\begin{align}\n",
    "x(t) & = \\bar{x} + \\left[x(t_0) - \\bar{x}\\right] e^{a(t-t_0)}\n",
    "\\end{align}\n",
    "\n",
    "which is an exact, analytical solution.\n",
    "\n",
    "#### Apply to Model Equation\n",
    "\n",
    "We now apply this textbook solution to the model equation. Comparing equations, we make the following identifications\n",
    "\n",
    "\\begin{align}\n",
    "T_1' \\sim x \\\\\n",
    "-\\frac{U_a}{C_p} \\sim a \\\\\n",
    "\\frac{P_1}{C_p} \\sim b \\\\\n",
    "u_1 \\sim u\n",
    "\\end{align}\n",
    "\n",
    "Substituting these terms into the standard solution we confirm the steady-state solution found above, and provides a solution for the transient response of the deviation variables.\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{x} = -\\frac{b}{a}\\bar{u} \\qquad & \\Rightarrow \\qquad \\bar{T}_{1}' = \\frac{P_1}{U_a}\\bar{u}_{1} \\\\\n",
    "x(t) = \\bar{x} + \\left[x(t_0) - \\bar{x}\\right] e^{a(t-t_0)} \\qquad & \\Rightarrow \\qquad\n",
    "T_1'(t) = \\frac{P_1}{U_a}\\bar{u}_{1} + \\left[T_1'(t_0) - \\frac{P_1}{U_a}\\bar{u}_{1}\\right]e^{-\\frac{U_a}{C_p}(t-t_0)}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.3.2.3 Plotting the Analytical Solution](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2.3-Plotting-the-Analytical-Solution)",
     "section": "2.3.2.3 Plotting the Analytical Solution"
    }
   },
   "source": [
    "### 2.3.2.3 Plotting the Analytical Solution\n",
    "\n",
    "The following cell demonstrates use of these results to plot the transient response for a particular choice of model parameters. \n",
    "\n",
    "The steady state analysis provided an estimate for the gross heat transfer coefficient $U_a$. Rerun this cell for different values of gross heat capacity $C_p$. Try to find a value that at least mimics the experimental response shown above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.3.2.3 Plotting the Analytical Solution](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.2.3-Plotting-the-Analytical-Solution)",
     "section": "2.3.2.3 Plotting the Analytical Solution"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd0HPW5//H3oy5L7kW2cccUm2aQ6YFYBEiAEJJDCpAQbkJiEn7JJZXATe6FFFK5IclJbioQUkD0EkIChtimFxvLDeFu3G1JtmzLstru8/tjRo4wsrQrabQr7ed1zhzNfjXlY3m1j+Y7M98xd0dERDJXVqoDiIhIaqkQiIhkOBUCEZEMp0IgIpLhVAhERDKcCoGISIZTIRARyXAqBCIiGU6FQEQkw+WkOkAiRowY4ZMmTerSuvv27aOoqKhnA/UA5UqOciVHuZKTrrmge9kWLlxY7e4jO13Q3dN+Ki0t9a6aO3dul9eNknIlR7mSo1zJSddc7t3LBizwBD5jI+0aMrP1ZrbUzCrMbEHYdrOZbQ7bKszswigziIhIx3qja6jM3asParvN3W/thX2LiEgndLJYRCTDRV0IHHjKzBaa2ew27V8wsyVmdoeZDY04g4iIdMA8wucRmNlYd99iZqOAOcAXgRVANUGR+C4wxt0/3c66s4HZACUlJaXl5eVdylBXV0dxcXEX/wXRUa7kKFdylCs56ZoLupetrKxsobvP7HTBRM4o98QE3Ax87aC2ScCyztbVVUO9R7mSo1zJUa7k9emrhsysyMwGts4D5wPLzGxMm8U+BCyLKoOIiHQuyquGSoCHzax1P3e7+z/N7M9mNoOga2g9cE2EGUREEubutMSdlpjTEo/TEnOa43FiB9qcWDxOc8yDtvB1LA4trcvFnVjMiXmwTCzuxD1YP+ZOPN7ma9yJebDfWJv2uHNgvcOa45H/uyMrBO6+FjihnfYro9qniPQ98bizvzkWTE0xGltiNDTHD3xtaI7R2BK8bnzb66CtqSUeTLGgrfV1cyz4wG6KtX0dtrXEaQk/0JtjcRqbWog/9QQt8fR6hrsZfPmk/Mj30yeGmBCR9OAefGjXNbSwp6GFvQ3N1DW2sK+xhbrGWPg1eL2vsYX6phj1TTH2NQXz22v2k71wHvub2n7wd/0v3iyD/Jxs8nKygik7i/zW+ZwscrODtkGFueRlG7nZWeRkZ5GbbeRlZ5ETtm3dvJkpkyYE38sycrKzyMkycrLtQFt26+us4HvZWcG6WVl24HVOlr3tdXaWkW3/ns9qb94My+LAclkHvoKZMW/evJ77DzwEFQKRDNTQHKO2vpld9U3sqm9id30zu8LXe/Y3szuc9jSEX/e3sKehmb0NLcQS+Ks5J8soys+hKC+bwrxsivJzKMzNZki+MW7MIApzsxkQfq8wN5zysinIySY/N4uC3OxgysmiMC+b/JxsCnKzyM/JJj8ni/zcrPCDvGdOc86bV8WsWUf3yLb6IhUCkX7A3dnT0ELV3kaq9jby8tYW1jy/jqq9jdTUNbJzXxM1+5qCr3WN7GuKHXJb+TlZDC7MZVBhLoMLcxk1sICpI3MYVJjLwIIcBhbkUpyfE87nUJwfvC7Oz6EoP/jQz8/JIjw/+Dbz5s1j1qyTovxRSBeoEIikuXjc2bG3kc21+9m6ez/bdjcE054Gtu9p/dpI08FdLIvfIDfbGF6Uz7CiPIYX5zFx+IBgviiPYUX5DBmQy5ABuQwdkHfga0Fudmr+oZIyKgQiKebuVO1t5K2d9WyoqeetnfVs2lXPltr9bK4NPvibY2/vjinIzWL0oAJKBhVw0oShlAwqYNTAfEaG0/rKJVx0zlkMKsxp9y9zkbZUCER6yc59TaytqmNt1T7WVAdf36rZx4ad9TS0uUTQDEYPKuCwIYWcNGEoY4cUclg4jR1SyOhBBZ1+wDdtzGLwgNze+GdJP6BCINLDauoaWbF9Lyu37WXF9jpWbt/Lmqo6auubDyyTl53FxOEDmDSiiLOPGMmE4QOYMCyYDhtaSH6Oumek96gQiHRRLO5sqYvzyKLNLNu8mze27mHl9r1U1zUdWGbIgFyOLBnIhceNYcqIIg4fWcyUkUWMGzqA7Cx12Uh6UCEQSYC7s76mnkUbdrF4Yy3LtuzhjS172N8cAyrIy8li2uiBvOfoEo4cPZCjSgZy5OhiRhbnq49e0p4KgUg76hpbqNhQy6INu1i0Mfi6K+zaKcrL5pixg/nYyePJ2buFD7/nVA4fWUxuD13TLtLbVAhEgL0NzSxYv4uX19bw8rqdLNu8+8CNU0eMKua86SWcNGEoJ04YytRRxQe6debNq+Lo0YNSGV2k21QIJCM1tsRYuH4X81dW8dLaGpZt3k3cITfbmDF+CNfOOpyZk4YxY/wQBhfq6hvp31QIJGNsqKln/sodzF9ZxYtraqhvipGbbZw4fihfKJvKaVOGc+KEoRTm6YodySwqBNJvxePOoo21PLV8G0+9sZ111fsAGD+skEtPGse7jxzJ6YcPpyhfvwaS2fQbIP1KcyzOy2treHL5Np5avp0dexvJzTZOmzKcq06fyLuPGsWk4QN0JY9IGyoE0ufF484r63byaMVmnli6lT0NLRTmZjPrqJG895jRlB09Sv38Ih1QIZA+a8W2vTy8aDOPVWxmy+4GivKyOf+Y0Vxw7GjOPnKkBk8TSZAKgfQpu+ubefD1Tdy/cBOVW/eQnWWcfcQIvnHB0Zw3vYQBeXpLiyRLvzWS9tyd1bUx/nbfYh5fsoXGljgnjBvMtz9wDBcdP4YRxdE/yk+kP1MhkLRV19jCI4s289dXNlC5tYGivK18uHQcV5w6gWPGDk51PJF+Q4VA0s72PQ3c8cI67n55A3sbW5g2ZhBXTc/j6x8ro1iXeor0uEh/q8xsPbAXiAEt7j7TzIYB9wKTgPXAR919V5Q5pG9YtX0vv3t2LY9UbCYWdy44dgxXnzWZE8cPYf78+SoCIhHpjd+sMnevbvP6BuAZd/+hmd0Qvv5GL+SQNPXa+p38Zt4annlzBwW5WVx+ygSuftdkJg4vSnU0kYyQij+xLgFmhfN3AfNQIchIFRtr+d+nVvDcqmqGDsjlS+cewSdPn8SworxURxPJKFEXAgeeMjMHfuvuvwNK3H0rgLtvNbNREWeQNFO5dQ//+9RKnq7czrCiPL554TQ+cdpEjfEjkiLm7p0v1dWNm4119y3hh/0c4IvAY+4+pM0yu9x9aDvrzgZmA5SUlJSWl5d3KUNdXR3FxcVdWjdKmZhr2744D69q4tVtMQpy4ILJuZw3MZfCnM6He8jEn1d3KFdy0jUXdC9bWVnZQnef2emC7t4rE3Az8DVgBTAmbBsDrOhs3dLSUu+quXPndnndKGVSrtp9TX7To8t8yo1/92n//Q//8T8rfde+xpTn6gnKlRzlSl53sgELPIHP58i6hsysCMhy973h/PnAd4DHgKuAH4ZfH40qg6RWLO7ct2AjP3lyBbX1TVxx6gS+dO6RugFMJM1EeY6gBHg4HOUxB7jb3f9pZq8B95nZ1cAG4CMRZpAUWfjWTm56bDnLNu/hlMnDuPniY5g+Vk/yEklHkRUCd18LnNBOew3wnqj2K6m1Y28DP3jiTR5etJnRgwr4xeUncvHxYzTss0ga0x060iPcnUcrtnDz35ZT3xjjC2VTubbscA0CJ9IH6LdUum37nga++fBSnq7cwUkThvDjD5/A1FHpeQWGiLyTCoF0mbvz4Oub+c7fltPYEudbF03jU2dOJjtL3UAifYkKgXTJtt0N3PjQEuauqOLkSUP50aXHM2WkjgJE+iIVAkna3Dd38JX7KtjfHOOmi6dz1emTyNJRgEifpUIgCWuOxbn1qRX8dv5ajh49kF99/CQO11GASJ+nQiAJ2VK7ny/es4iFb+3iilMn8D/vn65nAov0EyoE0qlnKrfz1fsX09wS5+eXzeCSGYelOpKI9CAVAjmkWNz5yZMr+M38NUwbM4hfXXGiTgiL9EMqBNKufY0tXFdewdOV27n8lAncdLG6gkT6KxUCeYetu/dz9R8X8Oa2Pdx88XT+48zJqY4kIhFSIZC3Wbc7xvW/fIH6phi3/8fJlB2l5waJ9HcqBHLAP5Zu5QevNDByUCEPfv5Ujho9MNWRRKQXqBAI7s6v56/hx/9cweGDsyj/f2cycqCeGSCSKVQIMpy78/0nKvn9c+u4+ISxXDyqVkVAJMNkpTqApE487nzrkWX8/rl1XHX6RH7+sRnkZWuoCJFMoyOCDNUSi3P9g0t46PXNXPPuKdzwvqP18BiRDKVCkIGaWuJ8+d4K/r50K18570i+eM5UFQGRDKZCkGEammNc+9fX+debO/jWRdP4zFlTUh1JRFJMhSCD7G+K8Zk/vcYLq2v43geP5ROnTUx1JBFJAyoEGaKpJc7n/rKQl9bU8L8fOYFLS8elOpKIpInIrxoys2wzW2Rmj4ev/2hm68ysIpxmRJ0h08XizlfvX8z8lVXc8qHjVARE5G1644jgOqASGNSm7evu/kAv7DvjuTs3P7acvy3ewjfedzSXnzIh1ZFEJM1EekRgZuOAi4A/RLkfObTb5qzkzy+/xTVnT+Hzsw5PdRwRSUPm7tFt3OwB4AfAQOBr7v5+M/sjcDrQCDwD3ODuje2sOxuYDVBSUlJaXl7epQx1dXUUF6ffGPq9keup9c3c/WYTZ4/L4VPH5CV0iWgm/7y6QrmSo1zJ6062srKyhe4+s9MF3T2SCXg/8H/h/Czg8XB+DGBAPnAX8D+dbau0tNS7au7cuV1eN0pR53pgwUaf+I3H/Zo/LfDmlljC62Xqz6urlCs5ypW87mQDFngCn9dRdg2dCXzAzNYD5cA5ZvYXd98aZmwE7gROiTBDRpr75g6uf3AJZ04dzs8vn0FOtkYSEZFDO+QnhJn9wMyuaaf9y2Z2S2cbdvcb3X2cu08CLgP+5e6fMLMx4XYM+CCwrMvp5R1W79jLF+9ZxNGjB/LbK2eSn6OniolIxzq6augS4Nh22n8BLAa+2cV9/tXMRhJ0D1UAn+viduQgtfVNfOauBRTkZvP7T86kOF+3iYhI5zr6pIi7e/zgRnePWZID07j7PGBeOH9OMutKYlpicb5w9yI21+6nfPZpjB1SmOpIItJHdNR53GBm77jeMGxriC6SdMX3n3iT51dXc8sHj6N04rBUxxGRPqSjI4KbgCfM7LvAwrBtJkGX0NeiDiaJu++1jdzxwjo+deYkPnry+FTHEZE+5pCFwN3/bmabgOuBr4fNy4DL3L2iN8JJ5xa+tZNvPrKUd00dwTcvnJbqOCLSB3V4NtHdFwMf76UskqQttfu55s+vM3ZIIb+84kRdJioiXXLIQmBmo4Avhi9/5e7beieSJKI5Fufav75OQ3OMez57KkMG5KU6koj0UR39CfkX4AHg/nBe0shtc1ZSsbGWH156HEeUDEx1HBHpwzrqGmoERgEOtPROHEnEi6ur+fX8NXxs5njef/zYVMcRkT6uoyOCK4AjgWkEdwZLGti5r4kv3VvB5BFF3PSB6amOIyL9QEdXDe0FftWLWaQT7s71Dyymtr6ZOz91MgPydOewiHSfLjPpQ/700ls8XbmDGy44mmPGDk51HBHpJ1QI+ojKrXu45YlKyo4ayafOnJTqOCLSj3Q0+ujJvRlEDm1/U4z/vGcRgwtz+clHTkjoATMiIonq6IjgmvCh838xs0+EI4ZKCtzyxBus2lHHTz96AiOK81MdR0T6mY5OFn8GwMyOBS4A7jGzAuBfwD+Bl9sbnVR61ourq/nLyxv47FmTOesI1WIR6XmdniNw92Xu/hN3Pxd4L7AAuDL8KhFqaI5x48NLmTR8AF89/6hUxxGRfiqp6w/dfR/wWDhJxH729Creqqnn7s+eSkGunjQmItHQVUNpatnm3fz+ubVcdvJ4zjh8RKrjiEg/pkKQhlpicW54aAnDivK48QINLS0i0UqoEJjZODMrC+fzzawo2liZ7fbn17Fs8x6+84FjGDwgN9VxRKSf67QQmNmnCc4J/CFsmgg8GmWoTLa+eh8/nbOS86eX8L5jR6c6johkgESOCP4TOA3YA+DuKwlGJU2ImWWH9yM8Hr6ebGavmNkqM7vXzDSQfsjd+a+Hl5KXncV3LjlWN46JSK9IpBA0uHtT6wszywaS+YS6Dqhs8/pHwG3ufgSwC7g6iW31a/cv3MSLa2q44cKjGT24INVxRCRDJFIIXjCz64GC8DzBvcDjiWzczMYBFxF2K1nwJ+45BA+8AbgL+GCyofujmrpGbvl7JadMGsblJ09IdRwRySCJFILrgb3AmwR/3T8DfDPB7f8sXL/1DuThQK27tz7oZhNwWMJp+7Hbnl5JXWMLt3zoWLKy1CUkIr3H3P3Q3wy6ge5w96uS3rDZ+4EL3f1aM5sFfA34FPCSu08NlxkPPOHux7Wz/mxgNkBJSUlpeXl5shEAqKuro7i4uEvrRqltrs1743zrhf2cMyGHK6endiyhvvDzSifKlRzlSl53spWVlS1095mdLujuHU7AU0BuZ8u1s94PCP7iXw9sA+qBvwLVQE64zOnAk51tq7S01Ltq7ty5XV43Sq254vG4f+IPL/txN/3Td9Y1pjaUp//PK90oV3KUK3ndyQYs8AQ+rxMZYmIt8JyZPQrsa1NAftFJgbkRuBGg9YjA3T9uZvcDHwbKgavI8EtR562o4rlV1XzromkMLdIFVCLS+xI5R1AFzAEGACPbTF31DeArZraa4JzB7d3YVp/WHIvzvb+/weQRRXzy9EmpjiMiGarTIwJ3/+/u7sTd5wHzwvm1wCnd3WZ/cPcrG1hTtY/ff3ImeTka7UNEUqPTQmBmc4B3nFF29/MjSZQh9jU7tz27kjMOH8650xK+P09EpMclco7gW23mC4BLgcZo4mSOR1c3sXt/C9+6aLruIBaRlEqka+iVg5rmm9n8iPJkhLVVdTyzoYWPzRzP9LGDUh1HRDJcIl1DbT+psoBSYExkiTLA9594k9ws9NQxEUkLiXQNLSc4R2BAC7AO+GyUofqzl9bU8HTldj58ZC4jB+pB9CKSeokUginu3ty2wcySesSlBNyd2+aspGRQPudP1KMnRSQ9JHLN4sHnCABe7ekgmeClNTW8un4n186aSl62ThCLSHo45F/2ZjaK4FxAoZkdx7+Hnh5EcHOZJMHd+dnTqxg9qICPnTyel19Yn+pIIiJAx11DFwGfBsYB/9emfS/Q7ZvMMs2L4dHAdy45hoJcdQuJSPo4ZCFw9zuBO83so+5+Xy9m6neCo4GVjB5UwEdnjk91HBGRt0nkPoL7zOy9wDEEN5S1tn8/ymD9yQura3ht/S6+q6MBEUlDidxH8H/AEOBs4E6CO4tfjjhXv9F6NDBmcAEfPVlHAyKSfhK5auhd7n4FUBMOQHcqwXkDScALq2tY8NYuri2bSn6OjgZEJP0k9PD61q9mNjp8PSmyRP2Iu3Nb69HATNVOEUlPiRSCJ8xsCHArUEHwxLEHOlxDAHh+dTULdTQgImmuw3MEZpYF/MPda4H7zexxoNDdd/ZKuj6s9S7isToaEJE01+ERgbvHgZ+3eb1fRSAxz62q5vUNtToaEJG0l0jX0BwzuyTyJP3Mb+avCc8N6EohEUlviQwe9wVgsJk1AvsJhppwdx8WabI+bPmW3by4poYbLjhaj6AUkbSXSCEYEXmKfuaO59czIC+by0+ekOooIiKd6vTPVXePAR8BvhHOjwFmRB2sr9qxp4HHFm/mozPHM3hAbqrjiIh0qtNCYGa/BMqAK8OmeuA3CaxXYGavmtliM1tuZt8O2/9oZuvMrCKc+lVR+dNLb9ESdz515qRURxERSUgiXUNnuPtJZrYIwN13mlleAus1Aue4e52Z5QLPm9k/wu993d373b0I+5ti/PWVtzhvWgkThxelOo6ISEISKQTN4f0EDmBmw4F4Zyu5uwN14cvccPIu5uwTHlq0iV31zVz9rsmpjiIikrBELmn5FfAgMDLs3nke+FEiGzezbDOrAHYAc9y99Wlnt5jZEjO7zcz6xYN743Hn9ufXcdxhgzllsi6oEpG+w4I/3DtZyOwY4Nzw5TPuviypnQRDVDwMfBGoAbYBecDvgDXu/p121pkNzAYoKSkpLS8vT2aXB9TV1VFcXNyldZOxuKqF2xY2Mvv4fM4Y2/mBVm/lSpZyJUe5kqNcyetOtrKysoXuPrPTBd290wk4HrgW+DxwfCLrtLONm4CvHdQ2C3i8s3VLS0u9q+bOndvldZNxxe9f8lNvedobm2MJLd9buZKlXMlRruQoV/K6kw1Y4Al8Pidy1dA3gXuAsQTDT99tZjcmsN7I8EgAMyskOKJ408zGhG0GfBBI6ugiHb2xZQ8vrK7hqjMm6QYyEelzEjlZ/Amg1N3rAczsFmAh8INO1hsD3GVm2QTnIu5z98fN7F9mNpLgDuUK4HNdTp8m7nhhHYW52Vxxim4gE5G+J5FC8NZBy+UAaztbyd2XACe2035Owun6gB17Gni0YjOXnzJBN5CJSJ+USCGoB5ab2ZMEl3+eT3BPwE8B3P0rEeZLe39+ufUGMl0yKiJ9UyKF4O/h1ErPKw41x+Lc8+pGyo4axeQRuoFMRPqmTguBu9/eG0H6omcqd1Bd16hzAyLSpyVy1dD7zOw1M9thZjvNbJeZ6eE0wL2vbaBkUD6zjhqZ6igiIl2WSNfQL4GPAktJYGiJTLGldj/zV1Zx7ayp5GTrklER6bsSKQSbgAoPHlspofsXbCLu8LGT9QQyEenbEikE1wN/M7N5BCOKAuDuv4gqVLqLxZ37FmzkXVNHMH7YgFTHERHplkT6NL4NxIAhwMg2U8Z6blUVm2v3c9kpOhoQkb4vkSOCUe5eGnmSPuTe1zYydEAu500vSXUUEZFuS+SI4Bkz61d3A3dH1d5G5ryxnUtPGkd+Tnaq44iIdFsiheCzwNNmVqfLR+Gh1zfREnd1C4lIv5FI19CIyFP0Ee7Ova9tZObEoUwdNTDVcUREekSnRwTuHgM+AnwjnB8D9KsHzifq1XU7WVu9T5eMiki/ksidxb8EyoArw6Z64DdRhkpX5a9tZGB+DhcdPybVUUREekwi5wjOcPdrgAYAd99J8JjJjLK7vpknlm7lAzPGMiAvkR41EZG+IZFC0GxmWQRDUGNmw8nAoSYeqdhMY0ucyzXAnIj0M4csBGbW+mfvr4AHgZFm9m3geeBHvZAtrZS/tpFjxg7i2MMGpzqKiEiP6qiP41XgJHf/k5ktJHjmsAEfcfc+/5zhZKzYtpfKrXu4+eLpqY4iItLjOioE1jrj7suB5dHHSU+PLd5MlsFFx49NdRQRkR7XUSEYaWaHfAylu/80gjxpx9352+KtnDl1BCMH5qc6johIj+uoEGQDxbQ5MshEFRtr2bCzni+cMzXVUUREItFRIdjq7t/p6obNrAB4FsgP9/OAu99kZpOBcmAY8Dpwpbs3dXU/UXts8RbysrN47zGjUx1FRCQSHV0+2t0jgUbgHHc/geBO5PeZ2WkEVxzd5u5HALuAq7u5n8jE4s7jS7Yy66iRDC7MTXUcEZFIdFQI3tOdDXugLnyZG04OnAM8ELbfBXywO/uJ0ivraqja28gHZugksYj0X+bu0W3cLBtYCEwluB/hJ8DL7j41/P544B/ufmw7684GZgOUlJSUlpeXdylDXV0dxcXFXVr3zmWNvLK1hZ+fM4D87J49VdKdXFFSruQoV3KUK3ndyVZWVrbQ3Wd2uqC7Rz4RPN1sLnAWsLpN+3hgaWfrl5aWelfNnTu3S+s1Nsf8+Juf9Ovueb3L++5IV3NFTbmSo1zJUa7kdScbsMAT+IxOZIiJbnP3WmAecBowpM1dy+OALb2RIVnPrapi9/5mdQuJSL8XWSEws5FmNiScLyS4M7mS4Mjgw+FiVwGPRpWhOx5bvIXBhbm8a2pGP55ZRDJAlMNojgHuCs8TZAH3ufvjZvYGUG5m3wMWAbdHmKFL9jfFmPPGdi6ZMZa8nF45aBIRSZnICoG7LwFObKd9LXBKVPvtCU9Xbqe+KcbFJ6hbSET6P/25247HFm9h1MB8Tp08PNVRREQip0JwkN37m5m/oor3Hz+W7KyMHl1DRDKECsFBnly2jaZYXFcLiUjGUCE4yGOLtzBx+ABOGKcH0IhIZlAhaKO6rpEX11Rz8fFjMVO3kIhkBhWCNp6p3E7c4cLjxqQ6iohIr1EhaGPOG9s5bEgh08YMTHUUEZFeo0IQqm9q4blV1Zw3vUTdQiKSUVQIQs+tqqaxJc5500tSHUVEpFepEISefmM7gwpyOGXysFRHERHpVSoEBE8i+9ebOyg7ehS52fqRiEhm0ace8PqGXdTsa+LcaeoWEpHMo0JA0C2Um23MOkpDTotI5lEhILhs9LQpwxlYoAfUi0jmyfhCsHpHHWur9+lqIRHJWBlfCJ6u3A6g8wMikrEyvhDMeWM7xx42iLFDClMdRUQkJTK6EFTtbeT1Dbs4b9roVEcREUmZjC4Ec9/cgTucO31UqqOIiKRMRheCp8JB5qaPGZTqKCIiKRNZITCz8WY218wqzWy5mV0Xtt9sZpvNrCKcLowqQ0f2N8V4fnWVBpkTkYyXE+G2W4CvuvvrZjYQWGhmc8Lv3ebut0a47049t6qKhua4rhYSkYwXWSFw963A1nB+r5lVAodFtb9kPV25nYEFOZw6RYPMiUhm65VzBGY2CTgReCVs+oKZLTGzO8xsaG9kaCsWd56p3EHZURpkTkTE3D3aHZgVA/OBW9z9ITMrAaoBB74LjHH3T7ez3mxgNkBJSUlpeXl5l/ZfV1dHcXHx29pW7YpxyysNfO6EfE4bE2XvWHK50oFyJUe5kqNcyetOtrKysoXuPrPTBd09sgnIBZ4EvnKI708ClnW2ndLSUu+quXPnvqPt1iff9Ck3/t1r65u6vN3uai9XOlCu5ChXcpQred3JBizwBD6ro7xqyIDbgUp3/2mb9rZPhv8QsCyqDIfy7KpqZowfwuBCDTInIhJlv8iZwJXAUjOrCNv+C7jczGYQdA2tB66JMMM71NY3sWRTLde954je3K2ISNqK8qqh54HtrazaAAAKeklEQVT2LtB/Iqp9JuLFNTW4w1lHjEhlDBGRtJFxl8w8t6qKgfk5nDBuSKqjiIikhYwqBO7OsyurOf3w4eToslERESDDCsH6mno21+7nrCP1SEoRkVYZVQieW1UFwNk6PyAickBGFYJnV1YzflghE4cXpTqKiEjayJhC0ByL8/LaGs46Qt1CIiJtZUwhqNhYS11ji7qFREQOkjGF4LmVVWQZnH64CoGISFuZUwhWV3OChpUQEXmHjCgEu+ubWbyxVucHRETakRGF4MU11cQ1rISISLsyohA8t7qa4vwcZozXsBIiIgfr94UgGFaiitOmDNfTyERE2tHvPxl31Dubdu3n7CPVLSQi0p5+XwiW18QAdKJYROQQ+n0hWFYd47AhhUwaPiDVUURE0lK/LgQtsTiVO2OcfeQIgidniojIwfp1IVi8qZb9LeoWEhHpSL8uBM+urMaAMw4fnuooIiJpq18XgrFDCnjXYTkMGZCX6igiImkrsofXp4OPnTyBkn1rUx1DRCStRXZEYGbjzWyumVWa2XIzuy5sH2Zmc8xsVfh1aFQZRESkc1F2DbUAX3X3acBpwP8zs+nADcAz7n4E8Ez4WkREUiSyQuDuW9399XB+L1AJHAZcAtwVLnYX8MGoMoiISOd65WSxmU0CTgReAUrcfSsExQIY1RsZRESkfebu0e7ArBiYD9zi7g+ZWa27D2nz/V3u/o7zBGY2G5gNUFJSUlpeXt6l/dfV1VFcXNy18BFSruQoV3KUKznpmgu6l62srGyhu8/sdEF3j2wCcoEnga+0aVsBjAnnxwArOttOaWmpd9XcuXO7vG6UlCs5ypUc5UpOuuZy7142YIEn8Fkd5VVDBtwOVLr7T9t86zHgqnD+KuDRqDKIiEjnoryP4EzgSmCpmVWEbf8F/BC4z8yuBjYAH4kwg4iIdCLycwQ9wcyqgLe6uPoIoLoH4/QU5UqOciVHuZKTrrmge9kmunung631iULQHWa2wBM5WdLLlCs5ypUc5UpOuuaC3snWr8caEhGRzqkQiIhkuEwoBL9LdYBDUK7kKFdylCs56ZoLeiFbvz9HICIiHcuEIwIREelAvy4EZvY+M1thZqvNrFdHOTWzO8xsh5kta9PW7hDcFvhFmHOJmZ0UYa6khgfvrWxmVmBmr5rZ4jDXt8P2yWb2SpjrXjPLC9vzw9erw+9PiiJXuK9sM1tkZo+nS6Zwf+vNbKmZVZjZgrAtHd5jQ8zsATN7M3yfnZ7qXGZ2VPhzap32mNmXUp0r3NeXw/f8MjO7J/xd6N33WCK3H/fFCcgG1gBTgDxgMTC9F/d/NnASsKxN24+BG8L5G4AfhfMXAv8AjGDI7lcizDUGOCmcHwisBKanOlu4/eJwPpdggMLTgPuAy8L23wCfD+evBX4Tzl8G3Bvhz+wrwN3A4+HrlGcK97EeGHFQWzq8x+4CPhPO5wFD0iFXm3zZwDZgYqpzEYzIvA4obPPe+o/efo9F+gNP5QScDjzZ5vWNwI29nGESby8E7Y6zBPwWuLy95Xoh46PAeemUDRgAvA6cSnAjTc7B/6cEY1idHs7nhMtZBFnGETw34xzg8fCDIaWZ2mRbzzsLQUr/H4FB4QebpVOug7KcD7yQDrkICsFGYFj4nnkceG9vv8f6c9dQ6w+41aawLZUONQR3SrJaYsOD91q2sAumAtgBzCE4oqt195Z29n0gV/j93cDwCGL9DLgeiIevh6dBplYOPGVmCy0YrRdS//84BagC7gy70/5gZkVpkKuty4B7wvmU5nL3zcCtBMPtbCV4zyykl99j/bkQWDtt6XqJVK9ntWB48AeBL7n7no4WbactkmzuHnP3GQR/hZ8CTOtg35HnMrP3AzvcfWHb5lRmOsiZ7n4ScAHBEwDP7mDZ3sqWQ9Al+mt3PxHYR8dPIezVn1nY1/4B4P7OFm2nrcdzheckLgEmA2OBIoL/z0PtO5Jc/bkQbALGt3k9DtiSoiyttpvZGIDw646wvVezmlkuQRH4q7s/lE7ZANy9FphH0Dc7xMxaB0dsu+8DucLvDwZ29nCUM4EPmNl6oJyge+hnKc50gLtvCb/uAB4mKJ6p/n/cBGxy91fC1w8QFIZU52p1AfC6u28PX6c617nAOnevcvdm4CHgDHr5PdafC8FrwBHh2fc8gsPBx1Kc6VBDcD8GfDK8UuE0YHfr4WpPM0t6ePBeyWZmI81sSDhfSPALUgnMBT58iFyteT8M/MvDjtOe4u43uvs4d59E8P75l7t/PJWZWplZkZkNbJ0n6PdeRor/H919G7DRzI4Km94DvJHqXG1czr+7hVr3n8pcG4DTzGxA+LvZ+vPq3fdYlCdlUj0RnPlfSdDX/M1e3vc9BH1+zQRV/GqCvrxngFXh12Hhsgb8Ksy5FJgZYa53ERxKLgEqwunCVGcDjgcWhbmWAf8Ttk8BXgVWExzO54ftBeHr1eH3p0T8/zmLf181lPJMYYbF4bS89f2d6v/HcF8zgAXh/+UjwNA0yTUAqAEGt2lLh1zfBt4M3/d/BvJ7+z2mO4tFRDJcf+4aEhGRBKgQiIhkOBUCEZEMp0IgIpLhVAhERDJcTueLiKQHM2u91A9gNBAjGM4AoN7dz4hovwsJxndpimL7XWVmde5enOoc0vepEEif4e41BNeoY2Y3A3XufmuU+wzHY9qcbkVApCepa0j6BTOrC7/OMrP5Znafma00sx+a2ccteNbBUjM7PFxupJk9aGavhdOZh9j0BcA/29nfD83sjXCs+ls72qaZFZvZneH+l5jZpWH75WHbMjP7Udt/i5ndYsGzGV42s5KwfbKZvRRu+7ttlh9jZs9aMM7+MjM7q2d+qpIpVAikPzoBuA44DrgSONLdTwH+AHwxXObnwG3ufjJwafi99ryPgwqBmQ0DPgQc4+7HA9/rZJv/TTBEwXHh8v8ys7HAjwjGL5oBnGxmHwyXLwJedvcTgGeBz7bZ/q/D7W9rE+kKgmGKZ4T/9ooEfkYiB6hrSPqj1zwcF8bM1gBPhe1LgbJw/lxgejC8CwCDzGygu+9tbQjHqBrn7msP2v4eoAH4g5n9nWAM+UNuM2y/rLXR3XeFI4XOc/eqcF9/JXiY0SNAU5ttLiR4XgQEg+BdGs7/maCQQDCu1h3hYIKPuLsKgSRFhUD6o8Y28/E2r+P8+z2fRXACeH8H2zkLeP7gRndvMbNTCAYIuwz4AsFf9u1uMxxM7OCxXNobTrhVs/977JcYb/89fceYMO7+bFhYLgL+bGY/cfc/dbB9kbdR15BkqqcIPsABMLMZ7SzzPoLHFb6NBc9yGOzuTwBfIjyB3cE2D24fSvAwoHeb2QgzyyYYFXN+J5lf4N9HFh9vs72JBM9N+D3ByLKRPV9X+icVAslU/wnMDE/evgF8rp1lZtH+h/NA4HEzWxJ+/8udbPN7wNDwRO5ioCzsurqRYLjhxQRj5D9Kx64jeADNawTj0LfNWWFmiwi6jn7eyXZE3kajj4q0w8zGAb939/aeFiXSr6gQiIhkOHUNiYhkOBUCEZEMp0IgIpLhVAhERDKcCoGISIZTIRARyXAqBCIiGe7/AxgJ3r3e89ATAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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",
    "\n",
    "# parameter values and units\n",
    "P1 = 4                 # watts\n",
    "Ua = 0.06              # watts/deg C\n",
    "Cp = 8                 # joules/deg C\n",
    "u1_steadystate = 0.5   # steady state value of u1 (fraction of total power)\n",
    "T_ambient = 21         # ambient temperature\n",
    "\n",
    "# initial conditions\n",
    "T1_deviation_initial = 0\n",
    "\n",
    "# steady state solution\n",
    "T1_deviation_steadystate = P1*u1_steadystate/Ua\n",
    "\n",
    "# compute the transient solution\n",
    "t = np.linspace(0,800)\n",
    "T1_deviation = T1_deviation_steadystate \\\n",
    "     + (T1_deviation_initial - T1_deviation_steadystate)*np.exp(-Ua*t/Cp)\n",
    "\n",
    "# plot\n",
    "plt.plot(t, T1_deviation + T_ambient)\n",
    "plt.xlabel('Time / seconds')\n",
    "plt.ylabel('Temperature / °C')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.3 Matching the Model to Experimental Data](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.3-Matching-the-Model-to-Experimental-Data)",
     "section": "2.3.3 Matching the Model to Experimental Data"
    }
   },
   "source": [
    "## 2.3.3 Matching the Model to Experimental Data\n",
    "\n",
    "The following cell provides an interactive tool for 'tuning' the model to fit the experimental data. Work with the sliders to find good choices for each of the parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.3 Matching the Model to Experimental Data](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.3-Matching-the-Model-to-Experimental-Data)",
     "section": "2.3.3 Matching the Model to Experimental Data"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0dee00c8602e4c47bf9b9461170fc321",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.05, description='Ua', max=0.1, step=0.001), FloatSlider(value=6.5, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.compare(Ua, Cp)>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "\n",
    "# known parameter values\n",
    "P1 = 4\n",
    "u1_steadystate = 0.5   # steady state value of u1 (fraction of total power)\n",
    "T_ambient = 21         # ambient temperature\n",
    "\n",
    "def compare(Ua, Cp):\n",
    "    T1_deviation_initial = 0\n",
    "    T1_deviation_steadystate = P1*u1_steadystate/Ua\n",
    "    T1_deviation = T1_deviation_steadystate \\\n",
    "         + (T1_deviation_initial - T1_deviation_steadystate)*np.exp(-Ua*t/Cp)\n",
    "    T1_model = T1_deviation + T_ambient\n",
    "    plt.plot(t, T1, t, T1_model)\n",
    "    plt.xlabel('time / seconds')\n",
    "    plt.ylabel('temperture / °C')\n",
    "    plt.grid()\n",
    "    plt.text(200,35,'Sum of errors = ' + str(round(sum(abs(T1_model-T1),2))))\n",
    "    plt.text(200,30,'Ua = ' + str(Ua))\n",
    "    plt.text(200,25,'Cp = ' + str(Cp))\n",
    "\n",
    "from ipywidgets import interact\n",
    "interact(compare, Ua=(0.0,0.10,0.001), Cp=(2.0,11.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[2.3.3 Matching the Model to Experimental Data](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.3-Matching-the-Model-to-Experimental-Data)",
     "section": "2.3.3 Matching the Model to Experimental Data"
    }
   },
   "source": [
    "**Exercise:** Determine values for $U_a$ and $C_p$. \n",
    "\n",
    "**Exercise:** The sum of absolute errors is shown on the chart. Try to find values of $U_a$ and $C_p$ that minimize this error criterion. In your opinion, is that the best choice of model parameters? Why or why not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.3.3.1 Does this solution make sense?](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.3.1-Does-this-solution-make-sense?)",
     "section": "2.3.3.1 Does this solution make sense?"
    }
   },
   "source": [
    "### 2.3.3.1 Does this solution make sense?\n",
    "\n",
    "The parameter values in the above plot were chosen to (at least roughly) reproduce the measured response of temperature control laboratory. The value used for the heat capacity was $C_p = 7$ watts/degC. Is this reasonable?\n",
    "\n",
    "The [specific heat capacity for solids](https://en.wikipedia.org/wiki/Heat_capacity) is typically has values in the range of 0.2 to 0.9 watts/degC/gram. Using a value of 0.9 that is typical of aluminum and plastics used for electronic products, the estimated mass of the heater/sensor pair would be\n",
    "\n",
    "$$ m \\approx \\frac{7 \\mbox{ watts/degC}}{0.9 \\mbox{ watts/degC/grams}} \\approx 8 \\mbox{ grams}$$\n",
    "\n",
    "Does that seem reasonable?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.3.3.1 Does this solution make sense?](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.3.1-Does-this-solution-make-sense?)",
     "section": "2.3.3.1 Does this solution make sense?"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.3.3.1 Does this solution make sense?](https://jckantor.github.io/CBE32338/02.03-First-Order-Model-for-a-Single-Heater.html#2.3.3.1-Does-this-solution-make-sense?)",
     "section": "2.3.3.1 Does this solution make sense?"
    }
   },
   "source": [
    "<!--NAVIGATION-->\n",
    "< [2.2 Fitting Step Test Data to Empirical Models](https://jckantor.github.io/CBE32338/02.02-Fitting-Step-Test-Data-to-Empirical-Models.html) | [Contents](toc.html) | [2.4 Two-Input, Two-Output Model](https://jckantor.github.io/CBE32338/02.04-Two-Input-Two-Output-Model.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE32338/blob/master/docs/02.03-First-Order-Model-for-a-Single-Heater.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.03-First-Order-Model-for-a-Single-Heater.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
}