notebooks/06.02-Maximizing-Concentration-of-an-Intermediate-in-a-Batch-Reactor.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A4n0I0WGFg2I"
   },
   "source": [
    "\n",
    "# Maximizing Concentration of an Intermediate in a Batch Reactor\n",
    "\n",
    "Keywords: ipopt usage, scipy.minimize_scalar, scipy.odeint, unconstrained optimization\n",
    "\n",
    "This notebook presents an example of the finding the time period required to achieve an optimal result. Because the period of operation is an unknown "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Statement\n",
    "\n",
    "A desired product $B$ is as an intermediate in a series reactions\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "A \\overset{k_A}{\\longrightarrow} B \\overset{k_B}{\\longrightarrow} C\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $A$ is a raw material and $C$ is an undesired by-product. The reaction operates isothermally with rate constants $k_A = 0.5\\ \\mbox{min}^{-1}$ and $k_B = 0.1\\ \\mbox{min}^{-1}$. The raw material is a solution with concentration $C_{A,f} = 2.0\\ \\mbox{moles/liter}$.\n",
    "\n",
    "A 100 liter tank is available for use as a batch reactor. How long should the reaction be operated to maximize the concentration of $B$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "qceXwUrFIKB8"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "from scipy.optimize import minimize_scalar\n",
    "\n",
    "import shutil\n",
    "import sys\n",
    "import os.path\n",
    "\n",
    "if not shutil.which(\"pyomo\"):\n",
    "    !pip install -q pyomo\n",
    "    assert(shutil.which(\"pyomo\"))\n",
    "\n",
    "if not (shutil.which(\"ipopt\") or os.path.isfile(\"ipopt\")):\n",
    "    if \"google.colab\" in sys.modules:\n",
    "        !wget -N -q \"https://ampl.com/dl/open/ipopt/ipopt-linux64.zip\"\n",
    "        !unzip -o -q ipopt-linux64\n",
    "    else:\n",
    "        try:\n",
    "            !conda install -c conda-forge ipopt \n",
    "        except:\n",
    "            pass\n",
    "\n",
    "assert(shutil.which(\"ipopt\") or os.path.isfile(\"ipopt\"))\n",
    "from pyomo.environ import *\n",
    "from pyomo.dae import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NdetctYWFg2L"
   },
   "source": [
    "## Mathematical model\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",
    "$$\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",
    "\n",
    "Eliminating the common factor $V$\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{dC_A}{dt} & = - k_A C_A \\\\\n",
    "\\frac{dC_B}{dt} & = k_A C_A - k_B C_B\n",
    "\\end{align}\n",
    "$$\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": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 295
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 6102,
     "status": "ok",
     "timestamp": 1559058898921,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 300
    },
    "id": "9H59ye_PFg2M",
    "outputId": "8e2e5fba-debe-418e-93d4-3134343cb664"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABEoElEQVR4nO3dd3yUVfb48c9JJ70SWuhNUHqxLYINLNi7Ypd11V1dv+6uuz93dbuuq6uuvWBXLFiwl5XYRUA6SJFmaKGTAElIcn5/3Cc4hEkyKZOZTM779XpeM/PUcxmdk+fe+9wrqooxxhhTXVSoAzDGGBOeLEEYY4zxyxKEMcYYvyxBGGOM8csShDHGGL8sQRhjjPHLEoQx9SAiq0Tk2FDHYUxzsARhWjzvR3uPiBSLyDYReUdE8gI8tquIqIjEBCGuS0Wkwotrp4jMFZGTm/o61a45WkQKgnkN03pYgjCRYryqJgPtgY3Af0McT5WvvbjSgQeBySKSHtKIahGMRGlaLksQJqKoagnwKtCvap2InCQis72/4n8Ukdt8DvnMe93u/aV/mHfMVSKyWESKRGSRiAzxOWaQiMwTkR0i8pKIJAQQVyXwLJAE9PKuES8i/xaRNSKyUUQeFpE23rYMEXlbRDZ5d0Vvi0gnnzJlisiTIrLO2/6GiCQB7wEdvLIUi0gH7zr3ePuu897He+cZLSIFIvI7EdkAPFnff3MTuSxBmIgiIonAucA3Pqt3ARfj/oo/CfiFiJzmbRvlvaararKqfi0iZwO3ecekAqcAW3zOdw4wDugGDAAuDSCuaOAyYC+w2lt9B9AbGAT0BDoCf/K2ReF+rLsAnYE9wP0+p3wWSAT6A22B/6jqLuAEYJ1XlmRVXQf8P+BQ7zoDgRHALT7nagdketeaWFdZTCuiqrbY0qIXYBVQDGwHyoF1wCG17H8P7gcVoCugQIzP9g+A62u51kU+n/8FPFzDvpd68WzHJYY9wDneNsElrh4++x8GrKzhXIOAbd779kAlkOFnv9FAQbV1PwAn+nweC6zy2b8MSAj192hL+C12B2EixWmqmg7EA9cBn4pIOwARGSki07zqmh3A1UB2LefKw/2o1mSDz/vdQHIt+37jxZUBTAV+5q3Pwd0BzBKR7SKyHXjfW4+IJIrIIyKyWkR24qrC0r07kTxgq6puq+W6vjrw010L3vsOPp83qauaM2Y/liBMRFHVClV9DagAjvRWv4D7cc5T1TTgYdxf8ODuHqr7EejRxHEVA9cAE0RkMLAZd0fRX1XTvSVNXYM2wP8BfYCRqprKT1Vh4sWXWUNjt7/yrMNVH1Xp7K2r7RhjLEGYyCLOqbi/2Bd7q1Nwf3GXiMgI4AKfQzbhqmu6+6x7HLhJRIZ65+spIr4/sA2iqlu8c/9JXaP1Y8B/RKStF3tHERnrE/MeXON5JnCrz3nW4xqjH/Qas2NFpCqBbASyRCTN59IvAreISI6IZOPaOZ5rbHlM5LMEYSLFWyJSDOwE/g5coqoLvW3XAH8RkSLcj+PLVQep6m5v/y+9qp5DVfUVb90LQBHwBq4RtyncA5woIgOA3wHLgW+8aqSPcXcNVfu1wd1pfIOrfvI1Adeu8T1QCNzgled7XEJY4ZWnA/A3YCYwD5gPfOetM6ZWomp3l8YYYw5kdxDGGGP8sgRhjDHGL0sQxhhj/LIEYYwxxq+IGpgrOztbu3bt2qBjd+3aRVJSUtMGFCKRUpZIKQdYWcJRpJQDGleWWbNmbVbVHH/bIipBdO3alZkzZzbo2Pz8fEaPHt20AYVIpJQlUsoBVpZwFCnlgMaVRURW17TNqpiMMcb4ZQnCGGOMX5YgjDHG+BVRbRDGGBMse/fupaCggJKS8Bv4Ni0tjcWLF9e6T0JCAp06dSI2Njbg81qCMMaYABQUFJCSkkLXrl0RkboPaEZFRUWkpKTUuF1V2bJlCwUFBXTr1i3g8watiklE8rwx+BeLyEIRud7PPiIi94nIcm8KxyE+28aJyBJv283BitMYYwJRUlJCVlZW2CWHQIgIWVlZ9b77CWYbRDnwf6p6EG66w2tFpF+1fU7Azc/bCzfV4UOwb3rGB7zt/YDz/RxrjDHNqiUmhyoNiT1oCUJV16vqd977ItzY/B2r7XYq8Iw63+BmzGqPmzN3uaquUNUyYLK3b5OrrFTu/2QZ8zeVB+P0xhjTYjVLG4SIdAUGA9OrbeqImx2rSoG3zt/6kTWceyLeROu5ubnk5+fXO74HPtnF8Bxt0LHhqLi4OCLKEinlACtLOKpvOdLS0igqKgpeQI1QUVERUGwlJSX1KnPQE4SIJANTgBtUdWf1zX4O0VrWH7hS9VHgUYBhw4ZpQ54m7Dz7M4oqdttTlWEmUsoBVpZwVN9yLF68uNaG4FCqq5G6SkJCAoMHDw74vEF9DkJEYnHJ4XlvnuDqCnATsFfphJsrt6b1QdEuLYFtJTZxkjGmZXjllVcYOXIkAwYMoGfPnvzzn/8MynWC2YtJgCeAxap6dw27TQUu9nozHQrs8ObbnQH0EpFuIhIHnOftGxTt0xLYagnCGNMCPP3009xxxx1MmTKFefPmMWfOHBITE4NyrWBWMR2Bmzd3vojM8db9AegMoKoPA+8CJ+Lm5d0NXOZtKxeR64APgGhgks/8wk2ufVobdpYppeUVxMdEB+syxhjTKDt37uTGG29kxowZdOrUCYDk5GSuv/6ApwiaRNAShKp+gf+2BN99FLi2hm3v4hJI0LVPSwCgcGcpeZnBycTGmMjx57cWsmhd9SbVxunXIZVbx/evdZ/XX3+dkSNH0r179ya9dk1sLCZcGwTA+h3h9wi9McZUWbhwIYMGDapx+7nnnstdd93VZNezoTb46Q5i/Y49IY7EGNMS1PWXfrAkJSWxZ4//36k333yTk08+mY8//rjJrmd3EED79DaA3UEYY8LbiSeeyCuvvMLGjRsBKC0t5bHHHqOkpIRXXnmFCRMmsGPHjia7nt1BAMnxMbSJgQ2WIIwxYWz48OHcdtttjB07loqKCsrLy7nooou49957KS4u5uqrr2bhwoXs2bOHNm3aNPp6liA8GQliVUzGmLA3YcIEJkyYsO/zmjVruOWWW3jjjTcA+POf/8y8efMYOdLv4BP1YgnCk5kQZVVMxpgWp3PnzjzwwAP7Pt96661Ndm5rg/BkJoglCGOM8WEJwpMRL2wuLqWsvDLUoRhjTFiwBOHJTBBUYeNOu4swxhiwBLFPZoJ76HuDJQhjjAEsQeyTmeD+KawdwhhjHEsQngzvDmL9duvqaowxYAlin8RYITk+xu4gjDHGYwnCR7u0BHua2hhjPJYgfLRPS7CnqY0xxmMJwodLEHYHYYwJb5Ew5egkESkUkQU1bP+NiMzxlgUiUiEimd62VSIy39s2M1gxVtcurQ2bikvZW2EPyxljwlNzTjkazDuIp4BxNW1U1TtVdZCqDgJ+D3yqqlt9dhnjbR8WxBj30z4tAVUoLCptrksaY0zAqqYcffnll1v8lKOfiUjXAHc/H3gxWLEEat/Mctv30DG98UPlGmMi1Hs3w4b5TXvOdofACbfXuktzTzlaa4IQkSjgUFX9KlgBiEgi7k7jOp/VCnwoIgo8oqqP1nL8RGAiQG5uLvn5+Q2Ko7i4mO1LXW3Y/77+juJVLXeg2+Li4gb/O4STSCkHWFnCUX3LkZaWRlFREQDxe8uIqihv0ngq95ZR6p2/JrNnz6Zfv3774qhSUVHBkCFDGDp0KIsXL+buu+/moIMOOuD4kpKSepW51l9BVa0UkbuAwwI+Y/2NB76sVr10hKquE5G2wEci8r2qflZDjI8CjwIMGzZMR48e3aAg8vPz+dnII7jlyw/J7Nid0aOaJ0MHQ35+Pg39dwgnkVIOsLKEo/qWY/HixaSkpLgPp9wdlJji6tiekZHBnj17forDJ7ZDDz2UBx98kLvvvptt27YdsA9AQkICgwcPDjieQNogPhSRM0VEAj5r/ZxHteolVV3nvRYCrwMjgnTt/aQmxJAYF8066+pqjAlDNU05OmfOHJYuXcrll1/OJ598wtixY5vkeoEkiBuBV4AyEdkpIkUisrMpLi4iacBRwJs+65JEJKXqPXA84LcnVFMTETqmt2HtNksQxpjw4zvl6CGHHMKgQYMoLCxk9uzZ3HXXXUyaNIn4+Hh27drVJNers6JdVQ+8TwmAiLwIjAayRaQAuBWI9c75sLfb6cCHqupbmlzgde+GJQZ4QVXfb0gMDZGXmciPliCMMWGq+pSjAMcffzxFRUVERUVx0EEHkZSU1CTXqjNBeFVLFwLdVPWvIpIHtFfVb2s7TlXPr+vcqvoUrjus77oVwMC6jg2WThltmLFyK6pK8GrVjDGm6UyZMsVvm0NjBVLF9CCukfoC73Mx8EDNu7dseRmJFJWWs2PP3lCHYowxIRVIghipqtcCJQCquo26G9tbrLxM9/zDj1utmskY07oFkiD2ikg07tkERCQHiNixKDpluEfWf9y2O8SRGGNMaAWSIO7DdTVtKyJ/B74AgjMyVBjIy/QSxFZLEMaY/alqqENosIbEHkgvpudFZBZwDCDAaaq6uP7htQxpbWJJTYixOwhjzH4SEhLYsmULWVlZLa4Di6qyZcsWEhIS6nVcIL2YnlXVCcD3ftZFpLzMRGuDMMbsp1OnThQUFLBp06ZQh3KAkpKSOn/8ExIS9g3wF6hABhzq7/vBa48YWq+rtDB5GYksK6x9TBRjTOsSGxtLt27dQh2GX/n5+fUaQiNQNbZBiMjvRaQIGODzBHURUIjPk8+RKC+zDQXb9rTo+kZjjGmsGhOEqv7Te4r6TlVNVdUUb8lS1d83Y4zNLi8zkdLySjbZvBDGmFasxiomEemrqt8Dr4jIkOrbVfW7oEYWQnk+XV3bptavUccYYyJFbW0Q/wdcBdzlZ5sCRwclojDg+7Dc0C4hDsYYY0KkxgShqld5r2OaL5zwsO9hOXsWwhjTitVWxXRGbQeq6mtNH054SIiNJicl3p6FMMa0arVVMY2vZZsCEZsgAPIy2tizEMaYVq22KqbLmjOQcJOXmcis1dtCHYYxxoRMIGMxtUqdMtqwfkcJ5RUROy6hMcbUKmgJQkQmiUihiPidLlRERovIDhGZ4y1/8tk2TkSWiMhyEbk5WDHWJi8jkYpKZf2OklBc3hhjQq7WBCEiUSJyeAPP/RQwro59PlfVQd7yF++a0bgJiU4A+gHni0i/BsbQYPtGdbWGamNMK1VrglDVSvw/B1EnVf0M2NqAQ0cAy1V1haqWAZOBUxsSQ2NUPSxXYA3VxphWKpDB+j4UkTOB17TpByc6TETmAuuAm1R1IdAR+NFnnwJgZE0nEJGJwESA3Nxc8vPzGxRIcXHxfseWVyoCfD5nMW13/dCgc4ZK9bK0VJFSDrCyhKNIKQcEryyBJIgbgSSgQkT24OaEUFVNbeS1vwO6qGqxiJwIvAH08s5fXY2JSVUfBR4FGDZsmI4ePbpBweTn51P92A7ffkJUSgajRzf9KInB5K8sLVGklAOsLOEoUsoBwStLnY3U3gB9Uaoa6zNoX2OTA6q6U1WLvffvArEiko27Y8jz2bUT7g6j2XXLTmLVFmuDMMa0TnUmCHEuEpE/ep/zRGREYy8sIu3Em5bJO18UsAWYAfQSkW4iEgecB0xt7PUaolt2Eis2Fduw38aYVimQKqYHgUrc4Hx/BYpxvYyG13aQiLwIjAayRaQAuBWIBVDVh4GzgF+ISDmwBzjPa+MoF5HrgA+AaGCS1zbR7LrnJFFUUs6WXWVkJ8eHIgRjjAmZQBLESFUdIiKzAVR1m/eXfa1U9fw6tt8P3F/DtneBdwOILai6ZScBsHLzLksQxphWJ5AH5fZ6zyYogIjk4O4oIl737GQAVmwqDnEkxhjT/AJJEPcBrwO5IvJ34AvgH0GNKkx0zGhDXHQUKzbvCnUoxhjT7OqsYlLV50VkFnCMt+o0VV0c3LDCQ3SU0CUrkZWbLEEYY1qfQNogABJxDcYKtAleOOGnW3YSK+0OwhjTCgXSzfVPwNNAJpANPCkitwQ7sHDRLSeJ1Vt2U1FpXV2NMa1LIHcQ5wODVbUEQERuxz0F/bdgBhYuemQnU1ZRydpte+iclRjqcIwxptkE0ki9Ckjw+RwPtKzBiRqhW47r6rpis/VkMsa0LoEkiFJgoYg8JSJPAguAYhG5T0TuC254oVf1LMQKa6g2xrQygVQxve4tVfKDE0p4ykqKIzUhJvQN1ZUVgECUTQJojGkegXRzfbo5AglXIkK3nOTmTRDbVsGiqbD6K9i0GHashcq9IFHQJhOyekC7AdDjaOj2M4hPab7YjDGtRqDdXFu17tlJTF+xJfgXWv0VfPovWDHNfc7qCR2HQr9TITYJKspg1ybYshzmvAAzHoOoWOh6JAy9BPqeDNGxwY/TGNMqWIIIQLfsJF6fvZY9ZRW0iYtu+gsUb4L3fwcLpkBSWzj6FhhwLqR3rvmY8lL4cTos/xgWvAavXArJ7WDYZcTsPbjpYzTGtDr1ShAiEgUkq+rOIMUTlrp7PZlWbdnFQe0bPRXG/tZ8Ay9NgJLtMPoPcPgvIS6A7rQx8dBtlFuOuRWWfQQzn4D8fzIyJgnib4SRV0N8ctPGa4xpNQJ5UO4FEUkVkSRgEbBERH4T/NDCh++ork1qwRR46mT3Iz7xUxj9u8CSQ3VR0dBnHFz4Cvz8c3ak9YNP/gr3DYI5L4LNZ2GMaYBAusT08+4YTsMNwd0ZmBDMoMLNT11dm/BZiPmvwpQrodNwuOoTyO3XNOdtP4AFh9wCV3wMGV3hjavhyRNh46KmOb8xptUIJEHEikgsLkG8qap7qWWO6CoiMklECkVkQQ3bLxSRed7ylYgM9Nm2SkTmi8gcEZkZYFmCJjEuhvZpCU03qusP0+C1idDlCLjoVWiT0TTn9ZU3HC7/EMbf53pCPfIzyL8dKvY2/bWMMREpkATxCO5p6iTgMxHpAgTSBvEUMK6W7SuBo1R1AG6mukerbR+jqoNUdVgA1wq6HjnJLC9sgjuIzcvg5Usgpw+c9wLEJTX+nDWJinK9m66bBf1Ph/x/wuPHQGGrGIzXGNNIdSYIVb1PVTuq6onqrAbGBHDcZ8DWWrZ/parbvI/fAJ0CDToUeuemsGxjMZWNGbRvbwm8chlEx8AFL0FCEzd41yQpC858HM55BnYUwCOj4NvHrG3CGFMr0Tp+JEQkFzdBUAdVPUFE+gGHqeoTdZ5cpCvwtqrW2u9SRG4C+qrqld7nlcA2XFXWI6pa/e7C99iJwESA3NzcoZMnT64rLL+Ki4tJTq65x8+nBXt5ckEZ/xrVhraJDXuaucfyx8kreIt5h/yRrVnBuzGqrSyxZdvp+/19ZG2dRWHO4Szpcx0VMUG8i2mEur6TlsTKEn4ipRzQuLKMGTNmVo01Napa6wK8B5wDzPU+xwDz6zrO27crsKCOfcYAi4Esn3UdvNe2wFxgVCDXGzp0qDbUtGnTat0+e8027fK7t/X9BesbdoHV36jemqr6zk0NO74e6iqLVlSofnGP6m0Zqv85RLVgVtBjaog6y9GCWFnCT6SUQ7VxZQFmag2/qYH8KZytqi/jzUOtquVARUMyVXUiMgB4HDhVVfc9qqyq67zXQtw4UCOa4nqN0auty85LNxTV/+CKvfD2DZDayT2zEGpRUXDE9XD5+26MpyeOh+mPWJWTMWY/gSSIXSKShddzSUQOBXY09sIi0hl4DZigqkt91ieJSErVe+B43AiyIZUUH0PnzES+39iABDH9EShcBCfeGV4PruWNgKs/h57Hwnu/hTevde0kxhhDYE9S3whMBXqIyJdADnBWXQeJyIvAaCBbRAqAW4FYAFV9GPgTkAU8KCIA5erqwXKB1711McALqvp+/YoVHL1zU+p/B7FnO3x2J/Q4BvqeGJS4GiUx0/Wm+vQO+PR22LQEzn0OUtuHOjJjTIgFMprrdyJyFNAHEGCJumch6jru/Dq2Xwlc6Wf9CmDggUeEXt92KUxbUkhpeQXxMQGOyfTlPW4YjWNvC2JkjRQVBWN+D7n94fWr4dHRcN7z0CksehgbY0KkxgQhImfUsKm3iKCqrwUpprDVu10KFZXKik0BjslUvAm+eRgOORvaDwh+gI3V7xQ3lPiL58OTJ8DJ98DgC0MdlTEmRGq7gxhfyzbFtR+0Kn3buXkXlmwoCixBTH8IykvgqN8FObImlNsfJua70WHfvAY2zIexf3fjPRljWpUaE4SqXtacgbQE3bKTiI0WlgTSUF2yE7593P1Vnt0r+ME1pcRMuOg1+PAWl+S2rYQznwivBnZjTNAFMpprmojcLSIzveUuEUlrjuDCTWx0FD1ykgNrqJ71JJTugCN/HfzAgiE6Bk64HU66yw0l/uQ4N7OdMabVCKSb6ySgCPew3Dm4cZieDGZQ4ax3bgrf15UgKitgxuPQ5UjoMLh5AguW4VfCBS/D1lVuHKf1c0MdkTGmmQSSIHqo6q2qusJb/gx0D3Zg4apPuxTWbt9DUUktHbmWfwzb18CIAzpptUy9joUrPgCJhkknwJL3Qh2RMaYZBJIg9ojIkVUfROQIYE/wQgpvfXJdQ/XSjbWM7DrjcUjOdXNER4rc/nDV/yCnN0y+wPXOMsZEtEASxC+AB7w5GlYD9wNXBzes8NWnXVWCqKGaafuPrs5+yCUQHduMkTWDlHZw6TvQ50Q3h/a7v4GK8lBHZYwJkkAelJsDDBSRVO9zq5qPurqO6W1IiotmSU3tEPNeAjRynx+IS4JznoWP/wRf/Re2rYKznrQeTsZEoDoThIikAxfjRmaN8YbAQFV/FczAwlVUlNC7XQrfb/CTJ1Vdguh8uJvuM1JFRcHxf4PM7vDOTa6H0wUvQ2qHUEdmjGlCgVQxvYtLDvOBWT5Lq9W/QyoL1+08cPKgdd/B5qUw8LzQBNbchl3u9XBaCY8d4x6qM8ZEjEASRIKq3qiqT6rq01VL0CMLY4d0TKOopJzVW3fvv2HuSxAdD/1ODU1godDrWDdsOMCkcbDs49DGY4xpMoEkiGdF5CoRaS8imVVL0CMLY4d0TAdgXsH2n1ZWVsLiqdDrOGiTHoqwQqfdIa6HU2Y3eOEcmDkp1BEZY5pAIAmiDLgT+JqfqpdmBjOocNcrN5n4mCgWrPWZFmPtTCha37ruHnyldoDL3oOex8Dbv4YP/+iSpjGmxQp0Poieqro52MG0FLHRUfTrkMq8Ap8EsehNiIqF3mNDF1ioxafAeS+6LrBf3ed6OJ3xKMS2CXVkxpgGCOQOYiGwu869WplDOqaxYO0O11CtCoumQo8xkNAqh6n6SXQMnPhvOP7vsPgteHq8G/bcGNPiBJIgKoA5IvKIiNxXtdR1kIhMEpFCEfE7Xag494nIchGZJyJDfLaNE5El3rabAy9O8zmkYxq7yipYsXmX672zYw0cVNsI6a2ICBx+HZz7LGxY4MZw2rS07uOMMWElkATxBvB34Cvq1831KWBcLdtPAHp5y0TgIQARiQYe8Lb3A84XkX4BXK9ZDeiUDuDaIZZ94Fb2asXVS/4cNN49eb13NzxxLKz6ItQRGWPqIZAnqRvUpVVVPxORrrXscirwjKoq8I2IpItIe9wzF8u9qUcRkcnevosaEkew9MhJIiE2inkFOzit8CM3amtKbqjDCj+dhsKVH8Pz58Azp8Gp97ee50SMaeECaaQOlo7Ajz6fC7x1/taPrOkkIjIRdwdCbm4u+fn5DQqmuLi43sd2SoK5CxegJTNY3eVsVjXw2k2tIWUJtpg+f6L/wtvJeP3nrJydz+ou57qqqFqEYzkaysoSfiKlHBC8soQyQfj7ddBa1vulqo8CjwIMGzZMR48e3aBg8vPzqe+x+TsXUjzzRSSqkq7HTaRrp2ENunZTa0hZmsWY4+HtG+g253m6pQHj74OYuBp3D9tyNICVJfxESjkgeGWpsQ1CRH4vIsGc7aYAyPP53AlYV8v6sDOgUxpH6HeUJ2RChyF1H9DaxcTBqQ/AmFtg7ovw3BmwZ1uoozLG1KC2RuqVwPUiMltEnhKRc0UkowmvPRW42OvNdCiwQ1XXAzOAXiLSTUTigPO8fcPOIR1SOSJqIeszR7oB7EzdROCo38AZj8GP0+GJ42HrilBHZYzxo8YqJlWdDEwG8O4kxgGveb2MPgbeV9VvazpeRF4ERgPZIlIA3ArEeud+GDcI4InActxzFpd528pF5DrgAyAamKSqCxtXzODoLmuJlu1MjRm43y2PCcCAcyC1I7x0ITx2NJzzDHQbFeqojDE+AmqDUNXZwGzgn968EMcBVwI1JghVPb+OcypwbQ3b3sUlkLAWvepzAN7Z1ZtTQhxLi9T1CLjyf/Di+fDs6XDCv2D4FaGOyhjjqXe9iKruVNUpqjoxGAG1KCs/ZVt8B6ZtSKRkb0Woo2mZsnrAlR9Bj6PhnRvd/BIVtcz3bYxpNlZx3lCVFbDqc3Z1OJyyikrm+w7cZ+onIQ3OnwyH/xJmPAbPnQm7t4Y6KmNaPUsQDbVhHpTsIL3/cQDMWGU/aI0SFe1mqTvtIVjzNTx+DIm7fqz7OGNM0ATUBiEiHYEuvvur6mfBCqpFWPMNAMm9fkbPtsuZucq6azaJQRdAZg946UKGfPdb6N3OzbFhjGl2dd5BiMgdwJfALcBvvOWmIMcV/tZ8DemdIa0jw7tmMHPV1gOnIDUN03kkXDWNPW1y4fmz4bN/uxFzjTHNKpAqptOAPqp6oqqO95bW3WlH1d1BdD4MgGFdMtlZUs7SwqIQBxZB0vOYPfh2OPgM+OSv8NJFULIz1FEZ06oEkiBW4D2/YDzbVkLxRshzQ0QN7+pmYJ1h1UxNqjI6Ac58Asb+E5a8556X2LQk1GEZ02oEkiB204D5ICKa1/5QdQeRl9mGtinxzLSG6qYnAoddA5dMhZLtLkksCssH642JOIEkiKnAX6n/fBCRa83XrmtmTl8ARIThXTOtoTqYuh4JEz91/+YvT4CPb3NdjY0xQRPQfBDemEi9vVVLVLV1P8lUMBM6jdhv/KVhXTN4Z/561m7fQ8d0m4M5KNI6wmXvwnu/gy/+A+tmw5mTICkr1JEZE5EC6cU0GliGm+XtQWCpiLTeQXNKi6BwMVQb2ruqHcKqmYIsJh7G3wOn/BdWfwWPjII100MdlTERKZAqpruA41X1KFUdBYwF/hPcsMLY+rmAHjC8d992KSTFRdsDc81lyMVw+QcQHQNPngBf3AOVlaGOypiIEkiCiFXVfV1HVHUprblX01qv+aXj/gkiJjqKIV0ymLHS2iGaTcch8PPP4KCT4eNb4YVzYNeWUEdlTMQIJEHMFJEnRGS0tzxGa26kXjsL0rtAUvYBmw7tnsWSjUUUFpWEILBWKiENzn4aTvw3rPwUHj4SVn8d6qiMiQiBJIhfAAuBXwHXA4uAq4MZVFhb+x10HOp306heOQB8uXxzc0ZkRGDEVXDFR66N4qmT4PO7rcrJmEaqM0Goaqmq3q2qZ6jq6ar6H1UtbY7gwk5xIez4scYE0b9DKhmJsXy+zBJESHQY5Kqc+p0C//szvHA2FG8KdVTGtFi1zUn9svc6X0TmVV8CObmIjBORJSKyXERu9rP9NyIyx1sWiEiFiGR621Z5154jIjMbWsAmtW62e63W/lAlKko4omc2ny/bjNrYQaGRkApnPQkn3Q0rP4eHDoOlH4Q6KmNapNqeg7jeez25ISf2piZ9ADf7XAEwQ0Smquqiqn1U9U7gTm//8cCvVdW3G9AYVQ2fP8fXe3mx3SE17jKqVw5vz1vPko1F9G2X2kyBmf2IuJnpOh8Gr13lGq+HXeGGE49LDHV0xrQYNd5BqOp67+01qrradwGuCeDcI4DlqrpCVctw81ufWsv+5wMvBhp4SGyYC5ndIT6lxl2O7OUarz9fGj55rdXK7QdXfQKHXQczn3DPTFTdBRpj6iR1VYWIyHeqOqTaunmqOqCO484Cxqnqld7nCcBIVb3Oz76JuLuMnlV3ECKyEtgGKPCIqj5aw3UmAhMBcnNzh06ePLnW8tSkuLiY5OTkWvcZ+c1EilJ6sqj/b2vd7w+f7yYzIYqbhic0KJbGCqQsLUFTliN921wOWnwvsXu3s6rrBazpfDpIdJOcOxCR8p1A5JQlUsoBjSvLmDFjZqnqML8bVdXvguu9NB/YBczzWVYCz9V0nM/xZwOP+3yeAPy3hn3PBd6qtq6D99oWmAuMquuaQ4cO1YaaNm1a7Tvs3qZ6a6rqZ/+u81y3TV2gvf/fu7qnrLzB8TRGnWVpIZq8HLu2qL58ifsenxinunVV056/FpHynahGTlkipRyqjSsLMFNr+E2trRfTC8B43GB9432Woap6UQCJqQDI8/ncCVhXw77nUa16SVXXea+FwOu4KqvQ2bjAvbYbWOeuo3rlUFpeaU9Vh5vETNeAffojsGE+PHQEzHzSJiMypga1tUHsUNVVqnq+unaHPbjqnmQR6RzAuWcAvUSkmzfY33m4ZLMfEUkDjgLe9FmXJCIpVe+B44EF9ShX0wuggbrKyO6ZxEaLdXcNRyIw8Dz4xZfQcTC8fQM8cypsWx3qyIwJO4EM1jdeRJbhqpY+BVYB79V1nKqWA9cBHwCLgZdVdaGIXC0ivg/anQ58qKq7fNblAl+IyFzgW+AdVX0/wDIFx4Z5kJwLKbl17poYF8OwLpl8ttT64IetjC5w8VQ4+T/u6fgHD4NvH7OH64zxEciT1H8DDgWWqmo34BjcHNV1UtV3VbW3qvZQ1b976x5W1Yd99nlKVc+rdtwKVR3oLf2rjg2pDQsCunuoMrpPDt9vKKJg2+4gBmUaRQSGXQ7XfOPmwX73Jnh6PGxdEerIjAkLgSSIvaq6BYgSkShVnQYMCm5YYaaiHDYvgdz+AR9yfP92AHy4cGOwojJNJT0PLnoNTrnftU08eDh885DdTZhWL5AEsV1EkoHPgOdF5F6gPLhhhZmtK6CiDNr2C/iQbtlJ9MlN4f2FG4IYmGkyIjBkAlz7DXQbBe/fDJOOdwnDmFYqkARxKm5e6l8D7wM/4HoztR6F3sPfbQ+q12Fj++cyc9VWthS3zqGrWqTUDnDBS3D6o7B1JTxyFHzw/6C0ONSRGdPsak0Q3nAZb6pqpaqWq+rTqnqfV+XUehQuAomC7N517+vj+P7tqFT4eLFVM7UoIjDwXLhuhrur+Pp+eGAELH471JEZ06xqTRCqWgHs9rqitl6Fi9wQG7H1m2u6f4dUOmW04QNrh2iZEjNh/L1uGPE2GfDShfDCebB9TagjM6ZZBFLFVALM9yYNuq9qCXZgYaVwcb2rlwBEhLH92/HFss0UlewNQmCmWeSNgImfusH+Vn4GD4x0U5yWl4U6MmOCKpAE8Q7wR1wj9SxvCY/ht5vD3j2ukboeDdS+xvZvR1lFJflL7JmIFi06Bg7/JVw7HXoc7aY4fegwWPphqCMzJmgCSRDpXtvDvgXICHZgYWPzUtDKBt1BAAztkkF2chwfWG+myJCeB+c9Dxe+CoiblOi5M2HTkjoPNaalCSRBXOJn3aVNHEf4KvzeveY0LEFERwnH9ctl2veFlOytaMLATEj1Og6u+RrG/hN+nOGexH7vZtizLdSRGdNkaptR7nwReQvoJiJTfZZpQOvpxbR5qRsWOrN7g08xfkAHdpVV8OEia6yOKNGxcNg18KvvYOgl8O0jcN8QmPG4e7jSmBauthnlvgLWA9nAXT7ri3DDfrcOm5dCZjeIiWvwKQ7tnkWHtASmzCrglIEdmjA4ExaSst2YTsMuh/d/D+/8H3z7OBx7G/Qe67rNGtMC1ZggvBFcVwOHNV84YWjzsno//1BdVJRw+pCOPJT/A4U7S2ibGpqJhEyQtTsELnkLFr8FH98GL57rpj099s+hjsyYBglkNNczRGSZiOwQkZ0iUiQiO5sjuJCrKIetP0B2r0af6owhnahUeGPO2iYIzIQtEeh3iuvtdNLdrgfcpOM5eP4/fmrPMqaFCKSR+l/AKaqapqqpqpqiqqnBDiwsbF/txmBq5B0EQI+cZAblpTNl1tqqGfNMJIuOheFXwK9mw9G3kL59vusW+8a1sKMg1NEZE5BAEsRGVV0c9EjC0eZl7rUJEgTAmUM7sWRjEQvXtY4bMAPEJcGo3zB95CNw6DUw/2XXkP3+H6DIOi2Y8BZIgpgpIi95vZrOqFoCObmIjBORJSKyXERu9rN9tFd1Ncdb/hTosc1i81L3mtWzSU43fkB74qKjmPKd/QXZ2uyNS4Wxf4dfzoJDzoLpD8O9AyxRmLAWSIJIxY3mejw/zUt9cl0HeQP9PQCcAPQDzhcRf48jf66qg7zlL/U8Nrg2L4WkHDcmTxNIT4zjmIPaMnXOOsrKba6BVim9M5z2oBsI8OAzfRLF76HIHqY04aXOBKGql/lZLg/g3COA5d7scGXAZNzQ4YFozLFNpwl6MFV3zrA8tuwqs3kiWrusHtUSxSNw70BLFCasSF0NpiLSG3gIyFXVg0VkAK7R+m91HHcWME5Vr/Q+TwBGqup1PvuMBqYABcA64CZv3uo6j/U5x0RgIkBubu7QyZMnB1Tw6oqLi0lOTt5v3eFfXszm7BEs7XPAZRusUpXff76HlDjhlkPrNzpsoPyVpSWKlHJA3WVps3s9nde8TLsN+VRGxbCh3TH8mHcaJW3aNWOUgYmU7yVSygGNK8uYMWNmqeowvxtVtdYF+BT3F/1sn3ULAjjubOBxn88TgP9W2ycVSPbenwgsC/RYf8vQoUO1oaZNm7b/ij3bVW9NVf38Pw0+Z00mfbFCu/zubZ2zZluTn1vVT1laqEgph2o9yrLlB9U3r1P9S7bqbemqL1+qum5OUGOrr0j5XiKlHKqNKwswU2v4TQ2kDSJRVb+tti6QcQQKgDyfz51wdwm+yWmnqhZ7798FYkUkO5Bjg27rSvfaiCE2anLW0E4kxUXz9FermvzcpoXL7A6n/Beun+dGj13+MTwyCp45DVbkg3WRNs0okASxWUR6AAr7qo7WB3DcDKCXiHQTkTjgPGCq7w4i0k7EjUMgIiO8eLYEcmzQbf3BvQYhQaQkxHL2sDzemreOwqKSJj+/iQCp7eG4v8CvF7ghOwoXwTOnwqOjYcEUqLD5RUzwBZIgrgUeAfqKyFrgBuAXdR2kquXAdcAHwGLgZXXtC1eLyNXebmcBC0RkLnAfcJ531+P32PoVrZG2rnCvmd2CcvqLD+vC3grlxek/BuX8JkIkpMGRv4Yb5sP4+6CsGF693DVof34X7Go942aa5lfbYH0AqOoK4FgRSQKiVLUo0JN71UbvVlv3sM/7+4H7Az22WW1dCSnt3YNOQdA9J5nRfXJ4bvpqfjG6B3ExgeRq02rFxLsRYwdfBMs+ct1j//cXyL8DBpwNI692Y0EZ04QCGYvpHyKSrqq7VLVIRDJEpNYeTBFhyw9BqV7ydcWR3dhUVMqrs+zBOROgqGjoMw4ufgOumQ6DL4QFr8HDR8KTJ8GiqTbUuGkygfzZeoKqbq/6oKrbcD2OItvWFUFPEEf2zGZw53QemLbcHpwz9de2rxtm/MZFcNxfYfsaeHkC3HMITPunjflkGi2QBBEtIvFVH0SkDRBfy/4tX2kR7CoMeoIQEW44tjdrt++xuwjTcG0y4IhfuYEBz30ecvvBp3e4RPHCubDkfai02QxN/dXZBgE8B/xPRJ7E9WS6HHg6qFGF2r4G6uAmCIBRvbIZlOfuIs4a2snaIkzDRcfAQSe7Zdsq+O5ZmP0sLH0fUjvCkIth8ARI6xjqSE0LEchQG/8C/g4cBPQH/uqti1zNmCDcXUQv1m7fY4P4maaT0RWO+SP8eiGc8yzk9IX82+Geg+HZM2D+q7B3T6ijNGEukDsIVPU94L0gxxI+tq1yr0Hq4lrdUb1zGJSXzv2fLOfMIXYXYZpQdKybwKjfKe6/69nPwdzJMOUKiE+F/qfBwAug86E2Nao5gM0o58+21dAmE+JTmuVyIsKvj3NtEc98vapZrmlaoYyucPQt7intS96CvifD/Cnw5Di4b7DrMlv1x5Ex2Ixy/m1fDRldmvWSR/XOYXSfHO79eBmbi0ub9dqmlYmKgm6j4PSH4KalcNpDkNYJ8v/hHsB7/Fj45iHYGciACSaS2Yxy/mxbDenNmyAAbjmpH3v2VvDvD5Y0+7VNKxWfDIMugEvfdk9rH3sb7C2B92+Guw+Cp06GmU/C7q2hjtSEQCBtEDNF5CXgDWDfn7aq+lqwggqpykrY8aPrCdLMerZN5tLDu/LElyu5cGQXDumU1uwxmFYsvbMb1uPIX8OmJe4BvAWvwts3wLs3QY+jod+pxJZFdgWC+UkgCcJ3RrkqCkRmgihaDxVlIbmDAPjlMb14ffZa/vzWQl65+jDEGg5NKOT0gTG/h9E3w4Z5boDABa/Dsg85nChYdwQcNB76nuSqp0xECmQspsuaI5CwsX21e23mNogqaW1i+c3YPtz82nxemVnAOcPz6j7ImGARgfYD3XLsn2H9HNa8/wBdds2H937rlg6DvWQxHnKadgZGE1qB9GLqJCKvi0ihiGwUkSkiErl/MmzzEkR615CFcM6wPEZ2y+Svby9i7Xbrq27ChAh0GMzK7hfBtdPhuplwzK2AuIEDHxgO9w+Hj26FVV/amFARIJBG6idxczF0ADoCb3nrItP21YBAeuj+co+KEu48ayAVqtw8ZV7VrHrGhJfsXvCzG2HiNPj1IjjhTkhpB1/fD0+dCP/qDq9cCnNehOJNoY7WNEAgbRA5quqbEJ4SkRuCFE/obVvthvmOCe1wU52zEvnDiQdxyxsLeH76Gi46NDRVXsYEJK0jjJzolpIdbva7ZR+6ockXvg4IdBwCvY6HXsdB+8Guu60Ja4EkiM0ichHwovf5fNysb5EpBM9A1OTCkZ35YOEG/vHuYn7WK5suWcGZm8KYJpWQBv1OdUtlpWvkXvahW/Jvh/x/QmIWdDsKuo92S5j8P2f2F0gKvxw4B9iAm2r0LG9dnURknIgsEZHlInKzn+0Xisg8b/lKRAb6bFslIvNFZI6IzAysOE0gRM9A+CMi3HHmAGKjo7j6ue/YU2YjcpoWJioKOgyCo34LV34Mv/kBTn8Ueh4Hq7+Ct34F9w6AewfBWzfAwjfsmYswEkgvpjXAKfU9sYhEAw8AxwEFwAwRmaqqi3x2WwkcparbROQE4FFgpM/2Maq6ub7XbrCKcihaF9L2h+o6pLfh3vMGcdlTM/jD6/O5+5yB1vXVtFxJWTDwXLeowualrjpqRb4bQHDWk4DXc6rbKOh6JOSNhDbpoY27laozQYjI08D1VZMGiUgGcJeq1nUXMQJY7k1ZiohMBk4F9iUIVf3KZ/9vgND2jipaD1oJaeGTIABG92nLr4/tzd0fLWVQXjqXHN411CEZ03gi7nmLnD4w8ufuD7R13/2UML55CL66DxBodzB0OQK6HO5ek7JDHHzrIHX1kBGR2ao6uK51fo47Cxinqld6nycAI1X1uhr2vwno67P/SmAb7qG8R1T10RqOmwhMBMjNzR06efLkWstTk+LiYjqWr2HwnN8zd8BtbMustXjNrlKV+74rZf7mCn4zPIG+mdE17ltcXExycnIzRhcckVIOsLI0RFRFKak7l5K2YyHp2xeSuvN7oivLANiV2Ikdaf3Znt6Pnal9KUnIrfdotPadOGPGjJmlqsP8bQukkTpKRDK8qUYRkcwAj/P3bfnNRiIyBrgCONJn9RGquk5E2gIficj3qvrZASd0ieNRgGHDhuno0aMDCO1A+fn5DM7Mhjkw8GcnheUDP0MP3cuZD33FA3NLeOnnw+nXwf+QB/n5+TT03yGcREo5wMrScGN/elteBuvnwOovSVr9FUlrvqbD+g/ctqQc6DQcOg2DTiPcw3vxtf9g2ndSt0B+6O8CvhKRV3E/8OfgJhCqSwHgW1fTCVhXfScRGQA8jpv7el/vKFVd570WisjruCqrAxJEk9rxo3sN0xm30trE8szlIzjzoa+4eNK3vPaLw+mclRjqsIxpHjFxkDfCLUf+2k2junEhFMyAgpnudcm7bl+Jgrb9IW+4lziGQ1ZPm/OingJppH7G60V0NO6u4IxqDc01mQH0EpFuwFrgPOAC3x1EpDNuTKcJqrrUZ30SEKWqRd7744G/BFimhttR4OaBiAvf7qQd0tvwzOUjOPuRr5kwaTqvXH0YbVMSQh2WMc0vKhraD3DL8Cvcut1bYe0s+PFblzDmvwozJ7ltCWk/DRvSfhBtdpe5brj2PEaNAp1RbhE+jcsBHlMuItcBHwDRwCRVXSgiV3vbHwb+BGQBD3o9c8q9urBc4HVvXQzwgqq+X5/rN8iOgrDqwVSTXrkpTLp0OBc+Np1zH/mG564cScf0NqEOy5jQS8x0D+L1Os59rqyEzUtcwlg321VRTX8EKspcd8m5v4N2A1xXXC9xkNXTkoYnoATRUKr6LvButXUP+7y/ErjSz3ErgIHV1wfdjoJmmYe6KQzpnMGzV4zgsqdmcPZDX/HclSPpnhMZDW7GNJmoKGh7kFuGXuLWVeyFwsV8n/8SfVNLYN0cmPE4lJe47XHJkNsf2vZzr7kHQ24/dwfSygQ1QbQ4Owpc3+sWYljXTF686lAunvQt5zzyNU9dNoKDO7a+/4iNqZfoWGg/gA3tt9K3qmG3otzdaayb4+4yNi6Eha95z2V40vIOTBxZPSE6cn9GI7dk9RRdvgtKd7a4se0P7pjGyz8/jIufmM5ZD3/FnWcNpHlm0jYmgkTHeD/6/WHwhW6dKuxc55JF4UL3unEhLP8YKr2RaqPjILOH6/WY3Qeye7v3Wb0gruV3ILEE4Uko8UabbGEJAtxMdG9edyS/eG4Wv3xxNid1i+Vno5ToKOuxYUyDibgejWkdobfPfGnlZe4J8MJFsHEBbF4GGxbA4rfcg7buYNeemV2VOHq5BwKz+7inyVsISxCenxJE+DdS+5OTEs8LVx3KbW8t5IXpa9j6+HTuOmcgHazx2pimFRPnnuxudzCu17+nvBS2/OCSx+albtrWzUvc3BjlPvO6JKS7tk7fJauHe03MCquuuJYgPPGlLfcOokpcTBT/OP0QEnZtYPLS7Yy95zP+dtrBnDooPJ/rMCaixMS7xuzcfvuvr5rnvippbF3hloIZrp1j310HEJ8Kmd18kkcPyOjq5gtP7eC69jZnkZr1amEsvnQLRMVAUttQh9JoozrFcskJh3Pjy3O5fvIc3pm3nj+N70enjJZfJ2pMixMV5YYzz+jyU/fbKuVlsH3NT0lj6w/udf1cWDQV1GcE56gY9wdseheXMDK6eO+7EFcanBkYLEF44ku3uImCIqT/c5esJF6aeCiPfb6Se/+3lGPv/pTrxvTkqlHdiY9p3r9CjDE1iImD7J5uqa5ir0se21e7aQiq3m9fA0s/gF2F+3YdHpMCY89s+vCa/IwtVHzpFncLF0FioqP4xegenDKoA397exH//nApL377I9cf04szhnQkJjoykqExESk61rVNZPXwv71st6u62r6GpbOn0z8IIdgvhGffHUQE6pjehocuGspzV4wkOzmO306Zx3H/+YwpswooK6+s+wTGmPATl+h6RvU6jk1tj6x7/wawBAGg6t1BRHZj7pG9snnj2iN47OJhxMdE8X+vzOXIOz7hgWnL2bqrLNThGWPCjFUxAZTuJLqyBFIj8w7Cl4hwXL9cjj2oLZ8t28wTX6zkzg+WcO/HyziuXy5nDevEqF459gyFMcYSBOCeloSIa4OojYhwVO8cjuqdw9KNRUz+9kden13AO/PXk50cx3H92jG2fy6H98gmLsZuNI1pjSxBwE8JIqX1JAhfvXNT+NP4ftx8Ql8++X4jb81bz9Q5a3nx2zWkxMcwpm9bRvXO4bAeWTZqrDGtiCUIaJV3EP7ExUQx7uD2jDu4PSV7K/jqh818sGAjHy/eyNS57t8oL7MNh3bLYmT3LAZ2SqN7TrJVRxkToSxBABStd68R2oupIRJiozm6by5H982lslJZWljE1z9s4esftvDhoo28MqsAgMS4aPp3SOXgjmkc1C6V7jlJ9MhJJiMpLsQlMMY0liUIgJ1rKYtNIy7GftT8iYoS+rZLpW+7VC47ohuVlcryTcXML9jB/LVuefHbNZTs/anLbEZiLD1ykumRk0znrEQ6pCfQIa0NHdLbkJuaYO0axrQAQU0QIjIOuBc3o9zjqnp7te3ibT8R2A1cqqrfBXJsk9q5ntL4LCw9BCYqSuidm0Lv3BTOHOrGrqqoVAq27WbFpl38sKnYW3bxv+83srl4/y60IpCTHE9uagKZSXFkJcWRmRRHhs/7rOQ41hZXsn7HHpLjY0iKiyHKqrKMaVZBSxAiEg08ABwHFAAzRGRqtfmsTwB6ectI4CFgZIDHNp2d6yiNz7J5FBohOkrokpVEl6wkxvTdfzyr3WXlrNtewvode1i/vYR1O/awbvseCotK2bqrjOWFxWzdVcaevRUHnviLT/a9TY6PcckiPprkhFhS4mNoExdNfEwUCbEHvibERhEf89NrTLQQExVFTJT89D5aiI0Woqutd+uE2OgooqOEaBGiRJAoiBIhStwr7P9ZxPUQMyYSBPMOYgSw3Js+FBGZDJzK/nNbnwo8o6oKfCMi6SLSHugawLFNp2gdpenDg3JqA4lxMfRsm0zPtrVPibqnrIItu1zS2LKrjOnfzaNz997sKi2nqLSc4pJyikv3UlxaTnFpBUUle9lcXEpZeSUleyso9Xktr9RmKt2BRECoShwuaahWEvvJ+4j3uXpS8Tl6v/NUX7v/ugP39T1VTYlq374NPNee3btJnJm/304tMSXu3r2bxFn5+z635MQetXcPVZPjNaVgJoiOwI8+nwtwdwl17dMxwGMBEJGJwESA3Nxc8vPz6xelVtI35RAK47uxrL7Hhqni4uL6/zuEGQEOTikhec8KtyLWW2q8zYv2lljAVXntrYSySthb4d5XVEKFKhVa9d4tlaqUV0Kl99l3v0qF8ko3uVhVC0ulgqKouvWKt2jVNp/13mtpWQUxseKzXlHctavsl9L0wLc1pTz1s0Hrs93PxWq7VnliJTExJTXu01KktakkJtrNQ93SyxIbVxGU/+eDmSD8pePq30NN+wRyrFup+ijwKMCwYcN0dEPS6Jij+T4/nwYdG4byI6QskVIOsLKEo0gpBwSvLMFMEAWA7/RsnYB1Ae4TF8CxxhhjgiiYfQ1nAL1EpJuIxAHnAVOr7TMVuFicQ4Edqro+wGONMcYEUdDuIFS1XESuAz7AVQ5PUtWFInK1t/1h4F1cF9fluG6ul9V2bLBiNcYYc6CgPgehqu/ikoDvuod93itwbaDHGmOMaT72OKsxxhi/LEEYY4zxyxKEMcYYvyxBGGOM8UvU32OWLZSIbAJWN/DwbGBzE4YTSpFSlkgpB1hZwlGklAMaV5Yuqprjb0NEJYjGEJGZqjos1HE0hUgpS6SUA6ws4ShSygHBK4tVMRljjPHLEoQxxhi/LEH85NFQB9CEIqUskVIOsLKEo0gpBwSpLNYGYYwxxi+7gzDGGOOXJQhjjDF+tfoEISLjRGSJiCwXkZtDHU9jiMgqEZkvInNEZGao46kPEZkkIoUissBnXaaIfCQiy7zXjFDGGKgaynKbiKz1vps5InJiKGMMhIjkicg0EVksIgtF5HpvfYv7XmopS0v8XhJE5FsRmeuV5c/e+ib/Xlp1G4SIRANLgeNwkxfNAM5X1eDMfR1kIrIKGKaqLe7hHxEZBRTj5ig/2Fv3L2Crqt7uJe8MVf1dKOMMRA1luQ0oVtV/hzK2+vDmh2+vqt+JSAowCzgNuJQW9r3UUpZzaHnfiwBJqlosIrHAF8D1wBk08ffS2u8gRgDLVXWFqpYBk4FTQxxTq6SqnwFbq60+FXjae/807n/osFdDWVocVV2vqt9574uAxbj54lvc91JLWVocdYq9j1WztStB+F5ae4LoCPzo87mAFvofjUeBD0VklohMDHUwTSDXm2EQ77VtiONprOtEZJ5XBRX21TK+RKQrMBiYTgv/XqqVBVrg9yIi0SIyBygEPlLVoHwvrT1BiJ91LbnO7QhVHQKcAFzrVXWY8PAQ0AMYBKwH7gppNPUgIsnAFOAGVd0Z6ngaw09ZWuT3oqoVqjoI6ASMEJGDg3Gd1p4gCoA8n8+dgHUhiqXRVHWd91oIvI6rQmvJNnp1x1V1yIUhjqfBVHWj9z91JfAYLeS78eq4pwDPq+pr3uoW+b34K0tL/V6qqOp2IB8YRxC+l9aeIGYAvUSkm4jEAecBU0McU4OISJLX+IaIJAHHAwtqPyrsTQUu8d5fArwZwlgapep/XM/ptIDvxmsMfQJYrKp3+2xqcd9LTWVpod9Ljoike+/bAMcC3xOE76VV92IC8Lq13QNEA5NU9e+hjahhRKQ77q4B3FzjL7SksojIi8Bo3LDFG4FbgTeAl4HOwBrgbFUN+8bfGsoyGleNocAq4OdV9cXhSkSOBD4H5gOV3uo/4OruW9T3UktZzqflfS8DcI3Q0bg/8l9W1b+ISBZN/L20+gRhjDHGv9ZexWSMMaYGliCMMcb4ZQnCGGOMX5YgjDHG+GUJwhhjjF+WIExEEpF0EbnG53MHEXk1SNeKFZFZDTx2mIjc14hrXyoiHRp6vDG1sQRhIlU6sC9BqOo6VT0rSNc6EviqIQeq6kxV/VUjrn0pYAnCBIUlCBOpbgd6eGP83ykiXavmZ/D+6n5DRN4SkZUicp2I3Cgis0XkGxHJ9PbrISLve4Mffi4ifWu41jjgveorRaRYRO7wjv9YREaISL6IrBCRU7x9RovI297727wB46r2+ZW3fl/s3uebvH3PAoYBz3vlbCMiQ0XkU++aH/gMvfArEVnkDUo3ucn+lU1EswRhItXNwA+qOkhVf+Nn+8HABbixd/4O7FbVwcDXwMXePo8Cv1TVocBNwIM1XGsMbjyc6pKAfO/4IuBvuLlHTgf+UsO5+gJjvbhu9cYP8ktVXwVmAhd6A7eVA/8FzvKuOckrG7h/j8GqOgC4uqZzGuMrJtQBGBMi07x5AYpEZAfwlrd+PjDAG/XzcOAVN4wPAPHVT+LV/29V1d1+rlEGvO9z3lJV3Ssi84GuNcT1jqqWAqUiUgjk1qNMfXCJ7yMv5mjcCKUA83B3Gm/ghjAxpk6WIExrVerzvtLncyXu/4soYLv3l3ltTgA+qGHbXv1pLJt911DVShGp6f8937gqvFjK2f9uP6GGYwVYqKqH+dl2EjAKOAX4o4j0V9XyGs5jDGBVTCZyFQEpDT3YmytgpYicDW40UBEZ6GdXv+0PTWwj0FZEskQkHjjZZ5tvOZcAOSJyGOzrXdVfRKKAPFWdBvwW14CfHOSYTQSwBGEikqpuAb4UkQUicmcDT3MhcIWIzAUWUm06WnFzmvdS1e8bF23tVHUvrs1iOvA2bmjnKk8BD3uzi0UDZwF3eDHPwVWTRQPPeVVbs4H/ePMIGFMrG83VmAbyhpC+SFWt0ddEJEsQxhhj/LIqJmOMMX5ZgjDGGOOXJQhjjDF+WYIwxhjjlyUIY4wxflmCMMYY49f/B9gCP4aVaJW7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "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$'])\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cLCET4hnFg2S"
   },
   "source": [
    "## Optimization with `scipy.minimize_scalar`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-rrDxccbFg2S"
   },
   "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": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "N3NI_5GWFg2T"
   },
   "outputs": [],
   "source": [
    "def CB(tf):\n",
    "    soln = odeint(batch, [CAf, 0], [0, tf])\n",
    "    return soln[-1][1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ui89VA8hFg2V"
   },
   "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": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 102
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 6072,
     "status": "ok",
     "timestamp": 1559058898925,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 300
    },
    "id": "diuqirNEFg2W",
    "outputId": "280f6670-fe70-4b88-a7bc-29fd74530584"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -1.3374806339222158\n",
       " message: '\\nOptimization terminated successfully;\\nThe returned value satisfies the termination criteria\\n(using xtol = 1.48e-08 )'\n",
       "    nfev: 23\n",
       "     nit: 19\n",
       " success: True\n",
       "       x: 4.023594924340666"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minimize_scalar(lambda t: -CB(t), bracket=[0, 50])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 6064,
     "status": "ok",
     "timestamp": 1559058898925,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 300
    },
    "id": "ijeny_k_Fg2Z",
    "outputId": "652ad20d-a98e-4223-d3a5-bd2d0ce24521"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Concentration c_B has maximum 1.3374806339222158 moles/liter at time 4.023594924340666 minutes.\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.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gvhPt9DxFg2c"
   },
   "source": [
    "## Solution using Pyomo\n",
    "\n",
    "The variable to be found is the time $t_f$ corresponding to the maximum concentration of $B$. For this purpose we introduce a scaled time\n",
    "\n",
    "$$\\tau = \\frac{t}{t_f}$$\n",
    "\n",
    "so that $\\tau=1$ as the desired solution. The problem then reads\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\max_{t_f} C_B(\\tau=1)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "subject to\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{dC_A}{d\\tau} & = - t_f k_A C_A \\\\\n",
    "\\frac{dC_B}{d\\tau} & = t_f(k_A C_A - k_B C_B)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The solution to this problem is implemented as a solution to the following Pyomo model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 7492,
     "status": "ok",
     "timestamp": 1559058900360,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 300
    },
    "id": "7S2xp-pHFg2d",
    "outputId": "1a9f8c2f-fdba-43dd-b457-d41afd4813b4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Concentration c_B has maximum 1.3374805810221082 moles/liter at time 4.023594178375689 minutes.\n"
     ]
    }
   ],
   "source": [
    "V   = 40    # liters\n",
    "kA  = 0.5   # 1/min\n",
    "kB  = 0.1   # l/min\n",
    "cAf = 2.0   # moles/liter\n",
    "\n",
    "m = ConcreteModel()\n",
    "\n",
    "m.tau = ContinuousSet(bounds=(0, 1))\n",
    "\n",
    "m.tf = Var(domain=NonNegativeReals)\n",
    "m.cA = Var(m.tau, domain=NonNegativeReals)\n",
    "m.cB = Var(m.tau, domain=NonNegativeReals)\n",
    "\n",
    "m.dcA = DerivativeVar(m.cA)\n",
    "m.dcB = DerivativeVar(m.cB)\n",
    "\n",
    "m.odeA = Constraint(m.tau, \n",
    "    rule=lambda m, tau: m.dcA[tau] == m.tf*(-kA*m.cA[tau]) if tau > 0 else Constraint.Skip)\n",
    "m.odeB = Constraint(m.tau,\n",
    "    rule=lambda m, tau: m.dcB[tau] == m.tf*(kA*m.cA[tau] - kB*m.cB[tau]) if tau > 0 else Constraint.Skip)\n",
    "\n",
    "m.ic = ConstraintList()\n",
    "m.ic.add(m.cA[0]  == cAf)\n",
    "m.ic.add(m.cB[0]  == 0)\n",
    "\n",
    "m.obj = Objective(expr=m.cB[1], sense=maximize)\n",
    "\n",
    "TransformationFactory('dae.collocation').apply_to(m)\n",
    "SolverFactory('ipopt').solve(m)\n",
    "print('Concentration c_B has maximum', m.cB[1](), 'moles/liter at time', m.tf(), 'minutes.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "NPwJqR_IFg2h"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "06.02-Maximizing-Concentration-of-an-Intermediate-in-a-Batch-Reactor.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}