{
 "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",
    "< [Optimization](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/06.00-Optimization.ipynb) | [Contents](toc.ipynb) | [Linear Production Model](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/06.02-Linear-Production-Model.ipynb) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE30338/blob/master/notebooks/06.01-Unconstrained-Scalar-Optimization.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/06.01-Unconstrained-Scalar-Optimization.ipynb\"><img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unconstrained Scalar Optimization\n",
    "\n",
    "One of the problems studied in introductory calculus courses is the minimization or maximization of a function of a single variable. That is, given a function $f(x)$, find values $x^*$ such that $f(x^*) \\leq f(x)$, or $f(x^*) \\geq f(x)$, for all $x$ in an interval containing $x^*$. Such points are called local optima. If the derivative exists at all points in a given interval, then the local optima are found by solving for values $x^*$ that satisfy\n",
    "\n",
    "\\begin{align}\n",
    "f'(x^*) = 0\n",
    "\\end{align}\n",
    "\n",
    "Let's see how we can put this to work in a process engineering context."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Reactor for a Series Reaction\n",
    "\n",
    "A desired product $B$ is produced as intermediate in a series reaction\n",
    "\n",
    "\\begin{align}\n",
    "A \\overset{k_A}{\\longrightarrow} B \\overset{k_B}{\\longrightarrow} C\n",
    "\\end{align}\n",
    "\n",
    "where $A$ is a raw material and $C$ is a undesired by-product. The reaction operates at temperature where the rate constants are $k_A = 0.5\\ \\mbox{min}^{-1}$ and $k_A = 0.1\\ \\mbox{min}^{-1}$. The raw material is available as solution with a concenration $$C_{A,f} = 2.0\\ \\mbox{moles/liter}$.\n",
    "\n",
    "A 100 liter tank is avialable to run the reaction. \n",
    "\n",
    "1. If the goal is obtain the maximum possible concentration of $B$, and the tank is operated as a continuous stirred tank reactor, what should be the flowrate? \n",
    "\n",
    "2. What is the production rate of $B$ at maximum concentration?\n",
    "\n",
    "3. [Bonus] Would it be better to operate the tank as a batch reactor?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continuous Stirred Tank Reactor\n",
    "\n",
    "The reaction dynamics for an isothermal continuous stirred tank reactor with a volume $V = 40$ liters and feed concentration $C_{A,f}$ are modeled as\n",
    "\n",
    "\\begin{align}\n",
    "V\\frac{dC_A}{dt} & = q(C_{A,f} - C_A) - V k_A C_A \\\\\n",
    "V\\frac{dC_B}{dt} & = - q C_B + V k_A C_A - V k_B C_B\n",
    "\\end{align}\n",
    "\n",
    "At steady-state the material balances become\n",
    "\n",
    "\\begin{align}\n",
    "0 & = q(C_{A,f} - \\bar{C}_A) - V k_A \\bar{C}_A \\\\\n",
    "0 & = - q \\bar{C}_B + V k_A \\bar{C}_A - V k_B \\bar{C}_B \n",
    "\\end{align}\n",
    "\n",
    "which can be solved for $C_A$\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{C}_A & = \\frac{qC_{A,f}}{q + Vk_A} \\\\\n",
    "\\end{align}\n",
    "\n",
    "and then for $C_B$\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{C}_B & = \\frac{q V k_A C_{A,f}}{(q + V k_A)(q + Vk_B)}\n",
    "\\end{align}\n",
    "\n",
    "The numerator is first-order in flowrate $q$, and the denominator is quadratic.  This is consistent with an intermediate value of $q$ corresponding to a maximum concentration $\\bar{C}_B$. \n",
    "\n",
    "The next cell plots $\\bar{C}_B$ as a function of flowrate $q$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl4VdW5+PHvm5nMIwEJEGZFRZQo4BjqUPDW2vbaWq221rZ20Ku31f5qb3uttb2temuH21pb2zq2aq1Vi/OARJxQBgUEBEIIEMbMEELm9/fH3sGTmORsQnbO9H6e5zzss6fzrnPCec9ea+21RFUxxhhjusWFOgBjjDHhxRKDMcaYHiwxGGOM6cESgzHGmB4sMRhjjOnBEoMxxpgeLDEYz0SkWERURBJCHUu0EJG1IlLqw3kLRWSJiOwXkTuG+vwmulliiGIicoWIrBGRZhHZLSJ3iUj2YRxfKSLnHMFrvz6YY8PRkbwXAee4T0R+GrhOVY9V1bIjCq5vVwE1QKaqXu/D+QckjmtF5H0ROSAiVSLyDxE53t1eJCL/FJEaEWl0/06vEJEzRKTJfRxwf4g0BTzGiUiZiLS4z2tE5HERGT3cZYxmlhiilIhcD9wGfBfIAuYA44GXRCQplLFFozC8ihoPrNNB3ME6RGX5DXAdcC2QC0wFngT+zd3+ILDdjTMP+CKwR1VfU9V0VU0HjnX3ze5ep6rb3HXXuPtMBtKBXwxBzKabqtojyh5AJtAEfK7X+nRgL3Cl+/w+4KcB20uBKnf5QaALOOie6/8BxYACCe4+WcBfgF3ADuCnQDxwDNACdLrHNvQTZy5wL7ATqAeeDNj2NaAcqAMWAkcFbFPgG8Am97g7Ael17HpgP7AOOMldfxTwT6Aa2AJcG3DMzcCjwAPucWuBEg/vxVeAbcASd99/ALuBRmAJcKy7/iqgHWhzz/GUu74SOMddTgZ+7b4fO93l5MDPBrje/Qx3AV/u5329r9drnePx3N9zY3+wj3NOAl4BanGuRP6G84Xd1+tPcT/7Uwb4G20CZgb5O+5+jxN6rS8Dvhrw/FvA2lD/v4umR8gDsIcPHyrMBzp6/4dyt90PPOwu30c/icF9fuhLy33e4z8qzi/APwJpwEjgHeDr7rYrgNeDxPkM8HcgB0gEznLXf8z98jnJ/UL7bfcXr7tdgaeBbGAczhf9fHfbZ3GS1MmA4PyiHI9zdbwCuAlIAiYCFcDH3eNuxklm5+Mkt58DSz28Fw+45R/hrr8SyODDL+L3Ao7p8X73Pi9wC7DUfS8LgDeBnwR8Nh3uPolunM1ATj/vbe/P1su5b3PjHtHH+SYD57rbC3CS3q/7ee1vAFuDfPYvA28AnwfG9bNP93vcb2LAudp4GfhXqP/fRdMj5AHYw4cPFS4Ddvez7VbgJXe595dHKR4TA1AItAZ+iQCXAIvd5SsYIDEAo3F+hX/kiw3nKuT2gOfpOL+Ai93nCpwesP1R4EZ3+QXguj7OORvY1mvd94F73eWbgZcDtk0HDnp4LyYOUMZsd5+svt7v3ucFNgPnB2z7OFAZ8NkcDPySxLlymNPPa/f+bIOduw1IOYy/sU8B7/az7QcEJNV+9slx/xbX4lxdvAec3GufQ39vvdaX4STFRnf7e/STXOwxuIe1MUSnGiC/n7ri0e72IzUe55frLhFpEJEGnKuHkR6PHwvUqWp9H9uOArZ2P1HVJpwqjDEB++wOWG7GSR7d593cT7xHdcfqxvtfOAmuv3OmeKhv3969ICLxInKriGwWkX04X/oA+UHO0a1Hud3lowKe16pqR68Y0/Em2LmrVbWlv4NFZKSIPCIiO9yy/ZX+y1WL83fWL1WtV9UbVfVYnM/gPeBJEREPZQGnGjALmIGTZIo8Hmc8sMQQnd7C+TX/mcCVIpIGLAAWuasOAKkBu4zqdZ6BGi63u6+Rr6rZ7iPT/Y8e7Nju43P76SW1E+eLPDDuPJwqomC249SH97V+S0Cs2aqaoarnezgn9F+ewPWXAhfi1Oln4fziBadKa6BzdOtRbpxqsp0e4wsm2LmDxfZzd58ZqpqJc1Xa35f4IqBIREq8BKaqNTiNx0fhtDt5pqprcNq27jyMpGKCsMQQhVS1Efgx8FsRmS8iiSJSjNMwWoXTmArOr7TzRSRXREYB/9nrVHtw6uL7eo1dwIvAHSKSKSJxIjJJRM4KOLaovx5Q7vHPAb8XkRw3xjPdzQ8BXxaRmSKSDPwMeFtVKz0U/8/ADSIyy+0yOVlExuO0f+wTke+JyAj31/1xInKyh3N2l6fP9yJABk6yrMVJuD87zHM8DPxQRApEJB+nPeSvHuML5kjPnYHbkUBExuD0duuTqm4Cfg88LCKlIpIkIiki8nkRuRFARG5z3/8EEckAvgmUq2rtIMp2P86V6icHcazpgyWGKKWqt+NUlfwC2Ae8jfOr+WxVbXV3exBYhVPl8SJOQ3Cgn+N8mTSIyA19vMwXcRpy1+H0DnqMD6sQXsGpP94tIv1VXV2O03bwAU59+X+6sS8C/hunB9EunCuAz3ss9z+A/8FJLvtxGshzVbUTuACYidMjqQYniWR5OS/B3wtwGqK34lzZrMNp7A30F2C6e44n+zj+p8ByYDWwBljprhsKR3ruH+N0BmjE6TTweJD9rwV+h9NjrAGneu/TwFPu9lTgCXdbBc7VzKC+2FW1Dfg/nL8ZMwTEbcwxxhhjALtiMMYY04slBmOMMT34mhhE5B4R2Ssi7/ezXUTk/0SkXERWi8hJfsZjjDEmOL+vGO7DuQu3Pwtwbp+fgjNkwF0+x2OMMSYIXwf+UtUlbjfJ/lwIPKBOC/hSEckWkdFuV8Z+5efna3HxQKft34EDB0hLSxvUseHGyhJ+oqUcYGUJV0dSlhUrVtSoakGw/UI9IuQYAu4cxeljPwani2IPInIVzlUFhYWF/OIXgxtMsampifR0rzeLhjcrS/iJlnKAlSVcHUlZ5s2btzX4XqFPDH3dqdhn/1lVvRu4G6CkpERLS0sH9YJlZWUM9thwY2UJP9FSDrCyhKvhKEuoeyVV4Yxt062IoRsCwBhjzCCEOjEsBL7o9k6aAzQGa18wxhjjL1+rkkTkYZwhffNFpAr4Ec6InKjqH4BnccaVL8cZKfLLfsZjjDEmOL97JV0SZLsCV/sZgzHGmMMT6qokY4wxYcYSgzHGmB5C3V3VRKDOLqWhuY3aA23U7G+l5kAbza0dtHZ00drRSVtHF22dSrwICfFCQpyQEB9HSmIcGSmJZKQkkJmSQEZKIvnpyeSkJoa6SMaYAJYYTL/2t7SzZkcjm/c2sbn6AFtqDlBR08SO+oN0DeFo7UnxcWQkKuPWvUFhRgojM5MpyhnB+Lw0xuelMi43ldQk+1M1ZrjY/zZzSE1TK2+U1/BGeQ3vbmugvLqJ7uk6UpPimViQxsyxOVx4whgKMpLJS08iL835Nz05geSEOJIT40lOiCMhTuhS6OjqoqNT6ehUWjo62d/Szr6WDva3dLDvYDvV+1vZu7+V1Zu2EpeUwObqJt7cXMO+lo4esRVkJFOcl8qkgnSmFmZw9KgMpo7KID89OQTvlDHRzRJDjCvfu59nVu/mhbW7WbdrHwBZIxKZNT6HC044ihlFWRw9KpPCzGQOd0rdeIH4uHiS3b+yLBIpzEzpc9+yEbspLZ196Hljcztb6w6wtbaZrbXd/zbzwtrdPLLsw1FU8tKSmDYq41CyOG5MFtNGZZAYb81nxgyWJYYY1NDcxmMrqnhsRRUf7N6PCMwal8N3Pz6NM6bkc+xRWcTHhXZe9azURGakZjOjKLvHelWluqmVjbub+GD3Pjbu2c+GPU38fdl2DrZ3ApCUEMcxozOZMSaLGUVZzCjKZvLI9JCXyZhIYYkhhqzd2chfXtvC02t20dbRxcyx2dx8wXQWHD+631/y4UZEGJmRwsiMFE6fkn9ofVeXsq2umTU7Glmzo5FV2xt4fGUVDy51xgwbkRjP8WOyOGl8DiXjc5g1PoectKRQFcOYsGaJIQas2FrPnYvLeeWDvaQnJ3BxyVgunT2OY0Znhjq0IRMXJxTnp1Gcn8YFJxwFOMmiouYAa3Y0sGp7I+9tb+Avr1fwh1edhpNJBWnMGp9DyfhcZhXnMDE/7bCry4yJRpYYolhFdRM/e3Y9L6/fS05qIjecN5XL5xaTNSI2uofGxQmTR6YzeWQ6nz6xCICW9k5WbW9gxbZ6VlTW8+K6PTy6vAqAnNREZo3PZc7EXOZMzGP66EzirPrJxCBLDFGo8WA7v3l5Ew+8VUlKYjz/b/40rji12Lp8AimJ8cyemMfsiXnAh1cVK7bWsbyynmWVdby8fg/gNMKfMsFJEnMn5nH0qAxLFCYm2DdFlFn8wV6+//ga9u5v4eKTx/Gdc6dSkGFdOvsTeFVx8cnjANjVeJC3K+pYWlHLWxW1vLTOSRTZqYmcUpzL3El5zJmYx7RCSxQmOlliiBL7Wtq55al1PLaiiqmF6fzx8tM4YWx28APNR4zOGsGnThzDp04cA8DOhoO8vaWWtzbXsrSijhfdRJGblsSpk/I4Y0o+p08pYEz2iFCGbcyQscQQBdbt3Mc3/7aCqvqDXD1vEteePYXkhPhQhxU1jsoewadPLDrUTrGj4SBvba7lzc01vL6phqdXO1OITMxPY0JqK20Fu5kzKY/MlNhoyzHRxxJDhHtsRRU/eGIN2amJ/P2qOZQU54Y6pKg3JnsEF80q4qJZRagqm/Y28dqmGl7fVM3r5QdY9OAK4uOEE4qyOH1KAWdMyWfm2Gy76c5EDEsMEaqzS/nb+lZe2rqKORNz+e0lJ1lbQgiICFMLnTuvv3L6BF5+ZTHpxTN4o7yG1zbV8LtXNvF/izaRnpzAnIm5nDY5nzOnFljXWBPWgiYGcf56i1R1e7B9zfBoae/k+kdX8dLWDr58WjE/OP8YEuzXaFhIiBPmTHQap68/bxqNze28VeEkidfLa3h5/V4AinJGcNbUAs6aWsCpk/NJT7bfaCZ8BP1rVFUVkSeBWcMQjwliX0s7Vz2wnKUVdVw8LYkfXXBsqEMyA8hKTWT+caOZf9xoALbXNfPqxmpe3VjNk+/u4G9vbyMxXigZn8tZ05xEcfSoDLuaMCHl9WfKUhE5WVWX+RqNGdC+lnYu//PbrN25j19fPJPsxk2hDskcprG5qVw2ZzyXzRlPW0cXK7bWH0oUtz73Abc+9wEjM5Kdq4lpBZw+OZ/sVBu6wwwvr4lhHvANEakEDgCCczExw6/ATE8HWju48t5lrN25j7sum8W50wspK7PEEMmSEuKYOymPuZPyuHHB0ezZ18KSjdWUbazmxXV7+MeKKuIEZo7N5qypIymdVsDxY7Ls3gnjO6+JYYGvUZgBtbR38tX7l7NyWz2/u/Qkzp1eGOqQjA8KM1P4bMlYPlsylo7OLlZVNR66mvj1oo386uWN5KYlccaUfM6aWsAZUwqsw4HxhafEoKpbReR0YIqq3isiBUC6v6EZcHofXfPQuyzdUsuvPjeT848fHeqQzDBIiI9jljsK7HfOnUrdgTZe2+QkiSUbq/nXezsBOG5MptuIPZITx1mXWDM0PCUGEfkRUAJMA+4FEoG/Aqf5F5oBuP35D3h5/R5+/MljD92Ja2JPbloSF84cw4Uzx9DVpazbte/Q1cQfXq3gzsWbyUhO4LTJ+ZROK+DMqQUcZXdim0HyWpX0aeBEYCWAqu4UkQzfojIAPLp8O39cUsHlc8bzpVOLQx2OCRNxccJxY7I4bkwWV8+bzL6Wdt4sr+HVjdWUbajm+bW7AZhamE7ptJGcNbWAkuIcuxveeOY1MbS53VYVQETSfIzJAO9sqeMHT6zh9Mn53HTB9FCHY8JYZsqHXWK778R+dYNzNXHfG5XcvaSCEYnxnDopj9JpTrXTuLzUUIdtwpjXxPCoiPwRyBaRrwFXAn/2L6zYVtPUytUPrWRsTip3XnqS1RsbzwLvxP7amRNpbuvgrc21h64mFn2wF1jLhPy0Q11i50zIY0SSXU2YD3ltfP6FiJwL7MNpZ7hJVV/yNbIY1dWlXP/oKhoPtvPgV04hK9UGYjODl5qUwNnHFHL2MU5PtsqaA5Rt2MurG6t5ZNk27nuzkqSEOGZPyOWsqQWUThvJpAIbriPWeW18vk1Vvwe81Mc6M4TueWMLr26s5icXHsvRo6Jn6k0THorz07gifwJXnDaBlvZOllXWUeZWO/30mfX89Jn1jMkewVnTCshr66CktcOG64hBXj/xc4HeSWBBH+vMEVhT1chtz3/AedMLuWzO+FCHY6JcSmI8Z0xx7of4b6CqvpklG2so27CXf727gwNtndy16kVKinM4a6rTiH3MaBuuIxYMmBhE5JvAt4CJIrI6YFMG8IafgcWalvZOrvv7u+SnJ3P7RTPsP58ZdkU5qVw6exyXzh5HW0cX9yxcTGPqGMo2VHPb8x9w2/M9h+s4Y3KBVXVGqWBXDA8BzwE/B24MWL9fVet8iyoG/b5sMxXVB7j/ylNsbBwTckkJcRydG09p6dF8b/6Hw3W82mu4jhPH5RwaJdaG64gewRKDqmqliFzde4OI5FpyGBob9+znrrJyPjXzKM6aWhDqcIz5iH6H69iwl1+9vJFfvuQM13HmlHznamJKAfnpNlxHpPJyxfAJYAWgOIPndVNgYrAXEJH5wG+AeODPqnprr+3jgPuBbHefG1X1Wa8FiHRdXcr3H19DWnIC//0Ju1/BhL/ew3XUNrXyennNoXsnnnSH6zh+TJZ730QBM8dm25whEWTAxKCqn3D/nTCYk4tIPHAnTuN1FbBMRBaq6rqA3X4IPKqqd4nIdOBZoHgwrxeJHnpnGyu21vOLz55Anv3CMhEoLz25x3Ada3fu49WNeynbUM2di8v57SvlZKQkcMaUfEqnjuTMqQWMykoJddhmAMEan08aaLuqrgxy/lOAclWtcM/3CHAhEJgYFOjul5kF7AxyzqjR0NzG7c9/wGmT8/j3k2wcJBP54uKE44uyOL4oi2s+NoXG5nbe2Pzh1cSza5zhOiaPTOeMKfmcMSWf2RPySLMusWFFVLX/jSKLBzhWVfVjA55c5CJgvqp+1X1+OTBbVa8J2Gc08CKQA6QB56jqij7OdRVwFUBhYeGsRx55ZKCX7ldTUxPp6eExMOzD61t5cWsHPzltBEUZh3+ZHU5lOVLRUpZoKQcMfVlUlaom5f2aTtbWdLKhvpP2LogXmJwdx7H58RyXF09xVhxxQ9wrzz4Xx7x581aoakmw/YJVJc0b1Kt/qK9Pt3cmugS4T1XvEJG5wIMicpyqdvWK5W7gboCSkhItLS0dVEBlZWUM9tihtLX2AK+89CoXnzyWyy4Y3HxH4VKWoRAtZYmWcoD/ZWlp72TF1nqWbKrm9U01PL5pH49vaidrRCKnTc7j9MkFnDEln7G5Rz6uk30uh8fv67cqYGzA8yI+WlX0FWA+gKq+JSIpQD6w1+fYQur25zeQEBfHd86dGupQjAmJlMR4Tpucz2mT82GBM0bYG+U1vL6phtfLaw5VOxXnpXL6lHzOmFLA3El5ZKbYvRN+8zsxLAOmiMgEYAfweeDSXvtsA84G7hORY4AUoNrnuEJqxdZ6nlmzi+vOnsLITGuEMwYgP6ARW1XZXN3Ea5ucRPH4yh38dek24uOEE4qy3Du28zlhrE1O5IegiUGcW3CLVHX74Z5cVTtE5BrgBZyuqPeo6loRuQVYrqoLgeuBP4nIt3Gqma7QgRo+Ipyq8vNn11OQkcxVZwbt7WtMTBIRJo/MYPLIDL582gTaOrp4d1s9r5fXsGRTDb99ZRO/WbSJ9OQE5kzM48yp+Zw6Kd8GABwiQRODOw/Dk8CswbyAe0/Cs73W3RSwvI4Ymgnuzc21LN9az08uPNZ6YhjjUVJCHLMn5jF7Yh7XnzeNxuZ23tzsJInXy6t5ef0eAAozkzl1Uj5zJ+Vx2uR8xtgsdoPi9ZtpqYicrKrLfI0mBvzulXJGZiTz2ZKxwXc2xvQpKzWRBcePZoE7B/rW2gO8ubmWNzfX8tqmap54dwcA4/NSOXVSPtmtHRzX1Gp3Y3vkNTHMA74hIpXAAZzeRqqqg+tOE6NWbK3nrYpafvhvx5CSaBOjGDNUxuelMT4vjUtOGYeqsnFPE2+U1/Dm5lqeXrWT/a0d3LXqZY4eleFcTUzK55SJudaQ3Q+viWGBr1HEiDsXl5OTmsils8eFOhRjopaIMG1UBtNGZXDl6RPo6OzigacWczBrPG9truWht7dx7xuVxMcJx4/J4lS32mnW+Bz7webyOoPbVhE5HZiiqveKSAEQHXeLDJP3dzTyygd7ueG8qaQmWduCMcMlIT6OidnxlJZO5up5k2lp7+TdbQ28udm5ovjjkgp+X7aZpIQ4Zo3L4dRJeZw6OY8ZRbHb48nrDG4/AkpwpvW8F0gE/koMNRofqd+XlZORnMDlc4tDHYoxMS0lMZ65k/KYOymP64Gm1g6Wbanjzc01vFFeyx0vbeSOl2BEYjwlxTnMnpDLnIlOokhKiI1E4fWn66eBE4GVAKq6U0QyfIsqymyva+a593fzzbMmkTXC6jSNCSfpyQnMO3ok844eCUDdgTberqjl7S11LK2o5RcvbgQgJdEZVXb2hDzmTMzjhLFZJCdEZ9WT18TQ5nZbVQARSfMxpqjz16VbiRPh8rk2Xacx4S43LalHj6e6A228s6WOt7fUsrSijl+9vBFVSE6I46RxOcyZmMfsibnMHJsdNW0UXhPDoyLyRyBbRL4GXAn8yb+wokdLeyd/X76d86YXMjrL+lQbE2ly05KYf9wo5h83CnBGRX5nSx1LK5xk8etFG9GXnXstThybfShRnDQuchuzvTY+/0JEzgX2AVOBm1T1JV8jixILV+2kobmdL1rbgjFRITs1ifOOHcV5xzqJorG5nWWVTrXT21vq3LuyISk+jpljs5kzMZfZE/M4cVx2xHQ8OZwo1wAjcIatWONPONFFVXngrUqmFqYzZ2JuqMMxxvggKzWRc6YXcs70QgD2tbSzvNK5olhaUcvvFpfzf6+UkxAnHDsmi5PH53DyhFxOLs4lNy0853f32ivpq8BNwCs4N7f9VkRuUdV7/Awu0r27vYH3d+zjJ586zsZvMSZGZKYk8rGjC/nY0U6i2N/SzvKt9SyvrGPZlnoeWLqVP7++BYBJBWmc4iaJk4tzKcoZERbfFV6vGL4LnKiqtQAikge8CVhiGMCDb20lPTmBT59os7MZE6syUhKZN20k86Y5vZ5aOzpZU9XIO5V1LK+s55nVu3j4HWeM0lGZKe7VRA4nF+cyrTCDuLjhTxReE0MVsD/g+X7gsEdbjSX1B9p4ZvUuLjllLOk2WJ4xxpWcEE9JcS4lxU71cleXsmHPfpZX1vFOZT3LttTx1Cpn2prMlAR33xxOKc7l+KKsYYnR6zfWDuBtEfkXThvDhcA7IvIdAFX9pU/xRaynV++krbOLz51sg+UZY/oXFyccMzqTY0ZncvncYmcK1PqDLKuscx/1vPKBM29ZUkIcXz8+kVKfY/KaGDa7j27/cv+1m9z68fi7Ozh6VAbTR2eGOhRjTAQREcbmpjI2N5XPnFQEQG1T66F2iiJ2+x6D1+6qP/Y7kGhSUd3Eu9sa+P6Co8OiIckYE9ny0pP5+LGj+Pixoygr83/W49gY+GOYPfnuDuIEPmWNzsaYCGSJYYh1dSmPv7uD0ybnU2jzORtjIpAlhiG2fGs9VfUH+cxJdrVgjIlMnhKDiNwuIpkikigii0SkRkQu8zu4SPT4yipSk+L5uHu7vDHGRBqvVwznqeo+4BM49zRMxbnpzQRoae/kmTW7WHDc6IgZE8UYY3rzmhi6JxE4H3hYVet8iieivVFew/6WDj4586hQh2KMMYPm9WftUyLyAXAQ+JY7tWeLf2FFpuff301GSgJzJ+aFOhRjjBk0T1cMqnojMBcoUdV2oBnn7mfj6ujs4uX1ezj76JExM/2fMSY6eW18TgWuBu5yVx2FMwe0cS2rrKe+ud0anY0xEc/rT9t7gTbgVPd5FfBTXyKKUC+s3U1yQhxnTSsIdSjGGHNEvCaGSap6O9AOoKoHceZlMDgT8ry4djdnTCmw3kjGmIjnNTG0iUj37G2IyCSg1beoIsyaHY3sbGw5NCesMcZEMq8/b38EPA+MFZG/AacBV/gVVKR5Ye1u4uOEc44ZGepQjDHmiHkdXfUlEVkJzMGpQrpOVWt8jSyCvLB2D7Mn5JKdGp7ztxpjzOEYMDGIyEm9Vu1y/x0nIuNUdaU/YUWOLTUHKN/bxOVzxoc6FGOMGRLBrhjuGGCbAh8bwlgi0pKN1QCH5nM1xphIN2BiUNV5wxVIpFqysZrxeamMy0sNdSjGGDMkvN7gligi14rIY+7jGhFJDH4kiMh8EdkgIuUicmM/+3xORNaJyFoReehwChBKbR1dvFVRy5lT7N4FY0z08Nor6S6cgfR+7z6/3F331YEOEpF44E7gXJyb4paJyEJVXRewzxTg+8BpqlovIhFTJ7Niaz3NbZ2cMSU/1KEYY8yQ8ZoYTlbVEwKevyIiqzwcdwpQrqoVACLyCM4YS+sC9vkacKeq1gOoqv8Tmg6R1zZVkxAnzJ1kg+YZY6KH18TQKSKTVHUzgIhMBDo9HDcG2B7wvAqY3Wufqe453wDigZtV9fneJxKRq4CrAAoLCykrK/MYek9NTU2DPra3Z1YeZGKWsGLpG0NyvsM1lGUJtWgpS7SUA6ws4WpYyqKqQR/A2cA2oAx4FagE5nk47rPAnwOeXw78ttc+TwNP4FRVTcBJHtkDnXfWrFk6WIsXLx70sYFq9rfo+O89rb9dtHFIzjcYQ1WWcBAtZYmWcqhaWcLVkZQFWK4evvO93uC2yG0LmIZzg9sHquplSIwqYGzA8yJgZx/7LFVnOO8tIrIBmAIs8xJbqLxe7tzfd4Y1PBtjoozXXknxwMeBUpyrh6tF5DseDl0GTBGRCSKSBHweWNhrnyeBee7r5ONULVV4ij6ElmysISc1kePGZIU6FGOMGVKeZ3DDmbFtDdDl9eSq2iEi1wAv4LQf3KOqa0XkFpxLmoXutvNEZB1+lFBEAAAZL0lEQVROu8V3VbX2cAox3FSV1zZVc9rkfOLjbJBZY0x08ZoYilR1xmBeQFWfBZ7tte6mgGUFvuM+IsLGPU3s3d9q9y8YY6KS12G3nxOR83yNJIK8s8W5oLFuqsaYaOT1imEp8ISIxOFM1iM4P/YzfYssjL1TWc+ozBSKckaEOhRjjBlyXhPDHcBcYI1b9ROzVJVlW+o4eUIuIta+YIyJPl6rkjYB78d6UgCoqj/I7n0tnFycE+pQjDHGF16vGHYBZSLyHAFTeqrqL32JKoy9s6UOgJOLc0MciTHG+MNrYtjiPpLcR8xaVllHZkoC0wozQh2KMcb4wuudzz/2O5BI8U5lHSXFucTZ/QvGmCjltY3BADVNrVRUH7BqJGNMVLPEcBiWVzrtC6dMsIZnY0z0GjAxiMglImJ3cbne2VJPckIcx4/JDnUoxhjjm2BtDOOBf7jTeC4CngPeidVuq8sq65g5NpukBLvQMsZErwG/4VT1VlX9GHA+sAq4ElgpIg+JyBdFpHA4ggwHTa0drN3ZyCkTrH3BGBPdvPZK2o8zmc4TACIyHVgAPIAzHHfUe29bA10KJdbwbIyJcl7vY+hBVdfhzNt8x9CGE75WVTUAMLPI2heMMdHNKss9WlPVSHFeKlmpiaEOxRhjfGWJwaPVVQ0cb1cLxpgYcFiJQURSRaRERGJqhprq/a3sbGzhhCKbxtMYE/2C3cfwSRGpFJGVInI+sBb4HbBGRL40LBGGgTU7nPaF421+Z2NMDAjW+PwT4DwgC1gMzFDVChEZiXNfw/0+xxcWVlc1IgLHWWIwxsSAYImhS1U3AojIFlWtAFDVvSLS4Xt0YWJ1VSOTC9JJSx5UJy5jjIkowb7p4kQkB6fKqctd7h5WNCYarlWV1VWNnDU1pppVjDExLFhiyAJW8GEyWBmwLSaGxdjV2EJNUyszrOHZGBMjBkwMqlrs5SQicqyqrh2SiMLM6qpGAEsMxpiYMVTVQQ8O0XnCzuqqBhLihGNGZ4Y6FGOMGRZDlRiidjqzNTsamTYqg5TE+FCHYowxw2KoEkNUtjd0NzxbNZIxJpbERM+iwdpW10zjwXabmMcYE1OGKjG0DdF5wsqaHdbwbIyJPcGGxPi4iFzUx/oviMi53c9VdY4fwYXa+l37SIgTphSmhzoUY4wZNsGuGH4MvNrH+kXALUMfTnhZv2s/kwrSSU6whmdjTOwIlhhSVbW690pV3Q2k+RNS+Fi/ax/HjM4IdRjGGDOsgiWGFBH5yE1wIpIIjPAnpPDQ0NzGrsYWu3/BGBNzgiWGx4E/icihqwN3+Q/utqBEZL6IbBCRchG5cYD9LhIRFZESL+f12/pd+wE42hKDMSbGBEsMPwT2AFtFZIWIrAAqgWp324BEJB64E1gATAcuEZHpfeyXAVwLvH1Y0fto/a59AFaVZIyJOcHGSuoAbhSRHwOT3dXlqnrQ4/lPcfevABCRR4ALgXW99vsJcDtwg9fA/bZ+1z7y05MYmZES6lCMMWZYiap/Ny27XV3nq+pX3eeXA7NV9ZqAfU4Efqiq/y4iZcANqrq8j3NdBVwFUFhYOOuRRx4ZVExNTU2kpwfvfvqjNw+SngjfPTl8m1K8liUSREtZoqUcYGUJV0dSlnnz5q1Q1aDV9X7PPNPXGEqHMpGIxAG/Aq4IdiJVvRu4G6CkpERLS0sHFVBZWRnBju3sUna//DyXzxlPaelHar7ChpeyRIpoKUu0lAOsLOFqOMri95AYVcDYgOdFwM6A5xnAcUCZiFQCc4CFoW6A3l7XTGtHF1MLrX3BGBN7PF8xiMgYYHzgMaq6JMhhy4ApIjIB2AF8Hrg04PhGID/gNcropyppOG3c4/RIsjuejTGxyFNiEJHbgItxGo073dUKDJgYVLVDRK4BXgDigXtUda2I3AIsV9WFg47cRx8mBrtiMMbEHq9XDJ8Cpqlq6+G+gKo+Czzba91N/exberjn98PGPU2MyR5BerLfTTDGGBN+vLYxVACJfgYSTjbu2c9Uq0YyxsQorz+Jm4H3RGQRcOiqQVWv9SWqEOro7KKi+gBnTSsIdSjGGBMSXhPDQvcR9Sprm2nr7GLqSGtfMMbEJk+JQVXvF5EkYKq7aoOqtvsXVuhschuerauqMSZWee2VVArcjzNOkgBjReRLHrqrRpzyvU0ATBoZ9aOKG2NMn7xWJd0BnKeqGwBEZCrwMDDLr8BCZXO10yMpNcl6JBljYpPXXkmJ3UkBQFU3EqW9lDZXH2BigV0tGGNil9fEsFxE/iIipe7jT8AKPwMLBVVlc3UTkwqsq6oxJnZ5rS/5JnA1zpwJgnPH8+/9CipUdu9robmtk0kjLTEYY2KX115JrcAv3UfU2rz3AACTrCrJGBPDBkwMIvKoqn5ORNYQMFx2N1Wd4VtkIbC52umRNNmqkowxMSzYFcN17r+f8DuQcLC5uomM5AQKMpJDHYoxxoTMgI3PqrrLXfyWqm4NfADf8j+84bW5uomJI9MR6Wt+IWOMiQ1eeyWd28e6BUMZSDjYvPcAk/KtfcEYE9uCtTF8E+fKYKKIrA7YlAG84Wdgw+1gWye797UwwRKDMSbGBWtjeAh4Dvg5cGPA+v2qWudbVCFQWev0SCq2xGCMiXEDJgZ36s1G4BIAERkJpADpIpKuqtv8D3F4VNY4icGuGIwxsc5TG4OIXCAim4AtwKs4g+k952Ncw26LXTEYYwzgvfH5p8AcYKOqTgDOJsraGCprDpCfnmzTeRpjYp7XxNCuqrVAnIjEqepiYKaPcQ27yppmJuSnhjoMY4wJOa8/jxtEJB1njKS/icheoMO/sIbfltoDlE616TyNMcbrFcOFOPM+fxt4HtgMXOBXUMOtqbWD6v2t1r5gjDF4uGIQkXjgX6p6DtCFM5NbVLEeScYY86GgVwyq2gk0i0jWMMQTEltrmwEozrPEYIwxXtsYWoA1IvIScKB7pape60tUw6z75rbxedb4bIwxXhPDM+4j0EeG4Y5U2+uayU9PJs26qhpjjOfEkK2qvwlcISLX9bdzpNlW18y43BGhDsMYY8KC115JX+pj3RVDGEdIOYnBqpGMMQaCj656CXApMEFEFgZsygBq/QxsuLR3drGz4SDjThwT6lCMMSYsBKtKehPYBeQDdwSs3w+s7vOICLOz4SBdCmPtisEYY4Dgo6tuBbYCc4cnnOG3rc7pqmpVScYY4/A6uupnRGSTiDSKyD4R2S8i+/wObjgcSgzWVdUYYwDvjc+3A59U1SxVzVTVDFXN9HKgiMwXkQ0iUi4iN/ax/Tsisk5EVovIIhEZfzgFOFLb6ppJio+jMCNlOF/WGGPCltfEsEdV1x/uyd3hNO7EmR96OnCJiEzvtdu7QImqzgAew0lCw2Z7XTNFOSOIi5PhfFljjAlbXu9jWC4ifweeBFq7V6rq40GOOwUoV9UKABF5BGdAvnUB51gcsP9S4DKPMQ2JbXXN1vBsjDEBvCaGTJzRVc8LWKdAsMQwBtge8LwKmD3A/l9hmGeG21bbzMyx2cP5ksYYE9Y8JQZV/fIgz99X/UyfQ2mIyGVACXBWP9uvAq4CKCwspKysbFABNTU1HTq2uV3Z19JBW/0uysoi77aMwLJEumgpS7SUA6ws4WpYyqKqQR/AVGAR8L77fAbwQw/HzQVeCHj+feD7fex3DrAeGOklnlmzZulgLV68+NDyup2NOv57T+tTq3YM+nyhFFiWSBctZYmWcqhaWcLVkZQFWK4evmO9Nj7/yf1Sb3eTyWrg8x6OWwZMEZEJIpLkHhN4BzUiciLwR5xeT3s9xjMkquoPAlCUY20MxhjTzWtiSFXVd3qtCzq1p6p2ANcAL+BcETyqqmtF5BYR+aS72/8C6cA/ROS9XkNv+GpHvXMPw5hsG0DPGGO6eW18rhGRSbjtAyJyEc5QGUGp6rPAs73W3RSwfI7HGIbcjoaDJCfEkZ+eFKoQjDEm7HhNDFcDdwNHi8gOYAvD3K3UDzsaDjImewQidg+DMcZ089orqQI4R0TSgDhV3e9vWMNjR/1BxuRYNZIxxgTyOlbSz0QkW1UPqOp+EckRkZ/6HZzfquoPUmSJwRhjevDa+LxAVRu6n6hqPXC+PyENj4NtndQeaLOGZ2OM6cVrYogXkeTuJyIyAkgeYP+wt6PB6apqVUnGGNOT18bnvwKLRORenJ5JVwL3+xbVMDiUGLLtHgZjjAnktfH5dhFZA5yNM8zFT1T1BV8j81mVew+DtTEYY0xPXq8YUNXnGOYB7vy0s+Eg8XFCYabNw2CMMYFidga3XQ0tFGYkE2/zMBhjTA9erxhuBy7QQUzWE652NbYw2nokGWPMR/g6g1s429V4kNFZVo1kjDG9+T2DW1hSVXY1tnDu9MJQh2KMMWHH7xncwlJ9czutHV2MzrKqJGOM6c3vGdzC0k73HgarSjLGmI/y2iupSESeEJG9IrJHRP4pIkV+B+eX3Y0tANb4bIwxffDa+HwvzsxrRwFjgKfcdRFpV6NzxXCUXTEYY8xHeE0MBap6r6p2uI/7gAIf4/LVzsYWEuKEvPSIHu7JGGN84TUx1IjIZSIS7z4uA2r9DMxPuxtbKMxMsZvbjDGmD14Tw5XA54DdOFN6XuSui0g7Gw5yVLZVIxljTF+89kraBnzS51iGze59Lcwoyg51GMYYE5a89kq6X0SyA57niMg9/oXln+6b26zh2Rhj+ua1KmlGHzO4nehPSP5qaoe2ji4bVdUYY/rhNTHEiUhO9xMRyeUwhuwOJw2tCsAou2Iwxpg+ef1yvwN4U0QewxkK43PA//gWlY/qW7oAKMy0rqrGGNMXr43PD4jIcuBjODO4fUZV1/kamU/q3SsGq0oyxpi+Hc4MbuuAiEwGgRpanMRQkGFXDMYY0xevbQxRo75VyU1LIjkhPtShGGNMWIq5xNDQooy0qwVjjOlX7CWGVrUeScYYM4CYTAyFGZYYjDGmPzGVGDo6u2hsVeuqaowxA4ipxFDT1IYChVaVZIwx/YqpxLBnnzNzm1UlGWNM/3xPDCIyX0Q2iEi5iNzYx/ZkEfm7u/1tESn2K5bd3YnBbm4zxph++ZoYRCQeuBNYAEwHLhGR6b12+wpQr6qTgV8Bt/kVz95DicHaGIwxpj9+XzGcApSraoWqtgGPABf22udC4H53+THgbBHxZWq1PftaiRNsSk9jjBmA3yOkjgG2BzyvAmb3t4+qdohII5AH1ATuJCJXAVcBFBYWUlZWdtjBtNS0U1KgvLbk1cM+Nhw1NTUN6n0IR9FSlmgpB1hZwtVwlMXvxNDXL38dxD6o6t3A3QAlJSVaWlp62MGUAmVlZQzm2HBkZQk/0VIOsLKEq+Eoi99VSVXA2IDnRcDO/vYRkQQgC6jzOS5jjDH98DsxLAOmiMgEEUkCPg8s7LXPQuBL7vJFwCuq+pErBmOMMcPD16okt83gGuAFIB64R1XXisgtwHJVXQj8BXhQRMpxrhQ+72dMxhhjBub79Jyq+izwbK91NwUstwCf9TsOY4wx3sTUnc/GGGOCs8RgjDGmB0sMxhhjerDEYIwxpgeJxJ6hIlINbB3k4fn0uqs6gllZwk+0lAOsLOHqSMoyXlULgu0UkYnhSIjIclUtCXUcQ8HKEn6ipRxgZQlXw1EWq0oyxhjTgyUGY4wxPcRiYrg71AEMIStL+ImWcoCVJVz5XpaYa2MwxhgzsFi8YjDGGDMASwzGGGN6iKnEICLzRWSDiJSLyI2hjudIiEiliKwRkfdEZHmo4/FKRO4Rkb0i8n7AulwReUlENrn/5oQyRq/6KcvNIrLD/VzeE5HzQxmjVyIyVkQWi8h6EVkrIte56yPqsxmgHBH3uYhIioi8IyKr3LL82F0/QUTedj+Tv7tTGgzta8dKG4OIxAMbgXNxJgdaBlyiqutCGtggiUglUKKqEXXTjoicCTQBD6jqce6624E6Vb3VTdg5qvq9UMbpRT9luRloUtVfhDK2wyUio4HRqrpSRDKAFcCngCuIoM9mgHJ8jgj7XEREgDRVbRKRROB14DrgO8DjqvqIiPwBWKWqdw3la8fSFcMpQLmqVqhqG/AIcGGIY4o5qrqEj87QdyFwv7t8P85/5LDXT1kikqruUtWV7vJ+YD3OfOwR9dkMUI6Io44m92mi+1DgY8Bj7npfPpNYSgxjgO0Bz6uI0D8YlwIvisgKEbkq1MEcoUJV3QXOf2xgZIjjOVLXiMhqt6oprKte+iIixcCJwNtE8GfTqxwQgZ+LiMSLyHvAXuAlYDPQoKod7i6+fI/FUmKQPtZFcj3aaap6ErAAuNqt1jChdxcwCZgJ7ALuCG04h0dE0oF/Av+pqvtCHc9g9VGOiPxcVLVTVWcCRTi1Hsf0tdtQv24sJYYqYGzA8yJgZ4hiOWKqutP9dy/wBM4fTaTa49YNd9cR7w1xPIOmqnvc/8xdwJ+IoM/Frcf+J/A3VX3cXR1xn01f5YjkzwVAVRuAMmAOkC0i3bNv+vI9FkuJYRkwxW3RT8KZW3phiGMaFBFJcxvWEJE04Dzg/YGPCmsLgS+5y18C/hXCWI5I95eo69NEyOfiNnT+BVivqr8M2BRRn01/5YjEz0VECkQk210eAZyD02ayGLjI3c2XzyRmeiUBuF3Ufg3EA/eo6v+EOKRBEZGJOFcJ4Mzb/VCklEVEHgZKcYYO3gP8CHgSeBQYB2wDPquqYd+o209ZSnGqKxSoBL7eXUcfzkTkdOA1YA3Q5a7+L5z6+Yj5bAYoxyVE2OciIjNwGpfjcX7EP6qqt7j//x8BcoF3gctUtXVIXzuWEoMxxpjgYqkqyRhjjAeWGIwxxvRgicEYY0wPlhiMMcb0YInBGGNMD5YYzKCIyLXuCJZ/E5ErROR3w/S6xSJyqY/nHy0iL/axvsn99ygRecxdnhkJo3R6JSLfEJEvHsHx/zWU8ZjQscRgButbwPmq+oWhPnHAXZ19KQZ8SwzAfOCF/jaq6k5V7b65aCZwWIkhSNmGxGBfQ1X/oKoPHMFLW2KIEpYYzGFzh/qdCCwUkW/32jZeRBa5g5UtEpFx7kBgFeLIFpGu7rGdROQ1EZnsjpd/t/tr/QH3yuA1EVnpPk51X+JW4Ax3TP1vu+f+XxFZ5r7m1/uJ+QfizMXxsog8LCI39FO8+cBzA5S9WETed++evwW42I3lYveO9HvcWN4VkQvdY64QkX+IyFM4Ax+OFpEl7nHvi8gZfbxOpYjcJs54/O+IyGR3fYGI/NN9jWUicpq7vsf71+tcpSLyqog8KiIbReRWEfmCe941IjIp4Bw3uMtlAa+/sTvG3leHIvK0e/5bgRFumf7mbrvMPf49EfmjOEPfm0igqvawx2E/cO4ezXeXrwB+5y4/BXzJXb4SeNJdfh44FvgEzvAkPwCSgS3u9ptxxs4f4T5PBVLc5SnAcne5FHg6II6rgB+6y8nAcmBCr1hn4dwJmwpkAuXADX2UKR54r5/yNrn/FgPv9y63+/xnOHehAmTjzP+R5u5XBeS6264HfhDwmhn9vL/d+3yxu8zAQ8Dp7vI4nKEfPvL+9TpXKdAAjHbfox3Aj91t1wG/DjjHDe5yGXCHu3w+8HI/ZX4aKA18j9zlY9y/hUT3+e+BL4b679Ye3h6+X9aamDMX+Iy7/CBwu7v8GnAmMAH4OfA14FWcJNFtoaoedJcTgd+JyEygE5jaz+udB8wQke7qnSycRLIlYJ8zgCdUtRlARPobI2s2Hw7RPBjnAZ8MuBpJwfnyBnhJPxxKYhlwjziDvT2pqu/1c76HA/79lbt8DjBd5NBgwZnijptFz/evt2XqDgEhIpuB7naUNcC8fo7pHkhvBU5CPBxn4yTkZW6sI4iAAfiMwxKD8Vv3mCuvAd8AjgJuAr6L80t2ScC+BwKWv40z/tAJOFWeLf2cX4D/UNV+2wV6xTGQBThXNoMlwL+r6oYeK0VmE1A2VV3iVqX9G/CgiPyv9l23r30sxwFzeycA98s38P3rLXAsna6A5130/z3QvU9nwD4d9KyCTunnWAHuV9XvDxCTCVPWxmCG2ps4I9cCfAFnOkJwfomfCnSpagvwHvB1nITRlyxglzrDJF+OU+UCsB/ICNjvBeCb7q9vRGSqOCPOBloCfFpERri/ri/o5zXPBhYFL+IhfcXyH+J+S4vIiX0dJCLjgb2q+ieckUBP6uf8Fwf8+5a7/CJwTcC5Zh5GvEOhEpgpInEiMpaew1e3d38OOO/jRSIyEg7NHT1+eEM1g2VXDGaoXYtTTfJdoBr4MoCqtorIdmCpu99rOCNerunnPL8H/ikin8UZZrj71/BqoENEVgH3Ab/BqeZY6X4hV9NrqkN15v/9O04y2kofyUhECoAWPbzJaRYDN4ozw9bPgZ/gjN672o2lEqdNpbdS4Lsi0o4zZ3R/XUSTReRtnB9wl7jrrgXuFJHVOP9/l+BciQ2XN3Cq6dbgDF29MmDb3ThlX6mqXxCRH+I0tscB7cDVOO+/CXM2uqqJOSJyM70mhheRy4AiVb01ZIEFEJFKoERVa0Idi4k9dsVgDKCqfw11DMaEC7tiMMYY04M1PhtjjOnBEoMxxpgeLDEYY4zpwRKDMcaYHiwxGGOM6eH/Aw/2dv1Q5x2fAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11235b9b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "V = 40     # liters\n",
    "kA = 0.5   # 1/min\n",
    "kB = 0.1   # l/min\n",
    "CAf = 2.0  # moles/liter\n",
    "\n",
    "def cstr(q):\n",
    "    return q*V*kA*CAf/(q + V*kB)/(q + V*kA)\n",
    "\n",
    "q = np.linspace(0,30,200)\n",
    "plt.plot(q, cstr(q))\n",
    "plt.xlabel('flowrate q / liters per minute')\n",
    "plt.ylabel('concentration C_B / moles per liter')\n",
    "plt.title('Outlet concentration for a CSTR')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that, for the parameters given, there is an optimal flowrate somewhere between 5 and 10 liters per minute."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytical Solution using Calculus\n",
    "\n",
    "As it happens, this problem has an interesting analytical solution that can be found by hand, and which can be used to check the accuracy of numerical solutions. Setting the first derivative of $\\bar{C}_B$ to zero,\n",
    "\n",
    "\\begin{align}\n",
    "\\left.\\frac{d\\bar{C}_B}{dq}\\right|_{q^*} = \\frac{V k_A C_{A,f}}{(q^* + V k_A)(q^* + Vk_B)} - \\frac{q^* V k_A C_{A,f}}{(q^* + V k_A)^2(q^* + Vk_B)} - \\frac{q^* V k_A C_{A,f}}{(q^* + V k_A)(q^* + Vk_B)^2} = 0\n",
    "\\end{align}\n",
    "\n",
    "Clearing out the non-negative common factors yields\n",
    "\n",
    "\\begin{align}\n",
    "1 - \\frac{q^*}{(q^* + V k_A)} - \\frac{q^*}{(q^* + Vk_B)} = 0\n",
    "\\end{align}\n",
    "\n",
    "and multiplying by the non-negative denominators produces\n",
    "\n",
    "\\begin{align}\n",
    "{q^*}^2 + q^*V(k_A + k_B) + V^2k_Ak_B - q^*(q^* + Vk_B) - q^*(q^* + Vk_A) = 0\n",
    "\\end{align}\n",
    "\n",
    "Expanding these expressions followed by arithmetic cancelations gives the final result\n",
    "\n",
    "\\begin{align}\n",
    "q^*  =  V\\sqrt{k_Ak_B}\n",
    "\\end{align}\n",
    "\n",
    "which shows the optimal dilution rate, $\\frac{q^*}{V}$, is equal the geomtric mean of the rate constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flowrate at maximum CB =  8.94427191 liters per minute.\n",
      "Maximum CB = -0.954915028125 moles per liter.\n",
      "Productivity =  -8.5410196625 moles per minute.\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "V = 40     # liters\n",
    "kA = 0.5   # 1/min\n",
    "kB = 0.1   # l/min\n",
    "CAf = 2.0  # moles/liter\n",
    "\n",
    "qmax = V*np.sqrt(kA*kB)\n",
    "CBmax = cstr(qmax)\n",
    "print('Flowrate at maximum CB = ', qmax, 'liters per minute.')\n",
    "print('Maximum CB =', CBmax, 'moles per liter.')\n",
    "print('Productivity = ', qmax*CBmax, 'moles per minute.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Solution using `scipy.minimize_scalar`\n",
    "\n",
    "In many applications, a simple analysis leading to an analytical solution is either difficult or intractable. In these cases, numerical solution solution is the only choice for calculating a solution to the optimization problem.\n",
    "\n",
    "The following cell demonstrates the use of the `scipy` optimization library. The  value of $q$ that maximizes $\\bar{C}_B$ can be found using the function `scipy.minimize_scalar`. Note the change in sign needed to convert the  `cstr` function from a maximization to a minimization problem for use with `minimize_scalar`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -0.95491502812526285\n",
       "    nfev: 15\n",
       "     nit: 14\n",
       " success: True\n",
       "       x: 8.9442720169043834"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "\n",
    "V = 40     # liters\n",
    "kA = 0.5   # 1/min\n",
    "kB = 0.1   # l/min\n",
    "CAf = 2.0  # moles/liter\n",
    "\n",
    "def cstr(q):\n",
    "    return -q*V*kA*CAf/(q + V*kB)/(q + V*kA)\n",
    "\n",
    "results = minimize_scalar(cstr, bracket=[0,50])\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell shows how to access the components of the solution obtained from `scipy.minimize_scalar`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flowrate at maximum CB =  8.9442720169 liters per minute.\n",
      "Maximum CB = 0.954915028125 moles per liter.\n",
      "Productivity =  8.54101976458 moles per minute.\n"
     ]
    }
   ],
   "source": [
    "if results.success:\n",
    "    qmax = results.x\n",
    "    CBmax = -results.fun\n",
    "    print('Flowrate at maximum CB = ', qmax, 'liters per minute.')\n",
    "    print('Maximum CB =', CBmax, 'moles per liter.')\n",
    "    print('Productivity = ', qmax*CBmax, 'moles per minute.')\n",
    "else:\n",
    "    print('No solution found.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Bonus Discussion] Batch Reactor\n",
    "\n",
    "A material balance for an isothermal stirred batch reactor with a volume $V = 40$ liters and an initial concentration $C_{A,f}$ is given by\n",
    "\n",
    "\\begin{align}\n",
    "V\\frac{dC_A}{dt} & = - V k_A C_A \\\\\n",
    "V\\frac{dC_B}{dt} & = - V k_A C_A - V k_B C_B\n",
    "\\end{align}\n",
    "\n",
    "Eliminating the common factor $V$\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dC_A}{dt} & = - k_A C_A \\\\\n",
    "\\frac{dC_B}{dt} & = - k_A C_A - V k_B C_B\n",
    "\\end{align}\n",
    "\n",
    "With an initial concentration $C_{A,f}$. A numerical solution to these equations is shown in the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x112366710>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xd4VGXa+PHvPekkkE5LAoQuKCBGEHQVGyArYu+9sLbV1V3fdfe379rWV11XV117wS4qYsGurGIXCIpUKVIDSC+hBEhy//54TmSIKYdkJjNJ7s91nWtmTr2PI3PnOU8TVcUYY4ypTSDSARhjjGkcLGEYY4zxxRKGMcYYXyxhGGOM8cUShjHGGF8sYRhjjPHFEoYxPonIEhE5JtJxGBMpljBMo+b9iO8Qka0islFE3hWRPJ/HdhIRFZHYMMR1oYiUeXFtEZEfROT4UF+niutOEpFLw30d0zxZwjBNwUhVTQHaAauB/0Q4ngrfeHGlAQ8DL4tIWoRjqpY49ptgqmX/c5gmQ1VLgNeAXhXrROS3IvK991f+chG5OeiQz73XTV5JYJB3zGUiMldEikVkjoj0Dzqmn4jMEJHNIvKKiCT6iKsceB5IBroFxXaIiHwtIpu8EsiQoG0XBcWwSER+F3xOERklItO9+/pJRIaLyO3Ab4AHvft50Nt3sIhM9WKeKiKDg84zSURuF5GvgO1A59ruxzRjqmqLLY12AZYAx3jvWwDPAs8FbR8CHID746gPrgRyoretE6BAbND+pwErgIMBAboCHYOuNQVoD2QAc4HLq4nrQuBL730McBWwC2jtrcsB1gMjvNiO9T5ne9t/C3TxYjgC92Pe39s2ANjsHRPwztXT2zYJuDQojgxgI3AeEAuc5X3ODNp/GdDb2x4X6e/UluhdQv7s1pgIeFNESoEUYA0wrGKDqk4K2m+GiIzF/QC/Wc25LgX+qapTvc8LK21/QFVXAojI20C/GuI6REQ24UoWpcC5qrrG23Yu8J6qvud9/lhECnEJ5FlVfTfoPJ+JyEe40sN3wCXAGFX92Nu+ooYYfgssUNXnvc9jReQaYCTwjLfuGVWdXcM5jAHskZRpGk5U1TQgAbga9wPbFkBEBorIpyKyVkQ2A5cDWTWcKw/4qYbtPwe9345LUtX51osrHZiA+8Gv0BE4zXsctclLLIfh6mEQkeNE5FsR2eBtGxEUd20xBmsPLK20bimuVFJhuc9zmWbOEoZpMlS1TFVfB8pwP74AL+F+rPNUNRV4FPeYB9zjqMqW4x4FhTKurcCVwHkicmDQdZ5X1bSgJVlV7xSRBGA88C+gjZd03guKu6YYK9/TSlxyCtaBvUslNmS18cUShmkyvFY+o3B/0c/1VrcENqhqiYgMAM4OOmQtUM7eFb1PAn8SkYO883UVkco/uPtMVdd75/67t+oFYKSIDBORGBFJFJEhIpILxONKS2uBUhE5DhgadLqngItE5GgRCYhIjoj09LatrnQ/7wHdReRsEYkVkTNwjQLeqe89mebHEoZpCt4Wka3AFuB24IKgZ/JXAreKSDHux/rVioNUdbu3/1feY6FDVHWct+4loBhX15ERojjvA0aISB9VXQ6MAv6KSwzLgRuAgKoWA9d4sW7EJbkJQXFPAS4C/o2r/P6MPaWI+4FTvT4pD3iJ6njgj7hK9f8BjlfVdSG6J9OMiKqVRo0xxtTOShjGGGN8sYRhjDHGF0sYxhhjfLGEYYwxxpcm1dM7KytLO3XqFOkwjDGm0Zg2bdo6Vc32s2+TShidOnWisLAw0mEYY0yjISKVRwKolj2SMsYY44slDGOMMb5YwjDGGONLk6rDMMaYcNm9ezdFRUWUlJREOpQ6SUxMJDc3l7i4uDqfwxKGMcb4UFRURMuWLenUqRMiUvsBUURVWb9+PUVFReTn59f5PGF7JCUied48BHNFZLaIXFvFPiIiD4jIQm/ay/5B2y4QkQXeckG44jTGGD9KSkrIzMxsdMkCQETIzMysd+konCWMUuCPqvqdiLQEponIx6o6J2if43BzHHcDBgKPAANFJAO4CSjAjdU/TUQmqOrGMMZrjDE1aozJokIoYg9bCUNVV6nqd977Ytz8BDmVdhuFm39ZVfVbIE1E2uGm2PxYVTd4SeJjYHg44iwvVx78ZAGfzV8bjtMbY0yT0SCtpESkE3AgMLnSphz2nh6yyFtX3fqqzj1aRApFpHDt2n3/0Q8EhMc+X8Qnc1fv87HGGNOchD1hiEgKbrrJP6jqlsqbqzhEa1j/65Wqj6tqgaoWZGf76t3+K+1Tk1i5uXG2fDDGmIYS1oQhInG4ZPGiN9dyZUW4Ce0r5OLmIK5ufVi0TU3kZ0sYxphGYty4cQwcOJA+ffrQtWtXbrnllga5bjhbSQlu7uG5qnpvNbtNAM73WksdAmxW1VXAh8BQEUkXkXTcfMYfhivWdqmJrLKEYYxpBJ599lnuuusuxo8fz4wZM5g+fTotWrRokGuHs5XUocB5wEwRme6t+yvQAUBVH8VNUD8CWAhsx81TjKpuEJHbgKnecbeq6oZwBdouNYl1W3eys7SMhNiYcF3GGGPqZcuWLVx//fVMnTqV3NxcAFJSUrjhhhsa5PphSxiq+iVV10UE76PAVdVsGwOMCUNov9IuNRGANVt2kpfRMJnaGNN43fL2bOasrFwlWz+92rfippG9a9znjTfeYODAgXTu3Dmk1/bLxpLC1WEA9ljKGBPVZs+eTb9+/ardfsYZZ3DPPfeE7fo2NAh7ShirNu+IcCTGmMagtpJAuCQnJ7NjR9W/U2+99RbHH388EydODNv1rYQBtEtLAqyEYYyJbiNGjGDcuHGsXu36je3cuZMnnniCkpISxo0bx3nnncfmzZvDdn0rYQApCbG0TIi1prXGmKh28MEHc/PNNzNs2DDKysooLS3l3HPP5e6772br1q1cfvnlzJ49mx07dpCUlBTy61vC8LRNTbRHUsaYqHfeeedx3nnn/fJ52bJl3HLLLbz55psA3HLLLcyYMYOBAweG/NqWMDzt0pLskZQxptHp0KEDTz311C+fb7rpprBdy+owPO1aWec9Y4ypiSUMT9vURNZt3cmu0vJIh2KMMVHJEoanXWoiqrB6i5UyjDGmKpYwPBVNa3+2hGGMMVWyhOFpZ729jTGmRpYwPL8MD7LJmtYaY0xVLGF4WiXGkZIQayUMY4yphiWMIDaRkjHGVM8SRpB21tvbGGOqZQkjiM28Z4xpDJriFK1jRGSNiMyqZvsNIjLdW2aJSJmIZHjblojITG9bYbhirKxtahJrt+5kd5l13jPGRKdITtEazhLGM8Dw6jaq6t2q2k9V+wF/AT6rNA3rkd72gjDGuJeKzntrinc21CWNMca3iilaX3311SY3RevnItLJ5+5nAWPDFYtfwU1rc9JCPzSwMaaJeP9G+HlmaM/Z9gA47s4ad4nqKVpFJCAig8MZgIi0wJVExgetVuAjEZkmIqNrOX60iBSKSOHatWvrFUv7VJtIyRgTvWqaorWgoICrrrqKI444gtmzZ4fl+jWWMFS1XETuAQaF5erOSOCrSo+jDlXVlSLSGvhYRH5U1c+rifFx4HGAgoICrU8gFSUMa1prjKlRLSWBcKluitbly5czYMAAHnroIe69916Kioro3Tv008j6qcP4SEROEREJ+dWdM6n0OEpVV3qva4A3gAFhuvZeWiXG0iI+hpXWtNYYE4Wqm6J12rRpzJ8/n4svvphPPvmEYcOGheX6fuowrgeSgTIR2QEIoKraqr4XF5FU4Ajg3KB1yUBAVYu990OBW+t7LZ/xkJOWxIqNljCMMdGnuilaly1bxj333EPfvn055ZRT2LZtG8nJySG/fq0JQ1Vb1uXEIjIWGAJkiUgRcBMQ553zUW+3k4CPVHVb0KFtgDe8Ak0s8JKqflCXGOoiL6MFyy1hGGOiVOUpWsGVPNatW0cgEGC//fYLS7IAHwnDexR1DpCvqreJSB7QTlWn1HScqp5V27lV9Rlc89vgdYuAvrUdGy656UlMXbwBVSV8T+GMMSZ03nvvvQa5jp86jIdxld5ne5+3Ag+FLaIIy0tvQfHOUjbv2B3pUIwxJqr4SRgDVfUqoARAVTcC8WGNKoLyMlzT2uUb7LGUMcYE85MwdotIDK5vBCKSDTTZsTNy010X++Ubt0c4EmOMiS5+EsYDuKatrUXkduBL4I6wRhVBeRlewthgCcMYszfVenX1iqhQxO6nldSLIjINOBrXpPZEVZ1b7ytHqdSkOFolxloJwxizl8TERNavX09mZmajaxCjqqxfv57ExMR6ncdPK6nnVfU84Mcq1jVJeRktrA7DGLOX3NxcioqKqO8QRJGSmJj4y4CFdeWn495e/cu9+oyD6nXVKJeX3oIFa4ojHYYxJorExcWRn58f6TAiqto6DBH5i4gUA31EZIuIFHuf1wBvNViEEZCXkUTRxh2N+nmlMcaEWrUJQ1Xv8Hp5362qrVS1pbdkqupfGjDGBpeX0YKdpeWstXkxjDHmF9U+khKRnqr6IzBORPpX3q6q34U1sgjKC2pa27pV/SqJjDGmqaipDuOPwGXAPVVsU+CosEQUBYI77x3UMcLBGGNMlKg2YajqZd7rkQ0XTnT4pfOe9cUwxphf1PRI6uSaDlTV10MfTnRIjIshu2WC9cUwxpggNT2SGlnDNgWabMIAyEtPsr4YxhgTpKZHUhc1ZCDRJi+jBdOWbox0GMYYEzX8jCXVLOWmJ7FqcwmlZU12nEVjjNknYUsYIjJGRNaIyKxqtg8Rkc0iMt1b/h60bbiIzBORhSJyY7hirEleegvKypVVm0sicXljjIk6NSYMEQmIyOA6nvsZYHgt+3yhqv285VbvmjG4CZqOA3oBZ4lIrzrGUGe/jFprFd/GGAPUkjBUtZyq+2HUSlU/BzbU4dABwEJVXaSqu4CXgVF1iaE+KjrvFVnFtzHGAP4eSX0kIqdIeMbzHSQiP4jI+yJSMchhDrA8aJ8ib12VRGS0iBSKSGEoR5Fsl5ZIQKyEYYwxFfyMVns9kAyUicgO3JwYqqqt6nnt74COqrpVREYAbwLdvPNXVu0ogKr6OPA4QEFBQchGC4yLCdAuNYll1nnPGGMAHyUMb8DBgKrGBQ1CWN9kgapuUdWt3vv3gDgRycKVKPKCds0FVtb3enWRn5XMkvWWMIwxBnwkDHHOFZH/9T7niciA+l5YRNpWPObyzhcA1gNTgW4iki8i8cCZwIT6Xq8u8rOSWbR2qw1zbowx+Hsk9TBQjhts8DZgK64V08E1HSQiY4EhQJaIFAE3AXEAqvoocCpwhYiUAjuAM9X9MpeKyNXAh0AMMEZVZ+/7rdVf5+xkiktKWb9tF1kpCZEIwRhjooafhDFQVfuLyPcAqrrR+8u/Rqp6Vi3bHwQerGbbe8B7PmILq/ysZAAWr9tmCcMY0+z5aSW12+sboQAiko0rcTR5nbNSAFi0dmuEIzHGmMjzkzAeAN4A2ojI7cCXwP+FNaookZOeRHxMgEXrtkU6FGOMibhaH0mp6osiMg042lt1oqrODW9Y0SEmIHTMbMHitZYwjDHGTx0GQAtcBbQCSeELJ/rkZyWz2EoYxhjjq1nt34FngQwgC3haRP4W7sCiRX52MkvXb6es3JrWGmOaNz8ljLOAA1W1BEBE7sT10v5HOAOLFl2yUthVVs6KjTvokNki0uEYY0zE+Kn0XgIkBn1OAH4KSzRRKD/bNa1dtM5aShljmjc/CWMnMFtEnhGRp4FZwFYReUBEHghveJFX0RdjkVV8G2OaOT+PpN7wlgqTwhNKdMpMjqdVYmzkK77LywCBgE2SaIyJDD/Nap9tiECilYiQn53SsAlj4xKYMwGWfg1r58LmFVC+GyQASRmQ2QXa9oEuR0H+byChZcPFZoxptvw2q23WOmclM3nR+vBfaOnX8Nk/YdGn7nNmV8g5CHqNgrhkKNsF29bC+oUw/SWY+gQE4qDTYXDQBdDzeIiJC3+cxphmyRKGD/lZybzx/Qp27CojKT4m9BfYuhY++DPMGg/JreGov0GfMyCtQ/XHlO6E5ZNh4USY9TqMuxBS2kLBRTDwckhKC32cxphmbZ8ShogEgBRV3RKmeKJSZ6+l1JL129ivXb2nAtnbsm/hlfOgZBMM+SsM/j3E+2i+G5sA+Ye75eibYMHHUPgUTLoDvn0YBl/jEkdCSmjjNcY0W3467r0kIq1EJBmYA8wTkRvCH1r0CB61NqRmjYdnjnc/6qM/gyF/9pcsKgvEQI/hcM44+N0X0GEQfHIbPNAPpo8Fm8/DGBMCfprc9PJKFCfihhzvAJwX1qiizJ6mtSHsizHzNRh/KeQeDJd9Am16hea87frA2a/AJRMhvRO8eTk8PQJWzwnN+Y0xzZafhBEnInG4hPGWqu6mhjm2K4jIGBFZIyKzqtl+jojM8JavRaRv0LYlIjJTRKaLSKHfmwmXFvGxtEtNDN2otT99Cq+Pho6HwrmvQVJ6aM4bLO9guPgjGPmAa2n12G9g0p1Qtjv01zLGNAt+EsZjuN7eycDnItIR8FOH8QwwvIbti4EjVLUPbia/xyttP1JV+6lqgY9rhV2X7BQWrglBCWPdAnj1AsjuAWe+BPHJ9T9ndQIB13rq6mnQ+yRXv/Hk0bCmWQw2bIwJsVoThqo+oKo5qjpCnaXAkT6O+xzYUMP2r1V1o/fxWyDXb9CR0L1NSxas3kp5fQYh3F0C4y6CmFj32CgxxBXo1UnOhFOehNOfg81F8NjhMOUJq9swxuwTP5XebUTkKRF53/vcC7ggxHFcArwf9FmBj0RkmoiMriW+0SJSKCKFa9euDXFYe/Rom8KO3WUs37i97if57y2weiac+EjNTWbDpdcouHIydB4C7/0Jxl0AJZsbPg5jTKPk55HUM8CHQHvv83zgD6EKQESOxCWMPwetPlRV+wPHAVeJyOHVHa+qj6tqgaoWZGdnhyqsX+nR1pUGfvy5uG4nWDbZNXcdMBq6DwthZPsoJRvOegWOvRXmvgOP/gZWfBe5eIwxjYafhJGlqq/izeOtqqVAWSguLiJ9gCeBUar6S1dqVV3pva7BjWM1IBTXq49urV1/hvl1SRhlu+GdP0CrXNdnItICATj0Wrj4AzdG1VNDYfJj9ojKGFMjPwljm4hk4rWMEpFDgHo/xxCRDsDrwHmqOj9ofbKItKx4DwzFjZAbUckJsXTIaMGPq+uQMCY/BmvmwIi7o6sjXd4AuPwL6HoMvP8/8NZVrp7FGGOq4Ken9/XABKCLiHwFZAOn1naQiIwFhgBZIlIE3ATEAajqo8DfgUzgYREBKPVaRLUB3vDWxQIvqeoH+3Zb4dG9Tct9L2Hs2ASf3w1djoaeI8ITWH20yHCttT67Cz67E9bOgzNegFbtIh2ZMSbK+Bmt9jsROQLoAQgwz+uLUdtxZ9Wy/VLg0irWLwL6/vqIyOvZtiWfzlvDztIyEmJ9jin11X1u2I9jbg5naPUTCMCRf4E2veGNy+HxIXDmi5AbFS2ajTFRotqEISInV7Opu4igqq+HKaao1b1tS8rKlUVrfY4ptXUtfPsoHHCa64Ed7Xqd4IZOH3sWPH0cHH8fHHhOpKMyxkSJmkoYI2vYprj6h2alZ1s378S8n4v9JYzJj0BpCRzx59r3jRZtesPoSW7027euhJ9nwrDb3XhVxphmrdqEoaoXNWQgjUF+VjJxMcI8PxXfJVtgypPur/asbuEPLpRaZMC5r8NHf3NJb+NiOOWp6KqwN8Y0OD8d91JF5N6KznEico+IpDZEcNEmLiZAl+wUfxXf056GnZvhsOvCH1g4xMTCcXfCb+9xQ6c/PdzN/GeMabb8NKsdAxQDp3vLFuDpcAYVzbq3aVl7573yMpj6JHQ8DNof2DCBhcvBl8LZr8KGJW4cqlU/RDoiY0yE+EkYXVT1JlVd5C23AJ3DHVi06tG2JSs27aC4pIaGYgsnwqZlMOBXjcAap27HwCUfgsTAmONg3vu1H2OMaXL8JIwdInJYxQcRORTYEb6QoluPNq7ie/7qGkaunfokpLRxc2w3FW16w2X/hezu8PLZrvWXMaZZ8ZMwrgAe8uaoWAo8CFwe3rCiV4+2FQmjmsdSm5a7Z/79L4CYuAaMrAG0bAsXvgs9Rrg5yN+7AcpKIx2VMaaB+Om4Nx3oKyKtvM/Naj7vynLSkkiOj2FedfUYM14BtOn2X4hPhtOfh4l/h6//AxuXwKlPWwsqY5qBWhOGiKQB5wOdgFhvyA5U9ZqwRhalAgGhe9uW/PhzFXlT1SWMDoPd9KhNVSAAQ/8BGZ3h3T+5FlRnvwqt2td+rDGm0fLzSOo9XLKYCUwLWpqt3u1bMXvlll9PprTyO1g3H/qeGZnAGlrBxV4LqsXwxNGuk58xpsnykzASVfV6VX1aVZ+tWMIeWRQ7ICeV4pJSlm6oNJnSD69ATIKbqKi56HaMGyYdYMxwWDAxsvEYY8LGT8J4XkQuE5F2IpJRsYQ9sih2QE4aADOKNu1ZWV4OcydAt2MhKS1CkUVI2wNcC6qMfHjpdCgcE+mIjDFh4Cdh7ALuBr5hz+OownAGFe26tUkhITbArBVB04KsKITiVc2rdBGsVXu46H3oejS8cx189L8uiRpjmgy/82F0VdV14Q6msYiLCdCrfStmFAUljDlvQSAustOvRlpCSzhzrGty+/UDrgXVyY9DXFKkIzPGhICfEsZsYHutezUzB+SkMmvFZlfxrQpzJkCXIyGxWQ6ztUdMLIz4Fwy9Hea+Dc+OdMO8G2MaPT8JowyYLiKPicgDFYufk4vIGBFZIyJVTrEqzgMislBEZohI/6BtF4jIAm+5wN/tNJwDclLZtquMReu2udZBm5fBfjWNCN+MiMDgq+GM5+HnWW4MqrXzaz/OGBPV/CSMN4Hbga/Z92a1zwDDa9h+HNDNW0YDjwB4leo3AQOBAcBNIpLu85oNok+uq9ietWIzLPjQrezWjB9HVWW/ka5n+O7t8NQxsOTLSEdkjKkHPz2969yEVlU/F5FONewyCnhOVRX4VkTSRKQdbi7wj1V1A4CIfIxLPGPrGkuodclOJjEuwIyizZy45mM3Km3LNpEOK/rkHgSXToQXT4fnToRRDzaffirGNDF+ShjhlAMsD/pc5K2rbv2viMjoirk61q5tuGflsTEBerdPZfHyZVA0FboNbbBrNzrpndxotx0OgTd+B5PudPU+xphGJdIJQ6pYpzWs//VK1cdVtUBVC7Kzs0MaXG0OyEkl8+cvQcstYdQmKd3N4tfvHJh0B7x5BZTuinRUxph9UG3CEJG/iEi4Z/8pAvKCPucCK2tYH1X65KZyqH5HaWIGtO9f+wHNXWw8jHoIjvwb/DAWXjgZdmyMdFTGGJ9qKmEsBq4Vke9F5BkROSMMFc8TgPO91lKHAJtVdRXwITBURNK9aw711kWVA9q34tDAbFZlDHQD8pnaicARN8DJT8DyyfDUUNiwKNJRGWN8qLbSW1VfBl4G8Eoaw4HXRSQGmAh8oKpTajq5iIzFVWBniUgRruVTnHf+R3EDG44AFuL6elzkbdsgIrcBU71T3VpRAR5NOssKYmQTE2L77lUcMj70OR1a5cAr58ATR8Hpz0H+4ZGOyhhTAz89vVHV74HvgTu8eTGOBS4FakwYqnpWLdsVuKqabWNw84lHrZglXwDw7rbunBDhWBqlTofCpf+FsWfB8yfBcf+Egy+JdFTGmGrs83MUVd2iquNVdXQ4AmpUFn/GxoT2fPpzC0p2l0U6msYpswtc+jF0OQrevd7Nr1FWw3zpxpiIsQfvdVVeBku+YFv7wewqK2dm8ECEZt8kpsJZL8Pg38PUJ+CFU2B71D2BNKbZs4RRVz/PgJLNpPU+FoCpS+wHrl4CMW4WvxMfgWXfeMOJzIt0VMaYIL4ShojkiMhgETm8Ygl3YFFv2bcApHT7DV1bp1C4xJqHhkS/s+GCd2BnMTx5DCz4ONIRGWM8tSYMEbkL+Ar4G3CDt/wpzHFFv2XfQFoHSM3h4E7pFC7Z8OspW03ddBgIl30K6R3hxdPg839Zz3BjooCfEsaJQA9VHaGqI72leTcKUnUljA6DACjomMGWklLmrymOcGBNSFoeXPwh7H8yfHIbvHIulGyJdFTGNGt+EsYivL4TxrNxMWxdDXkDATi4k5uxdqo9lgqt+GQ45SkYdgfMe9/117B6DWMixk/C2E4d58Nosrz6i4oSRl5GEq1bJlBoFd+hJwKDroQLJkDJJpc05kyIdFTGNEt+EsYE4DbqNh9G07TsG9cUNLsnACLCwZ0yrOI7nDodBqM/c//NXz0PJt7smjYbYxqMr/kwRCQe6O6tmqeqzbtnVVEh5A7Ya/yogk7pvDtzFSs27SAnzeawDovUHLjoPXj/z/Dlv2Hl93DKGEjOjHRkxjQLflpJDQEWAA8BDwPzm3Wz2p3FsGYu5BbstbqiHsMeS4VZbAKMvA9O+A8s/RoeOxyWTY50VMY0C34eSd0DDFXVI1T1cGAY8O/whhXFVv0A6K+GM+/ZtiXJ8THWga+h9D/ftaKKiYWnj4Mv74Py8khHZUyT5idhxKnqL01TVHU+zbnV1Aqv+iZn74QRGxOgf8d0pi62eowGk9Mffvc57Hc8TLwJXjodtq2PdFTGNFl+EkahiDwlIkO85Qmac6X3immQ1hGSs3616ZDOmcxbXcya4pIIBNZMJabCac/CiH/B4s/g0cNg6TeRjsqYJslPwrgCmA1cA1wLzAEuD2dQUW3Fd5BzUJWbDu/mpoj9auG6hozIiMCAy+CSj10dxzO/hS/utUdUxoRYrQlDVXeq6r2qerKqnqSq/1bVnQ0RXNTZugY2L682YfRu34r0FnF8scASRkS07+ceUfU6Af57C7x0GmxdG+mojGkyaprT+1XvdaaIzKi8+Dm5iAwXkXkislBEbqxi+79FZLq3zBeRTUHbyoK2RUdPrZXfu9ecqufvDgSEQ7tm8cWCdaiNfRQZia3g1Kfht/fC4i/gkUEwP+pm9zWmUaqpH8a13uvxdTmxN5XrQ7jZ+YqAqSIyQVXnVOyjqtcF7f974MCgU+xQ1X51uXbYrPLyZNsDqt2cuwdbAAAgAElEQVTl8G7ZvDNjFfNWF9OzbasGCszsRcTN3NdhELx+masML7jEDZ8e3yLS0RnTaFVbwlDVVd7bK1V1afACXOnj3AOAhaq6SFV34eYHH1XD/mcBY/0GHhE//wAZnSGhZbW7HNbNVYZ/Md8eS0Vcm15w2Scw6GoofMr12agoJRpj9pmfSu9jq1h3nI/jcoDlQZ+LvHW/IiIdgXzgk6DViSJSKCLfisiJ1V1EREZ7+xWuXRvm59WrZkDbPjXu0j4tia6tU/h8gT07jwqxCTDsdjj/Ldi1zc2x8cU9NqyIMXVQUx3GFSIyE+hRqf5iMeCnDkOqWFfdg/0zgddUNfhfcQdVLQDOBu4TkS5VHaiqj6tqgaoWZGdn+wirjnZsgk1LoV3NCQPgN92ymLJ4g83zHU06D4ErvoL9RsJ/b4VnjoeNSyMdlTGNSk0ljJeAkbjBB0cGLQep6rk+zl0E5AV9zgVWVrPvmVR6HKWqK73XRcAk9q7faHirZ7nXtn1r3fXwbtnsLC23Xt/RpkWGqxA/6TH4eSY8cigUPm2TMxnjU011GJtVdYmqnuXVW+zAlRBSRKSDj3NPBbqJSL43eOGZuOSzFxHpAaQD3wStSxeRBO99FnAorv9H5Pio8K4wsHMGcTFizWujkQj0PdOVNnIOhHf+AM+NstKGMT74GXxwpIgsABYDnwFLgPdrO05VS4GrgQ+BucCrqjpbRG4VkeAZ+84CXta926Huh+th/gPwKXBncOuqiPh5BqS0gZZtat21RXwsBR0z+Hy+1WNErfSOcP4EOP7frvf+w4NgyhPW2c+YGtQ6vDnwD+AQYKKqHigiR+J+5Gulqu8B71Va9/dKn2+u4rivgdr/lG9IP8/yVbqoMKRHNne8/yNFG7eTm25NOaOSCBRcDF2Phbevgff+BLPfhFH/ca3hjDF78dNKareqrgcCIhJQ1U+B6OofEW5lpbBuHrTp7fuQob3bAvDR7NXhisqESloenPs6nPCgq9t4eDB8+4iVNoypxE/C2CQiKcDnwIsicj9QGt6wosyGRVC2C1r38n1IflYyPdq05IPZP4cxMBMyItD/PLjqW8g/HD64EcYMdQnEGAP4SxijcPN6Xwd8APyEay3VfKzxqk9a77dPhw3r3YbCJRtYv7V5Dr3VKLVqD2e/Aic9DhsWw2NHwIf/D3ZujXRkxkRcjQnDG97jLVUtV9VSVX1WVR/wHlE1H2vmgAQgq3vt+wYZ2rst5QoT59pjqUZFBPqeAVdPdaWObx6EhwbA3HciHZkxEVVjwvA60m0XkdQGiic6rZnjKkHj9m2u7t7tW5GbnsSHVo/ROLXIgJH3u2HTk9LhlXPgpTNh07JIR2ZMRPh5JFUCzPQmUXqgYgl3YFFlzdx9fhwFICIM692WLxeso7hkdxgCMw0ibwCM/swNXrj4c3hooJsStnRXpCMzpkH5SRjvAv+Lq/Se5i2F4Qwqquze4Sq996HCO9iw3m3ZVVbOpHnWJ6NRi4mFwb+HqyZDl6PclLCPDIL5H0U6MmMajJ+EkebVXfyy4HpmNw/r5oOW16mEAXBQx3SyUuL50FpLNQ1peXDmi3DOa4C4SZpeOAXWzqv1UGMaOz8J44Iq1l0Y4jii15of3Wt23RJGTEA4tlcbPv1xjQ1G2JR0Oxau/AaG3QHLp7qe4u/fCDs2RjoyY8KmptFqzxKRt4F8EZkQtHwKNJ9WUuvmg8TUq+fvyD7t2barjI/mWOV3kxITB4OuhGu+g4MugCmPwQP9YeqTrrOnMU1MTUODfA2sArKAe4LWF+NvePOmYd18yMiH2Pg6n+KQzpm0T01k/LQiTujbPoTBmaiQnOXGpCq4GD74C7z7R5jyJBxzM3Qf5prpGtME1DRa7VJVnaSqg1T1s6DlO29gweZh3YJ97n9RWSAgnNQ/hy8WrGXNlpIQBWaiTtsD4IK34fTn3cgAY8+Ap4+DZZMjHZkxIeFntNqTRWSBiGwWkS0iUiwiWxoiuIgrK4UNP0FWt3qf6uT+uZQrvDl9RQgCM1FLBHqd4FpT/fZe18JuzFAYe/ae+jBjGik/ld7/BE5Q1VRVbaWqLVW1VbgDiwqblrq/FOtZwgDokp1Cv7w0xk9bgdqEPU1fTBwcfAlc8z0c9TdY8oVrhvvmVbC5KNLRGVMnfhLGalWdG/ZIotG6Be41BAkD4JSDcpm3upjZK5tHAc0A8clw+A1wzXQ45EqY+aqrGP/gr1BsjSBM4+InYRSKyCteq6mTKxY/JxeR4SIyT0QWisiNVWy/UETWish0b7k0aNsF3qOwBSJSVdPe8Fs3371mdg3J6Ub2aUd8TIDx39lfmM1OciYMux1+Pw0OOBUmPwr397HEYRoVPwmjFW602qHsmdf7+NoO8gYufAg4DugFnCUiVXWXfkVV+3nLk96xGcBNwEBgAHCTiDR8Z8F18yE5240pFAJpLeI5er/WTJi+kl2lNtdCs5TWAU582A1suP8pQYnjL1BsnTtNdKs1YajqRVUsF/s49wBgoaouUtVdwMu4odL9GAZ8rKobVHUj8DEw3OexoROCFlKVnV6Qx/ptu2yejOYus0ulxPEY3N/XEoeJan5aSXUXkf+KyCzvcx8R+ZuPc+cAy4M+F3nrKjtFRGaIyGsikrePxyIio0WkUEQK164N8XhN6xe6f9ghdET3bDpltuCZrxaH9LymkaoucbxzvZuPw5go4ueR1BPAX4DdAKo6AzjTx3FV9Vaq3DzobaCTqvYBJgLP7sOxePE8rqoFqlqQnZ3tIyyfSjbD9nWQEdqEEQgIFwzuxHfLNvHD8k0hPbdpxCoSx+8Loc/p8P3z8J/+MO4iWPVDpKMzBvCXMFqo6pRK6/x03CsC8oI+5wIrg3dQ1fWqWjEd3RPAQX6PDbuKv+7qMSRIdU49KJfk+Bie/XpJyM9tGrmMznDCf+DaGW503IUT4bHD4bkTYdEksCbZJoL8JIx1ItIF7y98ETkVN2RIbaYC3UQkX0TicaWSCcE7iEi7oI8nABXNdz8EhopIulfZPdRb13A2/ORew5AwWibGcVpBHm/PWMmaYuv5barQqh0ceytcN8sNMbJmDjw3Ch4fArPGQ5nNr2Ianp+EcRXwGNBTRFYAfwCuqO0gb/iQq3E/9HOBV1V1tojcKiIneLtdIyKzReQH4Bq8UXBVdQNwGy7pTAVu9dY1nA2L3GtGflhOf/6gjuwuU8ZOXl77zqb5SkyFw66DP8yEkQ/Arq3w2sWunuOLe2Bb8xkH1ESe+O11LCLJQEBVi8MbUt0VFBRoYWGI5nZ680r46RP4Y/iGc7jw6SnMXrmFr/58FPGxfnK3afbKy2DBx6457qJPISYB+pwGAy93Y1kZs49EZJqqFvjZ108rqf8TkTRV3aaqxd5jon/UP8wot/6nsDyOCnbJYfmsLd7Ja9OsI5/xKRADPYbD+W/ClZPhwHNg1uvw6GHw9G9hzgQbWt2EjZ8/a49T1V+a83j9IkaEL6QosWFR2BPGYV2zOLBDGg99utA68pl917qnG1b9+jlw7G2waRm8eh7cdwB8eoeNWWVCzk/CiBGRhIoPIpIEJNSwf+O3sxi2rQl7whAR/nBMd1Zs2mGlDFN3Selw6DVuoMMzXoQ2veCzu1zieOkMmPeBe5RlTD3VNIFShReA/4rI07iWUhezp79E0/RLhXd4EwbA4d2y6JfnShmnHpRrdRmm7mJiYb/j3bJxCXz3vOvPMf8DaJUD/c+HA8+D1Cr7wBpTKz9Dg/wTuB3YD+gN3Oata7oaMGG4UkY3VmzaYYMSmtBJ7wRH/y9cN9tN6JTdEybdCfftD8+fDDNfg907Ih2laWT8lDBQ1feB98McS/TYuMS9hqlJbWVHdM+mX14aD36ykFP6WynDhFBMnJvQqdcJ7v/r71+AH16G8ZdAQivofSL0PRs6HGJTyZpa2Yx7Vdm4FJIyIKFlg1xORLjuWFeX8dw3SxrkmqYZSu/kJnO6doabSrbn8TBzPDw9HB44ECbdteePJWOqYDPuVWXTUkjv2KCXPKJ7NkN6ZHP/xAWs27qz9gOMqatAAPIPh5MegT/NhxMfgdRcmPR/rkPgk8fAt4/AFj8DOpjmxGbcq8rGpZDWsAkD4G+/7cWO3WX868N5DX5t00wlpEC/s+HCd1xv8mNuht0l8MGNcO9+8MzxUPg0bG/YgRZMdArrjHuNUnk5bF7e4CUMgK6tU7hwcCdeKVzOzKLNDX5908yldXDDkFzxJVw1BY74MxSvgnf+AP/qBi+e5upAbDiSZstPpXfwjHsVFHg9LBFFWvEqKNsVkRIGwO+P7sYb36/glrdnM+7yQYhVRJpIyO4BR/4FhtwIP89wAx7OegMWfAQSgI6Hwn4joedv3eMs0yzUmjBU9aKGCCRqbFrqXiNQwgBITYrjhmE9uPH1mYwrLOL0g/NqP8iYcBGBdn3dcswtsGo6zH0HfnwH3v8ft7Q/0EseIyE7tDNUmujip5VUroi8ISJrRGS1iIwXkab7J8VGL2GkdYpYCKcX5DEwP4Pb3pnDik3WVt5ECRGXHI7+X7hqMlxdCEffBAj891Z46GB48GD4+CZY8pWNadUE+anDeBo3j0V73DSpb3vrmqZNSwGBtMj9ZR8ICHef2pcyVW4cPwO/Iwob06CyusFvrofRn8J1c+C4u6FlW/jmQXhmBPyzM4y7EKaPha0hnj7ZRISfOoxsVQ1OEM+IyB/CFVDEbVwKLdtBbGSHy+qQ2YK/jtiPv705ixcnL+PcQyLziMwYX1JzYOBot5RsdrMDLvjIDcU++w1AIKc/dBsK3Y6Fdge65r2mUfGTMNaJyLnAWO/zWUDTbSYRgT4Y1TlnYAc+nP0z//feXH7TLYuOmcmRDsmY2iWmQq9Rbikvd5XmCz5yy6Q7YdId0CIT8o+AzkPcEiX/5kzNap1ASUQ6AA8Cg3Cto74GrlXVpbWeXGQ4cD8QAzypqndW2n49cClujvC1wMUV5xWRMmCmt+syVT2BWoRkAqV7e0Onw+Dkx+p3nhBZuWkHx93/Be3Tknj9isEkxcdEOiRj6m7bejdP+U+fuFLI1p/d+vT8Pckj/3BokRGxEJubfZlAyfeMe3UIIgaYDxwLFOGmWj1LVecE7XMkMFlVt4vIFcAQVT3D27ZVVVP25Zr1ThhlpfCPbPjNH90QClFi0rw1XPTMVE7sl8O9p/e1pramaVCFdfNd4lg0CRZ/AbuKAa9lVv7h7o+3vIGQlBbhYJuufUkYtT6SEpFncSWKTd7ndOAeVb24lkMHAAtVdZF33MvAKOCXhKGqnwbt/y1wrp+gw6Z4FWg5pEZXU9YhPVpz3THduffj+fTLS+OCwZ0iHZIx9Sfi+ntk94CBv3N/sK38bk8C+fYR+PoBQKDt/q7vR8fB7jU5K8LBN09+6jD6VJ5xT0QO9HFcDrA86HMRMLCG/S9h7xFxE0WkEPe46k5VfbOqg0RkNDAaoEOHDj7CqkHFDGVR2BHp6iO7MqNoE7e9M4cebVtySOfMSIdkTGjFxELeALcc8T9u+PWiQlj6NSz9CqY96+YyB8jq4SWPwZB7sBtY0UreYecnYQREJN2bmhURyfB5XFXfXpXPv7xK9QLgiKDVHVR1pYh0Bj4RkZmq+tOvTqj6OPA4uEdSPuKq3mYvv0VZCQNcU9t7Tu/HKY98zWXPFvLK7wbRq33THgPSNHNxSZD/G7cAlO5yHQeXfuWSyKzxMM1rwJmc7RJHbgHkDnD9RRL26Ym28cHPD/89wNci8hruB/903IRKtSkCgn95c4GVlXcSkWOA/wccoaq/DNOqqiu910UiMgk4EPhVwgipXxJGdM5IlpoUx3MXD+CUR77m/DFTeP2KwXTIbBHpsIxpGLHxe0ogh13npp1dPRuKprqSSNFUmPee21cC0Lo35B3sJZKDIbOrlULqyVelt4j0Ao7ClRr+G1xxXcMxsbhK76OBFbhK77NVdXbQPgcCrwHDVXVB0Pp0YLuq7hSRLOAbYFRt1613pfc718HsN+HPi+t+jgawYHUxpz32DalJcYy7fBCtWyZGOiRjosP2DbBiGiyf4hLIimmw05u+JzF1zzAn7fq5Ukh6frPvDxIVraS8QEYA9+Ga1Y5R1dtF5FagUFUniMhE4ACgYuD9Zap6gogMBh4DynG90e9T1adqu169E8aLp8HW1fC7z+t+jgby3bKNnPPEZNqmJvLCpQPJSUuKdEjGRJ/yclg3zyWQld+7R1qrZ7sBRsHNOti2D7TvtyeRZHZtVkkkahJGQ6t3wnh4kJvH+8wXQxdUGBUu2cBFz0ylZUIsL1w6kM7Z9szWmFqV7YY1c13yWPUDrJwOq2dBaYnbHp8CbXpD617utc3+0KaXK6E0QZYw6uqOPDeZzHF3hS6oMJu1YjPnj5lCQOCZiwawf07T/J/amLAqK3UlkZXT95RCVs9yw5xUSM37dSLJ7OpadzViIe2H0WyUbHbPOqOwSW1N9s9J5dXfDeL8pyZz6qNfc/epfRnZt32kwzKmcYmJ9ZJAbzjwHLdOFbasdMljzWwvicx2PdXLvZF4Y+Iho4sb1j2rB2R1d+8zu0F802uQYgmjwqaKFlKNK2GAm6nvrasP44oXpvH7sd8zZ9UW/jS0BzEBaxFiTJ2JuBaTqTnQPWj+uNJdrof6mjmuFLJuAfw8C+a+7Tr+uoPdiNdZFYmkm+ugmNUDkhtvHypLGBV+6bQXfX0w/MhumcBLlx3CzW/P5pFJPzF92SbuOb0v7a0y3JjQio13Pc/b7o/rZeAp3Qnrf3LJZN18WDvPPeZa8hWUBs1rk5jm6kqDl8wu7rVFZlQ3/bWEUWFz4y1hVIiPDfB/Jx1Av9w0bn57NsPu+5x/nLg/o/pFZ78SY5qU2ARXOd6m197ry8vd70tFEtmwyC1FU2H260GlElyrrYz8oGTSxfViT+sArdpDILKDj1rCqLBlJQRiIbl1pCOpt9MPzmNg5wyuf/UHrn15Ou/OWMXfR/YiN73pPVM1JuoFAm749vSObi6QYKW7YNOyPUlkw0/uddUPMGcCaFnQeWLdH7RpHV0CSe/ovfc+t2oX9luxhFGheJWbOKmJtL/umJnMK6MP4YkvFnP/f+dzzL2fcfWRXbns8M4kxNoQ6cZEhdh4yOrqlsrKdrtksmmpm9it4v2mZTD/Q9i2Zs++Senw5yXhDzfsV2gstqxwRb4mJDYmwBVDunBCv/b84505/Ouj+Yydspxrj+7Gyf1ziI1pGsnRmCYpJs7VbWR2qXr7ru3uUdemZbCzuEFCsl+MClu8EkYTlJOWxCPnHsQLlwwkKyWe/xk/g2P//TnjpxWxq7S89hMYY6JPfAvX8qrbsbD/yQ1ySUsYsKe9daumXTl8WLcs3rzqUJ44v4CE2AB/HPcDh931CQ99upAN23ZFOjxjTJSzR1LgOuzt3tYglUaRJiIc26sNx+zXms8XrOOpLxdz94fzuH/iAo7t1YZTC3I5vFu29eEwxvyKJQxwpQtocnUYNRERjuiezRHds5m/upiXpyznje+LeHfmKrJS4jm2V1uG9W7D4C5ZxMdaQdQYYwnDqUgYLZtPwgjWvU1L/j6yFzce15NPflzN2zNWMWH6CsZOWUbLhFiO7Nmaw7tnM6hLpo2Ka0wzZgkDmmUJoyrxsQGG79+O4fu3o2R3GV//tI4PZ61m4tzVTPjB/TfKy0jikPxMBnbOpG9uKp2zU+zxlTHNhCUMcH0woMm2kqqLxLgYjurZhqN6tqG8XJm/pphvflrPNz+t56M5qxk3zQ2l0iI+ht7tW7F/Tir7tW1F5+xkumSnkJ4cH+E7MMaEmiUMcH0wkrNdJxrzK4GA0LNtK3q2bcVFh+ZTXq4sXLuVmUWbmbnCLWOnLKNk954muukt4uiSnUKX7BQ6ZLagfVoi7VOTaJ+WRJtWiVYvYkwjFNaEISLDgftxM+49qap3VtqeADwHHASsB85Q1SXetr8AlwBlwDWq+mHYAm3CfTDCIRAQurdpSfc2LTnlIDf2Vlm5UrRxO4vWbuOntVu9ZRv//XE167bu3WRXBLJTEmjTKpGM5Hgyk+PJSI4nPeh9Zko8rRLjSEmMJSUhluT4WAL26MuYiApbwhCRGOAh4FigCJgqIhMqzct9CbBRVbuKyJnAXcAZ3hziZwK9gfbARBHprho8sEoIbVnZqAcdjAYxAaFjZjIdM5M5sufe43Ft31XKyk0lrNq8g1WbSli5eQcrN+1gTfFONmzbxcI1W9mwbRc7dtf89aYkeMkjIYaUxDhaJsSSFB9DQmyAxLhfvybGBUiI3fMaGyPEBgLEBmTP+xghLkaIqbTerRPiYgLEBIQYEQIiSAACIgTEvcLen0VcCzRjmqJwljAGAAtVdRGAiLwMjAKCE8Yo4Gbv/WvAg+L+tY0CXlbVncBiEVnone+bsERavBLyDg7LqQ20iI+la+sUuraueQrZHbvKWL/NJZH123ZRXFLK1pJStu0spXine79152627ixl684yikt2s27rTnaVllOyu4ydQa+l5ZGbSVIEhIpE4pJIRVKRSp8rtgcdvdd5Kq/de92v9w0+VXWJ65d963iuqsJtCimyMSf6jBbxvHr5oLBfJ5wJIwdYHvS5CBhY3T6qWioim4FMb/23lY6tshu2iIwGRgN06NBh36MsL4eux0DeIft+rAmppPgYcuNbhGRU3dKy8r0SSMnuMkrLldIypbS8nN1lSlm5Ulrmkktpebm3TdldVu5t01+2lZUr5QqqiiqUq/tc7k1xXF6+57N625TgdXv2UfacoywosQWnuL1nTtZfrdvrfVXbqzlXxb5Uu6/WcnzN+zZajfwmWiY2THV0OK9SVbqu/LVUt4+fY91K1ceBx8HN6b0vAQJudNqTH9/nw0x0i40JEBsTIDnB2nUYEyrhbKpSBARPX5cLrKxuHxGJBVKBDT6PNcYY04DCmTCmAt1EJF9E4nGV2BMq7TMBuMB7fyrwibpy7gTgTBFJEJF8oBswJYyxGmOMqUXYyutencTVwIe4ZrVjVHW2iNwKFKrqBOAp4HmvUnsDLqng7fcqroK8FLgqbC2kjDHG+CKqjby2J0hBQYEWFhZGOgxjjGk0RGSaqhb42de62xpjjPHFEoYxxhhfLGEYY4zxxRKGMcYYX5pUpbeIrAWW1vHwLGBdCMOJpKZyL03lPsDuJRo1lfuA+t1LR1XN9rNjk0oY9SEihX5bCkS7pnIvTeU+wO4lGjWV+4CGuxd7JGWMMcYXSxjGGGN8sYSxR1MagbCp3EtTuQ+we4lGTeU+oIHuxeowjDHG+GIlDGOMMb5YwjDGGONLs08YIjJcROaJyEIRuTHS8dSHiCwRkZkiMl1EGtUojCIyRkTWiMisoHUZIvKxiCzwXtMjGaNf1dzLzSKywvtupovIiEjG6IeI5InIpyIyV0Rmi8i13vpG973UcC+N8XtJFJEpIvKDdy+3eOvzRWSy97284k0rEdprN+c6DBGJAeYDx+ImbZoKnKWqc2o8MEqJyBKgQFUbXWckETkc2Ao8p6r7e+v+CWxQ1Tu9ZJ6uqn+OZJx+VHMvNwNbVfVfkYxtX4hIO6Cdqn4nIi2BacCJwIU0su+lhns5ncb3vQiQrKpbRSQO+BK4FrgeeF1VXxaRR4EfVPWRUF67uZcwBgALVXWRqu4CXgZGRTimZklVP8fNiRJsFPCs9/5Z3D/wqFfNvTQ6qrpKVb/z3hcDc4EcGuH3UsO9NDrqbPU+xnmLAkcBr3nrw/K9NPeEkQMsD/pcRCP9n8ijwEciMk1ERkc6mBBoo6qrwP2DB1pHOJ76ulpEZniPrKL+MU4wEekEHAhMppF/L5XuBRrh9yIiMSIyHVgDfAz8BGxS1VJvl7D8ljX3hCFVrGvMz+gOVdX+wHHAVd6jERMdHgG6AP2AVcA9kQ3HPxFJAcYDf1DVLZGOpz6quJdG+b2oapmq9gNycU9K9qtqt1Bft7knjCIgL+hzLrAyQrHUm6qu9F7XAG/g/kdqzFZ7z54rnkGviXA8daaqq71/5OXAEzSS78Z7Rj4eeFFVX/dWN8rvpap7aazfSwVV3QRMAg4B0kSkYtrtsPyWNfeEMRXo5rUuiMfNKT4hwjHViYgke5V5iEgyMBSYVfNRUW8CcIH3/gLgrQjGUi8VP7Cek2gE341XufoUMFdV7w3a1Oi+l+rupZF+L9kikua9TwKOwdXJfAqc6u0Wlu+lWbeSAvCa0d0HxABjVPX2CIdUJyLSGVeqAIgFXmpM9yIiY4EhuGGaVwM3AW8CrwIdgGXAaaoa9ZXJ1dzLENxjDwWWAL+rqAeIViJyGPAFMBMo91b/Fffsv1F9LzXcy1k0vu+lD65SOwb3R/+rqnqr9xvwMpABfA+cq6o7Q3rt5p4wjDHG+NPcH0kZY4zxyRKGMcYYXyxhGGOM8cUShjHGGF8sYRhjjPHFEoZpkkQkTUSuDPrcXkReq+mYelwrTkSm1fHYAhF5oB7XvlBE2tf1eGP2hSUM01SlAb8kDFVdqaqn1rB/fRwGfF2XA1W1UFWvqce1LwQsYZgGYQnDNFV3Al28OQ7uFpFOFfNTeH+Vvykib4vIYhG5WkSuF5HvReRbEcnw9usiIh94gzl+ISI9q7nWcOD9yitFZKuI3OUdP1FEBojIJBFZJCInePsMEZF3vPc3ewPgVexzjbf+l9i9z3/y9j0VKABe9O4zSUQOEpHPvGt+GDSExzUiMscbZO/lkP1XNs2KJQzTVN0I/KSq/VT1hiq27w+cjRs76HZgu6oeCHwDnO/t8zjwe1U9CPgT8HA11zoSN55PZcnAJO/4YuAfuLlXTgJureZcPYFhXlw3eeMfVUlVXwMKgXO8gehKgf8Ap3rXHOPdG7j/Hgeqah/g8urOaUxNYmvfxZgm6VNvXoRiEdkMvO2tnwn08UY1HQyMc8MQAZBQ+SRe/cEGVfFzr+MAAAGUSURBVN1exTV2AR8EnXenqu4Wkf/f3t3zQhBFYRz/P7agoBJUEo0otqDV6IXoVHS+g84HEIVWqfABFBSiUClUXjehUpOoVII9inuzhuwyu2jW8+smd2fmTLE5mXMn51wCYy3i2s/tHJ4k3QMjbTzTBCkRHuaYK6QOrAAXpDeRXVLLFbO2OWHYf1XssVMvHNdJ/4se0nyBqW+uMwsctFh7jvfeO417RES90FX0q7hecywvfKwG9LU4V0AtIqabrM0BM8ACsCapWpidYFaKS1LWrR6BgU5PzrMSbiUtQup2KmmyyU+b7l/8sjtgWNKgpF5gvrBWfM4bYEjSNDS+3qpK6gFGI+IIWCV9END/xzFbF3LCsK4UEQ/AsaQrSRsdXmYJWJF0DtT4NL5XaSb8eERc/yzar0XEM2nP4wTYA4r32wa28vS1Cqm99XqO+YxUVqsAO7kUdgps5jkKZm1xt1qzDuWW2csR4U1k+xecMMzMrBSXpMzMrBQnDDMzK8UJw8zMSnHCMDOzUpwwzMysFCcMMzMr5Q3YoHwhTkhExwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112360358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "V = 40     # liters\n",
    "kA = 0.5   # 1/min\n",
    "kB = 0.1   # l/min\n",
    "CAf = 2.0  # moles/liter\n",
    "\n",
    "def batch(X, t):\n",
    "    CA, CB = X\n",
    "    dCA_dt = -kA*CA\n",
    "    dCB_dt = kA*CA - kB*CB\n",
    "    return [dCA_dt, dCB_dt]\n",
    "\n",
    "t = np.linspace(0,30,200)\n",
    "soln = odeint(batch, [CAf,0], t)\n",
    "plt.plot(t, soln)\n",
    "plt.xlabel('time / minutes')\n",
    "plt.ylabel('concentration / moles per liter')\n",
    "plt.title('Batch Reactor')\n",
    "plt.legend(['$C_A$','$C_B$'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the maximum value, we first write a function to compute $C_B$ for any value of time $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CB(t):\n",
    "    soln = odeint(batch, [CAf, 0], [0, t])\n",
    "    return soln[-1][1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We gain use `minimize_scalar` to find the value of $t$ that minimizes the negative value of $C_B(t)$.|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -1.3374806339222158\n",
       "    nfev: 20\n",
       "     nit: 19\n",
       " success: True\n",
       "       x: 4.0235949243406663"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "minimize_scalar(lambda t: -CB(t), bracket=[0,50])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Concentration C_B has maximum 1.33748063392 moles/liter at time 4.02359492434 minutes.\n",
      "Maximum productivity of the batch reactor is 13.296374601 moles per minute.\n"
     ]
    }
   ],
   "source": [
    "tmax = minimize_scalar(lambda t: -CB(t), bracket=[0,50]).x\n",
    "\n",
    "print('Concentration C_B has maximum', CB(tmax), 'moles/liter at time', tmax, 'minutes.')\n",
    "print('Maximum productivity of the batch reactor is', V*CB(tmax)/tmax, 'moles per minute.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Optimization](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/06.00-Optimization.ipynb) | [Contents](toc.ipynb) | [Linear Production Model](http://nbviewer.jupyter.org/github/jckantor/CBE30338/blob/master/notebooks/06.02-Linear-Production-Model.ipynb) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE30338/blob/master/notebooks/06.01-Unconstrained-Scalar-Optimization.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/06.01-Unconstrained-Scalar-Optimization.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
}