{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NOTEBOOK_HEADER-->\n",
    "*This notebook contains course material from [CBE30338](https://jckantor.github.io/CBE30338)\n",
    "by Jeffrey Kantor (jeff at nd.edu); the content is available [on Github](https://github.com/jckantor/CBE30338.git).\n",
    "The text is released under the [CC-BY-NC-ND-4.0 license](https://creativecommons.org/licenses/by-nc-nd/4.0/legalcode),\n",
    "and code is released under the [MIT license](https://opensource.org/licenses/MIT).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Simulation and Optimal Control](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/07.00-Simulation-and-Optimal-Control.ipynb) | [Contents](toc.ipynb) | [Soft Landing a Rocket](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/07.02-Soft-Landing-a-Rocket.ipynb) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE30338/blob/master/notebooks/07.01-Simulation-and-Optimal-Control-in-Pharmacokinetics.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://raw.githubusercontent.com/jckantor/CBE30338/master/notebooks/07.01-Simulation-and-Optimal-Control-in-Pharmacokinetics.ipynb\"><img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "8C_vr6C0Jf8I"
   },
   "source": [
    "# Simulation and Optimal Control in Pharmacokinetics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "r1DkXPM5wif0"
   },
   "source": [
    "## Installations and Initialization\n",
    "\n",
    "Pyomo and ipopt are now preinstalled on [Google Colaboratory](https://colab.research.google.com/). On MacOS and Windows PC, a one-time installation of pyomo and ipopt is required. The following commands will perform the needed installation when using the Anaconda distribution of Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "fP3nFqWlidrB"
   },
   "outputs": [],
   "source": [
    "!wget -N -q \"https://ampl.com/dl/open/ipopt/ipopt-linux64.zip\"\n",
    "!unzip -o -q ipopt-linux64\n",
    "ipopt_executable = '/content/ipopt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ISnPZ93d5Zvr"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from pyomo.environ import *\n",
    "from pyomo.dae import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "GAe_MZ1QKG6j"
   },
   "source": [
    "On MacOS and Windows PC, a one-time installation of pyomo and ipopt is required. The following commands will perform the needed installation when using the Anaconda distribution of Python.\n",
    "\n",
    "    !conda install -c conda-forge pyomo\n",
    "    !conda install -c conda-forge pyomo.extras\n",
    "    !conda install -c conda-forge ipopt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "29eeNeG5HiH2"
   },
   "source": [
    "## A Pharmacokinetics Example\n",
    "\n",
    "Pharmacokinetics -- the study of the drugs and other substances administered to a living organism -- is a rich source of examples for simulation and optimization. Here we discuss a solution to a very simple problem.\n",
    "\n",
    "A patient is admitted to a clinic with a dangerous blood borne infection that needs to be treated with a specific antibiotic. The treatment requires a minimum therapeutic concentration of 20 mg/liter of blood, but, for safety reasons, should not exceed 40 mg/liter. The attending physician has established a target value of 25 mg/liter.  The antibiotic will be administered to the patient by continuous infusion over a 72 hour period using a pump capable of delivering up to 50 mg/hour.\n",
    "\n",
    "![Infusionspumpe](https://upload.wikimedia.org/wikipedia/commons/thumb/e/ed/Infusionspumpe.JPG/256px-Infusionspumpe.JPG)\n",
    "\n",
    "(photo from [Pflegewiki-User Würfel](https://commons.wikimedia.org/wiki/File:Infusionspumpe.JPG), [CC BY-SA 3.0](http://creativecommons.org/licenses/by-sa/3.0/))\n",
    "\n",
    "The research literature shows dynamics of antibiotic concentration in the body behaves according to a one-compartment pharmacokinetic model\n",
    "\n",
    "\\begin{align*}\n",
    "V\\frac{dC}{dt} & = u(t) - k C\n",
    "\\end{align*}\n",
    "\n",
    "where $V$ is compartment volume (e.g., blood is about 5 liters in an adult), and $kC$ is the rate of elimination due to all sources, including metabolism and excretion. Based on an observed half-life in the blood of 8 hours, $k$ has an estimated value of 0.625 liters/hour. Given this information, what administration policy would you recommend?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gPvUJdEuYqPr"
   },
   "source": [
    "## Simulation of Fixed Infusion Rate Strategy\n",
    "\n",
    "The first strategy we will investigate is constant infusion.  At steady-state, the pharmacokinetic model becomes\n",
    "\n",
    "\\begin{align*}\n",
    "\\bar{u} & = k \\bar{C} \\\\\n",
    "& = 0.625 \\times 25 \\\\\n",
    "& = 15.625 \\mbox{ mg/hr}\n",
    "\\end{align*}\n",
    "\n",
    "The following cell simulates the patient response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 294
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1061,
     "status": "ok",
     "timestamp": 1554231652418,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "QG8LbaTVJnkm",
    "outputId": "76b19a5c-76d3-4180-ef56-50e8d78fa7ef"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEVCAYAAAAM3jVmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcXGWV//FPdXd6y9adpDv7RoCT\nQNhB2QkKBEEHZRlmXFDBn6DgjNu4jMu46+i4oqKIisvICDIIiDM4oAQQkDWsySEha2ftJJ10dzrp\nrer3x3OrU2m6k0qna/++X696VdWte+9z6lbVqafOvfXcWCKRQERESkNZrgMQEZHsUdIXESkhSvoi\nIiVESV9EpIQo6YuIlBAlfRGRElKR6wCyzcwSwCtAD+FL7xXgWndfYWYLgJvc/dBhauvtwHvcfcFw\nrG8I7S8Elrj7mgNczoCJ7v6gmb0FeJO7XzlMMb0f+Cxwvbt/eTjWOcQ4/p+7/2QY13cT0OTun0tz\n/q8Cq939R/uZ79fAWYT30b0HGNNrgC+6+8IDWW4o9hWnmTUC/wnMTv1smdnNwEJgR8rsV7j74wOs\n/y3A14Fy4Bng3e7eamZTgR8BhwIx4LvufkO0TALwlNWsc/fXH+xzHYyZLQVOA54C3u7uD2eqrYNR\nckk/ssDdm6Dvw/dd4E25DSkjPgR8CTigpA+8hfDeeNDd7wDuGMaYLgE+5e4/HcZ1HhAzmwR8DBi2\npH+g3P2Tac76j8Dh7v7KENp4nJBUs2HAOM1sHLAI+B9g9gDLfdLdb97Xis1sNvBD4AxCJ+3bwBuB\n3wA/Bp509zeZ2RTgRTP7s7s7gLvPPahnlSYzmwHsdPetoc+Uv0o16af6M/B3/SeaWTXwHeBsIA78\nEfiYu/ea2dHADcB4YDfwcXe/18zKgO9F69tIeLMPyMw+DlxN+MXxB+Aj7p4ws38CriH8CnFCz6k5\n6hWtBk4FDgdeBi5y9w4zOwG4ERgNbADeBVwJvB6YZ2YfA+YBU4FjCB+W7wHXA+cAlcDD0TLnA58E\nusysHnie0Gs5J/oA/yhaRy/wC3f/9+j5JIArgA8Dk4Cvu/u3+z3nrwOnRDFNB762j228CvgZ8Dbg\n3NRfK9F2flXs7t7dr71RwK+AuUAVcD/wfuARYFrUMzua0Eu8AZgMdBJ6kU9G6/gM8HbCZ2VJtC22\nm9l44BbgMOAloANoMrNvADXufl20fD2wDpjh7ltSYrsZWO7uX4qe61eBq4DpwG/c/SNm9gDhfXBv\n9L74ISk9yGi5twOPRa/LGYSe8HOE98DxRL9c9/N+HrB9+okS20+AWUA34TX+Zf843f2PKYslgDdH\n2/ZVn7M0vR243d2XR/c/mPLYjwmvP+6+3sxWEt7rTpqi5/9N4N2Ez8j7CJ+d84Fm4A3u3hL9cr4J\naCd88fwHcLS7rwLOJby/kk40s/8AZgD/5e4fNrNZhPfeb4Hj3f2sdGMcTiVd0zezSsIb6q4BHv4g\n4QNwJOHDcwbwj1HC+S/g+1Ev4j3ALWY2mvAmOQ84gvBT98xB2j09Wu4YYD5wOnCpmZ0M/Avhl8hc\nQg/9qymLXgZcDswBGgg9cqJ4Pu3uhxN65d93988Qks3b3P230XwXABe4+3eiZc+I2p8HnABc7u53\nR+v47gAf/K8ALe5uUczvj55L0pHufhzhw/0VMytPXdjdPwY8Tkg2nxtsG6csMs3dbYDy1ICx82rv\nBLa7+zzCF2VP1NaVwJpoG/cAvwd+GW2/a4A7zawi+jK9DjiJkNyrovsAHwea3X02cC17etS3AJeZ\nWbJD9UbCL6a+hD+IMwlfiCcAHzCzaSllwQX9Eml/Cwm96LlRnC9G60q1v239qvYHaOdG4IHo9b8Q\n+J6ZzdpXnO7ekux1D+KtZvaEmb1kZv9qZrEB5jmG0An5PzN72cx+ZGa10frvdvcW6PtSOhx4Ormg\nmf06WveDZnbqPuKY7+7HA18kdBRuI3QGyoCLo/fyL4D3Ru+nw4CRKcu/KukTSj0nAtdFnRyACcDi\nXCV8KN2k/0DUy9tE+ED/fIB5LgRudPced99FqEmeR/hwTSIkWqIe4epoPWcC97h7e7TMrYO0f0E0\nX5u7dwELgP+O2vydu2+O5rspajPpHnff5u49hB74DDM7HJjg7v8TzfN9QgllIH9LJh93vx040d27\n3X038ARwyCDLpW6TH0bLb4tiTo3vV9H100A10JjG+gbaxkl/GGihA4h9M3CKmZ0HlLv7+9x9cb95\n5kZx/ixa918JvbtT3f0pYLq7t7p7nNBLS7ZzJtHrG/X0FkW3nwa2E3qKEL6gkl+6+/Ibd+919/WE\n9+X0/S2QopnQ0XgLUOvunxmg/r+/bb3P9s1sBCGxJV//1cBfgNcdQJz9LSJsm5MJX1xXAO8YYL66\nqO23AccROj3/2i++OuB24CspnYSfEH6NHEH4XNwdzTeQ30fXzwO73P0Bd08QvkCnEL5MqlI+Z9cT\n5c/oi+p04KGU9fXfnskv0REMb7n0gJVqeSe1pn8msMjMju83TwPQknK/hZAcGgi9x8QAj40D1veb\nPpAJqfO5e0cUS8MAy6cmztQdXr2En/ITUqdHXwg9g7S7LXkjauv66HnHCV9k3xlkuaSBtsmU/vFF\nJQOi+A50fanPdxsDSDd2d78tKkl9EZgb7Wz8cL/Z6oBaYElKLXYMMD7qTX472sEP4fW9J+V26uuR\n+jxuIfRgHyR8oaezE3yg1zYt7v64mX0A+ADwCzO7m1DGSrW/bb2/9scDMXfv/5z398W+r7hTO1tr\nzexGwi+jX/abdQfwaLIzZGY3AJ8APh3dn0TYZ3CXu38lZf3vTbl9q5l9mlAeHehXU1t03Uso35By\nvxyoZ+/tl/o5PZZQqutImdY6wDoAet099bGsK9Wefh93f5DQUz+930ObCG/0pPHRtE3AuH4/Q5OP\ntQBjU6Y3DNLsFkKyBsDMxkc14sHa3JctUTzJXseIqHa4P18m1GWPisoc9+xnfoYYXybWl3bs7v5j\nd38toSd8AqE3mWo90Oruc1MuUzzswP4g4Wf8CVFJ48aU5fb1Wt8CXBRd/uru29N4Tunon4zrkzfc\n/XfufjYwk/Al9i/9lj3Y124LEI/2UQx1HXsxs/lmVpUyqYLwuva3mr23dW90wczGAPcSynP/lrLu\nUfbqPaqDrT8drcColPuTUm6fA9w3xPVmXckn/ag8YsDSfg/9AbjKzMrNbCThZ+c9wCqgiaiGHNUJ\nJxFq1Y8CC82sNuolXjZIs3cBf2dm9VHt9/eEn7f3EOqHyQ/n1ew/GS+L4rk4un8Ve5JTN6EnO5BG\n4Hl37zSzYwj1x1H7We4PwHsBzGxC1GY6XxaDGWwb78++Yu9jZp8xsysB3H0dsJKwY7EbGBVt+9WE\nHbCXJp+Xmd0SxdMILHX3djObSSjLJdt5lGifipnNIaXTENWwXyHsqE6ntJOuDYT6NmZ2OaGEhpm9\nO9rhnCy7LY2eZ6qhbmui9fYQkuvVUZtzCCWug0l2NwL/FK2vnvCFPFBMtwKXm9m0qLZ+VUq7XwL+\n7P0OGiCUpx41s0Oj9Z9H6Gj9bYixLgNGpPzqu4Y927h/PT+vlWrSf8DMlkZ1/duAq939+X7zXA+s\nJdT0niR8aG6Lyjr/QNg5s4RwFMxl7r4TuBv4K+HIgUUM/DMSd38M+AawmHDkx9PALR4Osfsa8FAU\nWx3wqX09kSiey4BPmdky4K2Eow8Afgf8l5n1L2lAOFrhmug5XAt8BHiPmV0WPY9rzOx3/Zb5NFAf\nxfYg8DUf4JjqAzDgNk5juX3FnupXwDvMzKOYu6JpzxFKRxsJySH5eiaf1/3R6/kj4Cwz86jNDwOv\nN7MPEnawz7RwtMj1hP0bqW4BJgJ3prsx0vBF4MNm9gJhB/ZL0fQ7gRPMbFm0TY4AvtVv2aFu61TX\nAAui7XQH4ciytftawMzeFM3/S8I+qKVmlkyQVwDnm9nLhP0lvyFsN8zsq2Z2DfR9Xj5HOEpnKeFX\nx9eidVwNvDn5eY4u17j7EsIvtbui9j9PONptSKUVd+8kfK5uNrPFhKPn4oTEfzSh03fAzOyXZpbV\nw8VjGk9fZPiZ2d8Dl7r73+c6Fhl+0a+ldqCu336OvFeqPX2RjIlKex8n/AqUImHh0NLkocGXE/7t\nXlAJH5T0RYaVmb2RUIK42/P0b/gyZB8C/jUqR72f8D+QgqPyjohICVFPX0SkhOT1n7Oam9sO6mdI\nfX0tLS0d+58xxxTn8CqUOKFwYlWcwy+TsTY0jB5oOAugyHv6FRVp/6kxpxTn8CqUOKFwYlWcwy9X\nsRZ10hcRkb0p6YuIlBAlfRGREqKkLyJSQpT0RURKSEYP2bRwerwzona+Sjij0gnA1miWb7j7wYzS\nKCIiByBjSd/MziacguyUaKjgZwjno/2kuw94RiQREcmsTPb0H2TPcKPbCeeTLJyDaEUKXCKRoDee\nIB4P13vfjhNP0Hc/nrykLJNIJKfB6K0dtLR0EE8kiMfDuhPQN18iAfHouu+x1Pv9H09AggRE0wHi\nCaDvsTBfeB7JgeuTy+9ZNhEmh/tATU0lHR1dKcvtPV+0FlLuJifuuZ8yNE1ir7sp01MW7r+ePTcH\nmymoqh7B7t17zunSf0ScI2fXc+r8ya9e8CBlZewdM3svoczTSzjhSCXh/KXX7euE0T09vYlC+rOF\nSH+98QS7O3vY3dXD7q7e6HYvnV29dHb3RNfh0tUdp6u7N1x6wu3ulOvkpac3TndPLz29ib77fZee\nOD3xBL29cXp6Na5WITvykPF87dr+J/RL26D/yM140jeziwgnMT6PcGb4re6+2Mw+AUxz9+sGW/Zg\nh2FoaBhNc3Pb/mfMMcU5vIYzzt54nJ27emjb1c3OXd3s3N3Nzl09dHT20LG7m47d4fauvksvu7p6\n+pJ7V098WOJIVRaLUVERo6KsjIqKMirKY5SXxagoL6O8rIzy8hgVZTHKy8soL4v1XcqiaRVlMWKx\nGOXlMcpi0fSy5G0oS96OHiuLwejR1ezq6CKWMi0Wi66j+WMAMfqWTd6Opc4bzReLHk/eLouFZfse\nB+ibNzzvWHJd7FmW6H5ynvr6WrZv7yCWkvP2Wr5v4l5XfTPG+i2z56HYAPMPeHOvhVOn98/C48aP\nYtvW9n7r2XNn7KhKKsqHdqzNvoZhyPSO3IWEMz+dH407nXpKsbuAGzLZvshA4vEErR1dbG/vZHtb\nF9t3dtLa3sWOnV207uyitaOL1o5u2ju62Ll7sHPMD6xqRDnVVeXUVI+gfkw1NZXlVI0op6qynOrK\ncipHRPdHhNvj62vp2t3NiIoyKkeUh+uKMkYkL+UhsffdLi+jrGzQz3PGFNQXfnVeDynWp2FcLWW9\nvVlvN5M7cscSTgl4TnTeTszsduBf3H0FsAB4IVPtS+nq7ull47YOmrfvYsv2XWxp3c3WHbvZ1trJ\ntrbdbG/r6qsjDyQWg9E1I6gbVcX0xlGMrBnBqOgysnoEI6srqI2ua6oqqK0Ol5rKigNOyIWSTKV4\nZPIr8XLCiYhvTTkp/c+B35pZB+FUY+/OYPtSxBKJBNtaO9mwdSfrt3awaVsHm1rC9ba2zlftFINQ\naqgbXckhU8ZQN7qK+lFV1I2upG5kFWNHVTJ2ZCWjR1YyqnpETnrTItmQsaTv7jcSznbf3y8y1aYU\np12dPazd3M6aTW00Ne9kXXM7TVt20tn16p/GdaMqOWL2eOpHVtJQV82EuhrGj6lmwthq6kZVKZlL\nySuM4peUjO6eOGs2t7FifSsrN7SyckMbm7d17HXEW3lZjEnja5kyfiRTJoxk8vhaJo2rpbG+hurK\nCpVMRPZBSV9yaldnD8uaduBrW1jetIOVG9ro6d1zxEtNVQU2o44ZE0czc+JopjeOYtL42iEf1SBS\n6pT0Jat643FeWdfKCyu3sWTVNlZuaOvbqRqLwfTGURw6dSxzpoxl9pQxNNbXhMP/RGRYKOlLxrXv\n6ub5V7byzPItvLhyG7s6w2GQZbEYsyePZu7MemxGHXOmjKWmSm9JkUzSJ0wyYsfOLp72zTyxdDO+\ndnvf0TQTxlZz8hETmT97HHNn1ivJi2SZPnEybHZ19vD0y8088sJGlq5p6Uv0c6aM4djDJnDsYQ1M\nGV/b9+9GEck+JX05KIlEgmVNO1i0eD1PvbyZru6wE3bO1DGcNHciJ1oD48ZU5zhKEUlS0pch2dXZ\nw8PPbeCBxevYsLUDgMa6Gk6dP4mT50+isa4mxxGKyECU9OWAbN6+i/ueXMvDz21gd1cvFeUxXnvE\nRM46Zgo2o06lG5E8p6QvaWna3M7N9zoPLV5HIhH++XrByTM589gpjKmtzHV4IpImJX3Zp7Wb2/n9\nQyt4Zlk47cH0xlGc/9oZnDS3UX+QEilASvoyoE0tHfz+oZU8/tImEoQjcN76hnnMmqCjb0QKmZK+\n7GXn7m7ufHglf3l6Hb3xBDMmjuKSs+Ywf/Y4GhvHaEwbkQKnpC9AOLHIosXruOOhlbTv6qaxroaL\nzzqEE+c2ahgEkSKipC+s2dTGL/53KSs3tFFdWc5lC+ZwzonTGVGhmr1IsVHSL2Gd3b3c+fBK/vT4\nWuKJBCcfOZHLzz6UsaOqch2aiGSIkn6JWrmhlZ/c/RIbt3XQUFfNFQvncuTscbkOS0QyTEm/xPT0\nxvnDI6v4wyOriScSnHvidC4+6xCqRpTnOjQRyQIl/RKyrXU3P7rzRZav28G4MVVcdeERzJtZn+uw\nRCSLlPRLxPMrtvKTu1+ifVc3r5nXyBUL51JbrZdfpNToU1/k4okEd/91FXc+vJKK8hjvOO9wFhw3\nVX+wEilRSvpFrLOrl5/e8xJPejMTxlZz7VuOYuak0bkOS0RySEm/SG1r3c33bn+ONZvaOXx6He9/\ny3wNjCYiSvrFaF1zO9+69Vla2jo585jJvP080+BoIgIo6Red5U07+O7vnmXn7h4uWzCH8187Q/V7\nEemjpF9EnntlCz+84wV6ehNcdeE8Tjtqcq5DEpE8o6RfJJ5Z1swP73iB8rIY/3TpURw9Z0KuQxKR\nPKSkXwSefrmZG37/AhXlZXzwsqOxGfrDlYgMTEm/wD3lzfzozpDwP/T3x3D49LpchyQieUxJv4C9\nuGpbSPgVZXzoMiV8Edk/HcdXoFasb+X7tz9PLAb/dMnRSvgikpaM9vTN7OvAGVE7XwWeAH4FlAMb\ngHe4e2cmYyhG67fs5Du3PUtXTy/vf/NRGjRNRNKWsZ6+mZ0NzHf3U4Dzge8AXwB+4O5nAMuBKzPV\nfrHasbOLb9+6mPZd3bzz/LmcYA25DklECkgmyzsPApdFt7cDI4EFwF3RtLuBczLYftHp7unl+//9\nHFtbO3nzGbM585gpuQ5JRApMLJFIZLwRM3svocyz0N0bo2lzgF+5+6mDLdfT05uoqNDJPQASiQTf\n/M+nWfRMEwuOn8aH33q8/mkrIoMZNDlk/OgdM7sIuAo4D1iWTlBJLS0dB9V2Q8NompvbDmod2ZBO\nnHc/sopFzzQxZ8oY/vF1c9iypT1L0e1RTNszXxRKrIpz+GUy1oaGwUfTzejRO2a2EPgU8AZ33wG0\nm1lN9PBUYH0m2y8Wz6/Yyu8fXMH4MVVcd8nRjNCvHxEZokzuyB0LfAN4o7tviybfB1wS3b4E+N9M\ntV8stu7YzU/ufony8hjXXnwUY0dqeGQRGbpMlncuByYAt5pZcto7gZvM7GpgNfCLDLZf8Hp649xw\n5wu07+rmioXGrEljch2SiBS4jCV9d78RuHGAh87NVJvF5tY/L2fF+lZOOXISZx2rI3VE5ODpH7l5\navHyLdz3VBNTJ4zkioWmI3VEZFgo6eeh1p1d3PzHJVSUl3H1RUdSVakdtyIyPJT080wikeDm/1lK\na0c3l551CNMaRuU6JBEpIkr6eebBZ9ezePkW5s2s55yTpuc6HBEpMkr6eaR5+y5uuX8ZtVUVXHXh\nPMpUxxeRYaaknycSiQS/vNfp6o7z1nMPY9yY6lyHJCJFSEk/Tzz24iZeXLmN+bPHccqRk3IdjogU\nKSX9PLCjvZNb7l9G5YgyHZ4pIhmlpJ8Hbror/Ov24jMOYUJdzf4XEBEZIiX9HFuyahsPPNXErEmj\nOedEHa0jIpmlpJ9DvfE4v7l/GbEYXHG+UVamso6IZJaSfg498Mx61jXv5NzXzNRgaiKSFUr6OdLW\n0cUdD66gpqqCd7xhXq7DEZESoaSfI3c8tJKOzh4uOn02daOrch2OiJQIJf0caNrczqLF65g8vpbX\nHT811+GISAlR0s+B3y16hUQCLn/dYVSU6yUQkexRxsmyl9du57lXtmLT6zjqkHG5DkdESoySfhYl\nEgl+t+gVAC5dMEf/vBWRrFPSz6Jnl29ledMOjjtsAnOmjs11OCJSgpT0syQeT3D7g68Qi8HFZ83J\ndTgiUqKU9LPk8SWbWNe8k9PmT2bqhJG5DkdESpSSfhbEEwn+8OhqymIx3nTarFyHIyIlTEk/C572\nZtZv2ckp8yfSoFE0RSSHlPQzLJFIcPcjq4jF4MJTZuU6HBEpcUr6Gfbs8q2s3dzOa+ZNZNK42lyH\nIyIlTkk/g5K9fIALT5mZ22BERFDSz6iXVrWwckMrJxzewLSGUbkOR0REST+T7n18DQAXqJcvInlC\nST9DmprbeWHlNg6fXsfsyTpBiojkByX9DPm/J9YCsPAknfdWRPJHRSZXbmbzgTuBb7v7983sZuAE\nYGs0yzfc/Z5MxpALrTu7ePTFTTTW1XDMoRNyHY6ISJ+MJX0zGwlcD9zf76FPuvsfMtVuPvjLM+vo\n6Y1z7knTdbJzEckrmSzvdAIXAOsz2Ebe6e7p5c9PN1FbVcFpR03KdTgiInvJWE/f3XuAHjPr/9B1\nZvZhYDNwnbtvGWwd9fW1VFSUH1QcDQ2jD2r5A/V/f1tNW0c3l5x9KNOn1qe9XLbjHCrFOfwKJVbF\nOfxyEWtGa/oD+BWw1d0Xm9kngM8B1w02c0tLx0E11tAwmubmtoNax4G6Kxo++eS5jWm3nYs4h0Jx\nDr9CiVVxDr9MxrqvL5OsJn13T63v3wXckM32M23VxlZWbWzj2EMnMH5sda7DERF5lawesmlmt5vZ\nIdHdBcAL2Ww/0x54Juy+WHDclBxHIiIysLR6+mZ2nLs/cyArNrMTgG8Cs4BuM7uUcDTPb82sA2gH\n3n1g4eavXZ09/O2lTYwfU8382eNzHY6IyIDSLe98E3jdgazY3Z8i9Ob7u/1A1lMoHntxI53dvVx4\nykwdpikieSvdpL/GzB4AHgO6khPd/bOZCKrQJBIJ/vLMOsrLYpxx9ORchyMiMqh0k/7K6CIDeGV9\nK03NOzlxbiNjR1XlOhwRkUGllfTd/fNmNh6Y7e5PmlmZu8czHFvBeHBxtAP3WO3AFZH8ltbRO2b2\nD4TSzs3RpOvN7MpMBVVIOrt6ecI3M35MNXNnpv9nLBGRXEj3kM2PAMcAzdH9jwJXZySiAvP0smY6\nu3o5Zf4kymLagSsi+S3dpL/D3fv+Huvuu0jZoVvKHnl+AwCnzdc4OyKS/9LdkbvFzN4J1JjZ8cDl\n7On1l6xtrbt5aVULh04dy0Sd9FxECkC6Pf1rgJOA0cBNQA1wVaaCKhSPvriRBHCqRtMUkQKRbk//\nfHffa2A0M7sG+NHwh1QYEokEj7ywkYryMl4ztzHX4YiIpGWfSd/MjgOOBz5qZqn1ixHAZynhpL9y\nQxsbtnZw0txGaqtH5DocEZG07K+nvxuYCNQBZ6RMjwP/kqmgCsGjL2wE0IlSRKSg7DPpu/sSYImZ\n/dndH8tSTHkvHk/whG9mVM0Ijpg1LtfhiIikbX/lne+6+z8D3zCzRP/H3f3MjEWWx3ztdlp3dnHW\nsVOoKM/q6NQiIgdlf+Wdn0XXn850IIXkiSWbALQDV0QKzv6S/ngzO6AhlYtdbzzOk97MmNoR2AwN\nuyAihWV/Sf8z+3gsAfx5GGMpCEtXb6d9VzdnHz9V4+aLSMHZX9K/HbjX3ZdlI5hC8LhKOyJSwPaX\n9MuA75jZTOCvwJ+A+9x9R8Yjy0M9vXGefrmZsaMqOWx6Xa7DERE5YPs89MTdv+fuFxL+oPVbwlAM\n95vZX83s89kIMJ+8tKqFnbt7OMkaNaKmiBSkdE+i0kWo3/8ZwMwagXMzGFdeemJpVNqZNzHHkYiI\nDE1aSd/MHiLsuE3VY2anA19y93XDHlme6Y3HeXb5VupGVXLI1DG5DkdEZEjSHXDtPuBwwo7dXuAt\nwBqgBfg5cF5Gossjy5t20L6rmwXHTVVpR0QKVrpJ/3R3Ty3n3Glm97j7hWZ2USYCyzfPLNsCwHGH\nTchxJCIiQ5fuGAKNZtaX7cxsLDDLzOqAsRmJLI8kEgkWL9tCdWU5c/WHLBEpYOn29L9LGHhtNWGE\nzTnAl4E3Aj/OUGx5Y92WnWzevosT5zYyokJj7YhI4Uo36f8vMBWIAa8FHgNGufu3MhVYPlFpR0SK\nRbrd1v8BZhNOnvI0sCO6XRIWL2umLBbj6Dnjcx2KiMhBSbenv9Xdr8xoJHmqpa2TlRvamDeznpE6\nQ5aIFLh0k/4dZvY24FGgJznR3ddkJKo8snh5KO0cq9KOiBSBdJP+0cDbgK0p0xLAjGGPKM8sVj1f\nRIpIukn/ZKDe3TszGUy+6eruZemaFqZOGMmEsTW5DkdE5KClm/SfAKqBA0r6ZjYfuBP4trt/38ym\nA78CyoENwDvy+Yvk5bXb6e6Jc9Qh2oErIsUh3aQ/DVhlZkvYu6Y/6DlyzWwkcD1wf8rkLwA/cPfb\nzOwrwJXADQccdZY8v2IbAPMP0cnPRaQ4pJv0vzyEdXcCFwAfT5m2ALgmun038FHyOOm/sHIrVSPK\nOWyaxs4XkeKQ7tDKiw50xe7eQxiJM3XyyJRyzmZg8r7WUV9fS0VF+YE2vZeGhtFDWm7Ttg42bO3g\nNUdMYsrkzI80MdQ4s01xDr9CiVVxDr9cxJpuTz8T9jtUZUtLx0E10NAwmubmtiEtu+iZMFr04dPG\nDHkd6TqYOLNJcQ6/QolVcQ6/TMa6ry+TbA8k025mycNgpgLrs9x+2l5YEY5Ona+duCJSRLKd9O8D\nLoluX0IY0yfv9PTGeWl1CxOibstpAAAMpUlEQVTra2is06GaIlI8MlbeMbMTgG8Cs4BuM7uU8Aev\nm83samA18ItMtX8wljftoLOrl/lHqZcvIsUlY0nf3Z8iHK3TX96fW/f5laG0c5QO1RSRIqPB4Qfw\nwoptVJTHsOk6YYqIFBcl/X7aOrpYu7mdw6bVUVV5cIeLiojkGyX9fnzNdgDmzlQvX0SKj5J+P0tW\ntwAwT0lfRIqQkn4/S1a3UFVZzqxJhfOvPhGRdCnpp2hp62Tjtg5seh0V5do0IlJ8lNlSLF0TSjtz\nZ6i0IyLFSUk/her5IlLslPRTLF3dwsjqCqY3jsp1KCIiGaGkH2nevostO3ZjM+opK9vvAKAiIgVJ\nST+i0o6IlAIl/cjSKOnrT1kiUsyU9IFEIsGSNS2MGVnJlPG1uQ5HRCRjlPQJ9fwd7V0cPr2OWEz1\nfBEpXkr6wLKmHQAcNi3z58IVEcklJX1gWVMYZO3waXU5jkREJLOU9Ak9/erKcqY1jsx1KCIiGVXy\nSb+1o4sNWzuYM3Us5WUlvzlEpMiVfJZ7RfV8ESkhJZ/09+zEVT1fRIpfySf9l5u2U14W45ApY3Id\niohIxpV00u/s7mX1xjZmThpN1QidD1dEil9JJ/2V61vpjSdUzxeRklHSST95fL7q+SJSKko86Yed\nuIeqpy8iJaJkk348nmD5uh1MHFfLmNrKXIcjIpIVJZv012/dye6uXg6dqqN2RKR0lGzSX7G+FYBD\npqi0IyKlQ0l/snr6IlI6SjrpV1aUaZA1ESkpFdlszMwWALcBL0aTnnf3D2QzBoDdXT2s29LOoRpk\nTURKTFaTfmSRu1+ag3b7rN7YRiKBhl4QkZJTkt1c7cQVkVKVi57+EWZ2FzAO+Ly7/99gM9bX11JR\ncXBj4jQ0jH7VtKatHQCcOH8yDfX5cSL0geLMR4pz+BVKrIpz+OUi1mwn/WXA54FbgUOAv5jZoe7e\nNdDMLS0dB9VYQ8NompvbXjV96aptjBlZCd09Az6ebYPFmW8U5/ArlFgV5/DLZKz7+jLJatJ393XA\nb6O7r5jZRmAqsDJbMbS0ddLS1smxh04gFotlq1kRkbyQ1Zq+mb3NzD4a3Z4ETATWZTOGFevDeDva\niSsipSjb5Z27gN+Y2UVAJfC+wUo7mbJiQ3InrpK+iJSebJd32oA3ZbPN/laubyUGzNY/cUWkBJXU\nIZvxeIKVG9uYPGEkNVW5OHBJRCS3Sirpb9jWQWdXL7MnFc4hXSIiw6mkkv6ajeHwqJlK+iJSokoq\n6a/eFJL+rEmq54tIaSqppL9qYxsxYHrjqFyHIiKSEyWT9OOJBGs2tTFpfC1VlQc3tIOISKEqmaTf\n3LKL3V29queLSEkrmaSfrOfPnKikLyKlq3SS/sbkTlwlfREpXSWT9FdFSX96o5K+iJSukkj6iWgn\nbmN9DbXV+ieuiJSukkj6W3fsZufuHtXzRaTklUTS3/OnLCV9ESltJZX0Zyjpi0iJK42kv7Ed0OGa\nIiJFn/QTiQSrN7Yyfkw1o2pG5DocEZGcKvqkv729i9aObv0TV0SEEkj6azeH0s4MDbImIlL8Sb+p\nOST9aUr6IiIlkPSjnv60hpE5jkREJPeKP+k3t1M1opwJdTW5DkVEJOeKOul398TZsLWDaQ0jKYvF\nch2OiEjOFXXSb9rcRm88oXq+iEikqJP+qg2tAExrUNIXEYFiT/rrk0lfO3FFRKDYk/7GKOmrvCMi\nAhR70l/fyrgxVYys1vALIiJQxEm/raOLba27Vc8XEUlRtEm/qXknoJ24IiKpijjpJ4df0E5cEZGk\n4k360fAL09XTFxHpk/WzhJvZt4GTgQTwz+7+RCbaaWpup6K8jInjajOxehGRgpTVnr6ZnQUc5u6n\nAFcB38tEO/F4gnXNO5kxcTQV5UX7Y0ZE5IBlOyO+Hvg9gLsvAerNbMxwN7Jlxy66euLMnKwTp4iI\npMp2eWcS8FTK/eZoWutAM9fX11JRUX7AjYweU8NpR09h4cmzaGgojMSvOIdXocQJhROr4hx+uYg1\n6zX9fvY59GVLS8eQV3zVBXNpaBhNc3PbkNeRLYpzeBVKnFA4sSrO4ZfJWPf1ZZLt8s56Qs8+aQqw\nIcsxiIiUrGwn/T8BlwKY2fHAencvjK9lEZEikNWk7+6PAE+Z2SOEI3euzWb7IiKlLus1fXf/RLbb\nFBGRQAexi4iUECV9EZESoqQvIlJClPRFREpILJFI5DoGERHJEvX0RURKiJK+iEgJUdIXESkhSvoi\nIiVESV9EpIQo6YuIlBAlfRGREpLrk6hkRLZOvj5UZjYfuBP4trt/38ymA78CygnnF3iHu3fmMkYA\nM/s6cAbhffJV4AnyLE4zqwVuBiYC1cAXgWfJsziTzKwGeIEQ5/3kYZxmtgC4DXgxmvQ88HXyM9a3\nAR8DeoDPAs+RZ3Ga2VXAO1ImnQicBtxAyFHPufv7shVP0fX0s3Xy9aEys5HA9YQPfNIXgB+4+xnA\ncuDKXMSWyszOBuZH2/F84DvkYZzAm4An3f0s4O+Bb5GfcSZ9GtgW3c7nOBe5+4Lo8gHyMFYzGw/8\nG3A68EbgIvIwTnf/aXJbEuL9BeHz9M/ufhow1szekK14ii7pk6WTrx+ETuACwlnEkhYAd0W37wbO\nyXJMA3kQuCy6vR0YSR7G6e6/dfevR3enA03kYZwAZjYXOAK4J5q0gDyMcxALyL9YzwHuc/c2d9/g\n7u8lP+NM9Vng34HZKRWIrMZZjOWdAzr5era5ew/QY2apk0em/ATdDEzOemD9uHsvsDO6exXwR2Bh\nvsWZFJ2YZxqhx3dfnsb5TeA64J3R/bx73VMcYWZ3AeOAz5Ofsc4CaqM464HPkZ9xAmBmJwFrCaWo\nlpSHshpnMfb0+9vnydfzUF7Fa2YXEZL+df0eyqs43f1U4O+AX7N3bHkRp5ldATzq7isHmSUv4ows\nIyT6iwhfUD9l7w5ivsQaA8YDFwPvAn5OHr72Kd5D2P/UX1bjLMakX4gnX2+PdvABTGXv0k/OmNlC\n4FPAG9x9B3kYp5mdEO0Ix90XE5JTW77FCVwIXGRmjxE+/J8hD7cngLuvi8pmCXd/BdhIKJPmW6yb\ngEfcvSeKs438fO2TFgCPEKoP41OmZzXOYkz6hXjy9fuAS6LblwD/m8NYADCzscA3gDe6e3LHY97F\nCZwJfATAzCYCo8jDON39cnc/yd1PBm4iHL2Td3FCOCLGzD4a3Z5EODLq5+RfrH8CXmdmZdFO3bx8\n7QHMbArQ7u5d7t4NLDWz06OHLyaLcRbl0Mpm9jVCMogD17r7szkOqY+ZnUCo7c4CuoF1wNsIP/uq\ngdXAu6M3Rs6Y2XsJNdKXUya/k5Cw8inOGkL5YTpQQyhLPAn8kjyKM5WZfQ5YBdxLHsZpZqOB3wB1\nQCVhmz5DfsZ6NaH8CPAlwmHF+RjnCcCX3P0N0f0jgB8TOt5/c/cPZyuWokz6IiIysGIs74iIyCCU\n9EVESoiSvohICVHSFxEpIUr6IiIlRElfioaZvT26nmRmt2Vg/WVmttTMYinTFpjZw8PdlkimFOPY\nO1KCzKycMJjVr919I3sGixtOJwLPuLuOc5aCpaQvxeJnwEwz+xPwXuBhd59mZjcDW4B5wJHAJwjD\nMR8dzfM+ADP7CmGM8xpgEfCxAZL7eYR/gfZXbmY3AMcRRlG90N3bzexK4BqggzBkwP9z91YzSwAj\n3L3HzN4FnOPubzezVcBvgUOAdxP+IFUPjADudvcvH+Q2ElF5R4rGvwHN7n7eAI9NdPcLCf8w/gFw\nLfAa4F1mVmdmlwFT3f0sd38NcChhtM7+Bkv684DPRUMsdAMLzWwG4Z+sr4/GUV8LfCiN57HM3S8D\nziV8MZwBnEoYp0efVzlo6ulLKfhrdN0ELHH37QBmthUYC5wNnGJmD0TzjQVmp64gOifDGHdfN8D6\nl7r7ppQ26oDjgadSxn16gNDr359HUmL+gpndShjW+iZ3j6exvMg+KelLKegZ5DaEYW07gRvd/T/2\nsY7XAX9JY/3JdfYvDQ00DcLYNqm6ANx9s5kdA5xCGOL4STM73t137SNGkf3Sz0UpFnFC7XsoHgYu\nNrMKADP7rJkd1m+ewUo7g3kKOCEavAzCmZEei263EgaIg/Ar41XM7DzCvoG/uvvHgHag8QDaFxmQ\nkr4Ui/XARjN7inBqxwPx34RyyiNm9ihhKOEV/eY5i7CDNy3u3kQYM/8+M3sQaCCcFxXga8CfzOyP\nhNE2B1wF8BEzeygqO/3J3Ven277IYDTKpohICVFPX0SkhCjpi4iUECV9EZESoqQvIlJClPRFREqI\nkr6ISAlR0hcRKSH/H3+fNE9WNHTIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "tags": []
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# parameters\n",
    "t_final = 72     # hours\n",
    "V = 5.0          # liters\n",
    "k = 0.625        # liters/hour\n",
    "u = 25*0.625     # mg/hour infusion rate\n",
    "\n",
    "# create a model\n",
    "model = ConcreteModel()\n",
    "\n",
    "# define time as a continous set with bounds\n",
    "model.t = ContinuousSet(bounds=(0, t_final))\n",
    "\n",
    "# define problem variable indexed by time\n",
    "model.C = Var(model.t, domain=NonNegativeReals)\n",
    "model.dCdt = DerivativeVar(model.C, wrt=model.t)\n",
    "\n",
    "# the differential equation is a constraint indexed by time\n",
    "model.ode = Constraint(model.t, rule=lambda model, t: V*model.dCdt[t] == u - k*model.C[t])\n",
    "\n",
    "# initial condition\n",
    "model.C[0].fix(0)\n",
    "\n",
    "# transform the model to a system of algebraic constraints\n",
    "TransformationFactory('dae.finite_difference').apply_to(model, nfe=50, method='backwards')\n",
    "\n",
    "# compute a solution using the Pyomo simulator capability\n",
    "tsim, profiles = Simulator(model, package='scipy').simulate()\n",
    "\n",
    "# plot results\n",
    "plt.plot(tsim, profiles)\n",
    "plt.xlabel('time / hours')\n",
    "plt.ylabel('mg/liter')\n",
    "plt.title('Blood conctration for a steady infusion of ' + str(u) + ' mg/hr.');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PPEH2YEfLvc2"
   },
   "source": [
    "## Optimizing the Infusion Strategy\n",
    "\n",
    "The simulation results show that a constant rate of infusion results in a prolonged period, more than twelve hours, where the antibiotic concentration is below the minimum therapuetic dose of 20 mg/liter. Can we find a better strategy?  In particular, can a time-varying strategy yield reduce the time needed to achieve the therapuetic level?\n",
    "\n",
    "For this case, we consider $u(t)$ to be a time-varying infusion rate bounded by the capabilities of the infusion pump.\n",
    "\n",
    "\\begin{align*}\n",
    "0 \\leq u(t) \\leq 50 \\mbox{ mg/hr}\n",
    "\\end{align*}\n",
    "\n",
    "Next we need to express the objective of maintaining the antibiotic concentration at or near a value $C_{SP}$ = 25 mg/liter.  For this purpose, we define the objective as minimizing the integral square error defined as\n",
    "\n",
    "\\begin{align*}\n",
    "\\mbox{ISE} & = \\int_0^{t_f} (C(t) - C_{SP})^2 \\, dt\n",
    "\\end{align*}\n",
    "\n",
    "The next cell demonstrates how to formulate the model in Pyomo, and compute a solution using ipopt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 441
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 2244,
     "status": "ok",
     "timestamp": 1554236085393,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "ee6k_bfnLuii",
    "outputId": "831cec23-29e0-425b-bc18-6541920996a4"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGoCAYAAABL+58oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XuYZFV97/93T/fcGZhhaC4CKkT9\neuGYCNGIAR0NFxENiYjmBAkKRvGnOblITPIYb6gx0XiJyjExRiX6iz9vMUIwhgOKqOgJosZo4CsY\nwcAA08z0DDPMMDN9+f2xdg9F2z1TPV27qrrq/XoeHqp2Ve21+juX+sxaa+81MDk5iSRJUi9Z1OkO\nSJIktZoBR5Ik9RwDjiRJ6jkGHEmS1HMMOJIkqecYcCRJUs8x4Eg1iojJiLg1Im6OiB9FxL9ExLHV\na+si4tYWtvXiiLi2Vedrh4g4LCJ+dT8/+9sNj2+OiMNa17NZ2/x/IuLuiHhdi84XEfH06vGvR8RH\nqsfXRsSLZ3j/kRHxgybO+0sR8cTq8asj4i2t6K+0kAx1ugNSH1iXmXcARMTbgb8CntfZLnWNZwKn\nAJfP5UMRcTjwWuBvATLzsa3v2ozOBl6XmX/XovP9OuXv4esy8/PA5/f25sy8EziuifO+FPg68P3M\n/MC8eyktQAYcqb2+DPzMiEVELAPeS/nCnwC+CLw2M8erf4l/EFgLPAD8UWb+a0QsAt5Xne9u4Kuz\nNRoRfwS8AhgD/hl4TWZORsT/Ai6ijOYm8LLMHImIjwG3A08DHgP8CDgrM7dHxAnAh4BVwF3ASzLz\nJxHx+KqfRwA7gZdm5rcjYh3wduBa4NeAZcBLgK3AB4ChiDgA+GPgeuBTwPGZ+YxqdOdtwBJgG3Bh\nZn6vet9REXEz8MSqvaMz8479+Zma+bWofoYTgcdFxNGZ+aZpn3sZ8BrK36t3Aedl5u0R8RLgTOA+\n4OTq1+Ac4FjgT4BdEbEG+A/gxZl5SnXK/xER/1bV80vVz3Q0cGtmDlW//m+hhC6AbwGvAs4Dfgv4\n1Yg4FDgQOCozX1aNHn4MeBgwCrwiM7+D1IOcopLaJCKWAC9m5tGK36N8eT0BOJ7yRfg/qy+x/w/4\nQDVK8TLgkxGxCng2cBrweOAZwNNnafek6nM/T/nX/0nACyLiqcAfUkaYHgv8lPIlPuUc4EXAzwHD\nlNEGqv78aWY+hjLi8IGqn/8E/H11/CLgCxEx9Y+oJwHfyszHAf+7+vx3KAHns5n5G9X7DgG+V4Wb\nIeAy4LczM4AvAH9Zve8C4KeZ+djM3NXws+7vz9Roxl+LzHwt8G+U4PmmaTU+tPpZTs3MRwO3Aq9v\neMtzgP9d1eYrwO9l5hVV/f4qM18zQz+eCawDgvLr+9xpr78QOAM4oerrauD3M/OvG/r57mmf+RDw\nycx8FCU4fnyGdqWeYMCR6ndtNdJwD/Bk4KMzvOdM4EOZOZaZO4D/lxJejgEOp4QKMvPblFGIJ1MC\nzZWZua36zKdnaf851fu2VmFgHfCPVZufzcwN1fs+XLU55crM3JSZY5TRhYdHxGOAQzLzX6r3fIAy\ngvBY4FDgI1U/vwGMUEZLALZm5heqx98BHj5LXxdTTdNU7R6amd+qXvsaZdRjb+b8M81yjpl+LWZV\ntXfg1FTkDH39z8y8sXq8t5+/0Wczc3s1wnQlZfRoej8vy8z7M3Oc8vtq1n5WI1PPBD5ZHfoC8EtN\n9ENakJyikurXuAbn6cBXI+L4ae8ZpkwZTBmlBIZhYHNmTs7w2sHA+mnHZ3JI4/umpmQiYniGzx/a\n8HxLw+NxYLA6157jVVAYi4jVwArgpoiYevlAyrTa6Cznmsl4Zt7X8Px/RcT5wFLK1Na+Ns/bn59p\npnPM9Gsxq4gYBC6pptQGKdN3P5pju9ONTPv8EfPs58GUf9RuAah+T21roh/SgmTAkdooM6+LiNsp\n00SNX2D3UMLAlLXVsXuAgyNioCHkTL02ChzU8JnhWZq9lxJMAIiIqXZma3Nv7q36sygzJyJiMXAk\nJVTcN9Ni32oNzpxFxNOAPwKekpm3RcSpVIuK92J/fqZWnONFlLVQT8/Me6srvM6dY7vTHdzweA2w\naZ793EgJiGuBeyNigDJV9+NpAVrqCU5RSW1UTfEEcPO0l/4ZuDAiBiNiJWWh6JXAbcAdlC/QqS/9\nwylrLL4JnB4RKyJiBWV9yUwupyw4XVOta/kn4PTq/M9vCDyvqI7tzS1Vf55fPb+Qsq7jduCOiHhB\n1c9DIuKT1c+yN7spa0dmciiwAfhp9fOdD6ysvph3Awc0rPGZsj8/03Sz/VrszaHAbVW4WUtZH3NA\nE23t7ed/fkQsq/pwBmXaa3o/X1z9+g9Rfi2m+vkz583MncBVlAXeUH4PfNFwo15lwJHqd211n5ab\ngc9Qrlz5j2nveT/w38APgW9Tvrw+U335/Abw6oi4iXLV1DmZeT9wBfANypVCX6Vc7fMzqjUs7wS+\nB/wnZQ3IJzPz34A/B75W9W01sNf7u1T9OQd4XUTcAvwm8Mpp/bwZuA64purn3lwFPCsibpjhtS9R\nRoZ+XL3vvZTplc8C36eMaNwdEXvWs+zPzzSDGX8t9vGZTwJrq/safRL4U+DoiHjXPj53BXBRRHx2\nhteupixIvql6/KVpr3+W8mt+I/CDqs/vq177PPAXETF9kfHLgOdFxH8Bb6X8+kk9aWBy0vAuSZJ6\niyM4kiSp5xhwJElSzzHgSJKknmPAkSRJPaer74MzMrK11hXQa9asYHR0+77f2OOsQ2EdCutQWIfC\nOhTWoejGOgwPrxqY6Xhfj+AMDTVzM9HeZx0K61BYh8I6FNahsA7FQqpDXwccSZLUmww4kiSp59S6\nBici3gGcXLXzdspeLSdQ9kQBeGdmzvU26pIkSXtVW8CJiGcCx2XmidXeLN8Fvgz8SWb+c13tSpIk\n1TmCcx1lQ0CAzcBKYOGsTtoPu8fG+cp317Nz11inuzInK1cu5f77d3a6Gx1nHQrrUFiHwjoU1qHY\nnzo86THDHDXczN6zrdWWvagi4uWUqapxyk7ISyi7BL86M++d7XNjY+OTC2nF9pVf/y/++vPT91CU\nJKl/PfOEo/iD3zyhziZmvEy89vvgRMRZwIXAacAvAhsz83sR8cfAm4BXz/bZuq+1Hx5excjI1pad\n78s3/JQB4FXP/x8sXbJwgtnqg5azecuOTnej46xDYR0K61BYh8I6FPtTh2OPOLCl37XTDQ+vmvF4\n3YuMTwdeBzw7M7cA1zS8fDnwwTrbb6fRrTu55Y4tPPro1Rz/mOFOd2dOWh30FirrUFiHwjoU1qGw\nDsVCqkNtl4lHxEHAO4HnZuam6tjnIuLY6i3rgB/U1X67ffvmDUwCT3ncoZ3uiiRJfa/OEZwXAYcA\nn46IqWMfBT4VEduBbcBLa2y/rW64eQMDA3BCGHAkSeq02gJOZn4I+NAML11WV5udsum+B7j1zi08\n7hFrOGjlkk53R5KkvuedjFvghps3APDkxzp6I0lSNzDgtMAP/qvcmPn4WFiLiyVJ6lUGnBa48977\nWXvgUg5c4fSUJEndwIAzT9sf2M3mbbs44pCVne6KJEmqGHDmaf3GcjPCh6014EiS1C0MOPO0/t77\nAXiYIziSJHUNA848GXAkSeo+Bpx5Wr+xCjhrV3S4J5IkaYoBZ57uuvd+Vh+whBXLFne6K5IkqWLA\nmYcdO8fYeN9Op6ckSeoyBpx5uHuTV1BJktSNDDjz4AJjSZK6kwFnHgw4kiR1JwPOPBhwJEnqTgac\neVi/8X4OXLGYA5Z7BZUkSd3EgLOfdu4e597NDzh6I0lSFzLg7Ke7N25nEtxkU5KkLmTA2U8P3sHY\ngCNJUrcx4OwnFxhLktS9huo8eUS8Azi5auftwA3Ax4FB4C7gvMzcWWcf6mLAkSSpe9U2ghMRzwSO\ny8wTgWcD7wUuAS7NzJOBW4EL6mq/bus3bmflsiEOXOEVVJIkdZs6p6iuA86pHm8GVgLrgMurY1cA\np9TYfm12j02wYXQ7DztkJQMDA53ujiRJmqa2KarMHAfur55eCHwROL1hSmoDcMTezrFmzQqGhgbr\n6iIAw8Or5vyZ2+66j8lJOPao1fv1+W7UKz/HfFmHwjoU1qGwDoV1KBZKHWpdgwMQEWdRAs5pwC0N\nL+1z6GN0dHtd3QLKL9LIyNY5f+6Ht2wA4OCVS/br891mf+vQa6xDYR0K61BYh8I6FN1Yh9kCV61X\nUUXE6cDrgDMycwuwLSKWVy8fCayvs/26uMBYkqTuVuci44OAdwLPzcxN1eGrgbOrx2cDX6qr/ToZ\ncCRJ6m51TlG9CDgE+HRETB07H/hwRLwCuB24rMb2a7N+43aWLx1k9QFLOt0VSZI0gzoXGX8I+NAM\nL51aV5vtMDY+wT2btvPIw1d5BZUkSV3KOxnP0YbRHYxPTLoHlSRJXcyAM0d71t+4B5UkSV3LgDNH\nd01tsnnIig73RJIkzcaAM0d3byr35jn8YAOOJEndyoAzR3dv2sHgogHWHrSs012RJEmzMODMweTk\nJPds2s6ha5YzuMjSSZLUrfyWnoOt23ezfeeY01OSJHU5A84cuP5GkqSFwYAzB1MB5zADjiRJXc2A\nMwf3OIIjSdKCYMCZA6eoJElaGAw4c3D3pu0sXzrEqhWLO90VSZK0FwacJk1MTLJhdAeHH7zCTTYl\nSepyBpwm3bulbLJ5+MHLO90VSZK0DwacJt29aQfg+htJkhYCA06TvERckqSFw4DTJC8RlyRp4TDg\nNGnPCM4aA44kSd3OgNOkuzdt5+ADl7J0yWCnuyJJkvZhqM6TR8RxwBeA92TmByLiY8AJwMbqLe/M\nzCvr7EMr7Nw1zujWnTzuEWs63RVJktSEpgJORDwpM787lxNHxErg/cA10176k8z857mcq9PuGXX9\njSRJC0mzU1Tv2o9z7wSeA6zfj892FbdokCRpYWl2iuqnEXEt8C1g19TBzHzDbB/IzDFgLCKmv/Tq\niPgDYAPw6sy8d7ZzrFmzgqGhete8DA+v2ud7tn2vZLTHHLO2qfcvRL36c82VdSisQ2EdCutQWIdi\nodSh2YDzk+q/+fo4sDEzvxcRfwy8CXj1bG8eraaG6jI8vIqRka37fN+P/3sUgOVDA029f6Fptg69\nzjoU1qGwDoV1KKxD0Y11mC1wNRVwMvPNEbEWOCYzvx0RizJzYq6dyMzG9TiXAx+c6zk64e5NOxga\nHOCQA5d1uiuSJKkJTa3BiYjfoExPfaw69P6IuGCujUXE5yLi2OrpOuAHcz1Hu01OTnL3pu0cumYF\nixa5yaYkSQtBs1NUrwF+Hpi6pPti4FrgI7N9ICJOoCxOfiSwOyJeQLmq6lMRsR3YBrx0v3rdRlu3\n72bHzjEe+/DVne6KJElqUrMBZ0tmbp9aMJyZOyJi194+kJk3UkZppvvcnHrYYV5BJUnSwtNswLk3\nIs4HlkfE8cCLgJH6utU9DDiSJC08zd4H5yLgycAq4MPAcuDCujrVTe5xF3FJkhacZkdwnp2ZD7mc\nOyIuAv669V3qLntGcNYacCRJWij2GnAi4knA8cDFEdH4Db8YeAN9EnBWLB1i1fLFne6KJElq0r5G\ncB4ADgNWAyc3HJ8A/rCuTnWL8YkJNozu4BGHr2JgwEvEJUlaKPYacDLzJuCmiPhyZn6rTX3qGhu3\nPMD4xCSHrXF6SpKkhWRfU1R/lZm/C7wzIianv56ZT6+tZ13gntEdABx28PIO90SSJM3Fvqaopm7k\n96d1d6QbbagCzqGrDTiSJC0k+wo4ayPiWW3pSRca2VwFHKeoJElaUPYVcF6/l9cmgS+3sC9dZ88I\nzhpHcCRJWkj2FXA+B/xrZt7Sjs50m5HNO1i+dIiVy5q9XZAkSeoG+/rmXgS8NyIeAXwDuAq4OjO3\n1N6zDpuYnGTD5h08bO1KLxGXJGmB2etWDZn5vsw8k3Kzv09Rtmu4JiK+ERFvbkcHO2XLtl3sHptg\n2OkpSZIWnKbmXjJzF2W9zZcBIuJQ4NQa+9VxG0arPagMOJIkLThNBZyI+BplUXGjsYg4CXhrZt7Z\n8p512IbqCqphLxGXJGnBaXb17NXAYyiLjseBXwd+CowCHwVOq6V3HbTnEnEDjiRJC06zAeekzGyc\nkvpCRFyZmWdGxFl1dKzTvERckqSFa6+LjBscGhGHTD2JiIOAR0bEauCgWnrWYRtGdzA0uIjVq5Z2\nuiuSJGmOmh3B+SvKppu3U3YS/zngbcBzgb+pqW8dNbJ5B8Orl7HIS8QlSVpwmg04XwKOBAaAXwK+\nBRyQme/e24ci4jjgC8B7MvMDEXE08HFgELgLOC8zd+5v5+uybcdu7n9gjEcd2ZODU5Ik9bxmp6j+\nBTgGWAx8B9hSPZ5VRKwE3g9c03D4EuDSzDwZuBW4YK4dboepBcbeA0eSpIWp2RGcjZk51zCyE3gO\n8EcNx9YBF1WPrwAuBj44x/PWzl3EJUla2JoNOJ+PiHOBbwJjUwcz86ezfSAzxyj3ymk8vLJhSmoD\ncMTeGl2zZgVDQ4NNdnH/DA+v+plj9+9eD8CjH7l2xtd7Ub/8nPtiHQrrUFiHwjoU1qFYKHVoNuA8\nETgX2NhwbBJ4+Dza3ufq3dHqbsJ1GR5excjI1p85ftsdZautJQOTM77ea2arQ7+xDoV1KKxDYR0K\n61B0Yx1mC1zNBpynAmtasCB4W0Qsz8wdlEXL6+d5vlps2LyDgQE45CCnqCRJWoiaXWR8A7CsBe1d\nDZxdPT6bcnVW19kwup2DVy1j8VCz5ZEkSd2k2RGco4DbIuImHroG5+mzfSAiTgDeBTwS2B0RL6BM\nc30sIl4B3A5ctp/9rs2u3eNs3raLxz1iTae7IkmS9lOzAedtcz1xZt5IuWpquq7ehXzETTYlSVrw\nmgo4mfnVujvSLaZ2EXcPKkmSFi4XmUwz4j1wJEla8Aw409zjCI4kSQueAWeaqREc1+BIkrRwGXCm\n2bB5B6tWLGb50mbXX0uSpG5jwGkwPjHBxi0PuP5GkqQFzoDTYNN9OxmfmHT9jSRJC5wBp8EG199I\nktQTDDgNvAeOJEm9wYDT4MF74KzocE8kSdJ8GHAaTI3gDDuCI0nSgmbAabBhdDtLlwxy4IrFne6K\nJEmaBwNOZXJykpHNDzB80HIGBgY63R1JkjQPBpzKth272bl7nOHVyzrdFUmSNE8GnMq9Wx4A4JCD\nXH8jSdJCZ8CpPBhwHMGRJGmhM+BU7t1SrqAy4EiStPAZcCp7RnC8i7EkSQueAaeysQo4aw90BEeS\npIVuqJ2NRcQ64DPAD6tD/5GZv9POPsxmZPMOVi4bYsWytpZEkiTVoBPf5l/NzBd0oN1ZTU5OsnHL\nAxy+1i0aJEnqBU5RAVu372bX2ATDXiIuSVJP6MQIzuMj4nLgYODNmfl/ZnvjmjUrGBoarLUzw8Or\nGN0xCsDRRxzI8PCqWtvrVv36c09nHQrrUFiHwjoU1qFYKHVod8C5BXgz8GngWOArEfGozNw105tH\nR7fX2pnh4VWMjGzllts2ArBi8SJGRrbW2mY3mqpDv7MOhXUorENhHQrrUHRjHWYLXG0NOJl5J/Cp\n6umPI+Ju4EjgJ+3sx3TexViSpN7S1jU4EXFuRFxcPT4cOAy4s519mIl3MZYkqbe0e4rqcuAfIuIs\nYAnwytmmp9pp6i7Gaw04kiT1hHZPUW0FntfONpuxccsDHLB8McuXeg8cSZJ6Qd9fJj45Ocm9Wx5w\n9EaSpB7S9wHnvvt3sXtswvU3kiT1kL4POFMLjL3JnyRJvcOAM7XJpiM4kiT1DANOdQWVU1SSJPUO\nA473wJEkqecYcLyLsSRJPceAs+UBVq1YzNIl9W7qKUmS2qevA87ExCQbtzzg9JQkST2mrwPO6NYH\nGBufYK3TU5Ik9ZS+DjgbNnkFlSRJvaivA849o9sBGDbgSJLUU/o64GzYVAKOU1SSJPWW/g441QiO\nU1SSJPWWvg4492ycGsEx4EiS1Ev6O+CMbufAlUtYuth74EiS1Ev6NuBMTE4yMrrd6SlJknpQ3wac\nLdt2MTY+acCRJKkH9W3AGdlc7oHj+htJknrPULsbjIj3AE8FJoHfzcwb2t0HgI1usilJUs9q6whO\nRDwDeHRmnghcCLyvne03undLGcHxJn+SJPWedk9R/QrwTwCZeROwJiIObHMfgLKLODhFJUlSL2r3\nFNXhwI0Nz0eqY/fN9OY1a1YwNFTPJdyHHXIAhxy0jMc9apjFNbWxkAwPr+p0F7qCdSisQ2EdCutQ\nWIdiodSh7WtwphnY24uj1Z2G63DGk4/ixc9+LJs23V9bGwvF8PAqRka2drobHWcdCutQWIfCOhTW\noejGOswWuNo9RbWeMmIz5WHAXW3uAwADAwMMDvbtRWSSJPW0dn/DXwW8ACAijgfWZ2Z3RUFJkrTg\ntTXgZOb1wI0RcT3lCqpXtbN9SZLUH9q+Bicz/7jdbUqSpP7iIhRJktRzBiYnJzvdB0mSpJZyBEeS\nJPUcA44kSeo5BhxJktRzDDiSJKnnGHAkSVLPMeBIkqSeY8CRJEk9p9O7iXdMRLwHeCowCfxuZt7Q\n4S61TUQcB3wBeE9mfiAijgY+DgxSNj89LzN3drKP7RAR7wBOpvw5eDtwA31Wh4hYAXwMOAxYBrwF\n+Hf6rA5TImI58ANKHa6hz+oQEeuAzwA/rA79B/AO+qwOABFxLvBaYAx4A/B9+qwOEXEhcF7DoV8E\nfhn4IOW78/uZ+cpO9K0ZfTmCExHPAB6dmScCF1L2xeoLEbESeD/lL+8plwCXZubJwK3ABZ3oWztF\nxDOB46rfA88G3ksf1gF4HvDtzHwG8ELg3fRnHab8KbCpetyvdfhqZq6r/vsd+rAOEbEWeCNwEvBc\n4Cz6sA6Z+XdTvxco9biM8nfl72bmLwMHRcQZnezj3vRlwAF+BfgngMy8CVgTEQd2tkttsxN4DrC+\n4dg64PLq8RXAKW3uUydcB5xTPd4MrKQP65CZn8rMd1RPjwbuoA/rABARjwUeD1xZHVpHH9ZhBuvo\nvzqcAlydmVsz867MfDn9WYdGbwD+AjimYcajq+vQr1NUhwM3NjwfqY7d15nutE9mjgFjEdF4eGXD\nUOsG4Ii2d6zNMnMcuL96eiHwReD0fqvDlIi4HjiK8q/Vq/u0Du8CXg2cXz3vuz8XlcdHxOXAwcCb\n6c86PBJYUdVhDfAm+rMOAETEk4H/pkzXjTa81NV16NcRnOkGOt2BLtJXtYiIsygB59XTXuqrOmTm\n04BfBT7BQ3/2vqhDRPwW8M3M/Mksb+mLOgC3UELNWZSg93c89B/C/VKHAWAt8HzgJcBH6cM/Fw1e\nRlmrN11X16FfA856yojNlIdRFo31q23V4kqAI3no9FXPiojTgdcBZ2TmFvqwDhFxQrXInMz8HuXL\nbGu/1QE4EzgrIr5F+cv89fTh74fMvLOatpzMzB8Dd1Om8PuqDsA9wPWZOVbVYSv9+ediyjrgesps\nx9qG411dh34NOFcBLwCIiOOB9Zm5tbNd6qirgbOrx2cDX+pgX9oiIg4C3gk8NzOnFpX2XR2ApwOv\nAYiIw4AD6MM6ZOaLMvPJmflU4MOUq6j6rg4RcW5EXFw9Ppxydd1H6bM6UL4jnhURi6oFx3355wIg\nIh4GbMvMXZm5G7g5Ik6qXn4+XVyHgcnJyU73oSMi4s8pf7lPAK/KzH/vcJfaIiJOoKw1eCSwG7gT\nOJcy/LgMuB14afUbuWdFxMsp8+o/ajh8PuXLrZ/qsJwyDXE0sJwyPfFt4O/pozo0iog3AbcB/0qf\n1SEiVgH/AKwGllB+P3yXPqsDQES8gjJ9DfBWym0k+rEOJwBvzcwzquePB/6GMkDyfzPzDzrZv73p\n24AjSZJ6V79OUUmSpB5mwJEkST3HgCNJknqOAUeSJPUcA44kSeo5BhxJ8xIRL67+f3hEfKaG8y+K\niJsjYqDh2LqI+Hqr25LUO/p1LypJLRARg5RN+D6RmXfz4AamrfSLwHcz03taSGqaAUfSfHwEeERE\nXAW8HPh6Zh4VER8D7gUeBzwB+GPgecATq/e8EiAi/gz4ZcpNBr8KvHaGIHMa5c6y0w1GxAeBJwE7\ngTMzc1tEXABcBGyn3HL/tzPzvoiYBBZn5lhEvAQ4JTNfHBG3AZ8CjgVeSrnR3RpgMXBFZr5tnjWS\n1AFOUUmajzcCI5l52gyvHZaZZ1LuGH0p8CrgKcBLImJ1RJwDHJmZz8jMpwCPouxmPt1sAedxwJuq\n7RV2A6dHxMMpd9/9lcxcR9kB+feb+DluycxzgFMpIehk4GmU/aj8e1JagBzBkVSXb1T/vwO4KTM3\nA0TERuAg4JnAiRFxbfW+g4BjGk8QEQcCB2bmnTOc/+bMvKehjdXA8cCNDXvLXUsZzdmX6xv6fElE\nfBr4IvDhzJxo4vOSuowBR1JdxmZ5DDBAmVb6UGb+5V7O8SzgK02cf+qc06e3ZjoGZZ+lRrsAMnND\nRPw8cCJwFvDtiDg+M3fspY+SupBDr5LmY4KyVmV/fB14fkQMAUTEGyLi0dPeM9v01GxuBE6oNo0E\nOAX4VvX4PsqmolBGj35GRJxGWcvzjcx8LbANOHQO7UvqEgYcSfOxHrg7Im4EVs7xs/9ImRK6PiK+\nCRwG/Ne09zyDsvi4KZl5B/B64OqIuA4YBt5bvfznwFUR8UXKbuEzngJ4TUR8rZo6uyozb2+2fUnd\nw93EJUlSz3EER5Ik9RwDjiRJ6jkGHEmS1HMMOJIkqecYcCRJUs8x4EiSpJ5jwJEkST3HgCNJknqO\nAUeSJPUcA44kSeo5BhxJktRzDDiSJKnnGHAkNS0ibouIk5p4359FxPqIeOl+tHFkRPxg/3o4PxHx\nSxHxxE60Lam1hjrdAUk96UXAeZl5zVw/mJl3Ase1vktNeSnwdeD7HWpfUosYcCTtl4i4FrgceD5w\nDHAd8JvAJ4CHAx+JiLcC5wIfzsxPNHzuw5n5ier1c4AB4A7gxcAS4NbMHIqIRcBbgLOrZr8FvCoz\n75+t/cycnNbPjwGbgFOqc10JfBT4haqtz2XmxRFxEfBbwK9GxKHAe4DXV/1fBvwT8AeZOd6K+kmq\nl1NUkubjecCpwGOAZwFPy8zcnypBAAAZsklEQVRzgTuBczPzb2f7YEQ8AXghcFxmPgb4PCWENHoh\ncAZwAvAEYDXw+3trf5bmfgV4SmZ+BnglsAp4LHA88JKIOCkz/xr4N+C1mfluSth6IfAU4Oeq/165\nr4JI6g4GHEnz8dnM3JGZ9wM/oozcNGszMAycGxFrMvP9mfn3095zJnBZZt5fjZx8FDhtP9q/JjMf\nAMjMdwFnZeZkZo4CPwSOneEzzwM+kplbMnMM+DBltEjSAuAUlaT52NLweBwYbPaDmXlnRDwfuBh4\nf0RcB1w07W3DwGjD81Hg0P1of9PUg4h4NPDuiHhs9ZmjKcFputXAxRHx8ur5EDCy1x9KUtcw4Eiq\n2/TgsWbqQWZ+BfhKRKwE/hL4c+B1De+9B1jb8HxtdWw+LgVuBH4tM8cj4huzvG89cHlmfmCe7Unq\nAKeoJNXtLuDnASLiRMp6GSLitIi4NCIWVVNM/w5MTvvsPwMvjogVETEEXEhZJDwfhwLfrcLNqcCj\ngQOq13ZTRm4AvgCcFxErqv6+IiLOn2fbktrEERxJdXs38MmIOAO4FriqOn4d8D+BH0XETmADJcA0\n+izwRMqIywDwFeB98+zPW4H3RMQbKFdGvRm4JCK+S1no/M6IOBZ4DWVh83ciAuDHM/RPUpcamJyc\n/g8mSZKkhc0pKkmS1HMMOJIkqecYcCRJUs8x4EiSpJ7T1VdRjYxsrXUF9Jo1Kxgd3V5nEwuCdSis\nQ2EdCutQWIfCOhTdWIfh4VUDMx3v6xGcoaGmb7ra06xDYR0K61BYh8I6FNahWEh16OuAI0mSepMB\nR5Ik9RwDjiRJ6jkGHEmS1HMMOJIkqecYcCRJUs8x4EiSpJ5jwJEkST3HgCNJknpObVs1RMQ64DPA\nD6tD/wG8A/g4MAjcBZyXmTvr6oMkSepPdY/gfDUz11X//Q5wCXBpZp4M3ApcUHP7kiSpD7V7s811\nwEXV4yuAi4EPtrkPD3Hv5h2897PfZ/sDu+d1nqMOPYDfO+fnWTQw455fkiSpjeoOOI+PiMuBg4E3\nAysbpqQ2AEfs7cNr1qyofWOvjdt3s/7e+zlw5RJWLl+8X+cYve8BfvBfm1i+chkHrlzS4h62x/Dw\nqk53oStYh8I6FNahsA6FdSgWSh3qDDi3UELNp4Fjga9Ma2+fQx11b8k+PLxqTxu/dtIxrHvSkft1\nnr+94j/55g/v5s71m9m5enkru9gWw8OrGBnZ2uludJx1KKxDYR0K61BYh6Ib6zBb4Kot4GTmncCn\nqqc/joi7gSdHxPLM3AEcCayvq/1mjU9MAjC4aP+nlpYtLaNMD+web0mfJEnS/NS2yDgizo2Ii6vH\nhwOHAR8Fzq7ecjbwpbrab9aegDM4j4CzpAo4uww4kiR1gzqnqC4H/iEizgKWAK8Evgv8fUS8Argd\nuKzG9psyPj4BwOCi/c96yxaXgLPTgCNJUleoc4pqK/C8GV46ta4298dYNYIzNK8RnFLGB3aNtaRP\nkiRpfvr+Tsbj41NrcPa/FEudopIkqasYcCaqKSrX4EiS1DMMOK24isopKkmSuooBZ7wVAccRHEmS\nuokBZ89l4vO4imqJV1FJktRNDDhTa3AcwZEkqWcYcFowRbXUNTiSJHUVA04Lp6jcqkGSpO7Q9wFn\nrLqT8dA8RnCWDC1iYMApKkmSukXfB5xWXCY+MDDAsiWDLjKWJKlLGHBaMEUF5V44rsGRJKk7GHDG\n538VFcDSxYNOUUmS1CUMOHtGcOYXcJyikiSpexhwWrAGB0rA2TU2see+OpIkqXMMOHumqOa/Bge8\nm7EkSd3AgNPCKSrwUnFJkrpB3wecsYlJBgZg0YABR5KkXtH3AWd8fJKheV4iDrDUgCNJUtcw4ExM\nzHuBMTSuwfFeOJIkddpQnSePiOXAD4C3ANcAHwcGgbuA8zJzZ53tN2N8YrJFAccRHEmSukXdIzh/\nCmyqHl8CXJqZJwO3AhfU3HZTxscn530XY2iYonLDTUmSOq62gBMRjwUeD1xZHVoHXF49vgI4pa62\n56J1U1SO4EiS1C3qnKJ6F/Bq4Pzq+cqGKakNwBH7OsGaNSsYGhqsqXvFJAMsWbyI4eFV8zrPYYfc\nD8DQ4sF5n6sTFmKf62AdCutQWIfCOhTWoVgodagl4ETEbwHfzMyfRMRMb2lqyGR0dHtL+zXd8PAq\ndu8eZ3DpECMjW+d1rp07dgGwcXT7vM/VbsPDqxZcn+tgHQrrUFiHwjoU1qHoxjrMFrjqGsE5Ezg2\nIp4LHAXsBLZFxPLM3AEcCayvqe05GRufYKgVU1RLnaKSJKlb1BJwMvNFU48j4k3AbcDTgLOBT1T/\n/1Idbc9Vq66iWrp4KuB4mbgkSZ3WzvvgvBE4PyK+BhwMXNbGtmc1PjE5720a4MH74DiCI0lS59V6\nHxyAzHxTw9NT625vrsbHJ+e90SZ4FZUkSd2kr+9kPDk5ycRki6aoDDiSJHWNvg44Y+Ot2Ukcymad\nSxcPstOAI0lSx/V1wBkfnwBoyRQVlGkqFxlLktR5fR1wxiaqEZwWTFFBmaZyqwZJkjqvrwPOnhGc\nFkxRwdQIjgFHkqRO6+uAM7ZniqpFAadagzMxOdmS80mSpP3T1wFnfGqRcavW4CwtV93vcppKkqSO\n6uuAMzZRRnCGWjhFBV4qLklSp/V1wNkzgjPYmjI8uF2DAUeSpE7q64DT8jU41XYN3gtHkqTO6uuA\n8+AanFZPUXkvHEmSOqmpgBMRZ9TdkU6YWoPTysvEwSkqSZI6rdkRnD+IiNo35my3ll9FZcCRJKkr\nNBtaNgP/GRHfAXZNHczM36qlV20yNlZdRdXiNThOUUmS1FnNBpx/rv7rKa2eopraUdxFxpIkdVaz\nAedrtfaiQ5yikiSpNzUbcK4BJoEBYAkwDPwQeFJN/WqLVl8mvtSAI0lSV2gq4GTmMY3PI+IJwIW1\n9KiNHrzRX4vX4LhVgyRJHbVfczOZ+UPghBb3pe32rMFp0QjOcu+DI0lSV2hqBCciLpl26Ghg9T4+\nswL4GHAYsAx4C/DvwMeBQeAu4LzM3Dm3LrfO+J4pqhZt1TAVcHY6giNJUic1+80+3vDfGCWoPGcf\nn3ke8O3MfAbwQuDdwCXApZl5MnArcMH+dLpVxlo+RVVdReUUlSRJHdXsGpw3R8RKICiLjTMzt+/j\nM59qeHo0cAewDrioOnYFcDHwwTn2uWXGW7zIeHDRIhYPLXKKSpKkDmt2iurXKEHkvymjPodHxG9n\n5r808dnrgaOA5wJXN0xJbQCO2Ntn16xZwdDQYDNd3C9jN20A4OA1KxgeXtWScy5fOsTYxGTLztcu\nC62/dbEOhXUorENhHQrrUCyUOjR7mfgfAk/MzBGAiHgY8FlgnwEnM58WEb8AfIJymfmUfQ6bjI7u\ndZBo3qZGcLZt28nIyNaWnHPJ0CK2bd/dsvO1w/DwqgXV37pYh8I6FNahsA6FdSi6sQ6zBa5m1+Ds\nmgo3AJm5Htjr4uCIOCEijq7e/z1KmNoaEcurtxwJrG+y/VpMrcFp1VYNUC4V9z44kiR1VrMjONsi\n4jXA/6menw7sK8I9HXgE8HsRcRhwAPAl4GzKaM7Z1fOO2bMGZ7A1V1FBWWi8c9c4k5OTDAy0LjhJ\nkqTmNfvNfiHwaOAyyqXfx7DvG/39NXBoRHwNuBJ4FfBG4Pzq2MHV+TpmbGJqq4ZWjuAMMjE5ye5q\nI09JktR+zV5FtYEHr35qSmbuAH5zhpdOnct56vTgCE7rAk7jdg1LFte3QFqSJM2u2auofhP4I8rN\n/fakgcx8eE39aosH1+C0dooKynYNB7bsrJIkaS6aXYPzRspN+e6osS9t1+rNNqFhP6qd3gtHkqRO\naTbg3JKZ36i1Jx0wVsMU1TJ3FJckqeP2GnAi4lnVw+9HxJ8B11K2agAgM79cX9fqt2c38RYvMga3\na5AkqZP2NYLz+mnPT2x4PAks6IDz4G7irVyDU01ROYIjSVLH7CvgfA7418y8pR2dabfxFm+2CbB0\n8dQUlWtwJEnqlH0FnEXAeyPiEcA3gKso+0ltqb1nbVDPImPX4EiS1Gl7nZvJzPdl5pnA8cCngCcD\n10TENyLize3oYJ0eXINTw2XiBhxJkjqm2Rv97aKst/kyQEQcShfdsG9/7VmD09KrqEpJdxpwJEnq\nmGZv9Pc1yqLiRmMRcRLw1sy8s+U9a4PxWqeoXIMjSVKnNHsfnKuBx1AWHY8Dvw78FBgFPgqcVkvv\najZWw2XiS52ikiSp45oNOCdlZuOU1Bci4srMPDMizqqjY+0wPj7B4KKBlu76vec+OAYcSZI6ptnV\ntYdGxCFTTyLiIOCREbEaOKiWnrXB2MRkS0dvoPE+OE5RSZLUKc2O4PwVcFNE3A5MAD8HvA14LvA3\nNfWtduPjEy1dYAwwNDjA4KIBp6gkSeqgZgPOl4AjKTuJ/xLwLeCAzHx3XR1rh7HxyZZeIg4wMDDA\nsiWDPOBWDZIkdUyz3+7/AhwDLAa+A2ypHi9odYzgQFmH88BOA44kSZ3S7AjOxsy8oNaedMDYxCRD\nLV6DA7B0yRD33b+r5eeVJEnNaTbgfD4izgW+yUN3E/9pLb1qk7GxCYZaPEUFZQRnw6iLjCVJ6pRm\nA84TgXOBjQ3HJoGHt7xHbTQ+McHSxc2WoHlLFw8yNj7J2PgEQ4OtD1CSJGnvmv12fyqwJjN3zuXk\nEfEO4OSqnbcDNwAfBwaBu4Dz5nrOViqLjOtZgwPlZn8HLDfgSJLUbs1++94ALJvLiSPimcBxmXki\n8GzgvcAlwKWZeTJwK9DRdT3lRn91TFF5LxxJkjqp2RGco4DbIuImHroG5+l7+cx1wL9VjzcDK4F1\nwEXVsSuAi4EPzqG/LTU2PlnbVVTgdg2SJHVKswHnbXM9cWaOA/dXTy8Evgic3jAltQE4Ym/nWLNm\nBUNDg3NtumnjExMsWzrE8PCqlp734NXLAVi+cmnLz12XhdLPulmHwjoU1qGwDoV1KBZKHZoKOJn5\n1f1toNqr6kLKhpy3NLy0z6GT0dHt+9vsPk1MTDI5CRPjE4yMbG3tucfKyM3d92xl7Yruv13Q8PCq\nltdgIbIOhXUorENhHQrrUHRjHWYLXLWugI2I04HXAWdk5hZgW0Qsr14+ElhfZ/t7Mz4xAcBgDVc5\nLXUNjiRJHVVbwKk25Hwn8NzM3FQdvho4u3p8NmULiI4YG58EqP0qKkmS1H6tvwnMg14EHAJ8OiKm\njp0PfDgiXgHcDlxWY/t7NT5hwJEkqVfVFnAy80PAh2Z46dS62pyLPQGnhimqqYCz0w03JUnqiL69\nC934eFmDU8deVN4HR5KkzurfgFPjFNXSxdUUlTuKS5LUEQacOm70t7QKOE5RSZLUEf0bcKopqnq3\najDgSJLUCf0bcOq8impqiso1OJIkdUSdl4l3tamAM1TDVVRLFi9iYAB+fOcW/uzjN7b8/K22ePEg\nu51Osw4V61BYh8I6FNah2J86POv4I3nqEw6vqUez69uAMzY1RVXDGpyBgQEe/8iDufn2UX5y130t\nP78kSQvFXRvXdKTdvg044zXeyRjgNS/6hVrOW4du3FukE6xDYR0K61BYh8I6FAupDq7BqSngSJKk\nzunjgFPfZpuSJKmz+vbbve4pKkmS1Dn9G3CcopIkqWf1bcAZc4pKkqSe1bff7k5RSZLUu/o34DhF\nJUlSzzLg1HCjP0mS1Fn9G3CqOxkP1bDZpiRJ6qy+/XZ3ikqSpN5lwHGKSpKknlPrXlQRcRzwBeA9\nmfmBiDga+DgwCNwFnJeZO+vsw2ympqgGnaKSJKnn1PbtHhErgfcD1zQcvgS4NDNPBm4FLqir/X1x\nikqSpN5V5/DFTuA5wPqGY+uAy6vHVwCn1Nj+XjlFJUlS76ptiiozx4CxiGg8vLJhSmoDcMTezrFm\nzQqGhgZr6d/iJeVHP2TtAQwPr6qljYXEGhTWobAOhXUorENhHYqFUoda1+Dswz6HTkZHt9fW+LZt\nJWdtvW8HIyOdLEPnDQ+vYmRka6e70XHWobAOhXUorENhHYpurMNsgavdK2y3RcTy6vGRPHT6qq1c\ngyNJUu9qd8C5Gji7enw28KU2t7/HuJttSpLUs2qbm4mIE4B3AY8EdkfEC4BzgY9FxCuA24HL6mp/\nX9xsU5Kk3lXnIuMbKVdNTXdqXW3OhVNUkiT1rr6dnxnbc5l435ZAkqSe1bff7g/eydgRHEmSek3f\nXh/9c0cexJbtu1mxtG9LIElSz+rbb/fnPPURnP+847ruen5JkjR/fTtFJUmSepcBR5Ik9RwDjiRJ\n6jkGHEmS1HMMOJIkqecYcCRJUs8x4EiSpJ5jwJEkST3HgCNJknqOAUeSJPUcA44kSeo5BhxJktRz\nDDiSJKnnGHAkSVLPMeBIkqSeM9TuBiPiPcBTgUngdzPzhnb3QZIk9ba2juBExDOAR2fmicCFwPva\n2b4kSeoP7Z6i+hXgnwAy8yZgTUQc2OY+SJKkHtfuKarDgRsbno9Ux+6b6c3Dw6sG6u7Q8PCquptY\nEKxDYR0K61BYh8I6FNahWCh16PQi49oDjCRJ6j/tDjjrKSM2Ux4G3NXmPkiSpB7X7oBzFfACgIg4\nHlifmVvb3AdJktTjBiYnJ9vaYET8OfB0YAJ4VWb+e1s7IEmSel7bA44kSVLdOr3IWJIkqeUMOJIk\nqee0fauGbtHPW0ZExHHAF4D3ZOYHIuJo4OPAIOWqtvMyc2cn+9gOEfEO4GTKn4O3AzfQZ3WIiBXA\nx4DDgGXAW4B/p8/qMCUilgM/oNThGvqsDhGxDvgM8MPq0H8A76DP6gAQEecCrwXGgDcA36fP6hAR\nFwLnNRz6ReCXgQ9Svju/n5mv7ETfmtGXIzj9vGVERKwE3k/5y3vKJcClmXkycCtwQSf61k4R8Uzg\nuOr3wLOB99KHdQCeB3w7M58BvBB4N/1Zhyl/CmyqHvdrHb6ameuq/36HPqxDRKwF3gicBDwXOIs+\nrENm/t3U7wVKPS6j/F35u5n5y8BBEXFGJ/u4N30ZcOjvLSN2As+h3JNoyjrg8urxFcApbe5TJ1wH\nnFM93gyspA/rkJmfysx3VE+PBu6gD+sAEBGPBR4PXFkdWkcf1mEG6+i/OpwCXJ2ZWzPzrsx8Of1Z\nh0ZvAP4COKZhxqOr69CvU1Rz2jKil2TmGDAWEY2HVzYMtW4Ajmh7x9osM8eB+6unFwJfBE7vtzpM\niYjrgaMo/1q9uk/r8C7g1cD51fO++3NReXxEXA4cDLyZ/qzDI4EVVR3WAG+iP+sAQEQ8GfhvynTd\naMNLXV2Hfh3Bmc4tIx7UV7WIiLMoAefV017qqzpk5tOAXwU+wUN/9r6oQ0T8FvDNzPzJLG/pizoA\nt1BCzVmUoPd3PPQfwv1ShwFgLfB84CXAR+nDPxcNXkZZqzddV9ehXwOOW0Y81LZqcSXAkTx0+qpn\nRcTpwOuAMzJzC31Yh4g4oVpkTmZ+j/JltrXf6gCcCZwVEd+i/GX+evrw90Nm3llNW05m5o+BuylT\n+H1VB+Ae4PrMHKvqsJX+/HMxZR1wPWW2Y23D8a6uQ78GHLeMeKirgbOrx2cDX+pgX9oiIg4C3gk8\nNzOnFpX2XR0odxV/DUBEHAYcQB/WITNflJlPzsynAh+mXEXVd3WIiHMj4uLq8eGUq+s+Sp/VgfId\n8ayIWFQtOO7LPxcAEfEwYFtm7srM3cDNEXFS9fLz6eI69O2djPt1y4iIOIGy1uCRwG7gTuBcyvDj\nMuB24KXVb+SeFREvp8yr/6jh8PmUL7d+qsNyyjTE0cByyvTEt4G/p4/q0Cgi3gTcBvwrfVaHiFgF\n/AOwGlhC+f3wXfqsDgAR8QrK9DXAWym3kejHOpwAvDUzz6iePx74G8oAyf/NzD/oZP/2pm8DjiRJ\n6l39OkUlSZJ6mAFHkiT1HAOOJEnqOQYcSZLUcww4kiSp5xhwJM1LRLy4+v/hEfGZGs6/KCJujoiB\nhmPrIuLrrW5LUu/o172oJLVARAxSNuH7RGbezYMbmLbSLwLfzUzvaSGpaQYcSfPxEeAREXEV8HLg\n65l5VER8DLgXeBzwBOCPgecBT6ze80qAiPgz4JcpNxn8KvDaGYLMaZQ7y043GBEfBJ4E7ATOzMxt\nEXEBcBGwnXLL/d/OzPsiYhJYnJljEfES4JTMfHFE3AZ8CjgWeCnlRndrgMXAFZn5tnnWSFIHOEUl\naT7eCIxk5mkzvHZYZp5JuWP0pcCrgKcAL4mI1RFxDnBkZj4jM58CPIqym/l0swWcxwFvqrZX2A2c\nHhEPp9x991cycx1lB+Tfb+LnuCUzzwFOpYSgk4GnUfaj8u9JaQFyBEdSXb5R/f8O4KbM3AwQERuB\ng4BnAidGxLXV+w4Cjmk8QUQcCByYmXfOcP6bM/OehjZWA8cDNzbsLXctZTRnX65v6PMlEfFp4IvA\nhzNzoonPS+oyBhxJdRmb5THAAGVa6UOZ+Zd7OcezgK80cf6pc06f3prpGJR9lhrtAsjMDRHx88CJ\nwFnAtyPi+MzcsZc+SupCDr1Kmo8JylqV/fF14PkRMQQQEW+IiEdPe89s01OzuRE4odo0EuAU4FvV\n4/som4pCGT36GRFxGmUtzzcy87XANuDQObQvqUsYcCTNx3rg7oi4EVg5x8/+I2VK6PqI+CZwGPBf\n097zDMri46Zk5h3A64GrI+I6YBh4b/XynwNXRcQXKbuFz3gK4DUR8bVq6uyqzLy92fYldQ93E5ck\nST3HERxJktRzDDiSJKnnGHAkSVLPMeBIkqSeY8CRJEk9x4AjSZJ6jgFHkiT1nP8fu5AqFEEGoNcA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "tags": []
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# parameters\n",
    "t_final = 72     # hours\n",
    "V = 5.0          # liters\n",
    "k = 0.625        # liters/hour\n",
    "u_max = 50       # mg/hour infusion rate\n",
    "Csp = 25         # setpoint\n",
    "\n",
    "# create a model\n",
    "model = ConcreteModel()\n",
    "\n",
    "# define time as a continous set with bounds\n",
    "model.t = ContinuousSet(bounds=(0, t_final))\n",
    "\n",
    "# define problem variable indexed by time\n",
    "model.u = Var(model.t, bounds=(0, u_max))\n",
    "model.C = Var(model.t, domain=NonNegativeReals)\n",
    "model.dCdt = DerivativeVar(model.C, wrt=model.t)\n",
    "\n",
    "# the differential equation is a constraint indexed by time\n",
    "model.ode = Constraint(model.t, rule=lambda model, t: V*model.dCdt[t] == model.u[t] - k*model.C[t])\n",
    "\n",
    "# initial condition\n",
    "model.C[0].fix(0)\n",
    "\n",
    "# objective\n",
    "model.ise = Integral(model.t, wrt=model.t, rule=lambda model, t: (model.C[t] - Csp)**2)\n",
    "model.obj = Objective(expr=model.ise, sense=minimize)\n",
    "\n",
    "# transform the model to a system of algebraic constraints\n",
    "TransformationFactory('dae.finite_difference').apply_to(model, nfe=200, method='backwards')\n",
    "\n",
    "# simulation\n",
    "SolverFactory('ipopt', executable=ipopt_executable).solve(model)\n",
    "\n",
    "# extract solutions from the model\n",
    "tsim = [t for t in model.t]\n",
    "Csim = [model.C[t]() for t in model.t]\n",
    "usim = [model.u[t]() for t in model.t]\n",
    "\n",
    "# plot results\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(tsim, Csim)\n",
    "plt.title('Blood concentration of antibiotic')\n",
    "plt.xlabel('time / hours')\n",
    "plt.ylabel('mg/liter')\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(tsim, usim)\n",
    "plt.ylim((0, 1.1*u_max))\n",
    "plt.title('Infusion rate')\n",
    "plt.xlabel('time / hours')\n",
    "plt.ylabel('mg/hour')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "GUKn6xfIm-TJ"
   },
   "source": [
    "The result of this optimization is a much better therapy for the patient. Therapuetic levels of antibiotic are reach with 3 hours thereby leading to an earlier outcome and the possibility of an earlier recovery.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gbz3hEHdpcGo"
   },
   "source": [
    "## Alternative Objectives\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 441
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1415,
     "status": "ok",
     "timestamp": 1554236107978,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "vo_Y2R-hmwHt",
    "outputId": "16f74e27-33b6-460b-8d07-1333a449f11d"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGoCAYAAABL+58oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmcXFWd9/FPdVfv6XQ6SSeBJMgi\n/NgGB1AUZQnKIgIyw6IzAwiCCj46jzMj48zzckPU0dFRHIXHGYdREGccXEYBUWRAWQR5FETZfwZk\nkSSkO0mnl6TT3dVdzx/nVlJpe6nurlvb/b5fr7xy69ate3/166V+fc6556Sy2SwiIiIitaSu3AGI\niIiIFJsKHBEREak5KnBERESk5qjAERERkZqjAkdERERqjgocERERqTkqcERiZGZZM3vazJ4ys9+a\n2Y/MbN/ouTVm9nQRr3W+md1VrPOVgpktN7M3z/G178zbfsrMlhcvsimv+b/M7CUz+2CRzmdmdly0\n/adm9tVo+y4zO3+S41ea2WMFnPfVZnZYtP1eM/t4MeIVqSbpcgcgkgBr3P1FADP7FPDPwBnlDali\nnACcCNw8mxeZ2QrgA8C/Abj7gcUPbVJnAx90938v0vn+lPB7+B53/x7wvekOdvd1wKEFnPftwM+A\nR9z96nlHKVKFVOCIlNZPgD9osTCzZuALhA/8ceCHwAfcfSz6S/zLwBJgB/B37v5jM6sDvhid7yXg\n7qkuamZ/B1wKZIAfAO9396yZ/W/gMkJrrgPvcPceM7sOeB54LXAA8FvgTHffbmZHAl8B2oENwEXu\n/qyZHRzFuQcwDLzd3R80szXAp4C7gD8BmoGLgAHgaiBtZguAvwfuB24EjnD346PWnU8CjcAgcIm7\n/zo6bpWZPQUcFl1vtbu/OJf3VMjXInoPRwMHmdlqd79iwuveAbyf8Ht1A3CBuz9vZhcBpwH9wLHR\n1+BcYF/g/wAjZtYJPAqc7+4nRqf8IzP7RZTP26L3tBp42t3T0df/44SiC+AB4D3ABcDbgDeb2TJg\nIbDK3d8RtR5eB+wJ9AKXuvuvEKlB6qISKREzawTOZ/LWir8ifHgdAhxB+CD88+hD7L+Aq6NWincA\n3zSzduCNwMnAwcDxwHFTXPeY6HWvIPz1fwxwjpm9BvhbQgvTgcALhA/xnHOBtwL7AV2E1gaieD7k\n7gcQWhyujuL8PvD1aP9lwE1mlvsj6nDgAXc/CPi/0et/RShwvuPufxYdtxT4dVTcpIHrgXe6uwE3\nAf8UHXcx8IK7H+juI3nvda7vKd+kXwt3/wDwC0LhecWEHC+L3stJ7r4/8DTw4bxD3gT83yg3PwX+\nyt1vifL3z+7+/kniOAFYAxjh63v6hOffApwKHBnFugj4a3f/l7w4Pz/hNV8BvunuLycUjjdMcl2R\nmqACRyR+d0UtDRuBVwFfm+SY04CvuHvG3YeA/yAUL/sAKwhFBe7+IKEV4lWEguZWdx+MXvOtKa7/\npui4gagYWAP8d3TN77h7d3TctdE1c2519y3uniG0LuxlZgcAS939R9ExVxNaEA4ElgFfjeK8D+gh\ntJYADLj7TdH2r4C9poi1gaibJrruMnd/IHruXkKrx3Rm/Z6mOMdkX4spRddbmOuKnCTWJ9z9oWh7\nuvef7zvuvj1qYbqV0Ho0Mc7r3X2bu48Rvq+mjDNqmToB+Ga06ybg1QXEIVKV1EUlEr/8MTjHAXeb\n2RETjukidBnk9BIKhi5gq7tnJ3luMbB+wv7JLM0/LtclY2Zdk7x+Wd7jvrztMaA+OtfO/VGhkDGz\nRUAr8KSZ5Z5eSOhW653iXJMZc/f+vMf/28wuBJoIXVszLZ43l/c02Tkm+1pMyczqgSujLrV6Qvfd\nb2d53Yl6Jrx+j3nGuZjwR20fQPQ9NVhAHCJVSQWOSAm5+z1m9jyhmyj/A2wjoRjIWRLt2wgsNrNU\nXpGTe64X6Mh7TdcUl91EKEwAMLPcdaa65nQ2RfHUufu4mTUAKwlFRf9kg32jMTizZmavBf4OOMrd\nnzOzk4gGFU9jLu+pGOd4K2Es1HHuvim6w+u8WV53osV5253AlnnGuZlQIC4BNplZitBV98yEAlqk\nJqiLSqSEoi4eA56a8NQPgEvMrN7M2ggDRW8FngNeJHyA5j70VxDGWPwcOMXMWs2slTC+ZDI3Ewac\ndkbjWr4PnBKd/6y8gufSaN901kbxnBU9voQwruN54EUzOyeKc6mZfTN6L9MZJYwdmcwyoBt4IXp/\nFwJt0QfzKLAgb4xPzlze00RTfS2mswx4LipulhDGxywo4FrTvf+zzKw5iuFUQrfXxDjPj77+acLX\nIhfnH5zX3YeB2wkDvCF8D/xQxY3UKhU4IvG7K5qn5Sng24Q7Vx6dcMyXgN8DjwMPEj68vh19+PwZ\n8F4ze5Jw19S57r4NuAW4j3Cn0N2Eu33+QDSG5bPAr4EnCGNAvunuvwA+DdwbxbYImHZ+lyiec4EP\nmtla4C+Ad0+I8yngHuDOKM7p3A683sx+OclztxFahp6JjvsCoXvlO8AjhBaNl8xs53iWubynSUz6\ntZjhNd8ElkTzGn0T+BCw2sw+N8PrbgEuM7PvTPLcHYQByU9G27dNeP47hK/5Q8BjUcxfjJ77HvCP\nZjZxkPE7gDPM7HfAJwhfP5GalMpmVbyLiIhIbVELjoiIiNQcFTgiIiJSc1TgiIiISM1RgSMiIiI1\np6LnwenpGYh1BHRnZyu9vdtnPrDGKQ+B8hAoD4HyECgPgfIQVGIeurraU5PtT3QLTjpdyGSitU95\nCJSHQHkIlIdAeQiUh6Ca8pDoAkdERERqkwocERERqTkqcERERKTmqMARERGRmqMCR0RERIoum83y\nkHfzYvdgWa6vAkdERESKamR0jGt/8CTXfO8xbvvFC2WJoaLnwREREZHqkc1mWftiH9+8cy3PvzTA\nvnsu5Ozj9ytLLCpwREREZM627xjl2Q0DPLOuj1881c36TdsAOOaP9uCCUw6goUxz56jAKZOx8XEG\nhzIMDo2SyYyTGR9nbCxLZmyczFiWsbFxMuNZstldkznnNrPs3Nj134TnshPmgJ74OF97ey8DAzt2\nHcssJpCO59BZy073BguMo31BMwODO3bfGdf7K0K8RTj1pBYsaGJwcHiS88YUcwXkeLJzT5WHOZx6\nFjGU/2dv4te5ra2Jbdsmz8NszSZvFfHzlHdwa1sT26fJw+x/TuP5Joorx9lsluGRMbJ1KXr7djA0\nnGHb0Cib+3ewbUdm53H1dSmOOmgZJxy+kgNWLyKVmnSS4ZJQgROzrYPDPPHcFl7YOMiLPYNsHRyh\nf9sI24ZGY/3AFxERiVNTQz2LFzax754drF62gP1WLmT/VYtY0NJQ7tAAFTixGM9m+fljL3HvIxtY\n+/utuxUybc1p2lsb2XNJK+1tjSxoaaAhXUe6vo50fYp0XR3pdB3puhT19XXkit+dNXC0Y9fjXf9N\nrJT/8JjJK+n29ubdWnDyLlN0sz3vVDFPcfC8Dm1f2MJA/9A8z1v4wXH+YTObv5omHrlwYQv9k+Vh\n1jHM6uhYzjvrFOe9oGNhC33T5KGU35tTHluCHHd0tNDXN/33Q0xhVEaOo/87FrXSt3WGNZhm+UMd\n39d6Nuct/MTNjfWs3KODoW3DNDfWk66v7PuUVOAU2QsbB/j6j53fre8nBey/qoPDD+hinz0WsnrZ\nAlqaKi/lXV3t9PQMlDuMslMeAuUhUB4C5SFQHoKuJW30jI+XO4yCVN6nbRV78Klu/uWmxxnPZjnq\noGWcs2Y/lna0lDssERGRxFGBUyTPvdTPtT94goaGOt7zJ4dy6L5Lyh2SiIhIYlV2B1qV2Do4zJe+\n+yijmXEuPeMQFTciIiJlpgKnCG74sdM7MMw5a/bjj/dfWu5wREREEk8Fzjz1Dgzz66c38bIV7bzx\n1XuVOxwRERFBBc683ffoBrJZOP4Ve5Z1QiMRERHZRQXOPIxns/zskQ00pus46qDl5Q5HREREIipw\n5sFf2Er31iFeeeAyWpt1Q5qIiEilUIEzD/c+sh6AYw/bo8yRiIiISD4VOHM0PDLGQ97D8s4WDli9\nqNzhiIiISB4VOHP0QvcAo5lxDttvqQYXi4iIVJhYB46Y2WeAY6PrfAr4JXADUA9sAC5w96nXn69g\nL2wcBGCv5QvKHImIiIhMFFsLjpmdABzq7kcDbwS+AFwJXOPuxwJPAxfHdf24vbAxLLr2suXtZY5E\nREREJoqzi+oe4NxoeyvQBqwBbo723QKcGOP1Y/VC9yDp+jpWLGktdygiIiIyQWxdVO4+BmyLHl4C\n/BA4Ja9LqhuY9vajzs5W0un6uEIEoKtr9i0wmbFx1vVsY+892tljRUcMUZXeXPJQi5SHQHkIlIdA\neQiUh6Ba8hD75C1mdiahwDkZWJv31Iwjc3t7t8cVFhC+SD09A7N+3e+7B8mMjbPnktY5vb7SzDUP\ntUZ5CJSHQHkIlIdAeQgqMQ9TFVyx3kVlZqcAHwROdfc+YNDMWqKnVwLr47x+XHLjb1Yvq44qVkRE\nJGniHGTcAXwWON3dt0S77wDOjrbPBm6L6/pxyt1BpQHGIiIilSnOLqq3AkuBb5lZbt+FwLVmdinw\nPHB9jNePze+7B0gBq5a1lTsUERERmUScg4y/AnxlkqdOiuuapZDNZnlh4yDLFrfS3Kj1p0RERCqR\nZjKepU19O9g+nOFlmuBPRESkYqnAmaVdMxhr/I2IiEilUoEzS+s2hQJnVZdacERERCqVCpxZ6u4d\nAmDF4pYZjhQREZFyUYEzS929Q9TXpVjS0VzuUERERGQKKnBmqbt3O0s6mqmvU+pEREQqlT6lZ2Fo\nOEP/9lGWdap7SkREpJKpwJmF3Pib5Yu0griIiEglU4EzC91bQ4GjFhwREZHKpgJnFrqj1c1V4IiI\niFQ2FTizsLFXLTgiIiLVQAXOLHT3DpFKwdIOFTgiIiKVTAXOLHT3bmfJwmYa0kqbiIhIJdMndYGG\nR8bYOjii7ikREZEqoAKnQD0776DSLeIiIiKVTgVOgXYOMF6kFhwREZFKpwKnQN1bwy3iy9VFJSIi\nUvFU4BSoW7eIi4iIVA0VOAXKFThd6qISERGpeCpwCtTdu53O9iYaG+rLHYqIiIjMIB3nyc3sUOAm\n4Cp3v9rMrgOOBDZHh3zW3W+NM4ZiGM2MsaV/GNtrUblDERERkQLEVuCYWRvwJeDOCU/9H3f/QVzX\njUPP1h1k0fgbERGRahFnF9Uw8CZgfYzXKIldA4w1B46IiEg1iK0Fx90zQMbMJj71XjP7G6AbeK+7\nb4orhmLZuYq4BhiLiIhUhVjH4EziBmCzu//azP4euAJ471QHd3a2kk7HO6i3q6t9xmP6d2QAOHC/\npQUdX41q9X3NlvIQKA+B8hAoD4HyEFRLHkpa4Lh7/nicm4EvT3d8b9RyEpeurnZ6egZmPO75Df0A\npLPjBR1fbQrNQ61THgLlIVAeAuUhUB6CSszDVAVXSW8TN7Pvmtm+0cM1wGOlvP5cdfdup6OtkebG\nUjd4iYiIyFwU9IltZoe7+8OzObGZHQl8DtgbGDWzcwh3Vd1oZtuBQeDtswu39DJj42zq28HLV3aU\nOxQREREpUKFNEp8DXj+bE7v7Q4RWmom+O5vzlNvmvh1ks7pFXEREpJoUWuC8YGZ3AQ8AI7md7v6R\nOIKqJBt1i7iIiEjVKbTAeTb6lzi5W8S1iriIiEj1KKjAcfePmdkSYB93f9DM6tx9PObYKoJWERcR\nEak+Bd1FZWZ/Ruieui7a9SUzuziuoCpJ99aowNEkfyIiIlWj0NvE3w+8AuiJHl8OXBpLRBVmY+8Q\nC1oaaG1uKHcoIiIiUqBCC5w+d9856567D5E32LhWjY2Ps2nrkMbfiIiIVJlCBxlvMrMLgRYzOwJ4\nK7tac2rWlv5hxsazGn8jIiJSZQptwbkMeBXQDlwLtACXxBVUpciNv+nS+BsREZGqUmgLzhvdfbdF\nMc3sMuBfih9S5dikAkdERKQqTVvgmNnhwBHA5WaWP9NdA/ARar3A6dsBwNKO5jJHIiIiIrMxUwvO\nDmA5sAg4Nm//OPC3cQVVKXqiFpylHWrBERERqSbTFjju/iTwpJn9xN0fKFFMFWNT3w7q61J0tjeV\nOxQRERGZhZm6qP7Z3d8HfNbMshOfd/fjYousAmzaOsSShc3U1aXKHYqIiIjMwkxdVF+N/v9Q3IFU\nmuGRMfq3j7Jq2YJyhyIiIiKzNFOBs8TMXl+SSCrMpj6NvxEREalWMxU4H57muSzwkyLGUlF6ojuo\nuhbpDioREZFqM1OB813gx+6+thTBVJJNuoNKRESkas1U4NQBXzCzlwH3AbcDd7h7X+yRldnOOXDU\ngiMiIlJ1pl2qwd2/6O6nESb7u5GwXMOdZnafmX2sFAGWS24OnC614IiIiFSdgpZqcPcRwnibnwCY\n2TLgpBjjKrtNfTtobKijvbWh3KGIiIjILBVU4JjZvYRBxfkyZnYM8Al3X1f0yMpsU98OujpaSKU0\nB46IiEi1KXSxzTuAAwiDjseAPwVeAHqBrwEnT/YiMzsUuAm4yt2vNrPVwA1APbABuMDdh+f1DmKw\nbccoQ8MZlq7qKHcoIiIiMgfTjsHJc4y7n+fu/+3uN7n7RcCR7n4V0DjZC8ysDfgScGfe7iuBa9z9\nWOBp4OK5hx6fTVtzA4w1/kZERKQaFVrgLDOzpbkHZtYB7G1mi4CpmjmGgTcB6/P2rQFujrZvAU6c\nVbQlsmuRTd1BJSIiUo0K7aL6Z8Kim88TVhLfD/gkcDrwr5O9wN0zhHE6+bvb8rqkuoE9prtoZ2cr\n6XR9gSHOTVdX+x/sG3psIwD77dU56fO1KCnvcybKQ6A8BMpDoDwEykNQLXkotMC5DVgJpIBXAw8A\nC9z98/O49oyjd3t7t8/j9DPr6mqnp2fgD/Y/t34rAI0pJn2+1kyVh6RRHgLlIVAeAuUhUB6CSszD\nVAVXoV1UPwL2ARqAXwF90fZsDZpZbmDLSnbvvqoYO8fgaA4cERGRqlRoC85mdy/GgOA7gLOBb0T/\n31aEcxbdpr4h2prTtDYXmh4RERGpJIV+gn/PzM4Dfg5kcjvd/YWpXmBmRwKfA/YGRs3sHOA84Doz\nuxR4Hrh+jnHHJpvNsqlvB3suaSt3KCIiIjJHhRY4hxGKk815+7LAXlO9wN0fItw1NVFFz4Dct22E\n0cy41qASERGpYoUWOK8BOitxUr5iy42/0RpUIiIi1avQQca/BBLRpNHTF82BoxYcERGRqlVoC84q\n4Dkze5Ldx+AcF0tUZbRp5yR/asERERGpVoUWOJ+MNYoK0tMXdVGpBUdERKRqFVTguPvdcQdSKTZp\nmQYREZGqV+gYnMTY1LeDjgWNNMS8RISIiIjERwVOnrHxcbb0D+sOKhERkSqnAidPb/8w49msuqdE\nRESqnAqcPLkBxrpFXEREpLqpwMmjW8RFRERqgwqcPDtvEVcXlYiISFVTgZNn085ZjNWCIyIiUs1U\n4OTZtHUHdakUixc2lTsUERERmQcVOHl6+oZYvLCJ+jqlRUREpJrpkzwymhmnb3CEJQs1/kZERKTa\nqcCJ9A6EAcZLNMBYRESk6qnAiWzuHwZgsVpwREREqp4KnMiW/miSP7XgiIiIVD0VOJHNUYGjO6hE\nRESqnwqcyOZokj8NMhYREal+KnAiuS6qxe0qcERERKpdupQXM7M1wLeBx6Ndj7r7X5Yyhqls7h9m\nQUsDTY315Q5FRERE5qmkBU7kbnc/pwzXnVI2m2VL/w5WLGktdygiIiJSBOqiAgaHRhnJjGv8jYiI\nSI0oRwvOwWZ2M7AY+Ji7/89UB3Z2tpJOx9tl1NXVTt+LWwFYtWIhXV3tsV6vUiX1fU+kPATKQ6A8\nBMpDoDwE1ZKHUhc4a4GPAd8C9gV+amYvd/eRyQ7u7d0eazBdXe309Azw9HNbAGhJ19HTMxDrNStR\nLg9JpzwEykOgPATKQ6A8BJWYh6kKrpIWOO6+DrgxeviMmb0ErASeLWUcE+XuoNIyDSIiIrWhpGNw\nzOw8M7s82l4BLAfWlTKGyWiSPxERkdpS6i6qm4H/NLMzgUbg3VN1T5XSzmUaNMhYRESkJpS6i2oA\nOKOU1yzE5v5h0vUp2tsayx2KiIiIFIFuEye04Cxub6YulSp3KCIiIlIEiS9wRjNj9G0b0fgbERGR\nGpL4AmfLwDCgRTZFRERqiQqcvtwdVCpwREREakXiC5zN/VELjubAERERqRmJL3B2TvKnFhwREZGa\nkfgCR5P8iYiI1B4VOP0agyMiIlJrVOD0D7OgpYGmhnhXLRcREZHSSXSBk81m2dK/Q+NvREREakyi\nC5z+bSOMZsZ1B5WIiEiNSXSB09M7BGiAsYiISK1JdoGzdTugW8RFRERqTaILnO6oBUcFjoiISG1J\ndIGzq4tKBY6IiEgtSXaBs7OLSmNwREREakmyC5zeIdL1dbS3NZY7FBERESmiZBc4W4dYvLCJulSq\n3KGIiIhIESW2wBnNjLF1YFgDjEVERGpQYgucLf3DgObAERERqUWJLXByi2yqBUdERKT2pEt9QTO7\nCngNkAXe5+6/LHUMoAJHRESklpW0BcfMjgf2d/ejgUuAL5by+vl2dlFpHSoREZGaU+ouqjcA3wdw\n9yeBTjNbWOIYALXgiIiI1LJSd1GtAB7Ke9wT7euf7ODOzlbS6fp4Alm6gKUdzRy431IaYrpGNenq\nai93CBVBeQiUh0B5CJSHQHkIqiUPJR+DM8G0E9D09m6P7cKnvmoV57/xQLZs2RbbNapFV1c7PT0D\n5Q6j7JSHQHkIlIdAeQiUh6AS8zBVwVXqLqr1hBabnD2BDSWOAYBUKkV9fWJvIhMREalppf6Evx04\nB8DMjgDWu3tllYIiIiJS9Upa4Lj7/cBDZnY/4Q6q95Ty+iIiIpIMJR+D4+5/X+prioiISLJoEIqI\niIjUnFQ2my13DCIiIiJFpRYcERERqTkqcERERKTmqMARERGRmqMCR0RERGqOChwRERGpOSpwRERE\npOaowBEREZGaU+7VxMvGzK4CXgNkgfe5+y/LHFLJmNmhwE3AVe5+tZmtBm4A6gmLn17g7sPljLEU\nzOwzwLGEn4NPAb8kYXkws1bgOmA50Ax8HPgNCctDjpm1AI8R8nAnCcuDma0Bvg08Hu16FPgMCcsD\ngJmdB3wAyAAfAR4hYXkws0uAC/J2vRJ4HfBlwmfnI+7+7nLEVohEtuCY2fHA/u5+NHAJYV2sRDCz\nNuBLhF/eOVcC17j7scDTwMXliK2UzOwE4NDoe+CNwBdIYB6AM4AH3f144C3A50lmHnI+BGyJtpOa\nh7vdfU307y9JYB7MbAnwUeAY4HTgTBKYB3f/99z3AiEf1xN+V77P3V8HdJjZqeWMcTqJLHCANwDf\nB3D3J4FOM1tY3pBKZhh4E7A+b98a4OZo+xbgxBLHVA73AOdG21uBNhKYB3e/0d0/Ez1cDbxIAvMA\nYGYHAgcDt0a71pDAPExiDcnLw4nAHe4+4O4b3P1dJDMP+T4C/COwT16PR0XnIaldVCuAh/Ie90T7\n+ssTTum4ewbImFn+7ra8ptZuYI+SB1Zi7j4GbIseXgL8EDglaXnIMbP7gVWEv1bvSGgePge8F7gw\nepy4n4vIwWZ2M7AY+BjJzMPeQGuUh07gCpKZBwDM7FXA7wnddb15T1V0HpLagjNRqtwBVJBE5cLM\nziQUOO+d8FSi8uDurwXeDHyD3d97IvJgZm8Dfu7uz05xSCLyAKwlFDVnEgq9f2f3P4STkocUsAQ4\nC7gI+BoJ/LnI8w7CWL2JKjoPSS1w1hNabHL2JAwaS6rBaHAlwEp2776qWWZ2CvBB4FR37yOBeTCz\nI6NB5rj7rwkfZgNJywNwGnCmmT1A+GX+YRL4/eDu66Juy6y7PwO8ROjCT1QegI3A/e6eifIwQDJ/\nLnLWAPcTejuW5O2v6DwktcC5HTgHwMyOANa7+0B5QyqrO4Czo+2zgdvKGEtJmFkH8FngdHfPDSpN\nXB6A44D3A5jZcmABCcyDu7/V3V/l7q8BriXcRZW4PJjZeWZ2ebS9gnB33ddIWB4InxGvN7O6aMBx\nIn8uAMxsT2DQ3UfcfRR4ysyOiZ4+iwrOQyqbzZY7hrIws08TfrmPA+9x99+UOaSSMLMjCWMN9gZG\ngXXAeYTmx2bgeeDt0TdyzTKzdxH61X+bt/tCwodbkvLQQuiGWA20ELonHgS+ToLykM/MrgCeA35M\nwvJgZu3AfwKLgEbC98PDJCwPAGZ2KaH7GuAThGkkkpiHI4FPuPup0eODgX8lNJD8P3f/m3LGN53E\nFjgiIiJSu5LaRSUiIiI1TAWOiIiI1BwVOCIiIlJzVOCIiIhIzVGBIyIiIjVHBY6IzIuZnR/9v8LM\nvh3D+evM7CkzS+XtW2NmPyv2tUSkdiR1LSoRKQIzqycswvcNd3+JXQuYFtMrgYfdXXNaiEjBVOCI\nyHx8FXiZmd0OvAv4mbuvMrPrgE3AQcAhwN8DZwCHRce8G8DM/gF4HWGSwbuBD0xSyJxMmFl2onoz\n+zJwODAMnObug2Z2MXAZsJ0w5f473b3fzLJAg7tnzOwi4ER3P9/MngNuBPYF3k6Y6K4TaABucfdP\nzjNHIlIG6qISkfn4KNDj7idP8txydz+NMGP0NcB7gKOAi8xskZmdC6x09+Pd/Sjg5YTVzCeaqsA5\nCLgiWl5hFDjFzPYizL77BndfQ1gB+a8LeB9r3f1c4CRCEXQs8FrCelT6PSlShdSCIyJxuS/6/0Xg\nSXffCmBmm4EO4ATgaDO7KzquA9gn/wRmthBY6O7rJjn/U+6+Me8ai4AjgIfy1pa7i9CaM5P782K+\n0sy+BfwQuNbdxwt4vYhUGBU4IhKXzBTbAClCt9JX3P2fpjnH64GfFnD+3Dkndm9Ntg/COkv5RgDc\nvdvMXgEcDZwJPGhmR7j70DQxikgFUtOriMzHOGGsylz8DDjLzNIAZvYRM9t/wjFTdU9N5SHgyGjR\nSIATgQei7X7CoqIQWo/+gJmdTBjLc5+7fwAYBJbN4voiUiFU4IjIfKwHXjKzh4C2Wb72vwldQveb\n2c+B5cDvJhxzPGHwcUHc/UX990IAAAAehUlEQVTgw8AdZnYP0AV8IXr608DtZvZDwmrhk54CeL+Z\n3Rt1nd3u7s8Xen0RqRxaTVxERERqjlpwREREpOaowBEREZGaowJHREREao4KHBEREak5KnBERESk\n5qjAERERkZqjAkdERERqjgocERERqTkqcERERKTmqMARERGRmqMCR0RERGqOChwRERGpOSpwRKRg\nZvacmR1TwHH/YGbrzeztc7jGSjN7bG4Rzo+ZvdrMDivHtUWkuNLlDkBEatJbgQvc/c7ZvtDd1wGH\nFj+kgrwd+BnwSJmuLyJFogJHRObEzO4CbgbOAvYB7gH+AvgGsBfwVTP7BHAecK27fyPvdde6+zei\n588FUsCLwPlAI/C0u6fNrA74OHB2dNkHgPe4+7apru/u2QlxXgdsAU6MznUr8DXgj6NrfdfdLzez\ny4C3AW82s2XAVcCHo/ibge8Df+PuY8XIn4jES11UIjIfZwAnAQcArwde6+7nAeuA89z936Z6oZkd\nArwFONTdDwC+RyhC8r0FOBU4EjgEWAT89XTXn+JybwCOcvdvA+8G2oEDgSOAi8zsGHf/F+AXwAfc\n/fOEYustwFHAftG/d8+UEBGpDCpwRGQ+vuPuQ+6+DfgtoeWmUFuBLuA8M+t09y+5+9cnHHMacL27\nb4taTr4GnDyH69/p7jsA3P1zwJnunnX3XuBxYN9JXnMG8FV373P3DHAtobVIRKqAuqhEZD768rbH\ngPpCX+ju68zsLOBy4Etmdg9w2YTDuoDevMe9wLI5XH9LbsPM9gc+b2YHRq9ZTSicJloEXG5m74oe\np4Gead+UiFQMFTgiEreJhUdnbsPdfwr81MzagH8CPg18MO/YjcCSvMdLon3zcQ3wEPAn7j5mZvdN\ncdx64GZ3v3qe1xORMlAXlYjEbQPwCgAzO5owXgYzO9nMrjGzuqiL6TdAdsJrfwCcb2atZpYGLiEM\nEp6PZcDDUXFzErA/sCB6bpTQcgNwE3CBmbVG8V5qZhfO89oiUiJqwRGRuH0e+KaZnQrcBdwe7b8H\n+HPgt2Y2DHQTCph83wEOI7S4pICfAl+cZzyfAK4ys48Q7oz6GHClmT1MGOj8WTPbF3g/YWDzr8wM\n4JlJ4hORCpXKZif+wSQiIiJS3dRFJSIiIjVHBY6IiIjUHBU4IiIiUnNU4IiIiEjNqei7qHp6BmId\nAd3Z2Upv7/Y4L1EVlIdAeQiUh0B5CJSHQHkIKjEPXV3tqcn2J7oFJ50ueNLVmqY8BMpDoDwEykOg\nPATKQ1BNeUh0gSMiIiK1SQWOiIiI1BwVOCIiIlJzVOCIiIhIzVGBIyIiIjVHBY6IiIjUHBU4IiIi\nUnNU4IiIiEjNUYEjIiIiNSe2pRrMbA3wbeDxaNejwGeAG4B6YANwgbsPxxWDiIiIJFPcLTh3u/ua\n6N9fAlcC17j7scDTwMUxX19EREQSqNSLba4BLou2bwEuB75c4hh207N1iC98+zcMDWcKfk1rcwN/\nde5hLO1oiTEyERERmau4C5yDzexmYDHwMaAtr0uqG9hjuhd3drbGvrDX5sFRNmzeTseCRlqbG2Y8\nfmhHhvWbtrGxb5iDXr4s1thKqaurvdwhVATlIVAeAuUhUB4C5SGoljzEWeCsJRQ13wL2BX464XqT\nLm+eL+4l2bu62tncuw2As47bl2MP23PG1/x67Sa++N1HeH59Hz17LYo1vlLp6mqnp2eg3GGUnfIQ\nKA+B8hAoD4HyEFRiHqYquGIrcNx9HXBj9PAZM3sJeJWZtbj7ELASWB/X9Qs1khkHoCFd2HCkzvYm\nAHoHNDZaRESkUsU2yNjMzjOzy6PtFcBy4GvA2dEhZwO3xXX9Qo2OjgHQWGBX2KKowNmqAkdERKRi\nxdlFdTPwn2Z2JtAIvBt4GPi6mV0KPA9cH+P1CzI6FlpwGgtswWlvbaC+LkXvoAocERGRShVnF9UA\ncMYkT50U1zXnYmR0dl1UdakUixY0qYtKRESkgiV+JuPRnWNwCr9bq7O9ib7BEcbHs3GFJSIiIvOg\nAiczuy4qCONwxrNZ+rePxBWWiIiIzEPiC5yRTBhk3NBQeCo6F+hOKhERkUqW+AJnZxdV/SwKHN1J\nJSIiUtESX+Dk5sFpbCh8DM6i9kYA3UklIiJSoRJf4IzOcqI/UBeViIhIpVOBkxuDM5sCR11UIiIi\nFS3xBc5IZpx0fYq61IxLY+20KNeCoy4qERGRipT4Amc0Mz6rOXAgjNdpa06ri0pERKRCJb7AGcmM\nz6p7KqezvYmtasERERGpSIkvcEYzY7Oa5C9nUXsTQ8Nj7BjJxBCViIiIzIcKnLm24OhOKhERkYqV\n+AJnPl1UoDupREREKlHiC5zR0XEaZznIGEIXFcAWFTgiIiIVJ9EFTmZsnPFsdl5dVBpoLCIiUnkS\nXeCMjIZJ/uYyyDjXRaUxOCIiIpUn4QXO7JdpyMlN9tc3OFLUmERERGT+El7g5JZpmP0YnNbmNADb\nh3WbuIiISKVJdoETrUPV2DD7NKTr62hM16nAERERqUDpOE9uZi3AY8DHgTuBG4B6YANwgbuXdQDL\nzi6q+rnVeS3NaYZU4IiIiFScuFtwPgRsibavBK5x92OBp4GLY772jHZ2Uc2hBQegtSnN9h0qcERE\nRCpNbAWOmR0IHAzcGu1aA9wcbd8CnBjXtQu1s4tqDmNwAFqaQgtONpstZlgiIiIyT3F2UX0OeC9w\nYfS4La9LqhvYY6YTdHa2kp5j8VGI5zdtB2BRRwtdXe2zfv2i9mZ+t76fjs42mhrii7MU5vL+a5Hy\nECgPgfIQKA+B8hBUSx5iKXDM7G3Az939WTOb7JBUIefp7d1e1LgmGo66qEaGR+npGZj16+ujd/HC\ni707bxuvRl1d7XN6/7VGeQiUh0B5CJSHQHkIKjEPUxVccbXgnAbsa2anA6uAYWDQzFrcfQhYCayP\n6doFm89EfxC6qACGhjNVXeCIiIjUmlgKHHd/a27bzK4AngNeC5wNfCP6/7Y4rj0b85noDzQXjoiI\nSKUq5Tw4HwUuNLN7gcXA9SW89qR2teDMfZAxwJDupBIREakosc6DA+DuV+Q9PCnu683GrpmM536b\nOKgFR0REpNIkfCbjeXZR5Y3BERERkcqR7AKnWF1Uw2NFi0lERETmTwUOc2/BaWkKhdH24dGixSQi\nIiLzl+wCZ75dVM0NAAztUAuOiIhIJUl2gTPveXByLTgagyMiIlJJEl3gDO9cbHNuY3A0yFhERKQy\nJbrAGY0m+ptrC05zU5oUasERERGpNIkucOY7yLgulaK5qV4tOCIiIhUm0QXO8OgYqRTU1xW09uek\nWprSbNdMxiIiIhUl0QXOSGaMxnQ9qdTcC5zWprRacERERCpMsguc0fE5d0/ltDSlGRrJMJ7NFikq\nERERma+CPt3N7NS4AymHkdGxohQ42SwMj2guHBERkUpR6Kf735hZ7AtzltrI6Nic76DKaW3WreIi\nIiKVptCiZSvwhJn9ChjJ7XT3t8USVYmMZMZpa55f3daSt6L44mIEJSIiIvNW6Kf7D6J/NSV0Uc1t\nkr+c3GR/upNKRESkchRa4NwbaxRlMJ7NMpoZn38XlWYzFhERqTiFFjh3AlkgBTQCXcDjwOExxRW7\nzDwX2sxpUYEjIiJScQoqcNx9n/zHZnYIcEksEZXIfFcSz8kfgyMiIiKVYU6f7u7+OHBkkWMpqdGo\nwGmc40KbObqLSkREpPIU1IJjZldO2LUaWDTDa1qB64DlQDPwceA3wA1APbABuMDdh2cXcnGMZKJ1\nqOrVgiMiIlJrCv10H8v7lyEUKm+a4TVnAA+6+/HAW4DPA1cC17j7scDTwMVzCboYci04DQ3FGoOj\nif5EREQqRaFjcD5mZm2AEQYbu7tvn+E1N+Y9XA28CKwBLov23QJcDnx5ljEXxc4uqiLdRbV9x+i8\nYxIREZHiKLSL6k8IhcjvCa0+K8zsne7+owJeez+wCjgduCOvS6ob2GO613Z2tpKe5zw1U9nYH8JY\ntLCFrq72OZ+nfWELAGNZ5nWecqvm2ItJeQiUh0B5CJSHQHkIqiUPhd4m/rfAYe7eA2BmewLfAWYs\ncNz9tWb2x8A3CLeZ58y4hHdv77SNRPPSs3kQgJGRDD09A3M+TzabpS6Vom9geF7nKaeurvaqjb2Y\nlIdAeQiUh0B5CJSHoBLzMFXBVWj/zEiuuAFw9/XAtIODzexIM1sdHf9rQjE1YGYt0SErgfUFXr/o\nRkeL00WVSqVobU5rkLGIiEgFKbQFZ9DM3g/8T/T4FGCmEu444GXAX5nZcmABcBtwNqE15+zocVmM\nFGkMDkBLU71uExcREakghX66XwLsD1xPuPV7H2ae6O9fgGVmdi9wK/Ae4KPAhdG+xdH5yiI3yDhd\nlAJHLTgiIiKVpNC7qLrZdfdTQdx9CPiLSZ46aTbnictoNA9OYxEGMbc2pRkeGWNsfJz6uvkXTCIi\nIjI/hd5F9RfA3xEm99s5ONjd94oprtgVs4tqQWsjAAPbR1m0oGne5xMREZH5KXQMzkcJk/K9GGMs\nJVWstagAFreHomZL/7AKHBERkQpQaIGz1t3vizWSEhstYoGzZGEzAFv6d7DvngvnfT4RERGZn2kL\nHDN7fbT5iJn9A3AXYakGANz9J/GFFq+dY3DmudgmwOKowNncv2Pe5xIREZH5m6kF58MTHh+dt50F\nqrbA2dlFNc/FNgGWdIRuKRU4IiIilWGmAue7wI/dfW0pgimlYi22CbC4PbTg9PaXZWF0ERERmWCm\nAqcO+IKZvQy4D7idsJ5UX+yRxWzXYpvz76Jqb20gXV+nFhwREZEKMW3zhbt/0d1PA44AbgReBdxp\nZveZ2cdKEWBcijnIOJVKsWRhE1tU4IiIiFSEQif6GyGMt/kJgJkto0Im7Jurow9ZwYqlbbQ1F3oj\n2fQWL2zmyed7Gc2M0RDTCugiIiJSmEIn+ruXMKg4X8bMjgE+4e7rih5ZzI60Lt54zL5FWxV18cJo\nLpyBYZZ3thblnCIiIjI3hTZf3AEcQBh0PAb8KfAC0At8DTg5luiqyM65cPp2qMAREREps0ILnGPc\nPb9L6iYzu9XdTzOzM+MIrNrsmgtHd1KJiIiUW6EjbJeZ2dLcAzPrAPY2s0VARyyRVZn82YxFRESk\nvAptwfln4Ekzex4YB/YDPgmcDvxrTLFVldwYHN0qLiIiUn6FFji3ASsJK4m/GngAWODun48rsGqT\nm+xvy4C6qERERMqt0C6qHwH7AA3Ar4C+aFsiTY31LGhpUBeViIhIBSi0BWezu18cayQ1YPHCJl7a\nsp1sNksqlSp3OCIiIolVaIHzPTM7D/g5u68m/kIsUVWpxe3NvLBxkG07MixoUQOXiIhIuRRa4BwG\nnAdsztuXBfYqekRVLP9OKhU4IiIi5VNogfMaoNPdZzWC1sw+AxwbXedTwC+BG4B6YANwwWzPWckW\nd+y6k2qv5e1ljkZERCS5Ch1k/EugeTYnNrMTgEPd/WjgjcAXgCuBa9z9WOBpoKbG9SztaAHgpS3b\nyxyJiIhIshXagrMKeM7MnmT3MTjHTfOae4BfRNtbgTZgDXBZtO8W4HLgy7OIt6K9fGWY8/CJ53o5\n9dUvK3M0IiIiyVVogfPJ2Z7Y3ceAbdHDS4AfAqfkdUl1A3tMd47OzlbSMa/M3dVVvK6krq529t5j\nIWt/v5WFi1ppaqieVcWLmYdqpjwEykOgPATKQ6A8BNWSh4IKHHe/e64XiNaquoSwIOfavKdmvI+6\ntzferp6urvairSaeY6s7eG5DPz9/+EUO2WdxUc8dlzjyUI2Uh0B5CJSHQHkIlIdgtnlY++JWOhY0\nsWxRS6wxTabQMThzYmanAB8ETnX3PmDQzHLvciWwPs7rl8Oh+ywB4LFnN89wpIiISO3qGxzmM//5\nMN+565myXD+2AidakPOzwOnuviXafQdwdrR9NmEJiJqy/6oOGtJ1PP7slpkPFhERqVG/eLKbsfEs\nB6wqz5rchY7BmYu3AkuBb5lZbt+FwLVmdinwPHB9jNcvi8aGemz1Ih57dgu9A8N0tjeVOyQREZGS\ne+CJl6hLpTjqoOVluX5sBY67fwX4yiRPnRTXNSvFIfss5rFnt/DEc1t43R9NO45aRESk5mzcsp1n\nNwxw6L6LWdjWWJYYYh2Dk1SHRoOLH/2dxuGIiEjy/PzxlwA4+uAVZYtBBU4M9lzaxvLOFh7yHtZv\n2jbzC0RERGpENpvlgSc20thQx+EHLC1bHCpwYpBKpXjLCS9nbDzLf925lmw2W+6QRERESuKZ9f10\n9w5x+P5dNDfGOdR3eipwYvLH+y/l4L07eezZLfzmGXVViYhI7RseGeO6Hz0FwHGHlXcMqgqcmKRS\nKf78xAOoS6X4rzvXMjScmflFIiIiVew//ue3rN+0jTccuYqD9i7vZLcqcGK0cmkbJ75yFd29Q3z6\nP37F1sGaWThdRERkNz99eB0/e3QDL1vRzltOeHm5w1GBE7e3nPBy1hy+kt93D/LJrz/EM+v6yh2S\niIhI0WzfkeHfbnmCG37stDSlefeZh9CQLn95Ub7RPwlRV5figpMPoHNBI9+791k+ecNDHLLPYk45\najW2urMivglERERma9PWIX726Abu+c16tg6OsPeKdt55xsEs62wtd2iACpySSKVSnPG6fThg9SJu\nvu85Hn92C48/u4XGhjoOWL2I1V0L2GNJG50Lm1jQ3EBbS5oFLQ00NdSTSs24JqmIiEjRjWezjIyO\nMTg0yuDQKFv6hxl6opunfreJ323oZ8PmsCB2U2M9b37d3pz+2r1J11fOH+2pSr6FuadnINbgyrU6\n7NPr+njwqW4e/d3mnd8gk6mvS5FO15GuS1Ffl6K+vm7n/+m6FLtqn7CRXwtNeIrUro3dnk+lIJ2u\nJ5MZL8I7y5n/l61Y35azOU1Duo7RyfJQtFiKcKIS5CWdrivo+6F4vzqK8P1ShCgmnqg+XcfYHH4u\nihVLpfxurq8v7PthJsX7GpXn+6W+vo6xsQl5KMqPdHEyU4xvl2wWMmPjjGTGGc2Mk5n4fvM0N9az\n354LefXBK3jlgeW9Hbyrq33SlgAVOGUocPINbB9hw+btbNi8jb5tIwwOjbJtaJTBoQzbdoySGRtn\nbDzL2Fh21/Z4ducPWu7Lt3uispM/t/PxH74ot5nN7l4ozVVR2p2K1HiVKvBEqdQ0vySKFksRzlG0\nRr3JT1SXgvECf/KKFUolNVTmWk1TqVTZi4xK+Fmsq6tjfHy8KMHE/K0b6yl25mHieSrom3e+oaQI\nf/A21NfR2FBHQ30dTY31tDWnaWtpYHF7M/vu1UlbOsUeS9qoq6uM9z5VgaMuqjJrb22kvbWRA1Yv\nKlsMlVDoVQLlIVAeAuUhUB4C5SGopjxUTmeZiIiISJGowBEREZGaowJHREREao4KHBEREak5KnBE\nRESk5qjAERERkZqjAkdERERqTqzz4JjZocBNwFXufrWZrQZuAOqBDcAF7q4ltkVERKSoYmvBMbM2\n4EvAnXm7rwSucfdjgaeBi+O6voiIiCRXnF1Uw8CbgPV5+9YAN0fbtwAnxnh9ERERSajYuqjcPQNk\nzCx/d1tel1Q3sMd05+jsbCWdro8pwqCrqz3W81cL5SFQHgLlIVAeAuUhUB6CaslDOdeimnGVrt7e\nqVfaLoZqWlMjTspDoDwEykOgPATKQ6A8BJWYh6kKrlLfRTVoZi3R9kp2774SERERKYpSFzh3AGdH\n22cDt5X4+iIiIpIAsXVRmdmRwOeAvYFRMzsHOA+4zswuBZ4Hro/r+iIiIpJccQ4yfohw19REJ8V1\nTRERERHQTMYiIiJSg1TgiIiISM1RgSMiIiI1RwWOiIiI1BwVOCIiIlJzVOCIiIhIzVGBIyIiIjVH\nBY6IiIjUHBU4IiIiUnNU4IiIiEjNUYEjIiIiNUcFjoiIiNQcFTgiIiJSc1TgiIiISM1RgSMiIiI1\nRwWOiIiI1BwVOCIiIlJzVOCIiIhIzVGBIyIiIjUnXeoLmtlVwGuALPA+d/9lqWMQERGR2lbSFhwz\nOx7Y392PBi4BvljK64uIiEgylLqL6g3A9wHc/Umg08wWljgGERERqXGl7qJaATyU97gn2tc/2cFd\nXe2puAPq6mqP+xJVQXkIlIdAeQiUh0B5CJSHoFryUO5BxrEXMCIiIpI8pS5w1hNabHL2BDaUOAYR\nERGpcaUucG4HzgEwsyOA9e4+UOIYREREpMalstlsSS9oZp8GjgPGgfe4+29KGoCIiIjUvJIXOCIi\nIiJxK/cgYxEREZGiU4EjIiIiNafkSzVUiiQvGWFmhwI3AVe5+9Vmthq4Aagn3NV2gbsPlzPGUjCz\nzwDHEn4OPgX8koTlwcxageuA5UAz8HHgNyQsDzlm1gI8RsjDnSQsD2a2Bvg28Hi061HgMyQsDwBm\ndh7wASADfAR4hITlwcwuAS7I2/VK4HXAlwmfnY+4+7vLEVshEtmCk+QlI8ysDfgS4Zd3zpXANe5+\nLPA0cHE5YislMzsBODT6Hngj8AUSmAfgDOBBdz8eeAvweZKZh5wPAVui7aTm4W53XxP9+0sSmAcz\nWwJ8FDgGOB04kwTmwd3/Pfe9QMjH9YTfle9z99cBHWZ2ajljnE4iCxySvWTEMPAmwpxEOWuAm6Pt\nW4ATSxxTOdwDnBttbwXaSGAe3P1Gd/9M9HA18CIJzAOAmR0IHAzcGu1aQwLzMIk1JC8PJwJ3uPuA\nu29w93eRzDzk+wjwj8A+eT0eFZ2HpHZRzWrJiFri7hkgY2b5u9vymlq7gT1KHliJufsYsC16eAnw\nQ+CUpOUhx8zuB1YR/lq9I6F5+BzwXuDC6HHifi4iB5vZzcBi4GMkMw97A61RHjqBK0hmHgAws1cB\nvyd01/XmPVXReUhqC85EWjJil0TlwszOJBQ4753wVKLy4O6vBd4MfIPd33si8mBmbwN+7u7PTnFI\nIvIArCUUNWcSCr1/Z/c/hJOShxSwBDgLuAj4Ggn8ucjzDsJYvYkqOg9JLXC0ZMTuBqPBlQAr2b37\nqmaZ2SnAB4FT3b2PBObBzI6MBpnj7r8mfJgNJC0PwGnAmWb2AOGX+YdJ4PeDu6+Lui2z7v4M8BKh\nCz9ReQA2Ave7eybKwwDJ/LnIWQPcT+jtWJK3v6LzkNQCR0tG7O4O4Oxo+2zgtjLGUhJm1gF8Fjjd\n3XODShOXB8Ks4u8HMLPlwAISmAd3f6u7v8rdXwNcS7iLKnF5MLPzzOzyaHsF4e66r5GwPBA+I15v\nZnXRgONE/lwAmNmewKC7j7j7KPCUmR0TPX0WFZyHxM5knNQlI8zsSMJYg72BUWAdcB6h+bEZeB54\ne/SNXLPM7F2EfvXf5u2+kPDhlqQ8tBC6IVYDLYTuiQeBr5OgPOQzsyuA54Afk7A8mFk78J/AIqCR\n8P3wMAnLA4CZXUrovgb4BGEaiSTm4UjgE+5+avT4YOBfCQ0k/8/d/6ac8U0nsQWOiIiI1K6kdlGJ\niIhIDVOBIyIiIjVHBY6IiIjUHBU4IiIiUnNU4IiIiEjNUYEjIvNiZudH/68ws2/HcP46M3vKzFJ5\n+9aY2c+KfS0RqR1JXYtKRIrAzOoJi/B9w91fYtcCpsX0SuBhd9ecFiJSMBU4IjIfXwVeZma3A+8C\nfubuq8zsOmATcBBwCPD3wBnAYdEx7wYws38AXkeYZPBu4AOTFDInE2aWnajezL4MHA4MA6e5+6CZ\nXQxcBmwnTLn/TnfvN7Ms0ODuGTO7CDjR3c83s+eAG4F9gbcTJrrrBBqAW9z9k/PMkYiUgbqoRGQ+\nPgr0uPvJkzy33N1PI8wYfQ3wHuAo4CIzW2Rm5wIr3f14dz8KeDlhNfOJpipwDgKuiJZXGAVOMbO9\nCLPvvsHd1xBWQP7rAt7HWnc/FziJUAQdC7yWsB6Vfk+KVCG14IhIXO6L/n8ReNLdtwKY2WagAzgB\nONrM7oqO6wD2yT+BmS0EFrr7uknO/5S7b8y7xiLgCOChvLXl7iK05szk/ryYrzSzbwE/BK519/EC\nXi8iFUYFjojEJTPFNkCK0K30FXf/p2nO8XrgpwWcP3fOid1bk+2DsM5SvhEAd+82s1cARwNnAg+a\n2RH+/9u7e5QIgiiKwicwM1ZXIHcBLsBEXIBbEcw0NnQlhgaDkf/JLOAlA4KIsQhmYlAdjY7MOElT\nnC9q6OquCi+vHryqzz/OKGmELL1KWscXrVflP+6AoyQbAEnOkuzOrVl0PbXIFNgbhkYCHABPw/M7\nbagotOrRD0kOab0891V1AnwA2yvsL2kkDDiS1vEKvCWZApsrfntJuxJ6SPII7ACzuTX7tObjpVTV\nC3AKXCe5AbaAi+H1OTBJckWbFv7rL4DjJLfD1dmkqp6X3V/SeDhNXJIkdccKjiRJ6o4BR5IkdceA\nI0mSumPAkSRJ3THgSJKk7hhwJElSdww4kiSpO98EdS7AxOifbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "tags": []
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# parameters\n",
    "t_final = 72     # hours\n",
    "V = 5.0          # liters\n",
    "k = 0.625        # liters/hour\n",
    "u_max = 50       # mg/hour infusion rate\n",
    "C_mtd = 20         # setpoint\n",
    "\n",
    "# create a model\n",
    "model = ConcreteModel()\n",
    "\n",
    "# define time as a continous set with bounds\n",
    "model.t = ContinuousSet(bounds=(0, t_final))\n",
    "\n",
    "# define problem variable indexed by time\n",
    "model.err = Var(model.t, domain=NonNegativeReals)\n",
    "model.u = Var(model.t, bounds=(0, u_max))\n",
    "model.C = Var(model.t, bounds=(0, 25))\n",
    "model.dCdt = DerivativeVar(model.C, wrt=model.t)\n",
    "\n",
    "# the differential equation is a constraint indexed by time\n",
    "model.ode = Constraint(model.t, rule=lambda model, t: V*model.dCdt[t] == model.u[t] - k*model.C[t])\n",
    "\n",
    "# initial condition\n",
    "model.C[0].fix(0)\n",
    "\n",
    "# objective\n",
    "model.edef = Constraint(model.t, rule=lambda model, t: model.err[t] >= C_mtd - model.C[t])\n",
    "model.iae = Integral(model.t, wrt=model.t, rule=lambda model, t: model.err[t]**2)\n",
    "model.obj = Objective(expr=model.iae, sense=minimize)\n",
    "\n",
    "# transform the model to a system of algebraic constraints\n",
    "TransformationFactory('dae.finite_difference').apply_to(model, nfe=200, method='backwards')\n",
    "\n",
    "# simulation\n",
    "SolverFactory('ipopt', executable=ipopt_executable).solve(model)\n",
    "\n",
    "# extract solutions from the model\n",
    "tsim = [t for t in model.t]\n",
    "Csim = [model.C[t]() for t in model.t]\n",
    "usim = [model.u[t]() for t in model.t]\n",
    "\n",
    "# plot results\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(tsim, Csim)\n",
    "plt.title('Blood concentration of antibiotic')\n",
    "plt.xlabel('time / hours')\n",
    "plt.ylabel('mg/liter')\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(tsim, usim)\n",
    "plt.ylim((0, 1.1*u_max))\n",
    "plt.title('Infusion rate')\n",
    "plt.xlabel('time / hours')\n",
    "plt.ylabel('mg/hour')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "5k6wBKOhqEBF"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Simulation and Optimal Control](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/07.00-Simulation-and-Optimal-Control.ipynb) | [Contents](toc.ipynb) | [Soft Landing a Rocket](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/07.02-Soft-Landing-a-Rocket.ipynb) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE30338/blob/master/notebooks/07.01-Simulation-and-Optimal-Control-in-Pharmacokinetics.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://raw.githubusercontent.com/jckantor/CBE30338/master/notebooks/07.01-Simulation-and-Optimal-Control-in-Pharmacokinetics.ipynb\"><img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Simulation-and-Optimal-Control-in-Pharmacokinetics.ipynb",
   "provenance": [
    {
     "file_id": "1N-5hZtszVX-0zNIxxh_IQFRIbCGNBzaf",
     "timestamp": 1554226556061
    },
    {
     "file_id": "1TWUL-fRHE-VHji5VHLwKEf-5hY25zpE7",
     "timestamp": 1554216937956
    },
    {
     "file_id": "1dZEqItjD20N1qGb0kQbEFA4Ylc1k0QWq",
     "timestamp": 1554216470815
    }
   ],
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}