reactor_design.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6fvcYnfEJUez"
   },
   "source": [
    "# Chemical Reactor Design\n",
    "**Prepared by:** Zhicheng Lu (zlu3@nd.edu, 2024) and Shammah Lilonfe (slilonfe@nd.edu, 2024) at the University of Notre Dame."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Review syntax for Pyomo and Pyomo.DAE\n",
    "* Demonstrate the power of Pyomo in solving and optimizing complex dynamic process systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Statement: Design of Chemical Reactors for Chemical Transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chemical reactors are often the most important unit operations in a chemical plant. Reactors come in many forms, however two of the most common idealizations are the continuously stirred tank reactor (CSTR) and the plug flow reactor. The CSTR is often used in modeling studies, and it can be effectively modeled as a lumped parameter system. In this example, we will consider the following reaction scheme known as the Van de Vusse reaction:\n",
    "\n",
    "$$\n",
    "\\text{A} \\xrightarrow{k_1} \\text{B} \\xrightarrow{k_2} \\text{C}\n",
    "$$\n",
    "\n",
    "$$\n",
    "2\\text{A} \\xrightarrow{k_3} \\text{D}\n",
    "$$\n",
    "\n",
    "A diagram of the system is shown in figure below, where F is the volumetric flowrate. The reactor is assumed to be filled to a constant volume, and the mixture is assumed to have constant density, so the volumetric flowrate into the reactor is equal to the volumetric flowrate out of the reactor. Since the reactor is assumed to be wellmixed, the concentrations in the reactor are equivalent to the concentrations of each component flowing out of the reactor, given by $C_A$, $C_B$, $C_C$, and $C_D$.\n",
    "\n",
    "![process scheme](https://github.com/zhichenglu123/optimization/raw/Project-1---Reactor-Design-Adapted-from-Pyomo---Zhicheng-Lu-%26-Shammah-Lilonfe/notebooks/contrib-dev/research.png)\n",
    "\n",
    "This example (both text and figure) was taken from Hart et al (2017)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install Packages and Load Modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing the libraries needed for the project\n",
    "import sys\n",
    "if \"google.colab\" in sys.modules:\n",
    "    !wget \"https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py\"\n",
    "    import helper\n",
    "    helper.easy_install()\n",
    "else:\n",
    "    sys.path.insert(0, '../')\n",
    "    import helper\n",
    "helper.set_plotting_style()\n",
    "\n",
    "# Importing packages to be used in this project\n",
    "import pyomo.environ as pyo\n",
    "import pyomo.dae as dae\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Original Problem from the Reference Material"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal is to produce product B from a feed containing reactant A. If we design a reactor that is too small, we will obtain insufficient conversion of A to the desired product B. However, given the reaction scheme, if the reactor is too large (e.g., too much reaction is allowed to occur), a significant amount of the desired product B will be further reacted to form the undesired product C. As a result, the goal in this exercise will be to solve for the optimal reactor volume producing the maximum outlet concentration for product B.\n",
    "\n",
    "The steady-state mole balances for each of the four components are given by,\n",
    "\n",
    "$$\\begin{align*}\n",
    "&0=\\frac{F}{V}C_{Af}-\\frac{F}{V}C_A-k_1C_A-2k_3C_A^2 \\\\\n",
    "& 0=-\\frac{F}{V}C_B+k_1C_A-k_2C_B\\\\\n",
    "& 0=-\\frac{F}{V}C_C+k_2C_B \\\\\n",
    "& 0=-\\frac{F}{V}C_D+k_3C_A^2\n",
    "\\end{align*}$$\n",
    "\n",
    "The known parameters for the system are,\n",
    "\n",
    "$$\\begin{align*}\n",
    "C_{Af}=10^4\\,\\mathrm{\\frac{gmol}{m^3}}\\quad{k_1=\\frac{5}{6}\\,\\mathrm{min^{-1}}}\\quad{k_2=\\frac{5}{3}\\,\\mathrm{min^{-1}}}\\quad{k_3=\\frac{1}{6000}\\,\\mathrm{\\frac{m^3}{gmol\\cdot min}}}\n",
    "\\end{align*}$$\n",
    "\n",
    "Since the volumetric flowrate F always appears as the numerator over the reactor volume V, it is common to consider this ratio as a single variable, called the space velocity (sv). Our optimization formulation will seek to find the space-velocity that maximizes the outlet concentration of the desired product B.\n",
    "\n",
    "This original problem was taken from Hart et al (2017)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimization Model Formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mathematically, we want to solve the following optimal reactor volume problem:\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\max_{S_v} \\quad & C_B \\\\ \n",
    "\\mathrm{s.t.} \\quad & 0=S_vC_{Af}-S_vC_A-k_1C_A-2k_3C_A^2 \\\\\n",
    "& 0=-S_vC_B+k_1C_A-k_2C_B \\\\\n",
    "& 0=-S_vC_C+k_2C_B \\\\\n",
    "& 0=-S_vC_D+k_3C_A^2 \\\\\n",
    "& C_A, C_B, C_C, C_D, S_v \\geq 0.\n",
    "\\end{align*}$$\n",
    "\n",
    "Here \"s.t.\" means \"subject to\", i.e., the constraints."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Degrees of Freedom Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the optimization model formulated above:\\\n",
    "Number of variables = 5 \\\n",
    "Number of equality constraints = 4 \\\n",
    "Therefore, the degrees of freedom = 5 - 4 = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pyomo Implementation of the Original Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ipopt 3.13.2: \n",
      "\n",
      "******************************************************************************\n",
      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
      "         For more information visit http://projects.coin-or.org/Ipopt\n",
      "\n",
      "This version of Ipopt was compiled from source code available at\n",
      "    https://github.com/IDAES/Ipopt as part of the Institute for the Design of\n",
      "    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE\n",
      "    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.\n",
      "\n",
      "This version of Ipopt was compiled using HSL, a collection of Fortran codes\n",
      "    for large-scale scientific computation.  All technical papers, sales and\n",
      "    publicity material resulting from use of the HSL codes within IPOPT must\n",
      "    contain the following acknowledgement:\n",
      "        HSL, a collection of Fortran codes for large-scale scientific\n",
      "        computation. See http://www.hsl.rl.ac.uk.\n",
      "******************************************************************************\n",
      "\n",
      "This is Ipopt version 3.13.2, running with linear solver ma27.\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:       11\n",
      "Number of nonzeros in inequality constraint Jacobian.:        0\n",
      "Number of nonzeros in Lagrangian Hessian.............:        5\n",
      "\n",
      "Total number of variables............................:        5\n",
      "                     variables with only lower bounds:        5\n",
      "                variables with lower and upper bounds:        0\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:        4\n",
      "Total number of inequality constraints...............:        0\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        0\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0 -2.0000000e+03 7.50e+03 6.25e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0\n",
      "   1 -1.0801475e+03 5.54e+02 2.46e+00  -1.0 1.39e+03    -  5.54e-01 1.00e+00h  1\n",
      "   2 -1.0763574e+03 8.86e+01 1.87e+02  -1.0 3.66e+02    -  9.31e-01 1.00e+00h  1\n",
      "   3 -1.0727252e+03 1.07e+01 5.26e+00  -1.0 9.35e+01    -  9.71e-01 1.00e+00h  1\n",
      "   4 -1.0726714e+03 4.01e+00 1.22e-01  -1.0 6.03e+01    -  9.68e-01 1.00e+00h  1\n",
      "   5 -1.0724371e+03 3.44e-04 2.93e-05  -2.5 4.19e-01    -  1.00e+00 1.00e+00h  1\n",
      "   6 -1.0724372e+03 1.63e-08 3.93e-09  -3.8 3.85e-03    -  1.00e+00 1.00e+00h  1\n",
      "   7 -1.0724372e+03 5.46e-11 2.70e-11  -5.7 2.24e-04    -  1.00e+00 1.00e+00h  1\n",
      "   8 -1.0724372e+03 4.55e-13 3.56e-14  -8.6 2.78e-06    -  1.00e+00 1.00e+00h  1\n",
      "\n",
      "Number of Iterations....: 8\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:  -1.0724372001086319e+03   -1.0724372001086319e+03\n",
      "Dual infeasibility......:   3.5606114692563649e-14    3.5606114692563649e-14\n",
      "Constraint violation....:   4.5474735088646414e-14    4.5474735088646412e-13\n",
      "Complementarity.........:   2.5059065225787338e-09    2.5059065225787338e-09\n",
      "Overall NLP error.......:   2.5059065225787338e-09    2.5059065225787338e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 9\n",
      "Number of objective gradient evaluations             = 9\n",
      "Number of equality constraint evaluations            = 9\n",
      "Number of inequality constraint evaluations          = 0\n",
      "Number of equality constraint Jacobian evaluations   = 9\n",
      "Number of inequality constraint Jacobian evaluations = 0\n",
      "Number of Lagrangian Hessian evaluations             = 8\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      0.008\n",
      "Total CPU secs in NLP function evaluations           =      0.000\n",
      "\n",
      "EXIT: Optimal Solution Found.\n",
      "The optimum space-velocity is: 1.34 min^-1\n",
      "The outlet concentration of species B is: 1.07 kmol/m^3\n",
      "The optimum reactor volume is: 0.74 m^3\n"
     ]
    }
   ],
   "source": [
    "# function that creates pyomo model for the reaction system\n",
    "def create_model(k1, k2, k3, caf):\n",
    "    \"\"\"Creates reactor model\n",
    "    Argument:\n",
    "        caf: feed concentration of A, gmol/m^3\n",
    "        k1: rate constant for A to B, min^-1\n",
    "        k2: rate constant for B to C, min^-1\n",
    "        k3: rate constant for 2A to D, m^3/(gmol.min)\n",
    "\n",
    "    Returns:\n",
    "        model: pyomo model\"\"\"\n",
    "    # Create a Pyomo model\n",
    "    model = pyo.ConcreteModel()\n",
    "    \n",
    "    # create the variables\n",
    "    model.sv = pyo.Var(initialize=1.0, within=pyo.PositiveReals) # space-velocity\n",
    "    model.ca = pyo.Var(initialize=5000.0, within=pyo.PositiveReals) # concentration of A\n",
    "    model.cb = pyo.Var(initialize=2000.0, within=pyo.PositiveReals) # concentration of B\n",
    "    model.cc = pyo.Var(initialize=2000.0, within=pyo.PositiveReals) # concentration of C\n",
    "    model.cd = pyo.Var(initialize=1000.0, within=pyo.PositiveReals) # concentration of D\n",
    "    \n",
    "    # create the objective\n",
    "    model.obj = pyo.Objective(expr = model.cb, sense=pyo.maximize)\n",
    "    \n",
    "    # create the constraints\n",
    "    model.ca_bal = pyo.Constraint(expr = (0 == model.sv * caf - model.sv * model.ca - k1 * model.ca - \n",
    "                                          2.0 * k3 * model.ca ** 2.0)) # species A mole balance\n",
    "    model.cb_bal = pyo.Constraint(expr=(0 == -model.sv * model.cb + k1 * model.ca - k2 * model.cb)) # species B mole balance\n",
    "    model.cc_bal = pyo.Constraint(expr=(0 == -model.sv * model.cc + k2 * model.cb)) # species C mole balance\n",
    "    model.cd_bal = pyo.Constraint(expr=(0 == -model.sv * model.cd + k3 * model.ca ** 2.0)) # species D mole balance\n",
    "    \n",
    "    return model\n",
    "\n",
    "# define data\n",
    "k1 = 5.0/6.0 # min^-1\n",
    "k2 = 5.0/3.0 # min^-1\n",
    "k3 = 1.0/6000.0 # m^3/(gmol.min)\n",
    "caf = 10000.0 # gmol/m^3\n",
    "\n",
    "# reactor model object\n",
    "m = create_model(k1, k2, k3, caf)\n",
    "\n",
    "# solve model\n",
    "solver = pyo.SolverFactory('ipopt')\n",
    "solver.solve(m, tee=True)\n",
    "\n",
    "# print the solution of sv and the outlet concentration of B\n",
    "optimal_sv = pyo.value(m.sv)\n",
    "cb_max = pyo.value(m.cb)*10**(-3) # kmol/m^3\n",
    "print(f\"The optimum space-velocity is: {round(optimal_sv, 2)} min^-1\")\n",
    "print(f\"The outlet concentration of species B is: {round(cb_max, 2)} kmol/m^3\")\n",
    "\n",
    "# calculate optimum reactor volume for a volumetric flow rate of 1 m^3/min\n",
    "F = 1\n",
    "optimum_reactor_volume = F / pyo.value(m.sv)\n",
    "print(f\"The optimum reactor volume is: {round(optimum_reactor_volume, 2)} m^3\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimal Solution Found for the Original Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal space-velocity is $\\mathrm{1.34\\,min^{-1}}$, giving an outlet concentration for B of $\\mathrm{1.07\\,kmol\\cdot m^{-3}}$. For a volumetric flow rate of $\\mathrm{1\\,m^{3}\\cdot min^{-1}}$, the optimum reactor volume is $\\mathrm{0.74\\,m^{3}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_0X0jKuPwcWn"
   },
   "source": [
    "## How this Work Goes Beyond the Original Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2x8jafDG2F1X"
   },
   "source": [
    "In this project, species C is the desired product, and we are considering the dynamic state of the system (most often, industrial operations are dynamic). The goal of this work is to determine the volume of the reactor that will maximize the concentration of species C in the effluent after 10 minutes. Furthermore, for economic reasons, the reactor volume should not exceed $\\mathrm{100\\,m^3}$.\n",
    "\n",
    "The mole balance for each of the species in this system under dynamic state is given by:\n",
    "\n",
    "$$\\begin{align*}\n",
    "&\\frac{dC_A}{dt}=\\frac{F}{V}C_{Af}-\\frac{F}{V}C_A-k_1C_A-2k_3C_A^2 \\\\\n",
    "& \\frac{dC_B}{dt}=-\\frac{F}{V}C_B+k_1C_A-k_2C_B\\\\\n",
    "& \\frac{dC_C}{dt}=-\\frac{F}{V}C_C+k_2C_B \\\\\n",
    "& \\frac{dC_D}{dt}=-\\frac{F}{V}C_D+k_3C_A^2\n",
    "\\end{align*}$$\n",
    "\n",
    "Where $C_A$, $C_B$, $C_C$ and $C_D$ are the concentrations of species A, B, C and D, respectively. $F$ is the feed flow rate, $C_{Af}$ is the concentration of species A in the feed, $k_1$, $k_2$ and $k_3$ are the rate constants of reaction 1 (species A to B), reaction 2 (species B to C) and reaction 3 (species A to D), respectively.\n",
    "\n",
    "The data for the system are:\n",
    "$$\\begin{align*}\n",
    "C_{Af}=10^4\\,\\mathrm{\\frac{mol}{m^3}},\\quad{k_1=\\frac{5}{6}\\,\\mathrm{min^{-1}}},\\quad{k_2=\\frac{5}{3}\\,\\mathrm{min^{-1}}},\\quad{k_3=\\frac{1}{6000}\\,\\mathrm{\\frac{m^3}{mol\\cdot min}}},\\quad{F=1\\,\\mathrm{\\frac{m^3}{min}}}\n",
    "\\end{align*}$$\n",
    "\n",
    "The data was obtained from Hart et al (2017)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "E3hzUN6FjrEC"
   },
   "source": [
    "### Optimization Model Formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9cGBk30DkjgR"
   },
   "source": [
    "Mathematically, we want to solve the following optimal reactor volume problem:\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\max_{V} \\quad & C_C(t=10) \\\\\n",
    "\\mathrm{s.t.} \\quad & \\frac{dC_A}{dt}=\\frac{F}{V}C_{Af}-\\frac{F}{V}C_A-k_1C_A-2k_3C_A^2 \\\\\n",
    "& \\frac{dC_B}{dt}=-\\frac{F}{V}C_B+k_1C_A-k_2C_B \\\\\n",
    "& \\frac{dC_C}{dt}=-\\frac{F}{V}C_C+k_2C_B \\\\\n",
    "& \\frac{dC_D}{dt}=-\\frac{F}{V}C_D+k_3C_A^2 \\\\\n",
    "& V \\leq 100 \\\\\n",
    "& C_A(t=0) = 10000, ~~ C_B(t=0) = 0 \\\\\n",
    "& C_C(t=0) = 0, ~~ C_D(t=0) = 0 \\\\\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Nj_h-L7NpMYK"
   },
   "source": [
    "### Degrees of Freedom Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the optimization model formulated above:\\\n",
    "Number of variables = 5 \\\n",
    "Number of equality constraints = 4 \\\n",
    "Therefore, the degrees of freedom = 5 - 4 = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iAhG6jxzu5wA"
   },
   "source": [
    "### Pyomo Implementation\n",
    "We can use `Pyomo.dae` to automatically formulate the collocation equations (i.e., add constraints that numerically integrate the ODE model)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "aYW5M30J1JlR",
    "outputId": "5cab20c0-f777-4c31-f9b5-f6ac32c054c1"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ipopt 3.13.2: \n",
      "\n",
      "******************************************************************************\n",
      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
      "         For more information visit http://projects.coin-or.org/Ipopt\n",
      "\n",
      "This version of Ipopt was compiled from source code available at\n",
      "    https://github.com/IDAES/Ipopt as part of the Institute for the Design of\n",
      "    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE\n",
      "    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.\n",
      "\n",
      "This version of Ipopt was compiled using HSL, a collection of Fortran codes\n",
      "    for large-scale scientific computation.  All technical papers, sales and\n",
      "    publicity material resulting from use of the HSL codes within IPOPT must\n",
      "    contain the following acknowledgement:\n",
      "        HSL, a collection of Fortran codes for large-scale scientific\n",
      "        computation. See http://www.hsl.rl.ac.uk.\n",
      "******************************************************************************\n",
      "\n",
      "This is Ipopt version 3.13.2, running with linear solver ma27.\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:    12619\n",
      "Number of nonzeros in inequality constraint Jacobian.:        0\n",
      "Number of nonzeros in Lagrangian Hessian.............:     1006\n",
      "\n",
      "Error in an AMPL evaluation. Run with \"halt_on_ampl_error yes\" to see details.\n",
      "Error evaluating Jacobian of equality constraints at user provided starting point.\n",
      "  No scaling factors for equality constraints computed!\n",
      "Total number of variables............................:     1609\n",
      "                     variables with only lower bounds:      804\n",
      "                variables with lower and upper bounds:        1\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:     1608\n",
      "Total number of inequality constraints...............:        0\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        0\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0 -9.9999900e-03 1.00e+06 1.56e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0\n",
      "   1 -9.9998726e-03 1.00e+06 8.10e+02  -1.0 1.01e+06    -  9.90e-07 7.52e-05h  1\n",
      "   2r-9.9998726e-03 1.00e+06 9.99e+02   6.0 0.00e+00    -  0.00e+00 3.10e-07R  3\n",
      "   3r-3.2550526e+00 8.36e+05 2.18e+09   6.0 1.09e+09    -  1.95e-04 1.95e-04f  1\n",
      "   4 -3.2756934e+00 7.01e+05 6.67e+06  -1.0 8.47e+05    -  5.03e-04 2.84e-01h  1\n",
      "   5 -3.3599969e+00 4.50e+05 4.28e+05  -1.0 2.51e+05    -  6.14e-02 9.85e-01h  1\n",
      "   6 -3.4374817e+00 2.50e+05 2.67e+05  -1.0 4.42e+05    -  2.83e-03 4.99e-01h  1\n",
      "   7 -3.6563863e+00 1.27e+05 2.54e+05  -1.0 2.24e+05    -  7.45e-02 9.90e-01h  1\n",
      "   8 -4.0376969e+00 6.23e+04 2.60e+05  -1.0 5.04e+04    -  4.10e-02 9.90e-01h  1\n",
      "   9 -4.2238169e+00 4.57e+04 2.56e+05  -1.0 1.87e+04    -  1.46e-02 3.26e-01h  1\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "  10 -4.8943726e+00 2.20e+04 2.14e+05  -1.0 1.48e+04    -  1.98e-02 9.98e-01h  1\n",
      "  11 -5.0999886e+00 1.69e+04 2.38e+06  -1.0 4.00e+03    -  1.90e-02 2.92e-01h  1\n",
      "  12 -5.8592444e+00 8.22e+03 8.43e+06  -1.0 3.31e+03    -  3.47e-02 1.00e+00h  1\n",
      "  13 -6.6724932e+00 7.12e+03 8.23e+06  -1.0 1.71e+03    -  3.91e-02 1.54e-01h  1\n",
      "  14 -2.5003886e+01 3.51e+03 8.41e+06  -1.0 1.55e+03    -  8.22e-02 1.00e+00h  1\n",
      "  15 -8.4301553e+01 3.08e+03 8.00e+06  -1.0 1.70e+03    -  4.83e-02 1.40e-01h  1\n",
      "  16 -1.2308262e+02 2.89e+03 7.06e+06  -1.0 1.70e+03    -  1.19e-01 6.44e-02h  1\n",
      "  17 -8.1386603e+02 1.44e+03 5.95e+06  -1.0 1.70e+03    -  1.57e-01 1.00e+00f  1\n",
      "  18 -2.0541698e+03 8.40e+02 3.71e+06  -1.0 1.77e+03    -  3.77e-01 7.14e-01f  1\n",
      "  19 -3.5049209e+03 4.20e+02 1.79e+06  -1.0 1.45e+03    -  5.17e-01 1.00e+00f  1\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "  20 -3.7251819e+03 2.85e+02 2.82e+04  -1.0 4.62e+02    -  9.84e-01 4.76e-01f  1\n",
      "  21 -4.0002085e+03 1.42e+02 2.93e+03  -1.0 2.75e+02    -  8.98e-01 1.00e+00f  1\n",
      "  22 -4.0087092e+03 1.01e+02 2.07e+03  -1.0 7.06e+01    -  1.00e+00 4.15e-01h  1\n",
      "  23 -4.0087453e+03 1.00e+02 2.09e+03  -1.0 1.00e+02    -  7.72e-01 2.93e-03h  1\n",
      "  24 -4.0087456e+03 1.00e+02 5.80e+03  -1.0 1.00e+02    -  4.72e-03 2.92e-05h  1\n",
      "  25 -4.0087598e+03 1.00e+02 2.12e+03  -1.0 3.24e+02    -  4.40e-07 4.37e-05f  1\n",
      "  26 -4.3343427e+03 4.70e+00 1.83e+05  -1.0 3.26e+02    -  4.36e-05 1.00e+00f  1\n",
      "  27 -4.3141499e+03 8.39e-03 9.09e+01  -1.0 2.02e+01    -  1.00e+00 1.00e+00h  1\n",
      "  28 -4.3141284e+03 2.41e-08 1.08e+05  -1.7 2.15e-02    -  9.87e-01 1.00e+00h  1\n",
      "  29 -4.3141354e+03 7.14e-08 2.90e+04  -1.7 7.02e-03    -  7.31e-01 1.00e+00f  1\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "  30 -4.3150872e+03 1.32e-03 2.82e+04  -1.7 9.52e-01    -  2.70e-02 1.00e+00f  1\n",
      "  31 -4.5157522e+03 2.76e+02 2.80e+04  -1.7 2.01e+02  -4.0 7.26e-03 1.00e+00f  1\n",
      "  32 -5.3180445e+03 1.11e+02 1.95e+04  -1.7 8.02e+02    -  3.05e-01 1.00e+00f  1\n",
      "  33 -5.0535790e+03 1.43e+00 1.26e+01  -1.7 2.64e+02  -4.5 1.00e+00 1.00e+00h  1\n",
      "  34 -5.0655151e+03 9.00e+00 2.23e+00  -1.7 6.01e+01    -  1.00e+00 1.00e+00h  1\n",
      "  35 -5.0601188e+03 1.92e-01 7.30e-02  -1.7 1.58e+01    -  1.00e+00 1.00e+00h  1\n",
      "  36 -5.0598090e+03 5.39e-03 4.07e-04  -2.5 1.48e+00    -  1.00e+00 1.00e+00h  1\n",
      "  37 -5.0598054e+03 1.71e-04 1.11e-06  -3.8 2.78e-01    -  1.00e+00 1.00e+00h  1\n",
      "  38 -5.0598053e+03 5.24e-07 2.29e-09  -5.7 1.55e-02    -  1.00e+00 1.00e+00h  1\n",
      "  39 -5.0598053e+03 1.94e-10 1.67e-13  -8.6 1.91e-04    -  1.00e+00 1.00e+00h  1\n",
      "\n",
      "Number of Iterations....: 39\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:  -5.0598052797570590e+03   -5.0598052797570590e+03\n",
      "Dual infeasibility......:   1.6745249521991457e-13    1.6745249521991457e-13\n",
      "Constraint violation....:   1.9372237147763371e-10    1.9372237147763371e-10\n",
      "Complementarity.........:   2.5063723623922813e-09    2.5063723623922813e-09\n",
      "Overall NLP error.......:   2.5063723623922813e-09    2.5063723623922813e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 43\n",
      "Number of objective gradient evaluations             = 40\n",
      "Number of equality constraint evaluations            = 43\n",
      "Number of inequality constraint evaluations          = 0\n",
      "Number of equality constraint Jacobian evaluations   = 41\n",
      "Number of inequality constraint Jacobian evaluations = 0\n",
      "Number of Lagrangian Hessian evaluations             = 39\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      0.358\n",
      "Total CPU secs in NLP function evaluations           =      0.005\n",
      "\n",
      "EXIT: Optimal Solution Found.\n",
      "Optimal reactor volume is: 11.35 m^3\n",
      "Maximum Concentration of C at optimal reactor volume is: 5.06 kmol/m^3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtYAAAIwCAYAAACiDSxZAAAAP3RFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMS5wb3N0MSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8kixA/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACxeElEQVR4nOzdd3RU1doG8OdMT51U0huQhN4RpAlItWDBLoL1E732gqCioHgRULjXhh30qqigoiAiAgKKFKWGEiCQRiCkZ1InU873x5BJzkx6piTk+a2VlZw25yUJ4WHnPXsLoiiKICIiIiKiVpG5uwAiIiIioksBgzURERERkQMwWBMREREROQCDNRERERGRAzBYExERERE5AIM1EREREZEDMFgTERERETkAgzURERERkQMwWBMREREROQCDNVETZGRkYP78+Rg7dizCw8Oh0Wig0WgQGRmJiRMnYtGiRcjIyHB3meQmgiBY32JjY91dTpNt27ZNUrvtm4+PD7p37477778fe/fudXe50Ol0eP7559GzZ094enpKaj148CDS0tIk+0aPHi25ft68eZLjK1eudMufoyVyc3OxcOFCjBs3DhEREdBoNPDx8UFcXBymTp2KTz75BHq93t1lOpzt17Q5b9V/Fxv7viByJIW7CyBqy/R6PWbNmoX33nsPRqPR7nhWVhaysrKwadMmLFq0CAUFBW6o8tIQGxuL9PR067Yoim6sxhI6x4wZY92eMWNGuwpijlBaWork5GQkJyfj008/xaxZs/D666+7rZ5rr70WO3bscNv93WXp0qV48cUXUVFRIdmv1+tRWlqKtLQ0fP/993jppZfw+eef48orr3RTpUTEYE1Uj8rKSowfPx5//vmnZL+Pjw8GDRoEb29v5OTk4NChQ6isrITZbHZTpeRuU6dOtX7cqVMnN1bSOp6enpg8eTIAoKSkBH///TcKCwsBWP6js2jRIiQkJODee+91eW3JycmSUK1QKDBmzBj4+voCAPz8/ODl5SX5WvTs2dPldTraE088gf/+97+SfSEhIRgwYAD0ej12796N8vJyAMC5c+cwYcIErF69GjfeeKM7ynU4269pte3btyMvL8+6PWjQIMTExEjOqf67eCl+X1AbJhJRne677z4RgPVNEATx5ZdfFisqKiTnlZeXi5988onYu3dvN1V6aYiJiZF8vt3t999/l9QzY8YMd5fkcLZ/xpiYGMnxvLw8sVevXpJzunfv7pZat23bJqnjtttua/ZrvPzyy5LXWLFiheMLdaDVq1dL6gUgvvjii6LBYLCek5+fL06ePFlyjqenp5ienu7Gyp3viiuuaFdfS+o43P+vF1EblJSUJMpkMskP7ldeeaXBayorK+vcv3nzZvGOO+4QO3fuLHp6eopqtVqMiooSb7jhBnHNmjWiyWSyu2bFihWSe7/88svi+fPnxccee0yMjY0VVSqVGBISIt59993i2bNn663p7Nmz4ksvvSQOGzZMDAwMFBUKhRgQECD27t1bfOihh8Tk5GS7awoLC8XFixeLo0aNsl7j7+8vDh8+XFy6dKlYWlpqd01qaqqk3iuuuEKsrKwUlyxZIvbp00fUaDSir6+vOHHiRHHXrl2Sa20DdX1v1WbMmCHZ//vvv4tbt24VJ02aJAYEBIiCIFj/kT106JA4a9YsccKECWLXrl3FgIAAUaFQiN7e3mJiYqI4ffp0cceOHZJ6bMNmfW+1g3ZD4bRaRkaGOGfOHHHQoEGin5+f9WsxbNgwccGCBWJubm6d19m+tslkEj/66CNxyJAhopeXl+jl5SWOGDFC3LBhQ53XN6SxYC2Kovjpp5/a/dlLSkpEUaz7615SUiI+//zzYkJCgqhWq+1es6SkRHzrrbfEK6+8UuzUqZOoVCpFX19fsXfv3uKjjz4qHjt2rMEa63qrvkdd9dTWlGBtNpvFdevWiTfffLMYExMjajQa0cPDQ0xISBBnzpwpHj9+vM7PZWlpqbhkyRJx5MiRYnBwsKhUKkUvLy8xOjpaHDFihPjEE0+I69ata9LXpVq3bt0k9d500011nldWViZGRUVJzp05c6YoiqJ47tw5UaFQWPcPHjy4zteYPXu25PqvvvpKctxRPxca+/5oquYE65Z8X+zevVu86qqrRD8/P9HHx0ccPXq0uHXrVus1q1atEocOHSp6enqKvr6+4qRJk8Q9e/bUW0NLPn/UPjFYE9XhxRdflPygDQ4Orjc410ev14u33npro6FgzJgxYmFhoeRa22B91VVXiYGBgfWGCtvrq1/D09OzwXvb/mP0xx9/iKGhoQ1eEx8fL544cUJyne0/XL169RIHDBhQ5/VqtVrcvXu39drWButp06bV++dasmRJk1573rx51td3RrD+8ssvRS8vrwZfLygoSNy8ebPdtbXPCQkJESdMmFDn9YIgiN9//73d9Q1pSrD++eef7e517ty5Or/uffv2FXv37l3vax48eFCMjY1t8POgUCjEN954o1lfD0cFa51OZzf6a/umVCrF999/X3JdZWWlOHDgwEbrHDhwYJO/NkePHrW7fufOnfWev2DBArvvlWrXX3+95Jjtfw7MZrMkmAcGBkp+3jnq50Jj3x/N4cxgfe2114pyudzuzyiXy8WffvpJfOKJJ+r8HGg0GnHv3r1292/p54/aJ/ZYE9Vh586dku0rr7wSarW6Wa/x8MMP45tvvrFuKxQKDBw4EGq1Gnv37kVlZSUA4Pfff8fNN9+M3377rd7X2rBhAwCgf//+8Pb2xl9//QWTyQQASE9Px3vvvYfnn3/eev7atWtx7733Sh4A9PHxQZ8+feDv74/jx4/j9OnTknucPn0aV199NXQ6nXVfr169EBsbi9TUVBw9ehQAcOrUKUyePBlJSUnw9PSss94jR44AsDyQGB8fjz179lhfV6/XY+7cudi0aRMA4KqrrkJOTg5++eUXa68ogDr7KuvyxRdfALD0TcbGxiIlJcXunK5duyIkJAT+/v4wm804d+4cDh8+bO2LnzdvHqZMmYL+/fsjODgYU6dORW5urqSnNyYmBoMGDbJuDx48uEn1bdu2DdOnT7d+vQAgLi4OCQkJSEpKwrlz5wAAeXl5uO6667Bv3z4kJibW+VoXLlzApk2bEBYWhl69euHAgQPWPlNRFPHcc8/hhhtuaFJdTbV//37JtlKpRGBgYJ3nHjp0CICl33nAgAEQRREXLlwAYPnzTZw40boNAIGBgRgwYACysrJw7NgxAIDRaMQzzzyD0NBQ3HnnnU36ejiqr/3222/HL7/8Yt0ODg7GwIEDodfrsXPnTlRVVcFgMOChhx5CdHS0tR/9+++/x759+6zXVfdAA5YHnFNTU1FSUtKsWvbs2SPZVqlUDX7PDR8+XLJ94cIFpKenIyYmBv/3f/+HtWvXWo/973//w2uvvWbd3rZtGzIzM63bM2bMsP68c+TPhYa+P9qSdevWwcPDA0OHDkVmZqb1Z4rJZMJtt92G8vJyaLVaXHbZZUhKSkJ2djYAy3M5L774In799Vfrazny80fthHtzPVHb1KNHD8lowuzZs5t1/bFjx0RBECSjcNu3b7ceT0pKErVareQeGzdutB63HbGGzYiM7fExY8ZYj5nNZrtRweuuu07Mz8+X1Lh3717J6IrtyO+qVask5//73/+WHK89qmg7IgRAvPfee0Wj0SiKoigmJyeLKpXKekylUolVVVWS129qj7XtiLVCoRDXrl0rOad6tC0jI0PMycmp83XWr18veZ3nnntOcrw5Pda1z7MdgRs6dKjk+EMPPWRt/6moqBCvvvpqyXHb3mHbz+ukSZPE8vJyURRFMTs7W+zUqZPkeHN6axsasS4pKRG//fZb0dfX1+7+1er6uo8fP17yG5Tqr4Vtq8GQIUMk57366quS4xEREZI2qaZ8PVozYr1582bJsSlTpoh6vd56/MSJE6K3t7f1eK9evazHXnvtNet+Hx8fsaysTHJfo9Eo7ty5s1l9wIsWLZLUExYW1uD5ycnJdl+L6r/fJpNJ8vcrOjpaNJvN1mvvueceyXW1W8Qc/XOhvu+P5nLmiLWXl5d4+PBhURQtf0dt22wiIiLEzMxMURRF8fz586JarbYeU6vVkp9trfn8UfvEeayJmkBs5tRv69evl1wzdepUjBo1yrrdq1cv/N///Z/kmnXr1tX7ekOGDMHdd99t3Z4yZYrkeFZWlvXj/fv3Iy0tzbqt1Wrx2WefISAgQHLN4MGDrSNgZrMZP/30k/WYSqXCmjVrcNNNN1nftm3b1uR6NRoN3njjDcjlcgBAYmKiZBS2qqpK8kR/a8yYMQPXXXedZF/1aFtUVBT27duHO++8E926dYOPjw/kcjkEQcA111wjuSY5Odkh9dSWk5MjGXlUqVRYuHAhZDLLj16NRoPFixdLrtmwYUODM8wsW7YMHh4eACwjo0OGDJEcr/290Fzp6emSOaxvueUWyUibh4cH/v3vf9d7vVwux4cffgg/Pz/rvuqvRe3vL8DyW4La582ePRvh4eGSP4ftaLkz/fDDD5LtvLw83HHHHdbv/+effx5KpdJ6/MiRI9a/Z7VnoygpKcHTTz+Nr776Cnv37kVhYSHkcjmGDRsm+TvsaA39jJLJZLjvvvus2xkZGdi+fTsAoKKiAt9995312KhRo6x/Vx39c6Gh74+25NZbb0Xv3r0BWP6ODhw4UHL8wQcfRGRkJAAgNDRUMsuIXq+3/mxz9OeP2ge2ghDVISQkxPqraQCSoNoUtudX/5CurW/fvpLt1NTUel/P9lfAWq1Wsl17YYgzZ85IjvXr18/ufFv5+fmSAFVVVSX5x7YuDdXbtWtX+Pv7N7nm1mhosYfHH38cb731VpNep7i42CH11Jaeni4JPNHR0Xafh+7du0OlUqGqqgqAZRGU/Px8BAcH272et7c3unXrJtnnrM+rrc6dO2PFihXo379/vefExsbWu0BOY38nFAoFevToYW2NASzfY7Xbb5zJ9vv5r7/+atI1sbGxmDp1Kt544w0cPHgQAPD+++/j/ffft54XFxeHq666Cs8880yTFxCy/frn5eXBYDBIwn1t1e0ItdVukbnvvvswf/58a0vS559/jtGjR+PHH3+U/N2v/R9+R/9caOj7oy2x/d708fGRbPfq1avB49V/Bx39+aP2gcGaqA7Dhw/H77//bt3esmUL9Hp9k0dXbEePBEFoVT22Pa3VI8HuVFZWVu+xunpwnVVz7VHO2v755x+7UB0fH4/ExESo1WqUl5dL+mmb+1uJpnD29wHg2M9r7XmsBUGAt7c3oqKiMGzYMIwfP77Re9X3tQAc/7loC6r/Dmg0Gvz111/4+OOPsXbtWuzbt0/yH7XU1FS8++67+Oqrr3DgwAG7+Zbrctlll0m2DQYD9u7da9dLXc32uZBOnTohOjrauh0eHo6rr77aOoL63Xff4d1338X//vc/6zkBAQG46aabGq2tIQ39XGjo+6MtqT2iDsD6G6ZqtoMGjtTQ54/aB7aCENXh1ltvlfwwzcvLs/uVva3aI4VxcXGSY0lJSXbnHz58WLJte01Lde7cWbJ98ODBRkdjAwMDJaMuvr6+0Ov1EC0zB9X55qhWjmotDVq2/+hV++OPPyTbDz30EE6ePIl169ZhzZo1mDt3rlPqqc12dC4jI0MyggVYWlCqR6sBy+hXfQ8HOltwcDDWrFmDNWvWYPXq1VixYgVeeeUVTJo0qUkBvr6vBdD43wmj0Sj5LVFd1ziT7b2+/vrrBr//RVGUtBN5eHjg0UcfxZYtW1BUVIT8/Hzs2bNHMgJcWFiIFStWNKmenj172j3EartQTLWKigp8+OGHkn033HCD3fdw7Vp0Oh0+/PBD60PEADB9+nTJ4IGjfy409P1xKWoLP1fJ9TrWdzlRE/Xq1cuuH/Lll1/G/PnzrbN5VKuoqMAnn3wiade4+uqrJf+offfdd5IRpWPHjtn9Q2jb89tSAwYMkIxUFRcXY8aMGXbLrR88eBB///03AMs/eLXvr9Pp8NRTT9m1FYiiiD179uCJJ56w60ltreq+4Wqt6RUGLCN8tdV+0r64uFgyi4qz6unUqZNk5FGv1+P555+39lDr9XrMnj1bcs1VV111SQYQ2+/v+fPnS/7Dt2TJEkkbSHh4uHVmDVewfW5h7ty5df5aPisrC++++y4effRR676DBw/igw8+kNQfEBCAyy67zG4EuK6WjfrMnz9fsr169WrMmzdPMsNMYWEhbrnlFmRkZFj3eXh4YM6cOXavN3nyZERFRVm3Z8+eDaPRaN22fe6jLfxcaM/4+euY2ApCVI933nkHJ0+etC5pLooi5s2bhzfffBODBw+2Lml+8OBBVFZWSnpde/TogenTp+Ozzz4DYAl5o0ePxuDBg6FSqbB3715UVFRYzx8zZgwmTZrkkLoFQcCbb76Jm2++2brvxx9/RGxsLPr27Qt/f3+cPHkSJ06cwIoVK6z/IZg3bx7WrVuH0tJSAMC7776LVatWoW/fvvDx8UFeXh6OHj1qDUP9+vVzSL3VunXrhuPHj1u3L7/8cvTv3x9KpRKXX345nn766Wa93tChQyXbb775Jnbs2IHAwEDrQ2UNiY+Ph0wms4bgzZs34/LLL0dERAQAYM6cOXYPNdVl4cKFGD9+vPV13n33Xfzyyy/W6fZqB3ZPT0+8/PLLzfpzthdPP/00VqxYgdzcXADArl270LVrV+t0e9XTjlWr/ZCnK0yYMAHjx4+3Tnt56tQpxMfHY8CAAQgLC0N5eTlSUlKsveJXXHGF9dq0tDTMnDkTDz30ELp06YK4uDh4eXmhoKDAbtq87t27N7mmW2+9FX/88Qfeffdd67758+fjgw8+sC5pvmvXLsk0lTKZDJ999lmd7SbVDzHOmzcPACSDBCNGjKizNnf/XGjv+PnreBisierh4eGBzZs349lnn8V7771nHSUqKSnB1q1b7c63DQEffPABysrKsGbNGgCWX3Xv2rXL7rpRo0ZZz3GUm266CR999BEee+wxa4AvKSmx/iehLgkJCVi/fj1uu+0266haQUGBpNe8NoXCsT8+7r//fsloTWZmpmRu3eYaNWoUbrzxRnz//ffWfdUj9HK5HK+//jqee+65eq/39/fH1KlTsXr1auu+3bt3Wz9u6gwPY8eOxcqVK/Hggw9avxZnzpyxe8g0ICAAq1atalbwak86deqEjRs34oYbbrCOrubl5UlaEQDL1+a1117D9OnTXV7jmjVrcMstt1jnITaZTNbvGVt1ff+LooiUlJQ651IHLL9Nuv/++5tV09tvv42YmBi89NJL1iCcnZ1tndu+ttDQUKxcuRITJ06s9/Xuu+8+vPrqq5JRb8B+tLqau38utHf8/HU8/AoSNUCtVuOtt96yjrZt27YNJ06cQGFhIURRRHBwMHr06IGxY8fijjvusLt29erV2LRpE1auXIndu3cjOzsbJpMJwcHBGDRokHU6L2eMzN1///2YOHEiPvzwQ2zevBknT56ETqeDj48PIiIiMGLECFx++eWSa6644gokJydjxYoVWL9+PZKSklBYWAiZTIagoCAkJCRg2LBhuPbaa+2meWutq666Ct988w3++9//4vDhw9YRntb45ptv8MYbb2DlypU4c+YMfHx8MGTIELzwwguIiIhoMFgDwKefforo6GisXbsWmZmZkl7o5rjrrrswatQoLF++HJs3b0ZKSgrKysrg4+ODxMRETJ48GTNnznTYQidt1YABA3DkyBF88skn+Omnn5CUlISioiJoNBrExMRg9OjReOihhyTTl7mSr68vNm7ciA0bNuDLL7/Enj17cP78eej1evj6+iI2NhYDBgzA+PHjce2111qvGzFiBN5//33s2rULBw4cQG5uLvLz8yGKIoKCgtCzZ09cd911uP/++6HRaJpVkyAIePbZZzF9+nR8/PHH2LJlC5KTk1FQUAC5XI7g4GD0798fV199NaZNm9bo60dGRmLy5MlYv369dZ+/v7/kN1y23Plz4VLAz1/HIojOeBSeiIiIiKiDufSekCEiIiIicgMGayIiIiIiB2CwJiIiIiJyAAZrIiIiIiIHYLAmIiIiInIABmsiIiIiIgfgPNZuZDabce7cOfj4+EiWvyYiIiKitkEURZSUlCA8PLzRdScYrN3o3LlziIqKcncZRERERNSIzMxMREZGNngOg7Ub+fj4ALB8oXx9fd1cDRERERHZ0ul0iIqKsua2hjBYu1F1+4evry+DNREREVEb1pS2XT68SERERETkAAzWREREREQOwGBNREREROQADNZERERERA7AYE1ERERE5ACcFYSIiIiojTAYDDCZTO4u45Inl8uhVCod/roM1kRERERuptPpkJeXB71e7+5SOgy1Wo2goCCHTnnMYE1ERETkRjqdDllZWfD29kZQUBCUSmWT5kymlhFFEQaDAcXFxcjKygIAh4VrBmsiIiIiN8rLy4O3tzciIyMZqF3Ew8MDPj4+OHv2LPLy8hwWrPnwIhEREZGbGAwG6PV6aLVahmoXEwQBWq0Wer0eBoPBIa/JYE1ERETkJtUPKjrjQTpqXPXn3VEPjDJYExEREbkZR6vdw9Gf9w4frL/44gs8+OCDGDRoENRqNQRBwMqVK+s9X6fT4amnnkJMTAzUajViY2Px7LPPorS01HVFExEREVGb0+EfXnzxxReRnp6OoKAghIWFIT09vd5zy8rKcMUVV+DgwYOYMGECbr/9dhw4cABvvPEGtm/fjh07dkCj0biweiIiIiJqKzr8iPXHH3+MtLQ05ObmYubMmQ2eu3jxYhw8eBDPPfccfv31V7z++uv49ddf8dxzz+Hvv//GsmXLXFQ1EREREbU1HT5Yjxs3DjExMY2eJ4oiPv74Y3h7e2Pu3LmSY3PnzoW3tzc+/vhjZ5XZaik5JdiQdB7f7z+LVXszsONkrrtLIiIiIrqkdPhg3VSnTp3CuXPnMHz4cHh5eUmOeXl5Yfjw4Thz5gwyMzPdVGHD1h06j4e/3I+nvj2EOd8nYdXeDHeXRERERCSxb98+3HfffYiPj4eXlxc8PDzQpUsX3HXXXfjtt9/cXV6jGKyb6NSpUwCA+Pj4Oo9X768+ry56vR46nU7y5ipqpfRLrTeaXXZvIiIiooaYzWY89dRTGDRoED7//HN07twZM2fOxOOPP46BAwfi559/xoQJE/Dqq6+6u9QGdfiHF5uquLgYAKDVaus8Xr1iT/V5dVm4cCHmz5/v+OKaQKOQS7YrDY6Zr5GIiIiotV588UUsW7YM/fr1w5o1a9ClSxfJ8YqKCrzzzjvIz893U4VNw2DtQnPmzMFTTz1l3dbpdIiKinLJvTliTURERG1RSkoKFi9ejMDAQGzcuBEhISF253h4eODZZ5+FXq93Q4VNx2DdRNUj1fWNSFe3ddQ3og0AarUaarXa8cU1AUesiYiI2g+zWURheZW7y2gyf08VZLKWLbaycuVKmEwmPPjgg3WG6trclaOaisG6iRrroW6sB9vdOGJNRETUfhSWV2Hggs3uLqPJ9r04DoHeLQu9O3fuBACMHTvWkSW5BR9ebKL4+HiEh4dj586dKCsrkxwrKyvDzp07ERcX57LWjubiiDURERG1RdnZ2QCAyMhIN1fSegzWTSQIAu6//36UlpbaPZH66quvorS0FA888ICbqmscR6yJiIjoUvPLL79ALpfj/vvvd3cpANgKgo8//hh//vknACApKcm6b9u2bQCAESNGWL9Ys2bNwo8//ohFixbhwIEDGDBgAPbv349NmzZh8ODBeOKJJ9zxR2gSjZIj1kRERNT2hIaGIjk5GVlZWUhMTGzWtStWrMDTTz+NDz/8EG+99RY8PT2dVGXTdPhg/eeff+Kzzz6T7Nu5c6e13weANVh7eXlh+/btmDdvHr777jv8/vvvCAsLw9NPP42XX34ZHh4eLq29OdQKjlgTERG1F/6eKux7cZy7y2gyf09Vi68dPnw4tm3bhi1btjSrz7qgoAAbNmxAamoqdu/eje+++w533XVXi+twBEEURdGtFXRgOp0OWq0WxcXF1nmwneXkhRJMWLZDsu/Mv69q8RO8RERE1HqVlZVITU1FXFwcNBqNu8txi5SUFHTr1g3+/v44duwYgoOD6z1Xr9dbZwZ555138PPPP+OXX37B+++/j2+//RZbt25t1r2b8vlvTl5jj3UHYTtiDQBVJo5aExERkXt17doVs2bNQl5eHiZPnozU1FS7cyorK7F06VLMmzfPum/FihW48847AQA333wz/vrrL6Slpbmo6rp1+FaQjsK2xxqw9FnXtZ+IiIjIlRYsWIDKykosW7YMiYmJGDt2LHr16gWlUonU1FRs3rwZ+fn5WLBgAQDg8OHDSE5OxvXXXw8ACAwMxPjx47Fy5UpJ+HY1BusOoq4R60oDR6yJiIjI/WQyGZYuXYo77rgDy5cvx44dO7Bjxw6YzWaEhYVh4sSJuOeeezBunKXv/NNPP0V5eTn8/Pysr2E2m3H06FG8/PLLEAT3tLoyWHcQdY1M642cGYSIiIjajkGDBuGTTz5p8ByDwYAvv/wS7733HkaOHGndbzKZrA9Cjhkzxtml1onBuoOoa8S6glPuERERUTuzbt06lJWVYcaMGXbT602aNAkrVqxwW7Dmw4sdhCAI8LAZta6oYrAmIiKi9mXFihWYMGFCnXNWT506Fd999x1KSkrcUBlHrDsUD5VcMkrNHmsiIiJqb9atW1fvsdtvvx233367C6uR4oh1B6KxaQfh6otEREREjsNg3YFoVDatIAzWRERERA7DYN2B2PZYc8SaiIiIyHEYrDsQ2yn3OGJNRERE5DgM1h0IZwUhIiIich4G6w7EdsRab+SsIERERESOwmDdgWiU0i83R6yJiIiIHIfBugOxawVhjzURERGRwzBYdyAeKs4KQkREROQsDNYdCGcFISIiInIeBusOxDZYc8SaiIiIyHEYrDsQ+wViOCsIERERkaMwWHcgnBWEiIiIyHkYrDsQzgpCREREbdm+fftw3333IT4+Hl5eXvDw8ECXLl1w11134bfffnN3eY1isO5AOCsIERERtUVmsxlPPfUUBg0ahM8//xydO3fGzJkz8fjjj2PgwIH4+eefMWHCBLz66qvuLrVBCncXQK6jVjBYExERUdvz4osvYtmyZejXrx/WrFmDLl26SI5XVFTgnXfeQX5+vpsqbBoG6w7EdsSarSBERETkbikpKVi8eDECAwOxceNGhISE2J3j4eGBZ599Fnq93g0VNh2DdQdi12PNhxeJiIjaJrMZqChwdxVN5xEAyFrWYbxy5UqYTCY8+OCDdYbq2tRqdYvu4SoM1h2I7awglUZOt0dERNQmVRQAS7o0fl5b8expwCuoRZfu3LkTADB27FhHVuQWfHixA7Edsa4ymmEyi26qhoiIiAjIzs4GAERGRrq5ktZjsO5AbFdeBPgAIxEREbU/d999NwRBsL4FBwfjxhtvRFpamlvrYrDuQBisiYiIqK0JDQ0FAGRlZTXrumuuuQbnz5/HuXPnsGHDBuTk5GDatGnOKLHJ2GPdgdjOCgJwZhAiIqI2ySPA0rfcXngEtPjS4cOHY9u2bdiyZUuz+qzVarU1lIeFheGRRx7Bfffd1+I6HIHBugPRKOx/QcERayIiojZIJmvxw4Dtzd13343XX38dH374IZ544gkEBwfXe65er69zZpCioiKsXr0aQ4YMcWapjWIrSAeikMuglAuSfZUGzgxCRERE7tO1a1fMmjULeXl5mDx5MlJTU+3OqaysxNKlSzFv3jzrvrVr18Lb2xteXl7w9/fH4cOHsXLlStcVXgeOWHcwGqUcBpPRus1WECIiInK3BQsWoLKyEsuWLUNiYiLGjh2LXr16QalUIjU1FZs3b0Z+fj4WLFhgvWb8+PF4++23AQC5ublYuHAhbrjhBuzZswcKhXsiLoN1B+OhlKOkslaw5iIxRERE5GYymQxLly7FHXfcgeXLl2PHjh3YsWMHzGYzwsLCMHHiRNxzzz0YN26c9RovLy907doVgGXU+6OPPkJoaCi2bt2KCRMmuOXPwWDdwdjODMIeayIiImorBg0ahE8++aRF18rlloxTWVnpyJKahcG6g7Fb1pzBmoiIiNohvV5vXVwmNzcXCxYsQEBAAIYNG+a2mhisOxiNiiPWRERE1P6tX78eYWFhAAB/f38MHDgQGzduRFCQ+2ZTYbDuYGyn3OOsIERERNTerFy50u0zgNSF0+11MLaLxLAVhIiIiMgxGKw7GLsea84KQkREROQQDNYdjN2sIEYGayIiIiJHYLDuYOyCNUesiYiIiByCwbqD4XR7RERERM7BYN3BaJScFYSIiIjIGRisOxiOWBMRERE5B4N1B2M73R4XiCEiIiJyDAbrDkZt+/AigzURERGRQzBYdzC2rSDlnBWEiIiIyCEYrDsYL9uVFxmsiYiIiByCwbqD8VQrJNuleqObKiEiIiK6tDBYdzDearaCEBERETkDg3UH46mSjliXVRkhiqKbqiEiIiKqsW/fPtx3332Ij4+Hl5cXPDw80KVLF9x111347bff3F1eoxisOxhvm1YQUeSoNREREbmX2WzGU089hUGDBuHzzz9H586dMXPmTDz++OMYOHAgfv75Z0yYMAGvvvqqu0ttkKLxU+hS4mnz8CJgGbX2UvNbgYiIiNzjxRdfxLJly9CvXz+sWbMGXbp0kRyvqKjAO++8g/z8fDdV2DRMUx1MXQG6TG8CfNxQDBEREXV4KSkpWLx4MQIDA7Fx40aEhITYnePh4YFnn30Wer3eDRU2HYN1B6NWyCCXCTCZa/qqyzgzCBERUZtiFs0o0he5u4wm81P7QSa0rMN45cqVMJlMePDBB+sM1bWp1eoW3cNVGKw7GEEQ4KWSQ1dZE6YZrImIiNqWIn0RrvjmCneX0WTbb92OAE1Ai67duXMnAGDs2LGOLMkt+PBiB2T7ACMfXiQiIiJ3yc7OBgBERka6uZLWY7DugLhIDBEREbV3Z8+exYMPPojo6Gio1WrExcVh2rRpSEpKcltNDNYdkO0DjGwFISIiIncJDQ0FAGRlZTX5mtOnT2PgwIE4evQo3n//fRw/fhzfffcdYmNjMXv2bGeV2ij2WHdAXjZT7pWxFYSIiKhN8VP7Yfut291dRpP5qf1afO3w4cOxbds2bNmypcl91g8//DCioqKwY8cOyGQ148QDBgxAUVFRi2tpLQbrDogj1kRERG2bTJC1+GHA9ubuu+/G66+/jg8//BBPPPEEgoOD6z1Xr9ejtLQUv/32G7788ktJqK7m5+fnxGobxlaQDsh+xJrBmoiIiNyja9eumDVrFvLy8jB58mSkpqbanVNZWYmlS5di3rx5SElJgSiK6NatmxuqbRhHrDsgjlgTERFRW7JgwQJUVlZi2bJlSExMxNixY9GrVy8olUqkpqZi8+bNyM/Px4IFC9xdaoMYrDsgu+n29OyxJiIiIveRyWRYunQp7rjjDixfvhw7duzAjh07YDabERYWhokTJ+Kee+7BuHHjkJ+fD0EQkJycjP79+7u7dAkG6w7IU8Xp9oiIiKjtGTRoED755JMGzwkMDMS4cePw5ptv4tZbb7Xrsy4qKnJbnzV7rDsgLzV7rImIiKj9evfdd5GRkYFRo0bhl19+wZkzZ3Dw4EG8/PLLuPPOO91WF4N1B2TfY81WECIiImo/4uPjsW/fPnTv3h0PPPAAunfvjuuvvx6pqalYvHix2+pyaStIcXExLly4gIKCAvj7+yM0NBRardaVJRD48CIRERG1f1FRUfjoo4/cXYaE00est27diocffhgJCQkICAhA9+7dMXz4cPTo0QMBAQFISEjAww8/jK1btzq7FIcRRRHff/89xowZg7CwMHh6eiIxMREPPvggzpw54+7yGmU73V45F4ghIiIiajWnBetPP/0UvXr1wvjx4/HBBx9Y5xy0fUtJScEHH3yA8ePHo1evXlixYgVEUXRWWQ7xzDPPYOrUqThx4gSuv/56PProo4iLi8NHH32Efv364ciRI+4usUG2I9Z8eJGIiIio9RwerP/8808MGDAADzzwAI4fP24NyYIg1Hl+9X5RFHHs2DHcf//9GDRoEP78809Hl+YQ2dnZ+M9//oOYmBgcP34cy5cvx6JFi7Bx40a8+eabKCkpwdKlS91dZoPsptvjw4tERERErebwHutRo0ZJwnK1Hj16oHfv3ggKCoKvry+Ki4uRl5eHpKQkHD9+XPIaBw4cwOjRo2E0tr3Al5aWBrPZjOHDh9v1h19zzTV46qmnkJub66bqmsbTphXEYBKhN5qgVsjruYKIiIiIGuOUhxdFUYRKpcJVV12FadOmYdy4cfD19a33/OLiYmzevBlffPEFNmzYAIPB0GbbQeLj46FSqbBz507odDrJn2v9+vUAgCuvvNJd5TWJ7Yg1YFkkhsGaiIiIqOUcHqy9vLzwyCOP4Omnn0ZQUFCTrtFqtZg6dSqmTp2K3NxcLFmyBMuXL3d0aQ4RGBiI119/HU8//TS6deuG6667Dr6+vjh06JD1Qc1HHnmkzmv1ej30er11W6fTuapsCc86gnWp3gh/L5UbqiEiIiK6NDg8WKelpSEwMLDF1wcHB2Px4sV47rnnHFiVYz355JOIiIjA/fffj/fff9+6f8SIEbjjjjugUNT9aV24cCHmz5/vqjLr5am0H5nmIjFERERErePwhxdbE6qd8TrO8Morr2DatGl4/vnnkZmZiZKSEvzxxx+orKzE6NGj8dNPP9V53Zw5c1BcXGx9y8zMdHHlFjKZYNdnzUViiIiIiFqHKy820+bNm/Hyyy/jkUcewezZsxEZGQlvb2+MGDEC69atg1KpxNNPP13ntWq1Gr6+vpI3d+EiMURERESO5dZgLZe3v4flfvnlFwDAmDFj7I6FhoaiW7duSElJQWlpqatLaxZOuUdERETkWG4N1m115o+GVFVVAUC9U+rl5uZCJpNBqVS6sqxms20FKWUrCBEREVGruDVY17doTFs2fPhwAMDSpUtRXFwsOfb+++/j7NmzuPzyy6FWq91RXpPZtoJwxJqIiIiodZwyj/Wl7Oabb8by5cuxY8cOJCQkYMqUKfDz88P+/fuxdetWeHh4tPmVFwHAy27EmsGaiIiIqDWaFazHjh3b4hsJgoAtW7a0+Pq2Qi6XY9OmTVi2bBm+/fZbfPXVV6iqqkJISIh1ppDu3bu7u8xG8eFFIiIiamv27duH9957Dzt27MC5c+dgNpsRHh6OYcOGYfr06Rg/fry7S2xQs4L1tm3bWtS+IYpiu2z7qI9arcbs2bMxe/Zsd5fSYl4q22DNHmsiIiJyD7PZjGeeeQbLli2DQqHA2LFjMWXKFCiVSpw5cwY///wzvvjiC7zyyiuYO3euu8utF1tBOiiOWBMREVFb8eKLL2LZsmXo168f1qxZgy5dukiOV1RU4J133kF+fr6bKmyaZgfr9jiTB9nzVkt7rMurOGJNRERErpeSkoLFixcjMDAQGzduREhIiN05Hh4eePbZZ6HX691QYdM1K1inpqY6qw5yMU+bEWs+vEhERNR2iGYzTEVF7i6jyeR+fhBkLZtsbuXKlTCZTHjwwQfrDNW1tfVZ15oVrGNiYpxVB7kYp9sjIiJqu0xFRTg1bLi7y2iy+L92QhEQ0KJrd+7cCaB1k2S0FQ7tsd6/f791VDsuLg4DBgxw5MuTA9lPt8dWECIiInK97OxsAEBkZKSbK2k9hywQs2rVKkRHR2Pw4MG45ZZbcMstt2Dw4MGIiorCl19+6YhbkIPx4UUiIiJqr+6++24IggBBEKBSqRAWFoarrroK3377rVvranWwXrp0KaZNm4azZ89CFEXJW1ZWFqZPn4433njDEbWSA9lOt8dWECIiInKH0NBQAEBWVlazrrvmmmtw/vx5pKamYv369Rg2bBjuvfde3Hfffc4os0la1QqSnJyM2bNnS+aprp41pPb2888/j6uvvrpdLJzSUXipufIiERFRWyX380P8XzvdXUaTyf38Wnzt8OHDsW3bNmzZsqVZfdZqtdoayiMiIjBw4EAMHToU48ePx2233eaWxWRaFayXL18Oo9EIQRAgiiLCw8PRo0cPAMCxY8dw7tw5AIDJZMLy5cvx1ltvtb5icghvm1aQSoMZJrMIuezSWciHiIiovRJkshY/DNje3H333Xj99dfx4Ycf4oknnkBwcHC95+r1+gZnBhk3bhz69OmDH374wS3BulWtINu2bbN+/OqrryI9PR2bNm3Cpk2bkJ6ejtdee63Oc8n9bKfbA4AytoMQERGRi3Xt2hWzZs1CXl4eJk+eXOf0zpWVlVi6dCnmzZvX6Ot1794daWlpji+0CVo1Yp2eng5BENCjRw+88MILkmNyuRxz5szBV199haNHjyI9Pb1VhZJjeavsv/TlehN8NUo3VENEREQd2YIFC1BZWYlly5YhMTERY8eORa9evaBUKpGamorNmzcjPz8fCxYsaPS13LmYYauCdUVFBQBLX0t9IiIicPToUVRWVjb6eocOHYJSqbS2k4iiiB9//BEnT55E586dcf3110Oh4CrsjuBp02MNsM+aiIiI3EMmk2Hp0qW44447sHz5cuzYsQM7duyA2WxGWFgYJk6ciHvuuQfjxo1r9LWOHz+OkSNHuqBqe61Kqf7+/sjJycH+/fuh0+ng6+srOa7T6bB//34AgF8DTe35+fm48sorkZSUBAAYM2YMfvzxR1x//fXYunWr9bwhQ4Zgy5Yt8PDwaE3ZBEApl0GlkKHKaLbu48wgRERE5E6DBg3CJ5980uLrf/vtNyQlJWHp0qUOrKrpWtVj3bt3bwCWYDx58mT89ttvOH/+PLKzs/Hbb7/hqquuQl5eHgRBsJ5bl1dffRWHDx+2TtP3+++/48Ybb8SWLVsk0/ft2bOHD0A6kP0iMQzWRERE1D7o9XpkZ2cjKysL+/btw4IFC3DDDTfg3nvvbdLItjO0KlhPmTLF+vHu3bsxadIkREZGIiIiApMmTcKuXbvqPNfWhg0bIAgCAgICcOONN8Lf3x+bN2+2BvL//Oc/6NevH0RRxPfff9+akqkW+0ViuPoiERERtQ/r169HWFgY4uLicPXVV+Ovv/7Cp59+2qoR79ZqVSvIvffei6VLlyIjIwOAfbN49VzW0dHRDU7WnZmZCQD44YcfMHLkSOzYsQOjR4+GIAj49NNPMXDgQIwaNQoDBgzAiRMnWlMy1cJFYoiIiKg9WrlyJVauXOnuMuy0asTay8sLa9euRXBwsHWRmNpvoigiODgYP/zwA7y8vOp9neoAPmTIEADA0KFDrceqW0iqH2gsLy9vTclUCxeJISIiInKcVi9p3rdvXxw9ehTPPfccunfvDg8PD2g0GnTv3h2zZs3CkSNH0K9fvwZfIygoCACQl5cHAFCpVNZj1R8XFhYCsDwwSY5h2wpSzlYQIiIiohZzyNx1gYGBWLhwIRYuXNii63v27ImsrCzs3LkTN998MwDAbDZLzjl06BAAIC4urnXFkpVtKwhHrImIiIhark1MCv3cc89h3LhxDY5G//TTT9BqtW5ZnvJSZTdizR5rIiIiohZzSLA+d+4ctm7diqysLOj1+nrPe+mll+rcP3r0aIwePbrBe7zzzjt45513WlMm2bDvsWYrCBEREVFLtTpYP/vss/jvf/8Lk6nxUFZfsCb3sJ9ujyPWRERE7uDOZbg7Mkd/3lsVrN966y28+eabkn3VM3zUVj1jSGu1dmScpLxtgnVJpcFNlRAREXVMcrnlt8cGg4ErS7uBwWDJPtVfh9ZqVbD+6KOPAMA6tR7gvP9xcWTc8bQeSsl2cQWDNRERkSsplUqo1WoUFxfDx8fHIQOR1DSiKKK4uBhqtRpKpbLxC5qgVcE6JSXF+g1w3XXX4corr4SPj49DCqvN1SPjHYWfp/SbqIjBmoiIyOWCgoKQlZWFs2fPQqvVQqlUMs84kSiKMBgMKC4uRmlpKSIiIhz22q0K1r6+vsjLy8OgQYPwww8/OKomO64cGe9I/DxUkm0dgzUREZHL+fr6ArCs55GVleXmajoOtVqNiIgI6+ffEVoVrMeNG4dVq1Y5rC+lPq4aGe9obFtBisoNHPUnIiJyA19fX/j6+sJgMDSp7ZVaRy6XO6z9o7ZWBet58+bhp59+wt69e/H5559j+vTpjqpLwlUj4x2NbSuI0SyivMpkN1sIERERuYZSqXRK4CPXaFWCio+Px8aNGzFx4kTcc889WLBgAXr27FnnQi+CIOCTTz5p0X1cNTLe0fh62P/FLaowMFgTERERtUCrEpTJZMJbb72F8vJyAJaWjdOnT9udV91e0NJg7aqR8Y7GR62ATADMtdrVi8qrEOHH6X6IiIiImqtVwXrJkiVYvXq103tyXTUy3tHIZAJ8PZQoKq95aJFT7hERERG1TKuC9YoVK6wfO3OWDleNjHdEfrbBupzBmoiIiKglWhWsMzIyIAgCBEHAM888g+HDh8PHxwcymcxR9QFw3ch4R6T1VAH55dZtjlgTERERtUyrgnVISAgyMzNxxRVXYNGiRY6qyY6rRsY7Irsp9xisiYiIiFqkVcF66tSpWLZsGQoKChxVT51cNTLeEflxWXMiIiIih2hVsH7ppZewfv16HD58GPPnz8fs2bOhVqsdVZuVq0bGO6K6FokhIiIiouZrVbC+4YYboFKpIIoiXnnlFSxduhRdunSBn5+f3bmCIGDLli0tuo+rRsY7IttFYoorqtxUCREREVH71qpgvW3bNmuLhiiKKCkpwcGDB+0eMmztMtmuGhnviGxHrNkKQkRERNQyDltiz5kzdrhqZLwjYisIERERkWO0Oli7YpYOV42Md0R+nirJNkesiYiIiFqmVcE6NTXVUXU0GYOzY9m1gnDEmoiIiKhFWhWsY2Jimnyu0Whsza04f7WT2D68WKI3wmgyQyHnVIZEREREzdHqhxdHjx7d6HkGgwE33XQTfvzxxxbdxx0j4x2F7Yg1AOgqjQjwUtVxNhERERHVp1XB+vrrr8fvv/+O/v3713uOwWDADTfcgF9++aXF92nOyDg1T13Buqi8isGaiIiIqJla9ft+nU6HyZMn49SpU3Ue1+v1uP7667Fhw4bW3IacSKOUQ62QfhvwAUYiIiKi5mt1I21ubi4mTJiAc+fOSfZXh+qNGze29hbkZLZ91kUM1kRERETN1qpgrVJZ2gXS09Mxfvx468qIlZWVuPbaa/Hrr79az42Li2vNrciJ/DykbR86BmsiIiKiZmtVsF69ejUUCgUEQcDx48dx1VVXIScnB1dffTU2b95snXe6W7du2LFjh6NqJgfjIjFERERErdeqYH3ttdfi66+/hlwuhyAI+Pvvv9GlSxfrgi6iKKJfv37YsWMHwsPDHVUzOZjWk8uaExEREbVWq3usb7jhBnz11VeQy+UAgLKyMoiiCFEUMXToUPz+++8ICgpqdaHkPByxJiIiImo9h6wCctNNN+F///sfZDKZdenxsWPHYvPmzdBqtY64BTmRn22wrqhyUyVERERE7Vez5rEeO3Zsg8e1Wi0KCgogCALKy8txzTXXWI8JgoAtW7a0rEpyKttZQfjwIhEREVHzNStYV/dON0YURezZs0ey3ZTryD3YCkJERETUeq1aebEuDNDtj9ZTOt0eH14kIiIiar5mB2tRFJ1RB7mR3Yg1gzURERFRszUrWKempjqrDnIj24cXiysMbN8hIiIiaqZmBeuYmBiH3pyj322D7Yh1ldGMSoMZHiq5myoiIiIian8cMt1eS5nNZnfeni6ynRUE4JR7RERERM3l8GCdnJzcpl6HGuejUcK264MPMBIRERE1j8ODde/evXHnnXfi6NGjLbr+6NGjuP3229G7d28HV0b1kcsE+KilXUGcco+IiIioeRwerE0mE77++mv06dMHl112Gd5+++1GR5+Tk5Px1ltv4bLLLkOfPn3wzTffsE3Exfw45R4RERFRqzh8HmulUgmDwRLK9u3bh3379gEA/Pz80L17dwQFBcHX1xc6nQ55eXk4duwYiouLrddXP9CoUqnsX5ycxvYBxmKOWBMRERE1i8OD9ZEjR/Dkk09iw4YN1inbRFFEYWEhdu3aZXd+dZCuPg8ArrrqKvznP/9xdGnUANsHGDliTURERNQ8Dm8FiY+Px/r167Fp0yZMmDChydeJoojx48dj06ZNWL9+Pbp27ero0qgBvnaLxHBWECIiIqLmcPiIdbVx48Zh3LhxOHXqFL777jts3rwZ//zzD3Q6nfUcHx8fDBw4EGPHjsVNN92Ebt26OascaoTtIjF8eJGIiIioeZwWrKvFx8dj9uzZmD17NgCgsrIShYWF8PPzg4eHh7NvT03EVhAiIiKi1nF6sLal0WgQFhbm6ttSI+weXmSwJiIiImoWt668SG2Hnwen2yMiIiJqDQZrAlDHw4vssSYiIiJqFgbrVvjhhx8wfvx4BAYGQqPRIC4uDrfffjsyMzPdXVqzsceaiIiIqHVc3mN9KRBFETNnzsSHH36ILl264LbbboOPjw/OnTuH7du3Iz09HVFRUe4us1lse6x1lQaYzCLkMsFNFRERERG1LwzWLfDWW2/hww8/xMMPP4y33noLcrlcctxoNLqpspazHbEWRaCk0mC31DkRERER1Y2tIM1UUVGB+fPno3Pnzvjvf/9rF6oBQKFof/9fsX14EWA7CBEREVFztL8E6GabNm1CYWEh7rnnHphMJvz00084efIk/Pz8MG7cuHa7YqRGKYNKLkOVyWzdV1RuQEygG4siIiIiakcYrJtp3759AAC5XI4+ffrg5MmT1mMymQxPPvkk3njjjTqv1ev10Ov11u3aq1C6myAI0HoqkVtSUx9HrImIiIiajq0gzZSTkwMAWLp0KbRaLfbu3YuSkhLs2LEDCQkJePPNN7F8+fI6r124cCG0Wq31ra094Gj7AGMRgzURERFRk7k0WJtMJly4cAEZGRn1vrV1ZrOlVUKlUmHt2rUYPHgwvL29MXLkSKxevRoymQxvvvlmndfOmTMHxcXF1re2Ni2fn+3qi+VVbqqEiIiIqP1xSSvI33//jblz52L79u2oqqo/rAmC0OZn1NBqtQCAQYMGITw8XHKsV69e6Ny5M1JSUlBUVAQ/Pz/JcbVaDbVa7apSm83fS/oAY14pgzURERFRUzk9WO/ZswdjxoyBXq+HKIrOvp3TJSYmAoBdaK5Wvb+ioqLec9qqTj7S0J9Tq9+aiIiIiBrm9FaQefPmobKyEoBlRFoQ7BccqWtfWzVmzBgAwPHjx+2OGQwGpKSkwMvLC8HBwa4urdU6+Wgk27kllW6qhIiIiKj9cfqI9a5duyAIAkRRhEqlQmJiIrRabZ3zP7cHXbp0wYQJE7Bp0yZ8/PHHuP/++63HXn/9dRQVFWHatGntci7rYI5YExEREbWY09NfdU+1VqvFnj17kJCQ4OxbOt17772HYcOG4YEHHsDatWvRrVs3HDhwAFu3bkVMTAyWLFni7hJbxK4VRMdgTURERNRUTm8FqQ7SQ4YMuSRCNWAZtf7nn39w9913Y9++fXjrrbdw6tQp/Otf/8LevXsRGhrq7hJbpJOvNFjnlephNrf/vngiIiIiV3D6iPW0adMwa9YsnDp1CqIotqt+6oZERUVhxYoV7i7DoWx7rI1mEQXlVQjybrszmRARERG1FU4fsX7ssccwYMAApKWl4cEHH0RxcbGzb0ktFOStgu3/e3LZZ01ERETUJE4fsZ40aRIMBgNEUcQnn3yCL7/8El27dkVgYKDduYIgYMuWLc4uieqhkMsQ6KWSzF+dU6JH9zA3FkVERETUTjg9WG/bts06zZ4oiqioqEBSUpJdS8il1CbSngX7aKTBWscp94iIiIiawqVLmlcHbAbotouLxBARERG1jEsmW74UVlzsKGznsmaPNREREVHTOD1Yp6amOvsW5ED2I9ZsBSEiIiJqCqcH65iYGGffghyIi8QQERERtYzL193ev3+/dRQ7Li4OAwYMcHUJ1IBOvtK5rNljTURERNQ0LgvWq1atwnPPPYesrCzJ/vDwcLz++uu48847XVUKNaCuVhDO2EJERETUOJfMCrJ06VJMmzYNZ8+ehSiKkresrCxMnz4db7zxhitKoUbYrr5YaTCjVG90UzVERERE7YfTg3VycjJmz55tHfWsPfJZe37r559/HsePH3d2OdSITr72y5ezHYSIiIiocU4P1suXL4fRaLQG6LCwMIwbNw7jxo1DWFiYdSo+k8mE5cuXO7scaoRGKYePRtohxAcYiYiIiBrn9GC9bds268evvvoq0tPTsWnTJmzatAnp6el47bXX6jyX3Md2LmtOuUdERETUOKcH6/T0dAiCgB49euCFF16AXC63HpPL5ZgzZw569uwJURSRnp7u7HKoCWwfYOQiMURERESNc3qwrqioAABERETUe071scpKjoy2BbYPMLLHmoiIiKhxTg/W/v7+EEUR+/fvh06nszuu0+mwf/9+AICfn5+zy6EmsF8khv/hISIiImqM04N17969AQD5+fmYPHkyfvvtN5w/fx7Z2dn47bffcNVVVyEvLw+CIFjPJfeynRmEI9ZEREREjXP6AjFTpkzBli1bAAC7d+/GpEmTGjyX3I+tIERERETN5/QR63vvvRcxMTHWbdsFYqpFR0fjvvvuc3Y51AR8eJGIiIio+ZwerL28vLB27VoEBwdLFompvThMcHAwfvjhB3h5eTm7HGoC21aQ4goDKg0mN1VDRERE1D64ZEnzvn374ujRo3juuefQvXt3eHh4QKPRoHv37pg1axaOHDmCfv36uaIUaoJgm1YQgKPWRERERI1xeo91tcDAQCxcuBALFy501S2phXw1CqgUMlQZzdZ9OSV6RAV4urEqIiIiorbNJSPW1L4IglBHnzWn3CMiIiJqCIM11cluLmu2ghARERE1yOGtIGPHjgUAjB49Gi+99JJ1uykEQbBOzUfuZTflno7BmoiIiKghDg/W27ZtgyAICAoKkmw3pnrGEGob7BeJYSsIERERUUPYCkJ1YisIERERUfM4ZVaQ2gu/1LVNbZ9tKwin2yMiIqK2RG/SQyEoIJfJ3V2KlcODdWpqKgBYF3up3qb2JdiuFYTBmoiIiFyjylSFC+UXcKHsArLLs5FdZnmz7ivLRqG+EN9P+R7x/vHuLtfK4cG69vLlAKx9015eXggMDHT07chJgr2lwTq/VA+TWYRcxj54IiIiajmzaEZeRR6yy7Jxvuy89f350vO4UG4JzfmV+U16rQvlFy7tYG0rNjYWgiBg6tSp+Pbbb+s8Z/bs2di6dSsEQcCePXucXRI1ge3Di2bREq47+dqvykhERERUrbSqtCYsXwzOtbcvlF+A0Wx0yL2yy7Id8jqO4rKVFxty+vRp/PPPP5wVpA0J9FJDJlgCdbWcEgZrIiKijsxgNiC3PLfe0Jxdmo0SQ4nL6rlQfsFl92qKNhGs9Xr277Y1cpmAIG+1pLfaMuWe1n1FERERkdOIoohifXG9ofl82XnkVeTBLJpdVpOHwgOhXqEI9QxFiFeI9eNQr1CEeIYg3DvcZbU0hVOCdUZGht2+8vLyOvdnZWVZ2z84Yt22dPK1CdZcJIaIiKjdMpgNuFB2AedKz+Fc2bma/ubS89YHBCuMFS6rRyFTINQzFGHeYQjzCkOI58Xg7FUTnH1Vvu0qHzolWFf3VVcTRRG//PIL4uLiGrxOq+VoaFtimXJPZ93mzCBERERtV5WpCufLziOrNAvnSy++LztvDdI55TkuHW0O0AQgzMsSmqvDcu3tQI9AyIRLa0kVp7aC1J6/uqG5rAVBgCAIGDhwoDPLoWayXySGqy8SERG5S4WxQhKYrQG6zPI+tyLXZbXUbtEI8w6zC80hniHQKDrec1ltosdaFEUoFArMmTPH3aVQLXbBmq0gRERETlNaVWpp0agjPJ8rO4eCygKX1CETZAj2CLYLy7W3tWptu2rRcBWnBOtRo0ZZP9nbt2+HIAgIDAxEz549JecJggAPDw9069YNM2bMQJ8+fZxRDrWQ7Qwg54pd13dFRER0qSnWF1vbMs6V1rxVB2hdla7xF3EAL6UXwrzCEO4dXmdoDvYMhlKmdEktlxqnBOtt27ZZP5bJLL0zo0ePrncea2qbogM8JdsZ+eVuqoSIiKjt05v0yCrNwtmSs8gqzUJWSZZlu/QsskqyXDYNna/KF+He4Qj3Cre8t/m4vT0Q2J44vRWkevTadrSa2r6YQGmw1lUaUVReBT9PlZsqIiIich+T2YSc8hycLT1bE55rBWlX9ThXPxRoF54vbnurvF1SB9lzerCuPXpN7Uu4nwfkMgGmWqvEZBSUM1gTEdElSRRFFOmLrEH5bOlZSXA+X3beYSsGNiTII0gami++j/COQKhXKDyVno2/CLmFSx9eLC8vR0pKCoqLi+udJWTUqFGuLIkaoJTLEO6nQWZBTW91en45+kT6ua8oIiKiVig3lONc6TlJaK7+OKskC+VG57Y9ChDQybMTIrwjEOYdZteuEeYdBrVc3fgLUZvkkmCdlZWFRx99FD///DOMxvr/pycIQoPHyfWiAzwlwTqjgH3WRETUdomiiNyKXGSWZEreskosI9CumFnDX+2PSJ9IRHhHWN58IhDpHYlI70iEeoVCKeeDgZcqpwfr0tJSjBgxAhkZGQ3OZU1tU3SAF3Yi37rNBxiJiMjdDGYDzpWeswvPZ0ssvc+VJueuu+Ch8ECEtyUsR/hE2H3spfRy6v2p7XJ6sF6+fDnS09OtT59Wv68O2bbb1LbYzQzCEWsiInKBMkMZzpacRWZJJjJKMiTh+XzZeaeuICgX5AjzCrOONEd4R0hGoAM0AZxVg+rk9GD9888/Wz8eNmwY/vrrLwiCgNjYWCQkJGDz5s0AgJtuugmenmzGb2sYrImIyBlEUUR+Zb41PNu+ObtlI8gjqKZVwzsCUT5R1raNEM8QKGRtYg09amec/l2TnJwMAAgMDMTvv/8OtdrSkD9o0CB8++23+OGHHzB16lRkZGRg+/btzi6Hmsl2yr3zxRWoMpqhUsjcVBEREbUXZtGM7LJsZJRkIEOXYReeK4zOW3hMLVcj0jsSUT5RiPSJRKRPpDU8h3uHw0Ph4bR7U8fl9GBdVFQEQRDQt29fKJU1zfrVrR833HADBgwYgL1792LhwoV4+eWXnV0SNUOUzYi1WQSyiioQF8T+MSIisoTnnPIcZOgykF6SbnmvS7cG6SpzldPu7avyRZRPlOQt0icS0T7RCPYMhkzgIBC5ltODtUKhgMFggI+PDwBApVLBYDCgoKDmVzyhoaEQRRGrVq1isG5jtB5K+HkqUVRusO5Lzy9jsCYi6kCqZ9qoDsy1A7QzHxasnprONjxXB2itWuuU+xK1lNODtb+/PyoqKlBWVgYA8PPzQ05ODvbs2YOsrCyo1Wrs2bMHAJCRkeHscqgFogM8UVRebN3OZJ81EdElp7rn2RqedenIKLG8d2bbhlKmtPY4275F+ERwTmdqV5werDt16oSsrCwUFhYCABISEpCTk4OKigrEx8dDJpOhosLyl9XDg/1ObVF0gCcOn60J1umcco+IqF0SRREFlQXWwFw7QGfoMpy2OIqnwhMxvjF1hudOnp0gl8mdcl8iV3N6sO7duzcOHDiAkydPQhRFjB8/Hn/++ScAoLKy5ldHgiBw1cU2ijODEBG1L+WGcqTp0pBWnGZ5r0uzBulSQ6lT7umh8EC0TzSifaMR4xuDaJ+L732jEagJ5PR01CE4PVgPGjQIa9asgclkQlJSEh5++GG88847yM3NlcxhrdFoMH/+fGeXQy3AYE1E1PaYzCacKz2HVF0q0nXpNSG6OA05FTlOuadGrkGUbxRifGKsATrKJwoxvjEI9ghmeKYOz+nB+pFHHsEjjzwi2bdjxw48/fTT+PPPP2EwGDB06FC8/vrr6NOnj7PLoRaIDrQP1qIo8gcoEZELFFUWIU2XhtTiVOvIc1pxGjJKMmAwGxp/gWZSyVSI9o22jjjXDtKdPDtxpg2iBrhl9vPExESsX7/eHbemFrAdsS6vMiG/rApB3nyghIjIEapMVcjQZVjbNmq3cBTrixt/gWZSypSI8omSBOho32jE+MQgxCuE4ZmohZwerDt37gwA8PLywv79+yVzWVP7EKb1gFIuwGCqWXY+Pb+cwZqIqBmqZ91ILU6tedOlIr04HefKzjl8iW4BAsK9wxHrG4tYbSxifGOsb6GeoXxgkMgJnB6sz58/j6qqKowZM4ahup2SywRE+nsiNa/Mui+zoBwDY/zdWBURUdtU3ft8pvgMUotTcab4jPVjXZXO4ffzVfkiVhtrCdAXQ3SsbyyifaM5VR2Rizk9WMfExODUqVPsx23nogKkwZpT7hFRR1dprES6Lt0uPKcVpzl8tUGFTIEonyhrcI7zjbOOQvur/flvLFEb4fRgPXXqVCxcuBD79++HTqeDr6+vs29JThDDmUGIqIMq1hfXjD4X1QTorNIsiBAbf4FmCPYIto44x/jGIE4bh1jfWIR7h0Mhc8tjUUTUDE7/W/r8889j/fr1OHLkCG699VZ8+OGHiIqKcvZtycHsp9wrq+dMIqL2RxRFZJdl19m+UVBZ4NB7qWQqy6izNs4anKuDtLfK26H3IiLXcnqwvvbaayGXyyGKIjZt2oQuXbogNjYW4eHhkMmkTx0LgoAtW7Y4uyRqgbqm3CMiam9EUcT5svNIKUrB6aLTNW/Fpx2+ZLevyhedtZ3R2a8zOms7W4N0uFc4HxwkukQ5PVhv27YNgiBAEASIogij0YiUlBScPn1ach7nRW7bbEesL+j0qDSYoFHyHwcianuqA3R1cK4O0meKzzh82e4wrzDEaeMk4bmztjMCNAH8d42og3FpwxZ/wLRftsEasMwMEh/i44ZqiIgsqls4rCPQxTWj0I4M0ApBgWjfaGl49uuMON84eCrtfz4SUcfkkmAtio59uINcz0utQJC3CnmlNU+6ZzBYE5GLVAfo6uBcu5XDkQHaS+mFON+LobnW6HOkTySUMk4ZS0QNc3qwTk1NdfYtyEWiAjwlwZpT7hGRMxTri3Gy8CROFZ7CqaJTOFV4CilFKSgzOO6haW+lN7r4dUFXv67o4tcFXbRd0NmvM0I8Q/jbVSJqMZfMY02XhpgATxzIKLJu8wFGImoNvUmP00WnLQH6YohOKUxBTkWOw+5RO0B31na2BulOnp0YoInI4ZwerDMyMgBYljQPDAx09u3Iieyn3GOwJqLGmUUzzpacxanCUzhZdNIapDNKMhy2jHd1gK4efWaAJiJ3cHqwjo2NhSAImDp1Kr799ts6z5k9eza2bt0KQRCwZ88eZ5dELRQd6CXZZrAmIlt5FXmSEehThadwpviMw6ay81Z6o7PfxZHniwGaLRxE1Fa0iWWcTp8+jX/++Yc/FNu4ukaszWYRMhm/bkQdjdFsRFpxGk4UnsCJghM4UXgCyQXJDltMRSlTorO2M+L94y1vfpb3DNBE1Ja1iWCt1+vdXUKrLFq0CLNnzwYA7Nq1C0OHDnVzRc4RY7NITJXRjGxdJcL9PNxUERG5QklViTU8nyw8ieSCZKQUpqDKXNX4xU0Q6R1ZE6D945Hgl4Bo32gu4U1E7Y5TfmpV91XXVl5eXuf+rKwsa/tHexyFOHLkCF5++WV4eXmhrOzSXuY72FsNL5UcZVUm676TF0oYrIkuEaIoIqs0q2YU+mKYzirNcsjr+6v97Uagu/p15TzQRHTJcEqwru6rriaKIn755RfExcU1eJ1Wq3VGOU5jMBgwY8YM9OvXD/Hx8fjiiy/cXZJTyWQCEkJ9JDODJGeXYHRiJ/cVRUQtUmmsxOmi09YWjhMFltHoUkNpq19bI9egi18Xa4Du6t8VCf4JCNQEtssBFCKipnLq79lqLwzT0CIx1UueDxw40JnlONxrr72Go0ePYv/+/Vi8eLG7y3GJbqG+kmB9IrvEfcUQUZOUG8qRXJCMY/nHcLzgOI7lH0NqcSpMoqnxixvRybMTugV0Q6J/IhIDEpHon4gonyjIZXIHVE5E1L60iQY2URShUCgwZ84cd5fSZPv378drr72GV155BT169HB3OS7TLVS60uLx8zo3VUJEdSmpKrGG6OognVacBhGtWwFXISjQ2a8zugV0Q4J/gvW9v8bfQZUTEbV/TgnWo0aNsv66b/v27RAEAYGBgejZs6fkPEEQ4OHhgW7dumHGjBno06ePM8pxOL1ej+nTp6Nfv36YNWtWs66r/aCmTtf+QqltsD6dWwqDyQylXOamiog6rmJ9sXUE+ni+5X1Gif2zLM2lVWuR6J9oDdCJAYnorO0MlVzlgKqJiC5dTgnW27Zts34sk1kC1+jRo+udx7q9eemll3Dq1Cns27cPcnnTf925cOFCzJ8/34mVOV+3UF/JtsEk4kxuGRJtAjcROVZhZaElPBccs45GO+KhwmifaGsLR3WI5pR2REQt4/RWkOrRa9vR6vZq165deOONNzBv3jz06tWrWdfOmTMHTz31lHVbp9MhKirK0SU6ldZTiTCtBueLK637krN1DNZEDlRuKMex/GM4kncESXlJOJJ3BOfKzrXqNeWCHF38uqB7QHf0COyB7oHdkeifyBk5iIgcyOnBuvbodXtnNBqtLSvV81Y3h1qthlqtdkJlrpUY6mMTrEtwnRvrIWrPDGYDThedtgbopLwknC463aqlvhUyBeL94tE9sDt6BFhCdIJ/AjQKjQMrJyIiWy59eNFkMiEvL6/BBWGio6NdWFHzlJaW4tSpUwAAlaruXsPLL78cAPDDDz/g+uuvd1VpLtUt1BfbTuRatzkzCFHTiKKIs6VnJSPRx/OPo9JU2fjF9VDKlEjwT7COQvcI7IF4v3j2QxMRuYFLgvXff/+NuXPnYvv27aiqqn+lLkEQYDQaXVFSi6jVatx33311HtuxYwdOnTqFKVOmIDg4GLGxsa4tzoVsH2BM5swgRHUqqSrB4dzDOJR7yBqki/RFLX49tVyNRP9EdA/sjp6BPdE9sDu6+HWBUqZ0XNFERNRiTg/We/bswZgxY6DX6xucy7o98PDwwMcff1znsbvvvhunTp3CnDlzLtklzat1C5MG63PFlSiuMEDrwX/cqeMSRRGZJZk4mHsQB3MO4mDuQaQUprR4mju5IEeCfwJ6BfVCr6Be6BnYE138unCZbyKiNszpP6HnzZuHyspK6yIwgP1iMYIgtPvQ3ZF0DvKGQibAaK75mp3ILsFlcQFurIrItfQmPY7lH7OE6ItBuqCyoMWvF+UThV5BvdA7qDd6B/VGYkAiPBQeDqyYiIiczenBeteuXdbgrFKpkJiYCK1W26xp6qhtUSlk6NrJG8m1eqtPZOsYrOmSlluei0O5h3Ag5wAO5h7EsfxjMJpb1rrmr/ZH7+De1iDdM7AnF1ohIroEOD1YV/dUa7Va7NmzBwkJCc6+pVusXLkSK1eudHcZLpMY6iMJ1sf5ACNdQqrbOv658A/2XdiHfRf2tXjOaKVMiZ6BPdEnuA96B1tGo8O9wjlPNBHRJcjpwTohIQFJSUkYMmTIJRuqO6Juob74ETXz6nJmEGrPRFFEanEq/rnwjyVMZ+9DTkVOi14rQBOA/p36o19wP/Tr1A/dA7tDLW//02wSEVHjnB6sp02bhlmzZuHUqVMQRZGjNJcI25lBTmSX8OtL7YZZNONU4SnJiHRL+qMFCIj3j7eG6H6d+iHSO5J/D4iIOiinB+vHHnsMX3/9NQ4cOIAHH3wQS5YsgVardfZtyclsZwYp1RtxtrACUQFcxY3aHlEUka5Lx+7zu7H7/G78nf03dFXNnybSS+mFvsF90S+4H/p26os+QX3grfJ2QsVERNQeOT1YT5o0CQaDAaIo4pNPPsGXX36Jrl27IjAw0O5cQRCwZcsWZ5dEDhDqq4GvRgFdZc3DW8nZJQzW1GbkVeRZgvQ5S5i+UH6h2a8RqAnEwJCBGBQ6CAM6DUBXv66Qy/jgNRER1c0lS5pXT7UniiIqKiqQlJRk96tSthG0L4IgoFuYL/am1vz6/ES2DuN7hLixKurISqtK8c+Ff7Dn/B7sPr8bKUUpzX6NTp6dMChkEAaFDsKgkEGI9Y3lzyUiImoyl640wH+gLi3dQn0kwZozg5ArmUUzjucfxx9Zf2Bn1k4k5SXBJJqa9RoR3hGWEemLYZr90URE1BouCdZc/OXS1C3UV7LNmUHI2Yr1xdh1bhf+yPoDf2b92ewHDoM9gjE0bCiGhA3BZaGXIcw7zEmVEhFRR+T0YJ2amursW5CbJNrMDJKaV4ZKgwkaJXtQyTFEUcSJwhP446wlSB/MPQizaG7y9V5KLwwOHYyhYUMxNGwoOms7c0SaiIicxunBOiYmxtm3IDexDdYms4iTF0rQJ9LPPQXRJaHSWIld53Zh29lt+PPsn82aT1ohU6BfcD/rqHSvoF5QyFza8UZERB2YW/7FqX6I0dOTM0i0Z95qBeKCvJCaV2bddzCziMGamq2osgjbz27H1oyt2HV+FyqMFU2+NsY3BiMjRmJY+DAMDBkITyV/rhARkXu4LFhnZ2dj4cKFWLduHTIzMyGKIoxGIz744ANcuHABCoUCzz//vKvKIQfpH+UnCdb70wsx/fJY9xVE7UZWaRZ+z/gdWzO3Yv+F/U1+8FAtV2Nw6GCMiBiBkREjEe0b7eRKiYiImsYlwXr37t2YMmUK8vPzrQ8yVvc5ZmVlYcGCBRAEAcOGDcPo0aNdURI5SP8Yf3x/IMu6vT+jyH3FUJuXocvAr2m/YlP6JiQXJDf5ukjvSIyMHIkRESMwOHQwPBQeTqySiIioZZwerAsKCnDjjTciLy9PMp91tZtuugkLFiwAAGzYsIHBup0ZEO0n2c4oKEdeqR5B3mr3FERtztmSs9iUvgkbUzfieMHxJl0jE2QY0GkAxkSNwcjIkZxPmoiI2gWnB+t33nkH2dnZdoG6Wp8+fRAUFIT8/Hzs3r3b2eWQgyWG+MBTJUd5Vc2v8Q9kFHGhmA4uuywbv6b9il/TfkVSXlKTrtHINRgWPgxjo8diVOQo+Gv8nVwlERGRYzk9WK9bt8768VdffYVvv/0Wa9eulZzTo0cP7NixAykpzV8pjdxLIZehT6QWu8/UzCe8P6OQwboDKqkqwaa0Tfjp9E/Yn7O/Sdf4qf1wReQVGBs9FpeHX84WDyIiatecHqxPnToFQRAwaNAg3Hbbbfjuu+/szvHz8wNgaRuh9mdAtL80WKcXurEaciWj2Yi/zv2FdafXYWvGVlSZqxq9JkATgHHR4zAxdiIGhAzgdHhERHTJcPq/aBUVlmmzQkLqH8EsKioCAMhkMmeXQ07QP1r6K/vDZ4thNJmhkPPreSmqXrTlp9M/YcOZDcivzG/0Gq1aaw3Tg0MHM0wTEdElyen/ugUEBODChQs4ceJEncdLSkqwf7/l18aBgYHOLoecoL/NA4wVBhOSs0vQK0LrnoLIKUqqSrDhzAasObWmSTN6+Kh8rGH6srDLoJQpXVAlERGR+zg9WPft2xebNm1CSkoKlixZAqPRaD124cIFPPzwwygtLYUgCOjXr5+zyyEnCPJWIybQE+n55dZ9BzIKGawvAaIo4mj+Uaw+uRq/pP7S6MItCpkCV0RegWu7XItREaOglDNMExFRx+H0YH399ddj06ZNAIDZs2dLjkVEREhmCrnhhhucXQ45yYBof5tgXYS7LndjQdQqpVWl2JC6AatPrm7S6HSfoD64tsu1mBQ7CX4aP+cXSERE1AY5PVjffffdWLZsmXXGD1EUrfPRms1m69zW8fHxuPPOO51dDjlJ/2g//CBZKIYPMLZHpwpP4cvjX2JD6oZGR6fDvMJwTedrcG2XaxGnjXNRhURERG2X04O1RqPB2rVrMWHCBGRlZUkWeaie2zo0NBTff/891GouKtJeDbB5gDEtvxz5pXoEcqGYNs8smvFn1p/437H/Yff5hueSlwtyjI4ajZsTbsbl4ZdDJvABVSIiomoueTS/e/fuSEpKwptvvom1a9ciNTUVABAXF4cpU6bg6aefRkBAgCtKISfpFuoDjVKGSoPZuu9gZhGu7M75rNuqckM51p1ehy+Of4E0XVqD54Z7hWNqwlTc0PUGBHsGu6ZAIiKidkYQ61oOkVxCp9NBq9WiuLgYvr6+7i6n1W75YBf2ptbMZ/3Q6C54blI3N1ZEdckpz8GXx7/EmpNroKvS1XueXJBjVOQo3JxwM4aFD4NcJndhlURERG1Dc/IaJ5MlhxkU4y8J1n+dbnx+Y3KdrNIsfJr0KX5I+QEGs6He8/zV/rgl8RbcnHAzQrz4GwciIqKmcnqw/vbbb/H+++8DAB566CHcfPPNkuNr1qzBe++9BwCYOXMmbrnlFmeXRE4yvGsQ3tt22rqddLYIxeUGaD055Zo7pevS8XHSx1h/ej2MorHe87r6dcVdPe7CVXFXQaPQuLBCIiKiS4PTg/XXX3+Nbdu2QaVS1bmc+fjx4zF9+nTo9Xr4+fkxWLdjA2P8oVbIoDda+qzNIrDrTD4m9Qp1c2UdU0phCj5M+hC/pv0Ks2iu97xRkaMwrfs0DA0bKnm4mIiIiJrH6cH6wIEDEAQBAwcOhL+/v91xrVaLAQMG4K+//sLBgwedXQ45kUYpx2VxAfjjVJ51386UPAZrFzuWfwwfHf4ImzM213uOUqbE9V2vx1097uJUeURERA7i9GCdm5sLAAgKCqr3nOoZQarPpfZreNcgu2BNrnGm+Aze3v92g4FaI9fgpoSbcHfPu9k/TURE5GBOD9Zms+VX0GlpafWek56eDgAwmUzOLoecbERX6X+gzuSVIauoAhF+Hm6q6NJ3oewClh9ajrUpa2ES6/475KnwxO3dbsddPe5CoEegiyskIiLqGJwerIODg5GZmYmjR49i7969uOyyyyTHd+/ejSNHjkAQBAQHc37c9q5HmC/8PJUoKq+ZdWJnSh5uGRTlxqouTcX6Ynx65FN8efxL6E36Os/xUflgWvdpuLP7ndCqtS6ukIiIqGNxerAeMGAAMjMzYTabcd1112HJkiUYNWoUAGD79u2YNWuWdZnz/v37O7sccjKZTMDwLkH4Oem8dR+DtWNVmarw1fGv8FHSR/XOQ61Va3F3z7txW+Jt8FZ5u7hCIiKijsnpwfqmm27Cjz/+CEEQcOHCBcyYMUNyvPb6NDfddJOzyyEXGN7VPlhX/+eJWk4URWw/ux1L/l6CjJKMOs/RyDW4q8dduKfXPfBR+bi4QiIioo7N6cH61ltvxeLFi63tHrYLPVaHrZ49e+K2225zdjnkArZ91nmlVThxoQTdQtv/6pLucqboDBb9vQh/nfurzuNyQY6p8VMxs+9MLjlORETkJk4P1gqFAt999x3GjRuHjIwMu1FLURQRFRWF77//HgoFF4K8FEQHeiIqwAOZBRXWfX+eymOwboGSqhK8d/A9rEpeVe+DiRNjJ+LR/o8ixjfGxdURERFRbTJX3KRr165ISkrC3Llz0bNnT3h4eMDDwwM9e/bEiy++iMOHD6Nr166uKIVcxHbU+k9Ou9dsm9M347q11+GL41/UGar7BvfFqqtX4Y0r3mCoJiIiagME0bY3g1xGp9NBq9WiuLgYvr6X1mju+sPn8MhXB6zbGqUMB+ZOgIdK7saq2ocLZRfw7z3/xtbMrXUe7+TRCU8OehJXx13NvnUiIiIna05eY+8FOcWIrkGQywSYzJb/t1UazPgzJQ/je3BRkvqYRTO+PfEt/rP/PygzlNkdV8lUmNFzBu7vfT88lZ5uqJCIiIga4tJgbTKZkJeXB72+7jl3ASA6OtqFFZGz+HmqMDjWH7vPFFj3/XYsm8G6HqcKT2H+rvk4lHuozuPDI4bjhSEvIMqH0xYSERG1VS4J1n///Tfmzp2L7du3o6qqqt7zBEGA0Wh0RUnkAuO6h0iC9ZbjOTCZRchl7m1fMFdVQSwvh7myEmJlJcx6PcSKCpgr9RD1lTBXVEI0VAEQAEGwvJPJLB9DAGSW/YIgAHI5ZBoNBLUGglpl/VimUUPQaCBTqyEolfXWYjAZ8GHSh/g46WMYzfbf+wGaADw3+DlMjpvMtg8iIqI2zunBes+ePRgzZgz0er3dVHt0aRvfIwQLfj5u3c4vq8LBzCIMjPF3yv3MVVUwZGXBcDYLhrOZMOTkwJRfAGNBfs37vHyYy+zbLJxKLofM0xMyb2/Ivb0h8/aGzMcb5Srg75KjKEUhpqgFVKiAcnXN22XxYzBjyKPw78RRaiIiovbA6cF63rx5qKyshFA9wgfUOZc1Q/elJybQCwkh3jh5odS677djF1odrEVRhDE7G5UnTkB/4iT0J06g8uQJVKWmAaa6p6RzK5MJ5pISmEtKYDsmPcj6UV3f/1uQgy3IASCoVJD7+TX9zd8Pcq3WMtJORERELuH0YL1r1y5rcFapVEhMTIRWq4VcztkhOoJx3UNsgnU2Zk/u1qzXEEURVampKN+7F+V796Js798w5XWs6fvEqioYc3JgzMlp+kVyOeQB/lAEBEIRGAB5fe8DA6EICIDMkw9EEhERtYbTg3V1T7VWq8WePXuQkJDg7FtSGzKuRwje23baun06twzJ2bpGF4sRzWZUHDgA3a+/omTTbzBmZzuvSEGw9ENrNDV90Wq15ZjZDEC0/EbFLAKi5U0UzYAIiEYDRH2VtVcbbekZAZMJptw8mHLzUP/jwjUEDw8oAmqCtjwwAIqAQMv7wEDIAwKgCAqGIjgIcj8/joYTERHZcHqwTkhIQFJSEoYMGcJQ3QH1i/RDmFaD88WV1n3rDp2rN1hXJiejaPUalGzaBGNubovuKddqoYyMhDI8DPKgoJpwGBAIRZAlIMq8vCHzsARpQal02IOBotEIUa+3PBB5MWzry0vw1b5P8NepzfDUA556wKMK8NSL8KwU0MsjDgmqKKCsHKbSEphLSmEuKYFJp7sY7F1DrKiw9KhnZTV+slwORWAgFEFBkAcHQREUZAndQUFQWLeDIA8KhszLkw9eEhFRh+D0YD1t2jTMmjULp06dgiiK/Ae2g5HJBFzTJwwf/ZFq3bfu0Hk8MyGxpufeZELptm0o+OxzlO/d2/TX9vaGOjERmsQEqBMSoU5MgLprV8h9fBz+52gqQaGAoFBA5uUFAEgtTsWzSa/ghHgC6Cod4Q33CseCEQswOHRwna8lms2WgF1UZPdmvPjeXFws2TYVFUMsL3f6nxMmU5NbUwQPD2vQrg7e8trb1hAeBJlK5fzaiYiInMTpKy9WVVVh2LBhOHDgAO677z4sWbIEWq3WmbdsNy7llRdrO3y2CFPe2SnZ9+O/hqNPqBeK1qxB/qcrYMjMbPR1FCEh8BxyGbwuuwyegwZBGRPTpv+jtiltE+bunItyo33QnRQ7CS9d/hJ8VI7/T4BZr7eE7Px8GC++2c2Okl8AY0EBTHl5EA0Gh9fQUjKt1i5w1wTxYOtouNzfn60oRETkEs3Ja04P1mPHjkV+fj6SkpIgCAI0Gg26du2KwMBA+2IEAVu2bHFmOW1KRwnWoihizBvbkJZfXr0DLwYWYsy2r1F15kyD16rj4+EzaSJ8J0yAqmvXNh2kqxnNRvxn33/w2bHP7I5p5BrMvmw2boy/sU38WURRhLms7GIIL4CpoI73efnWqQpNRUXuLtlCLrf0gdu2odRqRZEHWj6WeXu3ic81ERG1T20qWMtksjqn2bP9h666TcTUFqdLc5KOEqwB4M1NJ/D21hRElVzAg0k/YmDOyXrPlQcFwf+WW+B7zdVQd+7swipbL68iD89sfwb7LuyzO9bVryuWjFqCrv5d3VCZY4gGA4wFBTDm5sGYlwtTXh6MeXmW7fx8y75cyz6XzxdeD0Gtrmk3kQTxQEkvuCIoEDKNxt3lEhFRG9OcvObSJc05atRxXdM7FBkfrcC9R3+Gqo4VBgFA3aM7AmfMgM/kye2y1/Z4/nE89vtjyC6zn8Hkui7X4YWhL8BD4eGGyhxHUCqhDAmBMqTxpenN5eWWsH0xhBvz8ixB/GLwrv0GJ7ajiHp9kx/KlPn4WMJ2YODFEF7PSHhgAASFS398EhFRO+CSfxm4+EvHZriQA48X5mBm0l91Htf07IlOzz4LzyGXtdv/fP2W/hte+PMFVBgrJPuVMiXmDJmDm+Jvard/tpaSeXpC5ekJVVTDK0eKoghzcXFN0LYG71y7IG4qLLRMeegk5pISVJWUoCo1teETBQFyf/+6+8AlI+EXpybsYF97IqKOyunBOrWxf6DoklZx6BAyH3kEplz7BV0KPLRIeHEWgm64vt0+iCaKIj44/AHePfiu3bEwrzAsHb0UvYJ6uaGy9kMQBOuKkequDbfJWFpRCmu1oeRfDOO5MObnWdtQjHl5MJeWNvharSKKMBUUwFRQAP3J+tuaAABKpXVqwpqR8LqnJ6yeTYaIiNonpwfrmJgYZ9+C2qji9T/j/PPPQ7y4SFBtP3QZic+6T8K/Yy/D1HYaqiuNlZi7cy42pm20Ozag0wAsG7MMAZoAN1R26bK0onSCMqRTo+eaKyoutqLkWh7OtGtDqekHr+t71GEMBhizs5u0yBGnJiQiat9c3iS4f/9+6yh2XFwcBgwY4OoSyAWKvvse51980e7X9jpvfyzsczMOdrIsFvTV3gxMHRjpjhJbpaiyCI9sfQSHcg/ZHbsx/ka8OORFKOVKN1RG1WQeHlBFRkIV2fD3lyiKMJeWNvhApjHv4mh4fr5TF+0RKypgyMxs0vSTMq0WioCLq2IGBl58HwBFYJD1vWVBpEAu0kNE5CIuC9arVq3Cc889hyybB4jCw8Px+uuv484773RVKeRk9YVqz8uHIuOB53Dwp5olzvelF+JEdgkSQ923qEtznS89j5mbZ+JMsXSqQJkgwzODnsG07tMYYtoRQRAg9/GB3McH6s5xDZ4rmkyWBXmqe77t+sKrw7nzpyY0Fxejqri48X5wAIJGY5meMCjo4vvAWiuR1qxIqggKglyrhSCXO7V2IqJLlUuC9dKlS/Hss8/W+RBjVlYWpk+fjvPnz+OZZ55xRTnkRPWFav87bkfInDkIl8kR/PtZ5JborcdW7c3AvCk9XV1qi5wqPIWZm2cip1y64qC30htLrliCEREj3FQZuYJQvZR7YCCQmNjguWJVlWRqQsmsKNWtKRfbUcxOXi1TrKyE4dw5GM6da/xkmcwSsgMCaoJ3HaPiisAAyAMDIVOrnVo7EVF74vR5rJOTk9GnTx8YjUa7+axrbysUChw6dAjdu3d3ZjltyqU2j3XJli04+8ijdqE68IEHEPzUk9av95Jfk/Hu7zWj1j4aBfY+Pw4eqrY9Srbvwj48uvVRlFSVSPaHeIbg/XHvt+v5qcm9zGVltcJ2HuqbmrCtrZQJADJv74bbUS4GcEVgIGQ+PvxtDhG1O21qHuvly5dbQ7UoiggPD0ePHj0AAMeOHcO5iyMoJpMJy5cvx1tvveXsksgJ9GdScW7Wc42GagC4bXA03tt22npqSaUR6w+fw82DGp6WzZ22ZGzBrO2zUGWWPuTWWdsZH4z/AKFeoW6qjC4FMi8vqLy8oIqObvA8URRh1ums7SeSlTHzpatkGvPzIVZWOr12c2kpqkpLgfT0xk9WKqHw9784Iu4PuX+A5ePAgIsf+1taVfwt+2S+vgziRNSuOD1Yb9u2zfrxq6++itmzZ0N+sX/PZDJh8eLFeOGFF+zOpfbDVFqGs48+arfSXl2hGgCiAjwxMj4YO07mWvd9tTejzQbrb098i9f2vAazKH1orV9wP7xz5TvQqrVuqow6GkEQINdqIddqoe7SpdHzzWVllnaUvDyYCgpgzMu3W6beWFAAU14eTMXFzv8DGAww5uTAmJMDfeNnAwoF5P5+UFQH8ADLe9sALg8IgNzf39If3k5nGSKiS4PTW0H8/PxQUlKCHj16ICkpqc5zevfujaNHj8LHxwfFrvjh3kZcCq0goigi64knUfLrr5L9vlOuRfiiRfWONm08ko2ZX0iX/f7l8ZHoHta2Pg8rjqzA0n1L7faPjhqNJaOWQKPgEth0aaieI7w6eNcVwI35eTDlF8BYUODU1TJbTC6H3M/PMhoeEGgJ4P61w3hgTSgPCOCDmkTUJG2qFaSiwrISXURERL3nRERE4OjRo6h0wa8tybEKVqy0C9Xqbt0QNn9+g7/CvbJ7J3TyUSOn1kOMK3emYdFNfZxWa3N9cOgDvHPwHbv9U+On4sWhL0Ih45LWdOlozhzh1tUymzga7uyHM61MJpjyLW0xQErj58tklsWJagdwfz/Lqpr+/pZRcD9/yzn+/lD4+0Hw5NSFRFQ/pycDf39/5OTkYP/+/dDpdHZJX6fTYf/+/QAso9vUflSdPYvc//5Xsk+m1SLy7bcg8/Bo8FqlXIZbB0fh7a01//h9f+AsHh8Xj3C/hq91NlEU8e7Bd/HB4Q/sjs3sOxMP932Y/7BShyZZLbNz50bPtyzWUwBToWW1SsvIeIElgFd/XFi9rwCiq4K42WxdQbMKpxs/H4CgUlkCt7//xTYVafiW7L/4JtPwN1tEHYXTg3Xv3r2xZcsW5OfnY/LkyZg3bx569eoFQRCQlJSE+fPnIy8vD4IgoHfv3s4uhxzowuuvQ9TX6pQUBES8sQSqqKb1St91eQw+2HEGVUZL77LBJOLDHWfcOvWeKIpYtn8ZVhxZYXfs2UHPYnrP6W6oiqh9syzWEwFE1v+by9rMlZU1AbywAMb8iwG80BK8TfkFMBYWWEO57fMdziRWVcF44QKMFy40+RrBw8MyEu7nB4WfvySAS0bH/WsCOlfWJGqfnB6sp0yZgi1btgAAdu/ejUmTJjV4LrUPpX/8gdLNWyT7/O+4A94jRzb5NTr5aHDb4Ch8vqtmNoGv/87AI2O7Isjb9XPjNhSqXxjyAm7rdpvLayLqiGQaDWTh4VCGhzfpfLNebx3tbiiAGwst+8ylpU7+E0iJFRUwVlTAeO580x7aBCDz9KwVwGuPjvvZtKj4WR5o9fPjyDhRG+D0hxfLysrQq1cvZGRkAIDdIjHV0/DFxMTgyJEj8PLycmY5bUp7fXjRXFWF1GunoKrW9FrygAB02fgL5M38c5wtLMfoJdtgNNd8Xzw8ugtmTermsHqb6r2D72H5oeWSfQIEvHz5y5iaMNXl9RCRc5irqmCq1XpS3Q5iLCyEqbDIcqzI8t5YZNkHo9HdZTdK0GisIdvuvZ/ttp91hhmBo+NEDWpTDy96eXlh7dq1mDhxInJycux6U0VRRHBwMH744YcOFarbs4IVKyWhGgA6Pf1Us0M1AET6e+L6/hFYs++sdd//dqXjwSu6QOuhbHWtTfVx0sd2oVomyPDq8FcxpQt/k0J0KZGpVJCFhEAZEtKk80VRhLm01BK4L4ZuSQgvLISpqNCyr6jIsr+oCDCZnPsHsa2zshLGyspmtakAF0fH/fwg89NC4ecHmV0orx3UL46S+/hAUPABbiJbTh+xrpafn4833ngD69atQ1paGkRRRFxcHK655ho8/fTTCA4OdkUZbUp7HLE2nD+P01ddDfHibC8AoOnbB7GrVrV4/tjTuaUYt3S7ZG2ZZyYk4JGx8a0tt0n+d+x/WPz3Ysk+AQIWjFjAUE1ELSKazTCXlFhGvatDeFGRNIRLgvnFMO6af5IdQubrW8fIeAOj5FqtZfVNzjVO7Uxz8prLgjXZa4/BOuupp6Db8EvNDkFA7LffwqN3r1a97r++3I+fk85bt/09ldg5eyw8Vc4dEfn2xLd4dferdvtfuvwl3Jxws1PvTURUm2gywaTTXQzhhQ2PkFeHcZ2uXYVxyGSQ+fhYwravr+XNT2sJ6b7ai+0pvhdDux/k2ovnaLWc6pDcpk21gtCloyojA7pfNkr2+d18c6tDNQA8PKaLJFgXlhvw1Z4M3D+y8am8WurXtF+xYPcCu/2zL5vNUE1ELifI5VBcnCUEiGvSNaLZDLNOB1NxsXXUu+bjYul2rXNc/QCnldkMc3ExzMXFaPYSQwpFTRjXaiHTNh7GZb4Xj/HBTnIRhwfrzMxM9O/fH2azGdHR0fj777+hVNbdK1tVVYXBgwcjMzMTcrkcBw4cQGRkpKNLcqisrCysXr0aGzZsQHJyMrKzsxEQEIDhw4dj1qxZGDJkiLtLdJrCVV9LRkZkPj4IfvIJh7x2z3AtxnbrhK3JOdZ9y7edxi2Do+CrcXyv9Z7zezDnjzkQIR3peXLgk7iz+50Ovx8RkTMI1Yvc+PkBMTFNvk40GCyj4/WF8LoCeXGx6+YYr4vRaH3QtLkEtRpyX19LGNf61T1a7mcZRZcEdB8fPtxJzeLwYL1q1SoUFBRAEAQ8/vjj9YZqAFCpVHjsscfwwAMPQBAEfPXVV5g1a5ajS3Kot99+G4sWLUKXLl0wYcIEBAcH49SpU1i7di3Wrl2Lr776Crfeequ7y3Q4c0UFir77TrLP78YbL46sOMa/xnSRBOv8siq8szUFz1/V3WH3AIBj+cfw2NbHYDBLx0tm9p2Je3vd69B7ERG1RYJSCUVgIBSBgc26zqzXW0J4sTR8m21DeKE0pIuGZo9PO5So18OYmwvk5jb7WkGjgdzHxxLKfXwh8/WB3McXcl8fyKrf+/rWvc/bG0IDOYguPQ7vsZ40aRI2bdoEDw8PFBQUQK1ueD7iyspKBAQEQK/XY/z48di4cWOD57vb999/j8DAQFxxxRWS/X/88QeuvPJKeHt74/z5843+uYH21WNduHo1sue+JNnX5deNUDVjhKQp7v/sH2w+XvNEu1IuYNOTVyAuyDEzxmToMnDXL3ehoFI64nFr4q14YcgL7N8jInIwURQhVlRYe8JNxTqYdJZ2EMvHOpiKiy62tFwcRdfpLMfbWw95HQRPT8vouE91APexhnOZr49ltNzXx9J77utrfS/3sewT5HJ3/xE6PLf2WB89ehSCIKBv375NCpcajQb9+vXD7t27cezYMUeX43A33nhjnftHjhyJMWPGYNOmTUhKSsKgQYNcXJnziKKIwi+/kuzzGjXS4aEaAF64uju2n8yBwWT5QWowiXjt5+P4eEbrP5/5Ffn4v9/+zy5Uj48ZjzmXzWGoJiJyAkEQIHh6Qubp2eRFf6qJZrNlqkOdDqaiYph1xZbg3VggLy526YqcDRHLy2EsL4cxO7tF18u8vWtGySXh/GI7S+1Rctv33t6chcXFHB6s8/LyAACdOnVq8jXV5+a24Fc0bUl124viEpvbs2L/fuiTkyX7Au50Th9yXJAX7h4Wi4/+SLXu23z8Av48lYcR8UEtft1KYyUe2/oYskqzJPuHhA7B6yNfh1zGEQEiorZGkMms/dBo5jNYotEIU0mJpVVFp7OGc5Ou2LJdZAni1oBuDevFECsrnfQnaj5zaSnMpaUw4nzjJ9dB5uVlGQX38YbM2wcyH2/IvX3s9/n4WD729rKOlsu9vS3h/BLLNc7ktM9UaTOeOK4+tz2PGGZkZGDz5s0ICwtD79696zxHr9dDr69Z0Fan07mqvFYp/PJLybYyKgpezVi6vLkevTIe3+/PQn5ZlXXfq+uP4efHRkAhb/7/vM2iGc//+TwO5x2W7O8e0B3/GfMfqOR8MIWI6FIjKBS1ZllpHnNVlWXku6TEEsJLSiyj4SUlMOlKYC65OEpeooNZV2J3HtzcU16buawM5rIyGFs2YA7gYjuLt3dN2PZpPKDLfS6e5+3doRYUcvifMiAgAOfPn8c///yDqqoqqBp5mlav1+Pvv/8GAPg78EE4VzIYDLjrrrug1+uxaNEiyOvph1q4cCHmz5/v4upax3AhB7pNv0n2+d9xh1N/teSrUeLpCYl4/ock674TF0qw6u9M3DW0+e0n/9n3H/yWLv0zhHuF471x78Fb5d3qeomI6NIiU6kgCw6GogWL14miCLGysiaA13pvDeI6XYPHXL1qZ2Oq21mQk9P4yfUQPDxqQrm3dx0BvdaouW1AvzhyLmtCi7G7OTxYd+vWDefPn0dJSQkWLVqEuXPnNnj+okWLUFJSAkEQkJiY6OhynM5sNuPuu+/Gjh078MADD+Cuu+6q99w5c+bgqaeesm7rdDpERUW5oswWK/r2W8BotG4LGg38brzB6fe9dXAU/rc7HcfP14zqL910Alf3DkOAV9NHmL898S1WHF0h2eej9MF7495DkEfLW0uIiIjqIggCBA8PyDw8gJCmt8VWE0URYnm5zSh500bLq9/DbHbCn6x1xIoKGCsqWjQzSzVBqawJ2RfDuf/0u+A7frwDK20dhwfrK664Ar///jtEUcT8+fMhl8vx7LPP2k27ZzAYsHjxYrzyyivWfaNHj3Z0OU5lNptx77334quvvsK0adPw/vvvN3i+Wq1u0gOdbYVYVYXCb7+R7NNeey3kWq3T7y2XCXjpmh64/aPd1n2F5QbM/fEI3r1jQJNeY8/5Pfj3nn9L9ikEBZaOWYoufl0cWi8REZEjCIIAwcsLMi8vKENDm329KIqW9g+dDqaSUphLL4bv6o9LS2s+Lim1BPbSi8dLLh4vLW2b4dxgsK48Ws332mvcWJE9h0+3l5mZiS5dusBkMkEURQiCgODgYFx55ZWIjY2FKIpIS0vDli1bkJeXh+rbKxQKnD59us2P4FYzm82455578Pnnn+P222/H//73v3pbQOrT1qfb0/26CVmPPy7ZF7f2B2i6dXNZDQ99sQ+/HJE2hr19e39c27fhJ8szdZm4fcPtKNYXS/a/MuwV3BDv/BF3IiKi9soSzsthLi2xhO26AnpdobysDOaSEmtAd0VLS8SypfCdPNmp93DrdHtRUVF46KGH8Pbbb0MQBIiiiJycHHz99deS86oDdfUDiw8//HC7DNW33npri0J1e1D6+++SbY+BA10aqgFg/pSe+Ot0Poorah4EmfvjEQzpHIBOPnUvUVtaVYpHtj5iF6rv63UfQzUREVEjBEGA3NsLcm8voAWj5kCt+csvBvG6ArrdSHkdo+e121HrIvNuW89KOeURzSVLlmD//v3YuXOnNTjbDozX3j9y5EgsWbLEGaU4XHX7x+eff46bb74ZX3zxxSUZqkWzGaV//inZ5ztxosvr6OSrwSvX9cTjXx+07isqN+D574/go+kD7WaSMYtmzPljDs4Un5HsHxM1Bo8NeMwVJRMREXV4tecvb0mvOXAxnFdVWUbBS0thKi27OP1gdctKGdRduzq48tZxSrBWqVTYvHkz/vWvf2HlypUwX+zTqR2mRVGETCbDPffcg3feeafBpc/bkldeeQWfffYZvL29kZCQgAULFtidc/3116Nfv36uL86BKo8fh+ninOTVvEc5b4q9hkzpG45fj2ZjQ1JNS8jm4xfw3f4s3DRQOq/pR4c/wraz2yT7uvp1xcKRCyETOEk+ERFReyEIAgS12jIbSFD7mHDAaZMKqtVqfPzxx3jyySfxzTffYOfOnbhwwbJUdUhICIYPH45bb70VPXv2dFYJTpGWlgbAMvf2a6+9Vuc5sbGx7T5Yl/3xh2RbGR0NVWysW2oRBAGvXtcLe84USOa2nr/uKIbEBSAqwBMAsDNrJ949+K7kWj+1H94e+za8lI5ZEp2IiIioPg5/eJGari0/vJh2x52o2L/fuu1/550InfuiGysCfj2ajQf/t0+yr2e4L757aBjy9dm4df2tkr5qmSDDB+M/wNCwoa4ulYiIiC4Rzclr/N042TEVF6Pi4EHJPne1gdQ2sWcobuwfIdl39JwOz/+wH0/+/qTdw4qP9X+MoZqIiIhchsGa7JT99Zdk/kpBpYLnZZe5saIa86/ric7Btds6RGw4/y6OFxyXnHdl9JW4t9e9ri2OiIiIOjQGa7JTukPaX+152WWWFaTaAB+NEh9MGwhPlWUmFqXfXij9pO0hsb6xWDB8gd2MIURERETOxGBNEpZp9qTBui20gdQWH+KDxTf1gUyTBXXIT5JjGrkH/jPmP/BWta15LYmIiOjSx2BNEvrkZJhypdPseY1sW8EaAMZ290NolzUQZNJVnRQFt8JPEVnPVURERETOw2BNErZtIMqoKLdNs9eQhXsXosR8XrKvKn8Ezp/rhvs++wcVVc5fRpWIiIioNgZrkii1mb/ae+TINtervOHMBqxNWSvZZ6qIhD5nEgDgYGYRHl21H0aTuY6riYiIiJyDwZqsTDqd3TR7Xm2svzqzJBOv7n5Vss9T4QVtyT2ovd7R5uM5ePLbQzCZOU07ERERuQaDNVmV/fUXYKppoRBUKngNGeLGiqQMZgNm75iNUkOpZP/cy1/E59Mnw1cjXUh03aFzeGb1IY5cExERkUswWJOVbRuI5+DBbWaaPQBYfnA5Ducdluyb0mUKrul8DeJDfPDxjMFQK6Tf0j8cyMIjXx2A3sieayIiInIuBmuyqjycJNn2GjHCTZXY23N+Dz5O+liyL9onGs8Ped66fVlcAD6aPggqm3C98Wg2Hvh8Hx9oJCIiIqdisCYAgLmyEvozZyT7PPr2cVM1UoWVhZjzxxyIqOmXVsgUWDxqMbyUXpJzRyUE44O7BtqNXO84mYvbP9qN3BK9S2omIiKijofBmgAA+lMpkv5qCALUCYnuK+giURQxd+dc5FbkSvY/3v9x9AzqWec1YxI7YeU9l8Hr4uqM1Q5mFuH6d3fiRHaJ0+olIiKijovBmgAAlcePSbZV0dGQe3vVc7brrEpehe1nt0v2DQsfhuk9pzd43eVdAvHF/UOg9VBK9mcVVWDq8r/w+4kch9dKREREHRuDNQGwrLhYm7pHdzdVUiO1OBVL9y2V7AvQBOC1Ea9BJjT+rds/2h+rZ16OCD/pA5ileiPuW/k33v09BWZOx0dEREQOwmBNAIDKY8cl25pu7g3WRrMRL/75IvQmaU/0ayNeQ5BHUJNfJyHEB2v/NRz9o/0k+80isOTXE5ixYi/7romIiMghGKwJosmEypMnJfs0bh6xXnFkhd3Uerd3ux0jIpo/U0mwjxqrHhiKKX3D7Y79cSoPV731B/5KyWtxrURERERA7aXqqMOqSs+AWF4u2afp1q1lL1aUCZzZBuQcAy4cBfJOAvoSwFQFKDwAryDAKxjwjwXC+wHh/YHQ3oCqpp/7RMEJvHfoPcnLxvjG4MmBT7asJgAapRz/va0fEkK8sfS3k6jdAZJbosedn+zB/43qjCfHJUCjlNf/QkRERET1YLAm6JOlbSDy4CAogoOb/gJmE3DqN+CfTyzvUU/fsqkK0BcDBaeBzN3A4a8t+wUZED4ASJyMqvjxmLN3Poxmo/UymSDDayNeg4eidYvVCIKAR8bGY1BsAB7/+gAu6GpaQEQR+GD7GWw6egELb+yNoZ0DW3UvIiIi6njYCkKoPN6K/upTvwFv9QNW3Qqc2oR6Q3VDRDOQ9Q+w9VW89+21OFV0SnL43l73om9w3+a/bj2Gdg7EhsdG4ooE+/88pOaV4bYPd2PO90nQVRocdk8iIiK69DFYEyqPS2cE0XRvQrA2GYHN84AvbwKKMhxSR5JKhRVaX8m+BLMMD+lVQHmBQ+5RLdBbjRV3D8bsyd2glAt2x1ftzcCVb27HN39nwMSZQ4iIiKgJBFEUmRrcRKfTQavVori4GL6+vo1f4ASiKOLUiJEw5edb90UsWwrfyZPrv0h3DlhzH5DxV93H1b5A5yuATj2BkB6ATxggVwJV5UBZLlByHsg+Apw7AOQeB0QzqgDcEhGK0yqV9WUUooivs7KRaDAAchXQ8wZg6MOW3mwHOnmhBM99dxgHMorqPN49zBdzr+6OYV2bPhsJERERXRqak9cYrN2oLQRrQ04OUkZdIdnXZeMvUMXG1n1B+i7gm2lAeR2zaIT1BQbfD/SaKnkYsUFV5UDGLrx14C18VC5dUv2xgiI8UKyzvyZ6GDD0IaDb1YDMMQ8amswiPt+VhsUbT6DCYKrznLHdOuHxK+PRN8rPIfckIiKitq85eY0PL3Zwepv+apmnJ5TR0XWfXJhm6aWuLJbulymBCa8CQ2YCgn1bRYNUnjjuH45PK9Ilu7sbRdxdV6gGLCPlGX8BftGWe/afBmi0zbuvDblMwD3D4zCuewhe/ukotibbr8y4NTkHW5NzMCohGI+N7YpBsQGtuicRERFdWthj3cHZ9leru3WDIKvj28JQCXxzl32o1kYD9/5qGUFubqgGYDAZMHfnXJjEmlFihaDAq9evhnLad0DPGy1tIHUpygB+fR5Y2hP4ZTZQmF73ec0QFeCJT+8ejM/vvQwJId51nrPjZC5uen8XbvtwF/5KyQN/6UNEREQAR6w7PPsZQeqZv3rDM0C2dMEWdBkLTP0E8Gz5yO3HRz7GicITkn0P9HkAiYHdgcDuQNdxQFk+sO9TYO/HQGm2/YtUlQB7lgN7P7T0YQ9/zNKW0gqjEoKxoctIfPNPJpZuOon8siq7c3afKcDuM3vQO0KLGcNicU2fMM6BTURE1IGxx9qN2kKPdcrEiTCk18zqEbbgVfjddJP0pP2fAz89Kt0XlAA8sBVQ+7T43qcKT+GW9bdI5qyO94/HN1d/A6VcaX+BsQo4thbY9S5w/mDDLx53hSVgd7myRSPptZXqjfhydzo++uMM8krtA3a1QC8Vbh4UhVsHRyEuqIk95kRERNSm8eHFdsLdwdpUWoqTgwZL9sWuWQOPXj1rdpw/BHw8HjDVLKYCpZclVHdq4eqMAExmE6b/Ml2ybLlckOPLq79Ez8CeDVwJy2ouGbuB3e8Byest82DXJ6QXMOxRywOVdYX1ZqioMmHV3gx8sOO0ZHGZugyJC8DUAZGY1DsUvprW3ZeIiIjch8G6nXB3sC7ftw/pd06r2aFQIHHfP5Cp1ZbtSh3w/gigyKZ3eeonQG+bUe1mWpW8Cv/e82/Jvvt63YcnBj7RvBcqTAN2vw/s/wwwlNd/nm+EpQ98wAxA07rPtd5owup/zuL97adxtrCiwXNVChmu7NYJ1/ULxxUJneChYqsIERFRe8Jg3U64O1gXfPElLixYYN1WJySg808/1pywbRGwTRp+MWQmMHlRq+6bXZaN63+8HmWGMuu+GN8YfDflO6jl6pa9aHmBZUn1PR9Y5squj1oLDLrH8ufwDWvZvS4ymUVsOX4Bn+1Kw86U/EbP1yhluCIhGBN7huLKbiHQenIkm4iIqK3jdHvUJFXp0pFodUJCzYa+1PJAYG2RlwHjX23VPUVRxGt7XpOEagB4+fKXWx6qAcsDlKOeBS5/FDj8NfDX20B+iv15+mJg538sfdp9brW0ibSwpUUuEzChZygm9AzFqQsl+PrvTHy//ywKy+teCr3SYMavRy/g16MXoJAJGBDtjysSg3FFQjB6hPlCJmtdLzgRERG5F4N1B2Y4e1ayrYqOqtn451OgolB6wTXLAEU9U9810ZaMLdiWuU2y74auN2Bw6OA6z282pQYYeDfQfzpw8hdg51tA5m7788wG4OAXlrf4iZYHHWOGt/hBx/gQH8y9pgdmTUrE5mM5WHswC9tO5MBgqvsXQkaziL1pBdibVoAlv55AoJcKl3cJxGVxARgcG4DEEB8GbSIionaGrSBu5O5WkDPXXgv9qZpR3bDXXoPf1Bstc1b/tw9QeqHm5ITJwB1ft+p+uiodrl97PXIralo1AjQB+On6n6BVt26BlwZl7AH+egtI/hlAA9/u4QMsAbv7FIes6FhUXoUNSdlYd+gc9qYVwGRu+l81X40Cg2IDMCjWH5fFBqB3pBZqBfuziYiIXI2tINQoURRRdTZLsk8ZGWn54MD/pKEaAEY90+p7/nfffyWhGgCeG/ycc0M1AEQPAaK/BPJSgF1vAwdXSWc5qXZuP7D6bsA/Frj8EaDfnYDKs8W39fNU4Y4h0bhjSDQKy6qw+bilDeSPU7nQGxuYyQSArtJoXekRAFRyGRJDfdArQoveF98SQr0ZtomIiNoQjli7kTtHrI35+Tg1fIRkX9ctm6EM7QS8NQAorpnbGp1HA9N/RGvsv7AfMzbOkOwbETEC7135HoRWzjPdbKU5lsVk9n4EVBbVf55GawnXg+4FguIddvvyKiP2nCnA9pO52H4yF6l5ZY1fVAelXEDnIG/Eh3gjMcQH8SE+SAjxRkygF+RsIyEiInIIzgrSTrgzWFccOoS0W2+r2aFQoNuhgxAOfw38+LD05BnrgbiRLb5XlakKN6+7GWeKz1j3eSg88MN1PyDCO6LFr9tqVWXAgS+AXe9YlkdvSNwoYPD9QOJVrZ4P21ZGfjl2ns7D32kF+DutAJkFDU/h1xiVQobOQV6IDfRCTKAnogM9ERNg+ThMq4FCXseS9URERFQntoJQo6psHlxUhoVZntv7c6n0xKghQKx0ZLu5PjnyiSRUA8C/+v3LvaEaAFRewJAHgUH3Acd/tDzoWN+Kjqk7LG/eocDAGZb5sLWOqT860BPRgdG4/bJoAEB2cSX2phXgn7QC7E0twIkLJWjOf3+rjGYkZ5cgObvE7phCJiDS3wMxgV6I8PdAqK8Gob4ahGg11o99PRSu/y0CERHRJYDBuoMy2PVXRwDHfrSfom7kM61aEvxM8Rl8dPgjyb4egT1wZ/c7W/yaDidXWFZm7HkjkPaHJWCn/Fb3uaXZwPZFwI4lluXS+99pGcVWtGKqQBuhWg2m9A3HlL7hAIDiCgOOnivGkaxiJGXpcCSruMXtI0aziLT8cqTl17+YjodSjlCtBiG+amvo7uSjQaCXCv5eKgR6qRBw8U2jZI83ERFRNQbrDspuqr3ISGDfCulJoX2A+PEtvodZNGP+X/NhMP9/e3ceHUWV9g/8W71m6SwQlhCyQWQJm4hsGhgWJfjGGcRRBvTAJM7g4Ki/GZFBGJX1HUTmgKPjnFdHUcFRxuUFUcejBF4IEVlEWURA6ISdGEIMSSch6bV+f1R3p6uTdNKh05VOfz/n1Kmue29VPZ1mefrm1r0N8zqrBTWW37YcGlUH/KMnCNKQjz4/AyrOAN+8JQ0Vqato3FZ0SMl30XYgsgswdIY0HrvXzTf0RaQpcZFa3J7RDbdndHOXmeqtOFVajdNXqnG6tBqnr9TAWFaN8hrLDd+vzmrH2fLaViXv0To1uhp06BqlQ2ykVtoitIiN1Dj3WsRGaJoo0yJCq2LPOBERdSodMLuhYLBe9hoK0iMBOPeVvFHWH28oSdxi3IJDZYdkZXMGzUFmQmabrxk0XfsC2f8NTHpG6sk/uB649HXTbeuuOR+GfA3oNkBa7n3IfUBCRruFFxuhxah0ac5rTxW1Fpy+Uo2z5bU4/9N1XKhw7n+6jmqzLeBx1FrsqK2oa9O4cLVKQJRWjUidGlE6NSJ1GkS5XmsbyqK96l3tdWoVdBpp02vU0Dtfe5ZLdVIZk3giImpvfHhRQUo+vFiUPRXWCw0P7CX9v+mIu/o/DQ00kcBTZ9o83dzV61dxz9Z7UG1tGOfb29AbW6ZtQZS27VPYKar0GHDwDeDYh4ClpuX2vYZLSfbgXwZsPHZbiaKIiloLzldISfaFiuv4saoepVV1KDWZccVUj4raG+/t7shkCbdaBb22IQnXqlXQqASoVQI0agEaVcvHapUglakFaL2OpbYqqAVApRIgCAJUAqAWBKgEAYIAqAQBKpVz797Q0Fbl1bZR+xbaOr9HSN8npDrXVwtBECA46wRnHdDQ3lUPrzbOS8F5RXeZ60uL5zUhNHwvF7zu6Xm+5z0aXstjaI6v70r8IkVEgcKHF8kn0W6H9ccfZWW62u/ljfrdeUNzOK85uEaWVAPAkrFLQjepBoDEocAvXgSmrgJOfAIceVcak92cH49IW/6zQOptwMC7pfHY7diT3RxBEJBg0CPBoMeI1C5NtjHb7CgzmVFqqkdpVT2uOPelpnqU15hRUWtBRa0F165b/VrspqOw2B2w2B1AE1OYE7n4TNZbPLf5Fr7O9X3PFu7axnjbes+Wvq/w60zb8ctg2/z1/mHIGdpL6TDcmFiHIduVK4DVKivTVh0APGeRy5zW5uvvvrgb285tk5Xd3fduZPXOavM1OxRdNDD8AWmrOAsc/be06EyVjyn7LuyTtvxnpeEiA3OkJLv3SEDVMaa/02vUSOkahZSuvr/8OBwiTPVW/ORMtH+qseDadQuq662orrfBVGeFyb23wlRnc+6tqLXYg/RuiNrG1+9wW/w62S6/AA69L7FEwWS1+15wLdiYWIch76n2BL0Wao1HN55KC/TLbtO1a621+MuBv8jK4vRxWDhyYZuu1+F17QNMehqYsFgag33sf4HjHwHXy5s/p/wUsOcUsOdvQHQPoP9UqTe7z4Qb+i1BsKhUAuKjdIiP0iGju3/n2uwOVNfbUFVnRa3FhjqLHdctdly32Jx7e0OZtaG+zqNNnVUqs9gc0maX9mabHVY7kxAiIlIOE+sw5D3Vni5OI//1Xt8JQGR8m679j8P/QGltqazsTyP/hITIhDZdL2SoVEDqWGm763ng7G7g+y3AyU8Bc1Xz59WWSUvIH/4XoNZJ84b3nQj0nQQkDQdUnWs6O41ahS7Oafvag8Mhuod8uBNvZ/JttjpgsdthbqLc7hBhdUh7m12U9g4RNrsDNofvY9drae+Azd5w7BCdmwMNr0VpzLvdo1wUIR07X/s8zyG9dp/naOI8fr8gIlIEE+swZL10UXas1ZnkDQb+vE3XPXb1GN49+a6sbEziGNyTcU+brhey1Brgpjuk7e51wJldwA+fAae/AGqvNn+e3SKN2T73JbDzv6Ul1dPHOxPtiUDCTQGfyq+zUakERKjUYT+/tuuZdFGUBhKIougeUCCVie5RC+69s0z0PN9V732OxzU9z4fYcI7rjrIYPBJ+zzgaYvO8pu9vB20dsuH7vOYrWxrl0V7XDfY9ff/sWvhMfNYStY/0hGilQ5BhYh2GGq26GOn5NJcgDUvwk9VuxdK9S2X/mOvVeiy9bWl4P5ChjQAG/Je0ORzA5W+BU58Bpz4Hrv7g+9z6KuCH/0gbAER1k3rEU8ZI+143B3RhGuo83LN0uP/qhfHfQSKiIGJiHYYarbpo8HigLPU2wNDD72uuP7YeRZXyVRsfufkRpMamtinGTkmlAlJGSdudy4GfiqUE+9Tn0oONYgsP9l0vlyfaaj3Qe0RDop00Aojp2e5vg4iIiJrGxDoMNVp1Mdojocv8hd/XM14z4rVjr8nKBnQZgNzBuW2KL2wkZAC3Py5t9Sbg3B7gTIG0lZ9q+Xy7uWG2EdfaPoZEaWx2r5sbttjeHEJCREQUBEysw4zDbIatrExWpo32WJEv07/x1XaHHcv2LoPN0XANtaDGyqyV0Kq0Ps4kmYhYaQq+gTnSsakEOLNbGp99ZjdQU+r7fJeaUmks9+kvGsqiukkJdo9MaeueCXQfAOgNgX8fREREYYyJdZixXi5pVOYeCtJrOBDv39CNd06+g2Plx2RleYPzMChhUFtDJACITWqYK1sUgWtngQsHgIv7pf3Vk62/1vVyoPj/pM1TXArQfSDQY6CUbHfrLy3lHtWVPdxERERtwMQ6zFgvy4eBqHV2qLXOBw5vutOva100XcQ/Dv9DVpYem45Hbn7khmIkL4IgJbxd+0qJNgDUXQMufQNc2A9cPACUHG7dMuueqi5KW9F2ebk+Tpqf23VPz9eGnky6iYiImsHEOsx4L2Wu9RxfnTq21dcRRRHL9y1Hvb1eVr7i9hWI0ETcUIzUCpFdgH5TpA2QZhypKAZKjjiXUj8K/Pid7zm0m2OualiO3ZtaL/WmxyVLY7fjejv3HscR8Uy+iYgoLDGxDjO2cvmKgJpI11KgApA8qtXX2WzcjK9Lv5aVzRowCyN6jrjREKktVCqgWz9pGzZDKnM4gMpzUpJd+r00vV/ZSWlYidjGJWDtZun8a2ebb6ONlhJsQ09phpno7tJm6CGtNGnoLu2ju0vTERIREXUSTKzDjO2qfIESTYSzx7pHZqtXW7xSewXrvlknK+sV3QtP3PpEACKkgFGpGoZwDL63odxaD/xkBMp+kMZqXz0lJdyV5wGPh1DbzFoLlJ+WtpboYxsS76iuUm93ZBfpz2Jkl6aPI+KkRXiIiIg6GP7vFGYaJdauHuuUMa06XxRF/GX/X1BjlY/nXXbbMkRrO9bqR9QMbQSQOFTaPNlt0pjrijNSj3TFWel1xRng2jnAVt/k5W6I2SRtFcX+naePBXQGaWYTz72sLMZHnQHQREo/C02ktNAOh68QEdENYmIdZuxXvYaCuHqsWzm++otzX6DgUoGsbFrGNGT1zgpEeKQktcb5oGKfxnUOB1D9I2C6DFRdcu4vA6ZLzv1loOZK8GJ1JeTVgbqgAGginIm2c9NGyvfuemdCrtYDai2g1jk3rcfeq1yl9Wrj3c7ZRqUGVBpAUHm8du5VqkC9WSIiaidMrMPMjfRYX6m9glUHVsnKEiIS8NSopwIWH3VQKpU0bjquN5Ayuuk2NgtQXeJMtEuA2jKgpgyovdqwd712WIMbf4tEwFYnbR2ZO9H2TMA1Xkm4c/NMyD3rBFXDBngcC151gke50ES59zlC83WNyl1tnXvA4zcGnsd+1jVqBx91vs5rqc5XXPBRd4O/FQnIb1WUjoH3D+n7d0SpY4EuaUpH4cbEOoyIothEYm2XHjLrku7zXIfowDN7nkGV1ywTz4x9BnH6uECHSqFIo5P+HLXwZwmiCNRXAjXORLu2DKgtB+oqpWkE65177+P2GIoSahw2ADbA3mJLIqLwcN8bTKxJGQ6TCaJV3lOoiXBIvdUtfAvecHwDDpQekJVlp2VjStqUgMdJnZwgOB9I7AJ079/686x1UrJdXyntLTWAuVraW2oBcw1gqXbuXWXOeneZ87XIzJSIiAKPiXUY8e6tBgB1hL3F8dXHy4/j5UMvy8oSoxOx9LalAY2PyCdtpLTF9rqx64ii1PttrfPYm6VhINZ6j71HG1u9V10dYLcAdqvX3iI9BOp+bZW/dni1D8QsLERE1GEwsQ4j3nNYq3QOqNTwOb76uvU6Fn25CDaxIQEQIOC5cc9xCAiFJkFoSNKV5nDIk23RATjsUsItOvcOh8exq85ZJmvr8DrP2V52nl06VxSdmwOAc+8udzRs8GjnWQevdq25HryvI0plgPM1vI79rHO/ho86X+e1ts67HXzUNXFeW3nfr20XUTgG3j+k799RRcQrHYEME+sw0uQc1poIIHFYs+esObgG503nZWVzh87FqMTWLyZDRM1QqQCVXpruj4iIQh7nbwojtrImZgRJukV66KwJHxk/whbjFlnZ0G5D8fvhv2+3GImIiIhCFRPrMNJoOfMIe7O91QdLD2LlvpWysihNFJ4f/zy0Km27xUhEREQUqphYh5Em57DuObhRuwumC5hfMF82rhoAnh37LFJjU9s1RiIiIqJQxcQ6jNjK5CvjST3WQ2RlJosJj+98vNF81bmDcvGLjF+0e4xEREREoYqJdRixlV6WHWsiRaDHoIZ6hw0Ldy/E2aqzsnYTkidg/q3zgxIjERERUahiYh1GbOU/yY41PRPdU47V2+rxZMGT2FuyV9amf5f+WPOzNVCr1EGLk4iIiCgUMbFuo4MHDyInJwfx8fGIjo7G2LFj8cEHHygdVrMcZjMc182yMk3qAABAlbkKD+c/jF0Xd8nqu0Z0xcuTX0a0NjpocRIRERGFKs5j3Qa7du3C1KlTERERgVmzZiEmJgabN2/GzJkzcfHiRSxYsEDpEBuxV1Q0KtP0HYqSmhI8suORRsM/dCodXpr0EpIMScEKkYiIiCikCaIYkOWcwobNZsPAgQNx6dIl7N+/H8OHDwcAVFVVYfTo0Th37hxOnz6NtLS0Fq9lMpkQFxeHqqoqxMbGtmvcdd8fx7n7728oEER888p9eK20EJXmSlnbGG0MXpr8EheBISIiorDnT77GoSB+2rlzJ4qLi/Hggw+6k2oAiIuLw9NPPw2LxYKNGzcqF2AzLp/eIzuuihLw1/OfNEqqe0T1wMb/2sikmoiIiMhPTKz9VFBQAADIzs5uVDd16lQAwO7du4MZUqsc/UE+froyqnGbjLgMvJvzLvp16RekqIiIiIg6D46x9pPRaAQA9OvXOPlMTEyEwWBwt/FmNpthNjc8QGgymdonyCZEXpPPS10VLciOp6RNwbLbliFOHxe0mIiIiIg6EybWfqqqkhLUuLimE9DY2Fh3G2+rV6/GihUr2i02XyKvyZP4ammWPYzoMQLzb52P4T2GBz8oIiIiok6EiXUQ/fnPf8aTTz7pPjaZTEhJSQnKvTVp6ThVUgVDrQMx1x3o3b0n3sj+K0YljoIgCC1fgIiIiIh8YmLtJ1dPdXO90iaTCV26dGmyTq/XQ6/Xt1tsvoxf8m7DgcMO2C3uxWGIiIiI6Mbx4UU/ucZWNzWOurS0FDU1NU2Ov+5QVGom1UREREQBxsTaTxMmTAAA5OfnN6rbtm2brA0RERERhQ8uEOMnm82GAQMG4PLly80uEHPq1Cmkp6e3eK1gLhBDRERERP7zJ1/jGGs/aTQarF+/HlOnTsXPfvYz2ZLm58+fx9q1a1uVVBMRERFR58LEug0mTZqEPXv2YNmyZXj//fdhtVoxdOhQrFmzBjNnzlQ6PCIiIiJSAIeCKIhDQYiIiIg6Nn/yNT68SEREREQUAEysiYiIiIgCgIk1EREREVEAMLEmIiIiIgoAJtZERERERAHAxJqIiIiIKACYWBMRERERBQATayIiIiKiAGBiTUREREQUAEysiYiIiIgCgIk1EREREVEAaJQOIJyJoghAWoOeiIiIiDoeV57mytt8YWKtoOrqagBASkqKwpEQERERkS/V1dWIi4vz2UYQW5N+U7twOBwoKSlBTEwMBEEIyj1NJhNSUlJw8eJFxMbGBuWeFDj8/EIfP8PQx88w9PEzDG3B/vxEUUR1dTWSkpKgUvkeRc0eawWpVCokJycrcu/Y2Fj+YxLC+PmFPn6GoY+fYejjZxjagvn5tdRT7cKHF4mIiIiIAoCJNRERERFRADCxDjN6vR7Lli2DXq9XOhRqA35+oY+fYejjZxj6+BmGto78+fHhRSIiIiKiAGCPNRERERFRADCxJiIiIiIKACbWREREREQBwMSaiIiIiCgAmFiHiYMHDyInJwfx8fGIjo7G2LFj8cEHHygdFrXC5cuX8eKLLyI7OxupqanQ6XRITEzEfffdhwMHDigdHrXRmjVrIAgCBEHA/v37lQ6HWumjjz7ClClTkJCQgIiICPTp0wcPPPAALl68qHRo1AJRFLFlyxZMmjQJvXr1QlRUFAYMGIB58+bhzJkzSodHTu+88w7mzZuHkSNHQq/XQxAEbNiwodn2JpMJTz75JNLS0qDX65Geno6FCxeipqYmeEF74KwgYWDXrl2YOnUqIiIiMGvWLMTExGDz5s04f/481q5diwULFigdIvmwePFirFmzBhkZGZg4cSK6d+8Oo9GIrVu3QhRFbNq0CTNnzlQ6TPLD999/j5EjR0Kj0aC2thb79u3D2LFjlQ6LfBBFEY888ghee+01ZGRkYOrUqYiJiUFJSQl2796Nd999F+PGjVM6TPJhwYIFeOGFF9CrVy/cc889iI2NxdGjR5Gfnw+DwYC9e/diyJAhSocZ9tLT03H+/Hl069YN0dHROH/+PN566y3k5eU1altbW4tx48bhyJEjyM7Oxi233ILDhw8jPz8fo0aNQmFhISIiIoL7BkTq1KxWq5iRkSHq9Xrx8OHD7vLKykqxf//+ok6nE8+dO6dcgNSizZs3iwUFBY3KCwsLRa1WK3bp0kWsr69XIDJqC4vFIo4YMUIcM2aMOHv2bBGAuG/fPqXDoha8+OKLIgDx0UcfFW02W6N6q9WqQFTUWj/++KOoUqnEtLQ0sbKyUlb3wgsviADEhx56SKHoyNP27dvdecnq1atFAOJbb73VZNulS5eKAMRFixbJyhctWiQCEJ977rn2DrcRDgXp5Hbu3Ini4mI8+OCDGD58uLs8Li4OTz/9NCwWCzZu3KhcgNSiX/7yl5gwYUKj8vHjx2PSpEm4du0ajh07pkBk1BarVq3C8ePH8eabb0KtVisdDrVCXV0dVqxYgb59++Kll15q8nPTaDQKREatde7cOTgcDmRlZSEuLk5W9/Of/xwAcPXqVSVCIy933nkn0tLSWmwniiLWr18Pg8GAJUuWyOqWLFkCg8GA9evXt1eYzWJi3ckVFBQAALKzsxvVTZ06FQCwe/fuYIZEAaTVagHwP/VQcejQIaxatQrLli3DoEGDlA6HWik/Px/Xrl3D9OnTYbfbsWXLFjz//PN49dVXUVRUpHR41Ar9+vWDTqfDV199BZPJJKv7z3/+AwC44447lAiN2shoNKKkpARZWVmIjo6W1UVHRyMrKwtnzpwJ+vMP/N+4kzMajQCkf1S8JSYmwmAwuNtQaLlw4QJ27NiBXr16YejQoUqHQy0wm8349a9/jeHDh+Opp55SOhzyw7fffgsAUKvVGDZsGE6fPu2uU6lUmD9/PtauXatUeNQKCQkJeP7557FgwQIMHDhQNsZ6586dePTRR/H4448rHSb5wVd+4yrftm0bjEYjUlJSghYXE+tOrqqqCgAa/erLJTY21t2GQofVasWcOXNgNpuxZs0aDikIAUuXLoXRaMS3337LzyvElJWVAQBeeOEFjBgxAl9//TUyMzNx+PBh/O53v8O6deuQkZGB3//+9wpHSr7Mnz8fvXv3xty5c/Hqq6+6y8eNG4cHH3yQv/kLMa3JbzzbBQuHghCFGIfDgby8PBQWFuLhhx/GnDlzlA6JWrBv3z6sXbsWzz77LGcdCEEOhwMAoNPpsHXrVowaNQoGgwHjx4/Hhx9+CJVKhXXr1ikcJbVk5cqVmD17Np5++mlcvHgR1dXV+PLLL1FfX4+JEyfik08+UTpE6gSYWHdyrm9yzX1jM5lMzX7bo47H4XDgN7/5DTZt2oTZs2fLel2oY7LZbMjNzcWwYcOwePFipcOhNnD9Gzly5EgkJSXJ6oYMGYK+ffuiuLgYlZWVCkRHrbFjxw4sW7YMjz/+OBYvXozk5GQYDAaMGzcOn376KbRaLaeeDTGtyW882wULE+tOzjX2qKlx1KWlpaipqWl2fBJ1LA6HAw899BA2btyIBx54ABs2bIBKxb/CHV1NTQ2MRiOOHDkCnU7nXhRGEAT3jDy33XYbBEHA1q1blQ2WmjRgwAAAQHx8fJP1rvK6urogRUT++vzzzwEAkyZNalSXmJiIgQMHoqioSLFFRch/vvIbz/Jg5zgcUNTJTZgwAatXr0Z+fj5mzZolq9u2bZu7DXVsrqT67bffxsyZM/Gvf/2L43RDhF6vx29/+9sm6woLC2E0GjFt2jR0794d6enpwQ2OWsWVjJ08ebJRndVqRVFREaKjo9G9e/dgh0atZLFYADQ/pd7Vq1ehUqncMy1Rx9evXz8kJSXhq6++Qm1trWxmkNraWnz11Vfo06dPUB9cBMAFYjo7q9Uq9u3b1+cCMWfPnlUsPmqZ3W4Xc3NzRQDijBkzuBBFJ+L6XLlATMeXnZ0tAhBff/11WfnKlStFAOLs2bMVioxa49///rcIQBw8eHCjBWJeeeUVEYCYlZWlUHTUnFBcIIZLmocBLmke2pYvX44VK1bAYDDgj3/8Y5NPrk+fPl22ABCFhry8PGzcuJFLmoeA4uJi3H777SgrK8Pdd9+NgQMH4vDhw9i5cyfS0tKwf/9+JCYmKh0mNcNut2Py5MkoLCxEjx49MG3aNMTHx+PQoUPYuXMnIiMjUVBQgNGjRysdathbv3499uzZAwA4duwYDh06hKysLNx0000ApFlc5s6dC0Dqmc7KysLRo0eRnZ2NESNG4NChQ+4lzXfv3o3IyMjgvoGgp/KkiAMHDoh33XWXGBsbK0ZGRoqjR48W33vvPaXDolZw9Wr62pr7Nk8dG3usQ8uFCxfEvLw8MTExUdRqtWJKSor42GOPiVeuXFE6NGqF+vp6cfXq1eItt9wiRkVFiRqNRuzdu7c4e/Zs8cSJE0qHR04t/Z+Xm5sra19ZWSk+8cQTYkpKiqjVasXU1FRxwYIFoslkUiR+9lgTEREREQUApxQgIiIiIgoAJtZERERERAHAxJqIiIiIKACYWBMRERERBQATayIiIiKiAGBiTUREREQUAEysiYiIiIgCgIk1EREREVEAMLEmIiIiIgoAJtZERDcoPT0dgiC0aTt37lyT1+hsjEYj9Ho9BEGAWq3GiRMngnbvvLw82c+2oKAgYNcuLy+HwWBwX7uwsDBg1yai0MPEmoiI2t38+fNhsVgAADNmzMCgQYMUjigwunXrhscee8x9/Ic//AEOh0PBiIhISRqlAyAiCnU5OTkoKyuTlZ04cQInT550H6elpWHkyJGNzo2Ojm72Gp3Fjh078Nlnn7mPn3nmmaDef9SoUaipqXEfd+/ePaDX/9Of/oQXX3wRFosFR48exZtvvom5c+cG9B5EFBoEURRFpYMgIupsli9fjhUrVriPc3NzsWHDBuUCUtBdd92Fbdu2AQBuvfVWfPPNNwpHFHj33nsvtm7dCgDIzMwM6lAXIuo4OBSEiKgD8DXGuqCgQFaXl5eHy5cvY+7cuUhKSkJkZCSGDRuGf/7zn+5zDh48iGnTpqFr166IjIzEyJEjsWnTpmbvb7PZsGnTJkybNg3JycmIiIhATEwMhg4dioULF+LSpUttel+nT59Gfn6++3jOnDmtev92ux1///vfMWzYMERGRiIpKQnz5s1DeXk5AMBkMuGpp55Cnz59oNfrkZqaiieeeAImk6nRtVsaY+1Zl56eDofDgfXr12Ps2LEwGAwwGAwYP348Pv/882bf5+zZs92vT548iR07dvjzYyKizkIkIqKAW7ZsmQjAveXm5vpsn5aWJmvvadeuXbK622+/XezRo4eszLUtXLhQ/PDDD0WtVttk/csvv9zo3iUlJeLo0aObbO/aYmJixI8//tjvn8Nzzz0nu87x48db9f6nT5/eZBwZGRliUVGR2L9//ybrx4wZI1qtVtm1c3NzZW127dolq/es69mzp5idnd3ktQVBELds2dJk/JWVlaIgCO628+bN8/tnRUShjz3WREQhZu/evbh69SpGjx6NsWPHyurWrVvn7hUeP358o4cElyxZgrq6Ovex1WpFTk4Ovv76a3dZcnIycnJykJWVBZVK+m+iuroaM2fOxNGjR/2K1bPnNj4+HpmZma06b+vWrUhOTkZ2djYMBoO7vLi4GMOGDcPp06fRv39/3HHHHVCr1e76AwcO4MMPP/QrRk9XrlxBfn4+evXqhSlTpqBbt27uOlEUsWjRoibPi4uLk/2s2WNNFJ6YWBMRhaA33ngDBw4cwL59+zBjxgx3ucPhgNlsxhdffIHCwkJ89913GDFihLu+srISBw8edB+//fbbOHLkiPv40Ucfxfnz5/HZZ59hz549+PLLL91DU+rr6/Hss8/6FafneOrMzMxWTyU4ZcoUFBUVYdu2bXj//fdlddevX0deXh5++OEH7NixA2vXrpXVb9++3a8Yvd11110oLi5Gfn4+vv/+e/To0cNdZzQaceHChSbPGzp0qPt1cXExrl27dkNxEFHo4awgREQhJiMjAw899JD7OCsrS9ZLO2nSJEyePBkAoFarMXHiRBw6dMhdf/nyZffrjz76SHZto9GIX/3qV7IynU4Hs9kMQEpazWYz9Hp9i3HW1dXJxjx79v62ZMmSJe57ZGVlNapfuXKlO0m/4447ZHWe768t/va3vyEyMhIA0LNnT4wZMwaffvqp7PqpqamNzvN+f1euXEGXLl1uKBYiCi1MrImIQsyQIUNkxzExMX7Vu5JkADh79qysrqXeXrPZjJKSEvTp06fFOL17bGNjY1s8x8Wz99c7/tjYWKSkpDRb7/n+/GUwGDBw4EBZWVxcXKuu7/3+Kioq2hwHEYUmDgUhIgox8fHxsmPXOGiX9u4lra2tbVU77zibmrGjNecG8/0lJCQ0KvMcw+1LVVWV7Ji91UThh4k1EVEY8+553r9/P0RR9Ll594g3JyoqStab7Joqr7Pyfn+JiYkKRUJESmFiTUQUxqZNmyY7nj9/fpMrQBYVFWHNmjVYuXKlX9e/9dZb3a9PnDgBsROvSXbs2DH364yMDPZYE4UhJtZERGEsLy8PgwcPdh/v27cPqampyMrKwvTp0zF58mT07t0b/fr1w+LFi3HmzBm/ru96iBKQhkp4LvPemVRWVsre25133qlgNESkFCbWRERhTKfT4YsvvsDIkSPdZWazGXv37sXHH3+MXbt2oaSkxF2n0fj3zPv9998vO/ZchbEz2bFjh6w33vt9E1F4YGJNRBTmkpOTsX//frz33nu49957kZqaioiICGi1WnTr1g2jR4/GY489hk8++QSvvPKKX9fOzMzElClT3MfvvPNOoMPvEDzfV2ZmJnusicKUIHbmAW9ERKS47du3Izs723383XffyabTC3VlZWVITk6G1WoFALz++uuYO3euwlERkRLYY01ERO1qypQpyMnJcR+vWrVKwWgCb+3ate6k+uabb5Yt3kNE4YU91kRE1O6MRiOGDBkCi8UClUqFY8eOYdCgQUqHdcPKy8uRnp7untu7oKAAEyZMUDgqIlIKE2siIiIiogDgUBAiIiIiogBgYk1EREREFABMrImIiIiIAoCJNRERERFRADCxJiIiIiIKACbWREREREQBwMSaiIiIiCgAmFgTEREREQUAE2siIiIiogD4/2m0k8UUMsQeAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# function that creates pyomo model for the dynamic reaction system\n",
    "def create_dynamic_model(F, CAF, k1, k2, k3):\n",
    "    \"\"\"Creates reactor model\n",
    "    Argument:\n",
    "        F: feed volumetric flow rate in m^3/min\n",
    "        CAF: feed concentration of A, gmol/m^3\n",
    "        k1: rate constant for A to B, min^-1\n",
    "        k2: rate constant for B to C, min^-1\n",
    "        k3: rate constant for 2A to D, m^3/(gmol.min)\n",
    "\n",
    "    Returns:\n",
    "        model: pyomo model\"\"\"\n",
    "\n",
    "    # Create a Pyomo model\n",
    "    model = pyo.ConcreteModel()\n",
    "\n",
    "    # Define time as a continuous set\n",
    "    model.t = dae.ContinuousSet(bounds=(0, 10))  # time from 0 to 10 minutes\n",
    "\n",
    "    # Define variables for concentrations and reactor volume\n",
    "    model.CA = pyo.Var(model.t, within=pyo.NonNegativeReals) # concentration of A\n",
    "    model.CB = pyo.Var(model.t, within=pyo.NonNegativeReals) # concentration of B\n",
    "    model.CC = pyo.Var(model.t, within=pyo.NonNegativeReals) # concentration of C\n",
    "    model.CD = pyo.Var(model.t, within=pyo.NonNegativeReals) # concentration of D\n",
    "    model.V = pyo.Var(within=pyo.NonNegativeReals, bounds=(0, 100), initialize=0) # reactor volume\n",
    "\n",
    "    # Define derivatives of concentrations\n",
    "    model.dCA_dt = dae.DerivativeVar(model.CA, wrt=model.t)\n",
    "    model.dCB_dt = dae.DerivativeVar(model.CB, wrt=model.t)\n",
    "    model.dCC_dt = dae.DerivativeVar(model.CC, wrt=model.t)\n",
    "    model.dCD_dt = dae.DerivativeVar(model.CD, wrt=model.t)\n",
    "\n",
    "    # CA differential equation\n",
    "    def ca_rate_ode(m, t):\n",
    "        return m.dCA_dt[t] == (F/m.V)*CAF - (F/m.V)*m.CA[t] - k1*m.CA[t] - 2*k3*m.CA[t]**2\n",
    "    model.CA_ode = pyo.Constraint(model.t, rule=ca_rate_ode)\n",
    "\n",
    "    # CB differential equation\n",
    "    def cb_rate_ode(m, t):\n",
    "        return m.dCB_dt[t] == -(F/m.V)*m.CB[t] + k1*m.CA[t] - k2*m.CB[t]\n",
    "    model.CB_ode = pyo.Constraint(model.t, rule=cb_rate_ode)\n",
    "\n",
    "    # CC differential equation\n",
    "    def cc_rate_ode(m, t):\n",
    "        return m.dCC_dt[t] == -(F/m.V)*m.CC[t] + k2*m.CB[t]\n",
    "    model.CC_ode = pyo.Constraint(model.t, rule=cc_rate_ode)\n",
    "\n",
    "    # CD differential equation\n",
    "    def cd_rate_ode(m, t):\n",
    "        return m.dCD_dt[t] == -(F/m.V)*m.CD[t] + k3*m.CA[t]**2\n",
    "    model.CD_ode = pyo.Constraint(model.t, rule=cd_rate_ode)\n",
    "\n",
    "    # Initial conditions\n",
    "    def init_conditions(m):\n",
    "        yield m.CA[0] == CAF\n",
    "        yield m.CB[0] == 0\n",
    "        yield m.CC[0] == 0\n",
    "        yield m.CD[0] == 0\n",
    "    model.init_conditions = pyo.ConstraintList(rule=init_conditions)\n",
    "\n",
    "    # Objective: Maximize concentration of C at t = 10 minutes\n",
    "    def objective_function(m):\n",
    "        return m.CC[10]\n",
    "    model.objective = pyo.Objective(rule=objective_function, sense=pyo.maximize)\n",
    "\n",
    "    return model\n",
    "\n",
    "# Parameters from data\n",
    "conc_A_feed = 10000  # gmol/m^3\n",
    "k1 = 5/6  # min^-1\n",
    "k2 = 5/3  # min^-1\n",
    "k3 = 1/6000  # m^3/gmol/min\n",
    "feed_flowrate = 1 # m^3/min\n",
    "\n",
    "# reactor model object\n",
    "reactor_model = create_dynamic_model(feed_flowrate, conc_A_feed, k1, k2, k3)\n",
    "\n",
    "# Discretize the model using collocation\n",
    "discretizer = pyo.TransformationFactory('dae.collocation')\n",
    "discretizer.apply_to(reactor_model, nfe=20, ncp=10)\n",
    "\n",
    "# Solve the model\n",
    "solver = pyo.SolverFactory('ipopt')\n",
    "results = solver.solve(reactor_model, tee=True)\n",
    "\n",
    "# print optimal reactor volume\n",
    "optimal_reactor_volume = pyo.value(reactor_model.V) # m^3\n",
    "print(f\"Optimal reactor volume is: {round(optimal_reactor_volume, 2)} m^3\")\n",
    "\n",
    "# Extract the concentration of C at t = 10 minutes\n",
    "optimal_concentration_C = pyo.value(reactor_model.CC[10])*10**(-3) # kmol/m^3\n",
    "print(f\"Maximum Concentration of C at optimal reactor volume is: {round(optimal_concentration_C, 2)} kmol/m^3\")\n",
    "\n",
    "\n",
    "# Extract results for plotting\n",
    "time = [t for t in reactor_model.t] # minutes\n",
    "CA_values = [(pyo.value(reactor_model.CA[t])*10**(-3)) for t in reactor_model.t] # kmol/m^3\n",
    "CB_values = [(pyo.value(reactor_model.CB[t])*10**(-3)) for t in reactor_model.t] # kmol/m^3\n",
    "CC_values = [(pyo.value(reactor_model.CC[t])*10**(-3)) for t in reactor_model.t] # kmol/m^3\n",
    "CD_values = [(pyo.value(reactor_model.CD[t])*10**(-3)) for t in reactor_model.t] # kmol/m^3\n",
    "\n",
    "# Plot the results\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(time, CA_values, label='$\\mathrm{C_{A}}$')\n",
    "plt.plot(time, CB_values, label='$\\mathrm{C_{B}}$')\n",
    "plt.plot(time, CC_values, label='$\\mathrm{C_{C}}$')\n",
    "plt.plot(time, CD_values, label='$\\mathrm{C_{D}}$')\n",
    "plt.xlabel('Time (min)', fontsize=16,fontweight='bold')\n",
    "plt.ylabel('Concentration ($\\mathbf{\\dfrac{kmol}{m^3}}$)', fontsize=16,fontweight='bold')\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.title('Concentration Profiles Over Time', fontsize=16,fontweight='bold')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pyomo Enables Dynamic Chemical Reactor Design"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the concentration profile above, $C_{A}$ decreases with time. This makes sense beacause species A is the reactant and it is continuously been consumed to form the products. \n",
    "\n",
    "Also, it is observed that $C_{B}$ and $C_{D}$ increase at the early phase of the reaction, but then start to decrease at a later time. This is because at the early phase of the reaction, only species B and D are produced until a later time when species C begins to form, and since the rate constant of species C ($k_2$) is faster than the rate constants of species B and D ($k_1$ and $k_3$ respectively), $C_{C}$ continue to increase whereas $C_{B}$ and $C_{D}$ deplete.\n",
    "\n",
    "The optimization result shows that the optimal reactor volume required to maximize the desired product (species C)  is $\\mathrm{11.35\\,m^{3}}$ and the corresponding species C concentration is $\\mathrm{5.06\\,kmol\\cdot m^{-3}}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take Away Messages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Pyomo is a powerful tool in solving complex dynamic optimization problems\n",
    "* Proper modelling of the problem is key to obtaining quality results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AgUvMvTPAZlV"
   },
   "source": [
    "## Reference\n",
    "Hart, W. E., Laird, C. D., Watson, J. P., Woodruff, D. L., Hackebeil, G. A., Nicholson, B. L., & Siirola, J. D. (2017). Pyomo-optimization modeling in python (Vol. 67, p. 277). Berlin: Springer.\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}