notebooks/06.01-Unconstrained-Scalar-Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "8CY-WljpSRJl"
   },
   "source": [
    "# Unconstrained Scalar Optimization\n",
    "\n",
    "Introductory calculus courses introduce the minimization (or maximization) of a function of a single variable. 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. \n",
    "\n",
    "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",
    "$$\n",
    "\\begin{align}\n",
    "f'(x^*) = 0\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Let's see how this is put to work in the context of process engineering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\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 *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PBRBB1uWSRJp"
   },
   "source": [
    "## Application: Maximizing production of a reaction intermediate\n",
    "\n",
    "A desired product $B$ is produced as intermediate in a series reaction\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 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 a solution with concentration $C_{A,f} = 2.0\\ \\mbox{moles/liter}$.\n",
    "\n",
    "A 100 liter tank is available to run the reaction. Below we will answer the following questions:\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?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "S4Nhug1xSRJq"
   },
   "source": [
    "## Mathematical model for a 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",
    "$$\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",
    "\n",
    "At steady-state the material balances become\n",
    "\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",
    "\n",
    "which can be solved for $C_A$\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\bar{C}_A & = \\frac{qC_{A,f}}{q + Vk_A} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "and then for $C_B$\n",
    "\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",
    "\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": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 295
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 566,
     "status": "ok",
     "timestamp": 1559058282028,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 300
    },
    "id": "8Zy9i0xySRJr",
    "outputId": "658bce6c-ebb5-4e61-9bc6-0b6d08e1b893"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEWCAYAAACEz/viAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4qElEQVR4nO3deZxVdf348dd7dmD2hWEZZgEGFBVEBkQFHcwFzbL6WmlpmZmZmtY3+2Z9+1W2muX3W98yzUpTU0nTTMvdHMEFRZBFUGAYBhjW2RhmYZjt/fvjnMHLOMuZYc7c7f18PO6Dc8923597h/u+n+V8jqgqxhhjoltMsAMwxhgTfJYMjDHGWDIwxhhjycAYYwyWDIwxxmDJwBhjDJYMzCCISKGIqIjEBTuWSCEi60Wk1Ifz5orIUhFpFJHbhvv8JvJYMohgInK5iKwTkRYR2SMid4hI+iCOrxSRs47itV8ZyrGh6Gjei4Bz/FlEfhy4TlWPU9Wyowqud1cBNUCqqn7Dh/P3SxzXi8g7ItIsIlUi8oiInOBuzxORR0WkRkQa3L/Ty0VkoYg0uY9m98dHU8AjX0TKRKTVfV4jIo+JyPiRLmOksWQQoUTkG8DPgW8CacB8oAB4XkQSghlbJArB2lIBsEGHcFXpMJXl18ANwPVAJjANeBz4sLv9fmCHG2cW8Dlgr6ouU9VkVU0GjnP3Te9ep6rb3XXXuftMBZKBXw5DzNFNVe0RYQ8gFWgCPtVjfTKwD7jCff5n4McB20uBKnf5fqALOOie67+AQkCBOHefNOBPwG5gJ/BjIBY4FmgFOt1j9/cRZyZwD7ALqAceD9j2JaAcqAOeACYEbFPgamCze9ztgPQ49l2gEdgAnOSunwA8ClQDW4HrA475AfAwcJ973HqgxMN78UVgO7DU3fcRYA/QACwFjnPXXwW0A23uOZ5011cCZ7nLicCv3Pdjl7ucGPjZAN9wP8PdwBf6eF//3OO1zvJ47m+5sd/fyzmnAP8GanFqHA/gfEn39vrF7mc/r5+/0SbgxAH+jrvf47ge68uAKwOeXwOsD/b/u3B/BD0Ae/jwocJioKPnfyJ3273AQ+7yn+kjGbjPD39Ruc+P+M+J80vv98AYYCzwJvBld9vlwCsDxPkv4K9ABhAPnOGuP9P9wjnJ/RL7De6XrbtdgX8C6UA+zpf7YnfbJ3ES01xAcH45FuDUglcC3wMSgMlABXCue9wPcBLY+TgJ7WfAcg/vxX1u+Ue5668AUnj/y3d1wDFHvN89zwv8EFjuvpc5wGvAjwI+mw53n3g3zhYgo4/3tudn6+XcP3fjHtXL+aYCZ7vbc3AS3a/6eO2rgW0DfPYvAK8CFwP5fezT/R73mQxwahUvAP8I9v+7cH8EPQB7+PChwqXAnj623QI87y73/MIoxWMyAHKBQ4FfHMAlwEvu8uX0kwyA8Ti/tj/wZYZT27g14Hkyzi/dQve5AgsCtj8M3OQuPwvc0Ms5Twa291j3beAed/kHwAsB22YABz28F5P7KWO6u09ab+93z/MCW4DzA7adC1QGfDYHA78YcWoI8/t47Z6f7UDnbgOSBvE39jHg7T62/TcBibSPfTLcv8X1OLWI1cDcHvsc/nvrsb4MJxE2uNtX00dCsYf3h/UZRKYaILuPtt/x7vajVYDzC3W3iOwXkf04tYSxHo+fBNSpan0v2yYA27qfqGoTTvPExIB99gQst+AkjO7zbukj3gndsbrxfgcnqfV1ziQP7ec7uhdEJFZEbhGRLSJyAOeLHiB7gHN0O6Lc7vKEgOe1qtrRI8ZkvBno3NWq2trXwSIyVkSWiMhOt2x/oe9y1eL8nfVJVetV9SZVPQ7nM1gNPC4iMnBRAKeJLw2YiZNY8jweZ/pgySAyvY7zq/0TgStFZAxwHvCiu6oZGB2wy7ge5+mv83GH+xrZqpruPlLd/9wDHdt9fGYfo5t24Xx5B8adhdP8M5AdOO3bva3fGhBruqqmqOr5Hs4JfZcncP1ngAtx2ujTcH7ZgtNc1d85uh1RbpwmsF0e4xvIQOceKLafufvMVNVUnNpnX1/cLwJ5IlLiJTBVrcHpAJ6A04/kmaquw+mrun0QicT0wpJBBFLVBuBm4DcislhE4kWkEKdzswqnQxScX2Pni0imiIwDvtbjVHtx2tZ7e43dwHPAbSKSKiIxIjJFRM4IODavr5FL7vFPA78TkQw3xtPdzQ8CXxCRE0UkEfgp8IaqVnoo/h+BG0Vkjju8caqIFOD0ZxwQkW+JyCj3V/zxIjLXwzm7y9PrexEgBSdB1uIk2Z8O8hwPAd8VkRwRycbp3/iLx/gGcrTnTsEdDCAiE3FGqfVKVTcDvwMeEpFSEUkQkSQRuVhEbgIQkZ+773+ciKQAXwHKVbV2CGW7F6dG+tEhHGtclgwilKreitMM8kvgAPAGzq/jD6nqIXe3+4E1OM0Zz+F05gb6Gc4XyH4RubGXl/kcTmfsBpxRPX/j/eaBf+O0B+8Rkb6apS7D6Qt4D6f9+2tu7C8C/w9n5M9unF/6F3ss9yPAT3ASSiNOJ3emqnYCHwFOxBlJVIOTONK8nJeB3wtwOpO34dRgNuB02Ab6EzDDPcfjvRz/Y+AtYC2wDljlrhsOR3vum3E69BtwOv4fG2D/64Hf4oz02o/TdPdx4El3+2jg7+62Cpxay5C+zFW1Dfg/nL8ZM0TidsgYY4yJYlYzMMYYY8nAGGOMz8lARO4WkX0i8k4f20VE/k9EykVkrYic5Gc8xhhjeud3zeDPOFfD9uU8nEvXi3Eu17/D53iMMcb0wtfJtVR1qTuksS8XAvep04u9XETSRWS8O+ywT9nZ2VpY2N9p+9bc3MyYMWOGdGyosbKEpkgpS6SUA6ws3VauXFmjqjm9bQv2TIsTCbiCE2cM/ESc4YRHEJGrcGoP5Obm8stfDm2SwqamJpKTvV60GdqsLKEpUsoSKeUAK0u3RYsWbetrW7CTQW9XDPY61lVV7wLuAigpKdHS0tIhvWBZWRlDPTbUWFlCU6SUJVLKAVYWL4I9mqgKZy6ZbnkM3+X3xhhjPAp2MngC+Jw7qmg+0DBQf4Exxpjh52szkYg8hDM9braIVAHfx5npElW9E3gKZ172cpwZGL/gZzzGGGN65/dooksG2K7AtX7GYIwxZmDBbiYyxhgTAiwZGGOMCfrQUhOGurqUA63t1DQdoqapjdqmNpoOtXOoo4tD7V20dXZxqKOLWBHiYoW4GCE2RkiMjyU1KY6UpDhSkuJJSYoja0wiWWN6veWBMWYEWTIwfTrY1sn6XQ1s3tdERXUTW2uaqahuZkd9C+2dwzf1eVyMkBIP+e+8Qk5KErmpiUzMGEVB5hgKskaTnzWa1KT4YXs9Y8wHWTIwhzUcbOf1LTW8Wl7Lqu31vLenkc4u50s/MS6GouwxTB+XwjnHjWNsSiJZyQlkJzv/JifGkRQfS2JcDIlxscTHCqrQ3tVFZ5fS3qkcau/kQGsHTYc6aGxt58DBDmqaDrGvsZU1m7YRMzqBqvoWVm6ro76l/YjYMsckkJ85mik5yUwfl8y03BSmj0thXGoSdrdDY46eJYMoV1XfwtPr9vD0O7tZvWM/XQpjEmKZnZ/BV86YwqxJ6RwzLoWJ6aOIiRncl64IJMbEvr9iVDxjU3vftyxxD6Wl8w4/bzrUwfbaFrbVNrOtzvm3sqaFZZureXRV1eH9UpPimD4uhWm5KRwzLoXjJqYxY3wqSfGxvb2MMaYPlgyiUPOhDp5cs4u/vrWDt7fvB+C4Calct2gqC4pzmJ2fTnxscMcWJCfGMWNCKjMmfDB71De3sXFvI5v2NrJxj/N4Ys0uHnijA3CanablpjAzL40T8tKYOTGd6eNSSIiz8RLG9MWSQRTZVtvMH5dt5e9v76TpUAfTcpP5r8XT+fAJ4ynICp8ZHTPGJDB/chbzJ2cdXqeq7GpoZV1VA+t27mdtVQPPrN/DkhXOPIgJsTEcOyGVOfkZlBRmUFKQwdjUpGAVwZiQY8kgCmza28jtL5Xz5JpdxMXGcMHM8XxmXj5zCjIipr1dRJiYPoqJ6aNYfPw4wEkQVfUHWVO1n3VVDby9Yz8PvLGNu1/dCkBexihKCjKYU5hJSUEG03JTiB1kU5gxkcKSQQTbe6CVW5/ZyGNvVzEqPpYvLZzMFxcWMTYlOn4RiwiTMkczKXM0F8ycAEBbRxcbdh/grco6Vm6r59UttTy+2pkbMSUxjtkFGZxclMn8yVnMzEsLenOZMSPFkkEEam3v5K6lFdxRtoXOLuWq0ydz9elTyLDx/CTExXDipHROnJTOlQvfrz28tc1JDiu21vOLZzcCMDohlpLCTOZPzuSUyVmcMDGNOEsOJkJZMogwb2+v58ZH1rClupnzTxjHTYuPJT9rdLDDClmBtYePz84DoLbpEG9urWN5RS2vV9Ry6zNOchjjJodTpjj9FcdPSLXkYCKGJYMI0dreyf++sIk/LK1gXGoS910xj9On9Xp3OzOArOREzjthPOedMB6AGjc5vL6lluUVtdzy9HuA06w0f0oWC4uzOW1qNpOzx0RMH4yJPpYMIsCOuhaueWAV63Y2cPHcSXznw8faFbvDKDs5kfNPGM/5bnKobjzEG1trebW8llfKq3l+w14AJqQlsaA4m8z2Dk5oOkRWcmIwwzZmUCwZhLmyjfv42l9X09mp/P6yOZx73LhghxTxclISuWDmhMOd0ttqm3mlvIZXNtfwzDt7ONDawZ1rXuC4CaksmJrNguJs5hZm2oVwJqRZMghTqspTW9t45NkVTM9N4c5L51CYHT7XCkSSgqwxFGSN4bMnF9DZpdz7xL9pSc1n2eYa7n51K79fWkFiXAxzCzNZUJzNwuJsZoxPtSYlE1IGTAbi/MXmqeqOEYjHeNDVpfz4X+/y8MZ2PnzCeH75yVmMSrBfnaEgNkaYnB5LaWkx151ZTPOhDt6srOOVzU7N4Zan3+OWp53axenFOZwxPYeFU7NtpJcJugGTgaqqiDwOzPE/HDOQQx2dfOPhNfxz7W7OLojjN5fMHvScQWbkjEmMY9H0sSyaPhaAfQdaWbq5hpc3VfPie3t5dFUVMQKzJqVzxrQczpiWw8y8dLv4zYw4r81Ey0Vkrqqu8DUa06+2ji6u+csqXnxvHzeddwzTu7ZbIggzY1OTuGhOHhfNyaOzS1lbtZ+XN1Xz8qZqfv3iZn71wmbSR8ezsNhJDKdPy46aiwRNcHlNBouAq0WkEmgGBKfSMNOvwMyROjq7uGHJ27z43j5+9LHjuWx+AWVl1nIXzmJjhNn5GczOz+BrZ02jvrmNV8prKNvoJIcn1zhXRs8Yn8oZ053kMKcgw66KNr7wmgzO8zUK06/OLuXGR9bw9Dt7+O6Hj+Wy+QXBDsn4IGNMAh+ZNYGPzJpAV5fy7p4DTq1hYzV/cK8oT06M47SpWZwxbSynT8smL8MuKDTDw1MyUNVtIrIAKFbVe0QkB0j2NzTT7YdPrufx1bv45rnTuXLh5GCHY0ZATIxw3IQ0jpuQxjWlU2lsbee1LbWHk8Oz651rG6aOTT7c1zCvyIavmqHzlAxE5PtACTAduAeIB/4CnOZfaAbgvtcruff1bXxpYRHXLpoa7HBMkKQkxXPuceM497hxqCpbqpsP9zXcv3wbf3plK0nxMcyfnEXptBzOmD6WwqzRNnzVeOa1mejjwGxgFYCq7hKRFN+iMgAs21zNzU9u4EPHjOWm844NdjgmRIgIU8cmM3VsMl9cUMTBtk6Wb63l5Y3VLN1UzQ+e3ABPbiA/czRnTMuhdHoOp0zJYnSCXVZk+ub1r6PNHWKqACJiVzf5bEt1E9c8sIqpOcn8+pLZNtTQ9GlUQuwRw1e317bw8qZ9vLzJuUXo/cu3kRAbw9yiDDc5jKV4bLLVGswRvCaDh0Xk90C6iHwJuAL4o39hRbeDbZ185S8rSYiN4Y+fLyE50X7RGe/ys0Zz2SmFXHZKIYc6Onmrsv5wX8NPn3qPnz71HuPTkg7XGk6dmm1zWRnPHci/FJGzgQM4/QbfU9XnfY0siv3oXxvYtLeJ+66Yx6RMGy1ihi4xLpbTpjqzqn7n/GPZ3XCQl92hq/9au5slK3YQGyPMyc84PHx1xvhUu34lCnntQP65qn4LeL6XdWYYPb1uNw++sZ0vnz7ZpqA2w2582igunpfPxfPyae/sYvWO/ZRtdJqUfvHsRn7x7EaykxOdi906OpjV3GZTZUQJr+0PZwM9v/jP62WdOQpV9S1869G1zMpL4xvnTA92OCbCxcc6k+fNLczkm+ceQ3XjIZZtrqZsYzUvvbeP+pZ2fr/ueWbluVNlTM9hlk2VEbH6TQYi8hXgGmCyiKwN2JQCvOpnYNFGVfnmI2vpUvi/S2aTEGdXmZqRlZOSyCdOyuMTJ+Udnn21MTmfsk37+M2/N/PrF22qjEg2UM3gQeBp4GfATQHrG1W1zreootAjK6t4vaKWn3z8eAqybLCWCa7A2VdvOKv48FQZ3dc2dE+VcdyE1MMXvZ1kU2WEtYGSgapqpYhc23ODiGRaQhge1Y2H+Mm/3mVeYSaXzM0PdjjGfEBfU2WUbazmrqUV/K5sCymJcZw2NftwR/SE9FHBDtsMgpeawQXASkBxJqjrpsCAcyOIyGLg10As8EdVvaXH9jScq5nz3Xh+qar3eC1AJPjRPzdwsK2Tn37ieBvFYUJez6kyDrS281p591QZ+3hm/R4AiscmUzo9hzOmjWVuUQaJcTZVRijrNxmo6gXuv0VDObmIxAK343RAVwErROQJVd0QsNu1wAZV/Yg759FGEXlAVduG8prh5qWN+3hizS6+dlYxU8faRd0m/KQmxbP4+HEsPt6ZKqN8X9Ph5qR7X9vGH5ZtZVR8LKdMyaJ0eg6nF+dQYFNlhJyBOpBP6m+7qq4a4PzzgHJVrXDPtwS4EAhMBgqkuHdUSwbqgI4BzhsR2jq6uPmJ9UzJGcNXSqcEOxxjjpqIUJybQnFuClcunExLWwfLK5ypMso2VfPv9/YBMClzFAum5rCwOJtTp2SRPtqGrwabqGrfG0Ve6udYVdUz+z25yEXAYlW90n1+GXCyql4XsE8K8ARwDM4opU+r6r96OddVwFUAubm5c5YsWdLfS/epqamJ5OTQmHD1+cp2Hnivjf+ck8jMnMFfZRxKZTlaVpbQ40c59jR3sb62k3dqOnmvrpODHU7bc1FaDMdlxXJ8dixT0mOIG+bm0kj5TODoyrJo0aKVqlrS27aBmokWDekV39fbJ9oz+5wLrAbOBKYAz4vIMlU90COWu4C7AEpKSrS0tHRIAZWVlTHUY4dTw8F2vrb0JRZMzearF80bUpU5VMoyHKwsocfvcrR3drFmx36Wba7hlfIanqrcz5MV7YxOiGX+5CwWFmezsDibKTlHP49SpHwm4F9Z/J70pgqYFPA8D9jVY58vALeoU0UpF5GtOLWEN32OLah+91I5DQfb+fb5x1jbqYlK8bExlBRmUlKYydfPnsaB1nZe31LLK5trWLb5/SalcalJLCzOZkFxNgumZpOVnBjkyCOT38lgBVAsIkXATuBi4DM99tkOfAhYJiK5OHMfVfgcV1DtqGvhnlcr+cTsPI6bkBbscIwJCakB92wA5//JK+VOYnhuw14eWVkFONc2LCjOZuHUHEoKM+yGPsNkwGTgduzmqeqgb7irqh0ich3wLM7Q0rtVdb2IXO1uvxP4EfBnEVmH06z0LVWtGexrhZPbntuICNx47rRgh2JMyJqUOZpL5uVzybx8OruUd3Y2sGxzNcs213D3K1v5/csVJMbFMK8o0+2IzrZJ9o7CgMnAvY/B48CcobyAqj4FPNVj3Z0By7uAc4Zy7nBUUd3EP9bs4qqFkxmfZhflGONFbIwwa1I6syalc92ZxTQf6uCNrbVOf8PmGn761HsApI+O55TJWZw61RmlNDl7jDXDeuS1mWi5iMxV1RW+RhMF7ijbQkJsjN3L2JijMCYxjjOPyeXMY3IB2NPQymtbanhtSy2vldfw9DvOhW/jUpM4dUoWmR3tTNt/0K6K7ofXZLAIuFpEKoFmnOYcVdWZfgUWiarqW/j72zu5dH4BOSnWCWbMcBmXlnR4kj1VZVttC6+6yaFsUzV1zW38cd2/KcoewylTsjhtSjbzJ2daZ3QAr8ngPF+jiBJ3La1ABK463WoFxvhFRCjMHkNh9hg+e3IBXV3KA/96ibaMIl4rr+GJ1bt48I3tABw7PpVTp2Rx6pQs5hVlkhLFd3zzeqezbSKyAChW1XvcaSMi4wqOEbLvQCtLVuzgP07Ks6qqMSMoJkaYlBJD6YIivrigiI7OLtbubOD1LbW8Wl7D/cu38adXthIbI8zMS+O0KU5/w0kF0TVSyeudzr4PlOAM+7wHiMeZXO40/0KLLH98ZSsdnV1cfYZNO2FMMMXFxnBSfgYn5Wdw7aKptLZ3smpbPa9tqeXVLTXc8fIWfvtSOQmxMZw4KZ2TJ2cyf3IWJ+VnMCohcpOD12aijwOzgVXgjAByp5EwHjQd6uDBN7ZzwcwJFGbbvQqMCSVJ8bHO6KOp2dzIdBpb21lRWccbFXUsr6jl9pfK+c2/y4mPFWblpTN/chYnT85kTkEGoxP8vlRr5HgtSZs7xFQBRMS+0Qbh76uqaDrUwRdOKwx2KMaYAaQkxR8xUqmxtZ23ttWzvKKWNyrqDtcc4txmJSc5ZFFSkMGYxPBNDl4jf1hEfg+ki8iXgCuAP/gXVuRQVe57fRsnTEzjxEnpwQ7HGDNIKUnxLJo+lkXTxwJOTX/l4eRQe/jmPrExwgkTneQwf7IzzUZyGCUHrx3IvxSRs4EDwDTge6r6vK+RRYjXK2rZvK+JX1w00y5+MSYCJCfGHb7VJ0BLW2ByqONPr1Rw58tOcjh+YhrzizLdZqVM0kaF7milwaStdcAonFlH1/kTTuS5//VtpI+O5yOzJgQ7FGOMD0YnxLGwOIeFxe8nh1Xb9vPG1lqWV9Ry96tb+b07rHx6bgrzipxaw7zCTMalJQU5+vd5HU10JfA94N84F5z9RkR+qKp3+xlcuNvdcJDnNuzlyoVFUTVEzZhoNjohzplhtTgbgINtnby9o563KutZUVnH31ZWcd/r2wDnJj9zCzKZW5TJ3MJMpuQEb/oMrzWDbwKzVbUWQESygNcASwb9ePCN7XSpcunJBcEOxRgTJKMSYjl1ijORHkBHZxfv7m7kzco6VmytY+nmah57eycAmWMSKCnIYJ6bHGZMSCU+NmZE4vSaDKqAxoDnjcCgZzGNJh2dXTz05g7OnD6WSZmjgx2OMSZExMXGcEJeGifkpfHFBUWoKltrmllRWccKt/bw3Ia9AIxOiGV2fjpzC53kMDs/3b+4PO63E3hDRP6B02dwIfCmiPwngKr+j0/xha1lm2uoaTrEp+ZOGnhnY0zUEhEm5yQzOSeZT8/NB2DvgdbDzUpvbq3j1y9uRhXiYoSLp8dT6kMcXpPBFvfR7R/uv3bhWR8ee3sn6aPjDw9HM8YYr3JTk/jwzPF8eOZ4AA60trNqm5McMlp2+vKaXoeW3uzLq0eoA63tPLd+D58qmURC3Mi09xljIldqUjyl08dSOn0sZWV7fHkN+6bywTPr9nCoo4tPnDQx2KEYY4wnlgx88OiqKiZnj7Erjo0xYcOSwTCrqm/hja11fHz2RLvi2BgTNjwlAxG5VURSRSReRF4UkRoRudTv4MLR4+544Y/NtiYiY0z48FozOEdVDwAX4FxzMA3nQjQTQFV57O2dnFyUadcWGGPCitdk0D270vnAQ6pa51M8YW3T3iYqqpv56Ik2D5ExJrx4vc7gSRF5DzgIXOPe9rLVv7DC0zPv7EEEzp6RG+xQjDFmUDzVDFT1JuAUoERV24EWnKuQTYBn1+9hTn4GY1NCZyZCY4zxwmsH8mjgWuAOd9UEnHsiG9eOuhY27D7AuceNC3YoxhgzaF77DO4B2oBT3edVwI99iShMPbveuSrQkoExJhx5TQZTVPVWoB1AVQ/i3NfAuJ5dv4djxqWQn2WjiIwx4cdrMmgTke67nCEiU4BDvkUVZqobD/HWtnoWH2+1AmNMePI6muj7wDPAJBF5ADgNuNyvoMLNC+/uRdWaiIwx4cvrrKXPi8gqYD5O89ANqlrja2Rh5Nn1e8jPHM0x42xGb2NMeOo3GYjIST1W7Xb/zReRfFVd5U9Y4aP5UAevltdw+amFNheRMSZsDVQzuK2fbQqcOYyxhKXlFbW0d6rdxMYYE9b6TQaqumikAglXSzdVMyo+ljmFGcEOxRhjhszrRWfxInK9iPzNfVwnIvEDHwkislhENopIuYjc1Mc+pSKyWkTWi8jLgylAsC3bXMP8yZkkxsUGOxRjjBkyr6OJ7sCZrO537vPL3HVX9neQiMQCtwNn41yotkJEnlDVDQH7pLvnXayq20UkbNpbdtS1UFHTzKXzC4IdijHGHBWvyWCuqs4KeP5vEVnj4bh5QLmqVgCIyBKcOY02BOzzGeAxVd0OoKr7PMYUdMs2OwOqTp+WE+RIjDHm6HhNBp0iMkVVtwCIyGSg08NxE4EdAc+rgJN77DMNiBeRMiAF+LWq3tfzRCJyFXAVQG5uLmVlZR5DP1JTU9OQj+3p0bdbyUwSdqxfQdWGkR9JNJxlCTYrS+iJlHKAlcULr8ngm8BLIlKBc51BAfAFD8f19g2pvcQwB/gQMAp4XUSWq+qmIw5SvQu4C6CkpERLS0s9hn6ksrIyhnpsoI7OLr5a9jwfPiGPRYtmHvX5hmK4yhIKrCyhJ1LKAVYWL7xedPaiiBQD03G+4N9TVS/TUVQBkwKe5wG7etmnRlWbgWYRWQrMAjYRwtZU7aextYOFxdZEZIwJf15HE8UC5wKlOL/grxWR//Rw6AqgWESKRCQBuBh4osc+/wAWikicO1X2ycC7HuMPmqWbaogROG1qVrBDMcaYo+b5Tmc4dzZbB3R5PbmqdojIdcCzQCxwt6quF5Gr3e13quq7IvIMsNY99x9V9Z3BFCIYlm6uZmZeOumjE4IdijHGHDWvySBPVYfUMK6qTwFP9Vh3Z4/nvwB+MZTzB0PDwXbW7NjPdYumBjsUY4wZFl6nsH5aRM7xNZIwsmpbPV0K86dYE5ExJjJ4rRksB/4uIjE4N7gRQFU11bfIQtiblXXExwqzJ9kUFMaYyOA1GdwGnAKsU9WeQ0OjzoqtdRw/MY1RCTYFhTEmMnhtJtoMvGOJAFrbO1lb1cDcwsxgh2KMMcPGa81gN1AmIk8TcLtLVf0fX6IKYWt27Kets8uSgTEmonhNBlvdR4L7iForKusAKCmw/gJjTOTwegXyzX4HEi7erKxnWm4yGWOiOicaYyKM1z4DA3R2Kau21VsTkTEm4lgyGIR3dx+g6VAH84osGRhjIku/yUBELhERu7LK9eZWp7/AagbGmEgzUJ9BAfCIe4vLF4GngTejdYjpiso6JqaPYkL6qGCHYowxw6rfmoGq3qKqZwLnA2uAK4BVIvKgiHxORHJHIshQoKqsqKyzJiJjTETyOpqoEfi7+0BEZgDnAffhTG0d8bbVtlDT1EZJoQ0pNcZEHq/XGRzBvaH9BpxpKqLCmqr9AJw4KT2ocRhjjB9sNJFH66oaSIyLYVpuSrBDMcaYYWfJwKO1VQ3MmJBKfKy9ZcaYyDOobzYRGS0iJSISVTf+7exS3tnVwKy89GCHYowxvhjoOoOPikiliKwSkfOB9cBvgXUi8vkRiTAEbKluoqWtkxMmpgU7FGOM8cVAHcg/As4B0oCXgJmqWiEiY3GuO7jX5/hCwtqqBgBmTbJkYIyJTAMlgy5V3QQgIltVtQJAVfeJSIfv0YWItVX7GZMQS1F2crBDMcYYXwyUDGJEJAOnOanLXZbubb5GFkLWVjVw/MQ0YmNk4J2NMSYMDZQM0oCVvJ8AVgVsi4opKdo6utiw+wCfP6Ug2KEYY4xv+k0Gqlro5SQicpyqrh+WiELMpr2NtHV0MdNGEhljIthwNfXcP0znCTndnccz86zz2BgTuYYrGURsY/q6nftJGxVPfuboYIdijDG+Ga5kELH9B2urGpiZl4ZIxOY7Y4yJnhFBQ9Ha3snGPY12sZkxJuINVzJoG6bzhJT39jTS0aXWX2CMiXgDTUdxrohc1Mv6z4rI2d3PVXW+H8EF27u7DwAwY7wlA2NMZBuoZnAz8HIv618Efjj84YSWd3cfIDkxjrwMu82lMSayDZQMRqtqdc+VqroHGONPSKHj3d0HOGZcCjF25bExJsINlAySROQDF6aJSDwQ0T+XVZX3djdy7PjUYIdijDG+GygZPAb8QUQO1wLc5TvdbQMSkcUislFEykXkpn72mysinb31UQRDVf1BGg91cMx4u7OZMSbyDZQMvgvsBbaJyEoRWQlUAtXutn6JSCxwO3AeMAO4RERm9LHfz4FnBxW9j7o7j61mYIyJBgPNTdQB3CQiNwNT3dXlqnrQ4/nnuftXAIjIEuBCYEOP/b4KPArM9Rq4397d3YgIHDPOagbGmMgnqv5dPOw2+SxW1Svd55cBJ6vqdQH7TAQeBM4E/gT8U1X/1su5rgKuAsjNzZ2zZMmSIcXU1NREcvLA9yX4zdutVDV28fPTQ3caCq9lCQdWltATKeUAK0u3RYsWrVTVkt62DTSF9dHqbRhOz+zzK+BbqtrZ35QPqnoXcBdASUmJlpaWDimgsrIyvBz7w7fKmF2UQmnpnCG9zkjwWpZwYGUJPZFSDrCyeOF3MqgCJgU8zwN29dinBFjiJoJs4HwR6VDVx32OrU+t7Z1U1jZzwczxwQrBGGNGlOdk4DbnFAQeo6pLBzhsBVAsIkXATuBi4DOBO6hqUcBr/Bmnmehxr3H5oaK6mS6F4lzrLzDGRAdPyUBEfg58Gqfjt9NdrUC/yUBVO0TkOpxRQrHA3aq6XkSudrffOdTA/bRpbyMA063z2BgTJbzWDD4GTFfVQ4N9AVV9Cniqx7pek4CqXj7Y8/th095G4mKEwqyIv8jaGGMA77OWVgDxfgYSSjbtbaIoewwJcTbDtzEmOnitGbQAq0XkReBw7UBVr/clqiDbtLeRE2zaamNMFPGaDJ5wHxHvYFsnO+pb+I+T8oIdijHGjBhPyUBV7xWRBGCau2qjqrb7F1bwlO9rQhWm5UbGBSrGGOOF19FEpcC9OPMSCTBJRD7vYWhp2CmvdkYSFVsyMMZEEa/NRLcB56jqRgARmQY8BITu5blDtGVfM7ExQn6mjSQyxkQPr8Nl4rsTAYCqbiJCRxdtqW6iIHO0jSQyxkQVrzWDt0TkT8D97vPPAiv9CSm4tlQ3MTnHmoiMMdHF68/frwDrgeuBG3CuRL7ar6CCpaOzi8qaFqaMtSYiY0x08Tqa6BDwP+4jYlXVH6Sts4spVjMwxkSZfpOBiDysqp8SkXV8cOppVHWmb5EFwZbqJgBLBsaYqDNQzeAG998L/A4kFLyfDKyZyBgTXfrtM1DV3e7iNaq6LfABXON/eCNry75mspMTSB+dEOxQjDFmRHntQD67l3XnDWcgoWBLdROTs62JyBgTffpNBiLyFbe/YLqIrA14bAXWjkyII6eytpmibGsiMsZEn4H6DB4EngZ+BtwUsL5RVet8iyoIGlvbqWlqo9CSgTEmCvWbDFS1AWgALgEQkbFAEpAsIsmqut3/EEdGZU0LAEXZo4MciTHGjDxPfQYi8hER2QxsBV7GmbDuaR/jGnFba5sBrGZgjIlKXjuQfwzMBza5N7D/EPCqb1EFQWWNkwwKbII6Y0wU8poM2lW1FogRkRhVfQk40b+wRl5lTTPj05IYlRAb7FCMMWbEeZ2obr+IJANLgQdEZB/Q4V9YI29rbTOFWVYrMMZEJ681gwtx7oP8deAZYAvwEb+CCobKmmbrLzDGRK0BawYiEgv8Q1XPArpw7ngWURpa2qlvabeRRMaYqDVgzUBVO4EWEUkbgXiCorJ7JJE1ExljopTXPoNWYJ2IPA80d69U1et9iWqEVdqwUmNMlPOaDP7lPgJ9YErrcLWjzrngbFKGNRMZY6KT12SQrqq/DlwhIjf0tXO42V7XwtiURBtWaoyJWl5HE32+l3WXD2McQbW9roX8TKsVGGOi10B3OrsE+AxQJCJPBGxKAWr9DGwk7ag7yMlFmcEOwxhjgmagZqLXgN1ANnBbwPpGImQK67aOLnY1HGSS1QyMMVFsoFlLtwHbgFNGJpyRt3P/QVSxZiJjTFTzOmvpJ0Rks4g0iMgBEWkUkQN+BzcStrsjifKzLBkYY6KX1w7kW4GPqmqaqqaqaoqqpno5UEQWi8hGESkXkZt62f7ZgDuovSYiswZTgKN1OBlYzcAYE8W8JoO9qvruYE/uTmVxO879kmcAl4jIjB67bQXOUNWZwI+Auwb7OkdjR10LCXEx5CQnjuTLGmNMSPF6ncFbIvJX4HHgUPdKVX1sgOPmAeWqWgEgIktwJr3bEHCO1wL2Xw7keYxpWGyvbWFSxihiYmQkX9YYY0KK12SQijNr6TkB6xQYKBlMBHYEPK8CTu5n/y8ywndQ217XYiOJjDFRT1T9m1VCRD4JnKuqV7rPLwPmqepXe9l3EfA7YIF7I52e268CrgLIzc2ds2TJkiHF1NTURHJy8uHn17zQzPwJcXxuRvg1E/UsSzizsoSeSCkHWFm6LVq0aKWqlvS6UVUHfADTgBeBd9znM4HvejjuFODZgOffBr7dy34zce6RMM1LPHPmzNGheumllw4vNxxs04Jv/VPvKCsf8vmCKbAs4c7KEnoipRyqVpZuwFvax/eq1w7kP7hf5O1uAlkLXOzhuBVAsYgUiUiCe0zglcyISD5Oc9NlqrrJYzzDYmf9QQDyMkaN5MsaY0zI8dpnMFpV3xQ5opN1wNteqmqHiFwHPAvEAner6noRudrdfifwPSAL+J17/g7tqxozzLqTwcR0SwbGmOjmNRnUiMgU3GmrReQinGkqBqSqTwFP9Vh3Z8DylcCVHuMYVjv3u8nAagbGmCjnNRlcizP+/xgR2YlzbcClvkU1QnbuP0hCXAzZY8Kv89gYY4aTp2SgznUCZ4nIGCBGVRv9DWtk7Kw/SF66XWNgjDFe5yb6qYikq2qzqjaKSIaI/Njv4PxWVd9iTUTGGIP36SjOU9X93U9UtR4435eIRtDO/Qet89gYY/CeDGJF5HDDuoiMAsK6ob21vZOapjZLBsYYg/cO5L8AL4rIPTgjiq4A7vUtqhFgI4mMMeZ9XjuQbxWRdcCHAAF+pKrP+hqZz6oOX3Bm8xIZY4zXmgGq+jQjPImcn3a5NYMJ6UlBjsQYY4Ivau90tnv/QUQgN9WSgTHGeK0Z3Ap8RIdwg5tQtbuhlbEpicTHeu1DN8aYyOXrnc5C2e6GVsanWeexMcaA/3c6C1m7Gg4yPTcl2GEYY0xI8PtOZyFJVdnT0ErptLHBDsUYY0KC16GlX/A7kJF04GAHLW2djE+zzmNjjAHvo4nyROTvIrJPRPaKyKMiMqI3rh9Ouw84w0rH27BSY4wBvHcg34Nzh7IJODe5f9JdF5Z2728FsA5kY4xxeU0GOap6j6p2uI8/Azk+xuWrXQ1uzcCaiYwxBvCeDGpE5FIRiXUflwK1fgbmpz0NrcQIjE0J67n2jDFm2HhNBlcAnwL24Nzu8iJ3XVjatb+V3NQk4uyCM2OMAbyPJtoOfNTnWEbMngMHGWdNRMYYc5jX0UT3ikh6wPMMEbnbt6h8tnt/KxOs89gYYw7z2k4ys5c7nc32JSKfqSp7DrTaBHXGGBPAazKIEZGM7icikskgpr8OJa2d0NLWybg06zw2xphuXr/QbwNeE5G/4UxD8SngJ75F5aP6VgVs6mpjjAnktQP5PhF5CzgT505nn1DVDb5G5pP9hywZGGNMT4O509kGICwTQKD61i7AkoExxgSKuoH29W7NwC44M8aY90VdMtjfqqQkxjEmMSz7v40xxhfRlwwOKbl2wZkxxhwhOpNBqjURGWNMoKhLBvWtSm6K1QyMMSZQVCWDri61ZiJjjOlFVCWD+pY2OhVybSSRMcYcwfdkICKLRWSjiJSLyE29bBcR+T93+1oROcmvWPYccO5wZtcYGGPMkXxNBiISC9wOnAfMAC4RkRk9djsPKHYfVwF3+BXPvgOHABhrycAYY47gd81gHlCuqhWq2gYsAS7ssc+FwH3qWA6ki8h4P4LZ69YM7F4GxhhzJL+vvJoI7Ah4XgWc7GGfiTh3VDtMRK7CqTmQm5tLWVnZoIPZW93BCZnKu6uWszlGBn18qGlqahrS+xCKrCyhJ1LKAVYWL/xOBr194+oQ9kFV7wLuAigpKdHS0tJBB1MKzCorYyjHhqIyK0tIipSyREo5wMrihd/NRFXApIDnecCuIexjjDHGR34ngxVAsYgUiUgCcDHwRI99ngA+544qmg80qOrunicyxhjjH1+biVS1Q0SuA54FYoG7VXW9iFztbr8TeAo4HygHWoAv+BmTMcaYD/J96k5VfQrnCz9w3Z0Bywpc63ccxhhj+hZVVyAbY4zpnSUDY4wxlgyMMcZYMjDGGAOI038bXkSkGtg2xMOzgZphDCeYrCyhKVLKEinlACtLtwJVzeltQ1gmg6MhIm+pakmw4xgOVpbQFClliZRygJXFC2smMsYYY8nAGGNMdCaDu4IdwDCysoSmSClLpJQDrCwDiro+A2OMMR8UjTUDY4wxPVgyMMYYE13JQEQWi8hGESkXkZuCHc/REJFKEVknIqtF5K1gxzMYInK3iOwTkXcC1mWKyPMistn9NyOYMXrRRzl+ICI73c9ltYicH8wYvRKRSSLykoi8KyLrReQGd31YfS79lCPsPhcRSRKRN0VkjVuWm931vnwmUdNnICKxwCbgbJwb6qwALlHVDUENbIhEpBIoUdWwu5BGRE4HmnDufX28u+5WoE5Vb3ETdYaqfiuYcQ6kj3L8AGhS1V8GM7bBcu87Pl5VV4lICrAS+BhwOWH0ufRTjk8RZp+LiAgwRlWbRCQeeAW4AfgEPnwm0VQzmAeUq2qFqrYBS4ALgxxTVFLVpUBdj9UXAve6y/fi/AcOaX2UIyyp6m5VXeUuNwLv4tyLPKw+l37KEXbU0eQ+jXcfik+fSTQlg4nAjoDnVYTpH4lLgedEZKWIXBXsYIZBbvcd7tx/xwY5nqNxnYisdZuRQrpZpTciUgjMBt4gjD+XHuWAMPxcRCRWRFYD+4DnVdW3zySakoH0si6c28hOU9WTgPOAa90mCxN8dwBTgBOB3cBtQY1mkEQkGXgU+JqqHgh2PEPVSznC8nNR1U5VPRHn3vDzROR4v14rmpJBFTAp4HkesCtIsRw1Vd3l/rsP+DtOM1g42+u293a3++4LcjxDoqp73f/AXcAfCKPPxW2XfhR4QFUfc1eH3efSWznC+XMBUNX9QBmwGJ8+k2hKBiuAYhEpEpEE4GLgiSDHNCQiMsbtHENExgDnAO/0f1TIewL4vLv8eeAfQYxlyLr/k7o+Tph8Lm5n5Z+Ad1X1fwI2hdXn0lc5wvFzEZEcEUl3l0cBZwHv4dNnEjWjiQDc4WS/AmKBu1X1J8GNaGhEZDJObQCc+1g/GE5lEZGHgFKcqXj3At8HHgceBvKB7cAnVTWkO2f7KEcpTlOEApXAl7vbd0OZiCwAlgHrgC539Xdw2tvD5nPppxyXEGafi4jMxOkgjsX54f6wqv5QRLLw4TOJqmRgjDGmd9HUTGSMMaYPlgyMMcZYMjDGGGPJwBhjDJYMjDHGYMnADJGIXO/ODPmAiFwuIr8dodctFJHP+Hj+8SLyXC/rm9x/J4jI39zlE8Nh9kuvRORqEfncURz/neGMx4wsSwZmqK4BzlfVzw73iUUkrp/NhYBvyQDnCs9n+9qoqrtU9SL36YnAoJLBAGUbFkN9DVW9U1XvO4qXtmQQxiwZmEETkTuBycATIvL1HtsKRORFd0KwF0Uk351sq0Ic6SLS1T2XkogsE5Gp7nzzd7m/yu9zawDLRGSV+zjVfYlbgIXunPRfd8/9CxFZ4b7ml/uI+b/FuZfFCyLykIjc2EfxFgNP91P2QhF5x72K/YfAp91YPu1eGX63G8vbInKhe8zlIvKIiDyJM7ngeBFZ6h73jogs7OV1KkXk5+LMZ/+miEx11+eIyKPua6wQkdPc9Ue8fz3OVSoiL4vIwyKySURuEZHPuuddJyJTAs5xo7tcFvD6m7pj7FkLFJF/uue/BRjllukBd9ul7vGrReT34kwjb0KVqtrDHoN+4FzFme0uXw781l1+Evi8u3wF8Li7/AxwHHABztQg/w0kAlvd7T/AmXt+lPt8NJDkLhcDb7nLpcA/A+K4Cviuu5wIvAUU9Yh1Ds4VqaOBVKAcuLGXMsUCq/sob5P7byHwTs9yu89/ClzqLqfj3D9jjLtfFZDpbvsG8N8Br5nSx/vbvc/nussMPAgscJfzcaZd+MD71+NcpcB+YLz7Hu0Ebna33QD8KuAcN7rLZcBt7vL5wAt9lPmfQGnge+QuH+v+LcS7z38HfC7Yf7f26Pvhe5XVRJ1TcG6+AXA/cKu7vAw4HSgCfgZ8CXgZJzF0e0JVD7rL8cBvReREoBOY1sfrnQPMFJHupps0nOSxNWCfhcDfVbUFQET6mpPqZN6f7ngozgE+GlDrSML5wgZn+uHuKQNWAHeLM6Ha46q6uo/zPRTw7/+6y2cBM0QOT8KbKu48VRz5/vW0Qt3pF0RkC9DdL7IOWNTHMd2T1a3ESYKD8SGcJLzCjXUUYTDJXTSzZGD81j3fyTLgamAC8D3gmzi/WJcG7NscsPx1nPl+ZuE0Z7b2cX4Bvqqqfbbz94ijP+fh1GCGSoD/UNWNR6wUOZmAsqnqUreZ7MPA/SLyC+29rV57WY4BTun5pe9+4Qa+fz0dCljuCnjeRd/fA937dAbs08GRzctJfRwrwL2q+u1+YjIhxPoMzHB7DWdGWIDP4tyqD5xf3KcCXaraCqwGvoyTJHqTBuxWZ8rhy3CaUwAagZSA/Z4FvuL+ykZEpokzk2ugpcDHRWSU+yv6I3285oeAFwcs4ft6i+Wr4n4zi8js3g4SkQJgn6r+AWeGzZP6OP+nA/593V1+Drgu4FwnDiLe4VAJnCgiMSIyiSOngm7v/hxw3seLRGQsHL5vb8HIhmoGw2oGZrhdj9ME8k2gGvgCgKoeEpEdwHJ3v2U4M0mu6+M8vwMeFZFPAi/x/q/etUCHiKwB/gz8GqcJY5X7JVxNj9sAqnM/3L/iJKBt9JKARCQHaNXB3dDlJeAmce5E9TPgRziz4q51Y6nE6SPpqRT4poi049xDua/hnIki8gbOj7ZL3HXXA7eLyFqc/79LcWpcI+VVnCa4dTjTQK8K2HYXTtlXqepnReS7OB3mMUA7cC3O+29CkM1aaqKO9HLTehG5FMhT1VuCFlgAEakESlS1JtixmOhgNQNjAFX9S7BjMCaYrGZgjDHGOpCNMcZYMjDGGIMlA2OMMVgyMMYYgyUDY4wxwP8HLprwF39j7yQAAAAASUVORK5CYII=\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 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(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "AimarcucSRJy"
   },
   "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": {
    "colab_type": "text",
    "id": "xMzxrJOWSRJz"
   },
   "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",
    "$$\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",
    "\n",
    "Clearing out the non-negative common factors yields\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "1 - \\frac{q^*}{(q^* + V k_A)} - \\frac{q^*}{(q^* + Vk_B)} = 0\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "and multiplying by the non-negative denominators produces\n",
    "\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",
    "\n",
    "Expanding these expressions followed by arithmetic cancellations gives the final result\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "q^*  =  V\\sqrt{k_Ak_B}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "which shows the optimal dilution rate, $\\frac{q^*}{V}$, is equal the geometric mean of the rate constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 68
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 556,
     "status": "ok",
     "timestamp": 1559058282031,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 300
    },
    "id": "tvH0fKHXSRJ0",
    "outputId": "5f12a428-c7dc-4800-9143-ec2d993ecabc"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flowrate at maximum CB =  8.94427190999916 liters per minute.\n",
      "Maximum CB = 0.9549150281252629 moles per liter.\n",
      "Productivity =  8.541019662496845 moles per minute.\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",
    "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": {
    "colab_type": "text",
    "id": "C5krauY4SRJ4"
   },
   "source": [
    "## Numerical solution with Pyomo\n",
    "\n",
    "This problem can also be solved using Pyomo to create a model instance. First we make sure that Pyomo and ipopt are installed, then we proceed with the model specification and solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 68
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1812,
     "status": "ok",
     "timestamp": 1559058292199,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 300
    },
    "id": "BYbuwdy2SRJ5",
    "outputId": "624f95d9-9f6b-4759-f82c-bee6f23b0509"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flowrate at maximum CB =  8.944271964904416 liters per minute.\n",
      "Maximum CB = 0.954915028125263 moles per liter.\n",
      "Productivity =  8.541019714926701 moles per minute.\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",
    "# create a model instance\n",
    "m = ConcreteModel()\n",
    "\n",
    "# create the decision variable\n",
    "m.q = Var(domain=NonNegativeReals)\n",
    "\n",
    "# create the objective\n",
    "m.CBmax = Objective(expr=m.q*V*kA*CAf/(m.q + V*kB)/(m.q + V*kA), sense=maximize)\n",
    "\n",
    "# solve using the nonlinear solver ipopt\n",
    "SolverFactory('ipopt').solve(m)\n",
    "\n",
    "# print solution\n",
    "print('Flowrate at maximum CB = ', m.q(), 'liters per minute.')\n",
    "print('Maximum CB =', m.CBmax(), 'moles per liter.')\n",
    "print('Productivity = ', m.q()*m.CBmax(), 'moles per minute.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Fy4-evYVSRJ8"
   },
   "source": [
    "One advantage of using Pyomo for solving problems like these is that you can reduce the amount of algebra needed to prepare the problem for numerical solution. This not only minimizes your work, but also reduces possible sources of error in your solution.\n",
    "\n",
    "In this example, the steady-state equations are\n",
    "\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",
    "\n",
    "with unknowns $C_B$ and $C_A$. The modeling strategy is to introduce variables for the flowrate $q$ and these unknowns, and introduce the steady state equations as constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 68
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 273,
     "status": "ok",
     "timestamp": 1559058293753,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 300
    },
    "id": "Zie6HLy6SRJ9",
    "outputId": "aa26b431-d6c4-41e1-b069-c4ac4111e06a"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flowrate at maximum CB =  8.94427200287657 liters per minute.\n",
      "Maximum CB = 0.9549150281377385 moles per liter.\n",
      "Productivity =  8.541019751298467 moles per minute.\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",
    "# create a model instance\n",
    "m = ConcreteModel()\n",
    "\n",
    "# create the decision variable\n",
    "m.q  = Var(domain=NonNegativeReals)\n",
    "m.CA = Var(domain=NonNegativeReals)\n",
    "m.CB = Var(domain=NonNegativeReals)\n",
    "\n",
    "# equations as constraints\n",
    "m.eqn = ConstraintList()\n",
    "m.eqn.add(0 == m.q*(CAf - m.CA) - V*kA*m.CA)\n",
    "m.eqn.add(0 == -m.q*m.CB + V*kA*m.CA - V*kB*m.CB)\n",
    "\n",
    "# create the objective\n",
    "m.CBmax = Objective(expr=m.CB, sense=maximize)\n",
    "\n",
    "# solve using the nonlinear solver ipopt\n",
    "SolverFactory('ipopt').solve(m)\n",
    "\n",
    "# print solution\n",
    "print('Flowrate at maximum CB = ', m.q(), 'liters per minute.')\n",
    "print('Maximum CB =', m.CBmax(), 'moles per liter.')\n",
    "print('Productivity = ', m.q()*m.CBmax(), 'moles per minute.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ae2lkOJsSRKB"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "06.01-Unconstrained-Scalar-Optimization.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
}