{ "cells": [ { "cell_type": "markdown", "id": "4dae19b6", "metadata": {}, "source": [ "\n", "*This notebook contains material from [CBE60499](https://ndcbe.github.io/CBE60499);\n", "content is available [on Github](git@github.com:ndcbe/CBE60499.git).*\n" ] }, { "cell_type": "markdown", "id": "4ad3d2d1", "metadata": {}, "source": [ "\n", "| [Contents](toc.html) | [Tag Index](tag_index.html) | [1.1 60 Minutes to Pyomo: An Energy Storage Model Predictive Control Example](https://ndcbe.github.io/CBE60499/01.01-Pyomo-Nuts-and-Bolts.html) >

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 1, "link": "[1.0 Getting Started with Pyomo](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0-Getting-Started-with-Pyomo)", "section": "1.0 Getting Started with Pyomo" } }, "source": [ "# 1.0 Getting Started with Pyomo" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[1.0.1 Local Installation](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1-Local-Installation)", "section": "1.0.1 Local Installation" } }, "source": [ "## 1.0.1 Local Installation" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.1 Install anaconda](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.1-Install-anaconda)", "section": "1.0.1.1 Install anaconda" } }, "source": [ "### 1.0.1.1 Install anaconda" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.1 Install anaconda](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.1-Install-anaconda)", "section": "1.0.1.1 Install anaconda" } }, "source": [ "Anaconda is a Python distribution that supports the `conda` package manager and environments. Follow these steps if you do not have anaconda or miniconda already installed:\n", "* https://docs.anaconda.com/anaconda/install/\n", "* https://docs.conda.io/en/latest/miniconda.html" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.2 Create a new conda environment](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.2-Create-a-new-conda-environment)", "section": "1.0.1.2 Create a new conda environment" } }, "source": [ "### 1.0.1.2 Create a new conda environment" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.2 Create a new conda environment](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.2-Create-a-new-conda-environment)", "section": "1.0.1.2 Create a new conda environment" } }, "source": [ "`conda` envirorments allow us to maintain several difference Python installations on our computer. This allows for easy switching between Python versions or different versions of Python packages. Using an environment will ensure anything we install for this class does not disrupt your other Python workflows.\n", "\n", "After installing `conda`, run this command in the terminal:\n", "\n", "```\n", "conda create --name spring2021 python=3.8\n", "```" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.2 Create a new conda environment](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.2-Create-a-new-conda-environment)", "section": "1.0.1.2 Create a new conda environment" } }, "source": [ "This creates an environment called `spring2021`. Next we need to activate it:\n", "\n", "```\n", "conda activate spring2021\n", "```\n", "\n", "Whenever you open a new terminal, you will likely need to activate our new environment." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.3 Install IDAES-PSE](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.3-Install-IDAES-PSE)", "section": "1.0.1.3 Install IDAES-PSE" } }, "source": [ "### 1.0.1.3 Install IDAES-PSE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.3 Install IDAES-PSE](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.3-Install-IDAES-PSE)", "section": "1.0.1.3 Install IDAES-PSE" } }, "source": [ "The IDAES-PSE toolset (Institute for the Design of Advanced Energy Systems - Process Systems Engineering) is an open-source framework for multiscale process modeling, optimization, and advanced analytics. The IDAES toolset is built on top of Pyomo. We'll mainly use Pyomo in this class, although we'll highlight a few IDAES capabilities.\n", "\n", "To install `idaes`, run:\n", "\n", "```\n", " conda install -c idaes-pse -c conda-forge idaes-pse\n", "```\n", "\n", "This command will install `idaes` as well as all of its dependencies, including `pyomo`.\n", "\n", "After `conda` finishes, run:\n", "\n", "```\n", "idaes --version\n", "```\n", "\n", "To verify `idaes` is installed in your active environment." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.4 Install Ipopt](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.4-Install-Ipopt)", "section": "1.0.1.4 Install Ipopt" } }, "source": [ "### 1.0.1.4 Install Ipopt" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.4 Install Ipopt](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.4-Install-Ipopt)", "section": "1.0.1.4 Install Ipopt" } }, "source": [ "**Windows and Linux users** - you are in luck. IDAES comes with a nonlinear optimization solver Ipopt compiled with the HSL linear algebra routines; I may refer to this as \"good Ipopt\". Later in the semester, we'll talk about why the linear alegbra library is important. To install the \"good Ipopt\", run:\n", "\n", "```\n", "idaes get-extensions\n", "```\n", "\n", "Linux users: you may need to install some dependencies too: https://idaes-pse.readthedocs.io/en/stable/getting_started/index.html#linux" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.4 Install Ipopt](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.4-Install-Ipopt)", "section": "1.0.1.4 Install Ipopt" } }, "source": [ "This will install Ipopt in a hidden user folder. If you get the error `WARNING: Could not locate the 'ipopt' executable` add the line `import idaes` to the top of your notebook then restart the Python kernel. Our library `helper` takes care of adding Ipopt to the system path on Colab. If you use the instructions below, \"okay\" Ipopt will be installed." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.4 Install Ipopt](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.4-Install-Ipopt)", "section": "1.0.1.4 Install Ipopt" } }, "source": [ "**macOS users** - unfortunately, IDAES does not include binaries for \"good Ipopt\". You can either compile the solver yourself after obtaining an academic license for HSL or you can install the \"okay Ipopt\" using conda:\n", "\n", "```\n", "conda install -c conda-forge ipopt\n", "```" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.5 Install Additional Solvers](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.5-Install-Additional-Solvers)", "section": "1.0.1.5 Install Additional Solvers" } }, "source": [ "### 1.0.1.5 Install Additional Solvers" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.5 Install Additional Solvers](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.5-Install-Additional-Solvers)", "section": "1.0.1.5 Install Additional Solvers" } }, "source": [ "In additional to `ipopt`, we will also use `glpk`. To install it, run:\n", "\n", "```\n", "conda install -c conda-forge glpk\n", "```\n", "\n", "You may also wish to install other commercial solvers, especially if you want to consider large mixed integer optimization problems for research. These solvers offer free academic licenses:\n", "* CPLEX: https://www.ibm.com/products/ilog-cplex-optimization-studio\n", "* Gurobi: https://www.gurobi.com/downloads/end-user-license-agreement-academic/" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.1.6 Install Notebook Spellcheker](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.1.6-Install-Notebook-Spellcheker)", "section": "1.0.1.6 Install Notebook Spellcheker" } }, "source": [ "### 1.0.1.6 Install Notebook Spellcheker\n", "\n", "Here is a great Jupyter extension that adds a [spellchecker](https://jupyter-contrib-nbextensions.readthedocs.io/en/latest/nbextensions/spellchecker/README.html) to your notebooks. To install, run the following commands in the terminal (in your new conda environment):\n", "\n", "```\n", "conda install -c conda-forge jupyter_contrib_nbextensions\n", "jupyter contrib nbextension install --user\n", "jupyter nbextension enable spellchecker/main\n", "```\n" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[1.0.2 Cloud Computing with Google Colab](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.2-Cloud-Computing-with-Google-Colab)", "section": "1.0.2 Cloud Computing with Google Colab" } }, "source": [ "## 1.0.2 Cloud Computing with Google Colab" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "nbpages": { "level": 2, "link": "[1.0.2 Cloud Computing with Google Colab](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.2-Cloud-Computing-with-Google-Colab)", "section": "1.0.2 Cloud Computing with Google Colab" } }, "outputs": [], "source": [ "## Tip: Please put code like this at the top of your notebook.\n", "# We want all of the module/package installations to start up front\n", "\n", "import sys\n", "if \"google.colab\" in sys.modules:\n", " !wget \"https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py\"\n", " import helper\n", " helper.install_idaes()\n", " helper.install_ipopt()" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[1.0.3 Your First Optimization Problem](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3-Your-First-Optimization-Problem)", "section": "1.0.3 Your First Optimization Problem" } }, "source": [ "## 1.0.3 Your First Optimization Problem" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[1.0.3 Your First Optimization Problem](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3-Your-First-Optimization-Problem)", "section": "1.0.3 Your First Optimization Problem" } }, "source": [ "We are now ready to solve your first optimization problem in Pyomo." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.1 Mathematical Model](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.1-Mathematical-Model)", "section": "1.0.3.1 Mathematical Model" } }, "source": [ "### 1.0.3.1 Mathematical Model" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.1 Mathematical Model](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.1-Mathematical-Model)", "section": "1.0.3.1 Mathematical Model" } }, "source": [ "Let's start with a purely mathematical example:\n", "\n", "$$\\begin{align*} \\min_{x} \\quad & x_1^2 + 2 x_2^2 - x_3 \\\\\n", "\\mathrm{s.t.} \\quad & x_1 + x_2 = 1 \\\\\n", "& x_1 + 2 x_2 - x_3 = 5 \\\\\n", "& -10 \\leq x_1, x_2, x_3 \\leq 10 \n", "\\end{align*} $$\n", "\n", "We want to solve the constrained optimization problem numerically." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.2 Define the Model in Pyomo](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.2-Define-the-Model-in-Pyomo)", "section": "1.0.3.2 Define the Model in Pyomo" } }, "source": [ "### 1.0.3.2 Define the Model in Pyomo" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.2 Define the Model in Pyomo](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.2-Define-the-Model-in-Pyomo)", "section": "1.0.3.2 Define the Model in Pyomo" } }, "source": [ "

\n", "Activity: Fill in the missing constraint.\n", "
" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.2 Define the Model in Pyomo](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.2-Define-the-Model-in-Pyomo)", "section": "1.0.3.2 Define the Model in Pyomo" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3 Var Declarations\n", " x1 : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : -10 : 1 : 10 : False : False : Reals\n", " x2 : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : -10 : 1 : 10 : False : False : Reals\n", " x3 : Size=1, Index=None\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " None : -10 : 1 : 10 : False : False : Reals\n", "\n", "1 Objective Declarations\n", " OBJ : Size=1, Index=None, Active=True\n", " Key : Active : Sense : Expression\n", " None : True : minimize : x1**2 + 2*x2**2 - x3\n", "\n", "2 Constraint Declarations\n", " h1 : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 1.0 : x1 + x2 : 1.0 : True\n", " h2 : Size=1, Index=None, Active=True\n", " Key : Lower : Body : Upper : Active\n", " None : 5.0 : x1 + 2*x2 - x3 : 5.0 : True\n", "\n", "6 Declarations: x1 x2 x3 OBJ h1 h2\n" ] } ], "source": [ "import pyomo.environ as pyo\n", "\n", "# Create instance of concrete Pyomo model.\n", "# concrete means all of the sets and model data are specified at the time of model construction.\n", "# In this class, you'll use a concrete model.\n", "m = pyo.ConcreteModel()\n", "\n", "## Declare variables with initial values with bounds\n", "m.x1 = pyo.Var(initialize=1, bounds=(-10, 10))\n", "m.x2 = pyo.Var(initialize=1, bounds=(-10, 10))\n", "m.x3 = pyo.Var(initialize=1, bounds=(-10, 10))\n", "\n", "## Declare objective\n", "m.OBJ = pyo.Objective(expr=m.x1**2 + 2*m.x2**2 - m.x3, sense = pyo.minimize)\n", "\n", "## Declare equality constraints\n", "m.h1 = pyo.Constraint(expr= m.x1 + m.x2 == 1)\n", "\n", "# YOUR SOLUTION HERE\n", "\n", "## Display model\n", "m.pprint()" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.3 Solve using Ipopt](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.3-Solve-using-Ipopt)", "section": "1.0.3.3 Solve using Ipopt" } }, "source": [ "### 1.0.3.3 Solve using Ipopt\n", "\n", "Toward the end of the semester we will learn, in perhaps more detail than you care, what makes Ipopt work under the hood. For now, we'll use it as a computational tool." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.3 Solve using Ipopt](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.3-Solve-using-Ipopt)", "section": "1.0.3.3 Solve using Ipopt" } }, "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", "\n", "This is Ipopt version 3.13.2, running with linear solver ma27.\n", "\n", "Number of nonzeros in equality constraint Jacobian...: 5\n", "Number of nonzeros in inequality constraint Jacobian.: 0\n", "Number of nonzeros in Lagrangian Hessian.............: 2\n", "\n", "Total number of variables............................: 3\n", " variables with only lower bounds: 0\n", " variables with lower and upper bounds: 3\n", " variables with only upper bounds: 0\n", "Total number of equality constraints.................: 2\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+00 3.00e+00 3.33e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", " 1 4.3065612e+00 1.11e-16 2.89e-01 -1.0 4.36e+00 - 6.72e-01 1.00e+00h 1\n", " 2 4.2501103e+00 0.00e+00 5.55e-17 -1.0 1.31e-01 - 1.00e+00 1.00e+00f 1\n", " 3 4.2500000e+00 0.00e+00 4.32e-16 -2.5 6.02e-03 - 1.00e+00 1.00e+00f 1\n", " 4 4.2500000e+00 8.88e-16 1.14e-16 -3.8 4.41e-05 - 1.00e+00 1.00e+00f 1\n", " 5 4.2500000e+00 0.00e+00 3.13e-16 -5.7 1.98e-06 - 1.00e+00 1.00e+00h 1\n", " 6 4.2500000e+00 0.00e+00 1.55e-16 -8.6 2.45e-08 - 1.00e+00 1.00e+00f 1\n", "\n", "Number of Iterations....: 6\n", "\n", " (scaled) (unscaled)\n", "Objective...............: 4.2500000000000000e+00 4.2500000000000000e+00\n", "Dual infeasibility......: 1.5452628521041094e-16 1.5452628521041094e-16\n", "Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00\n", "Complementarity.........: 2.5059105039901454e-09 2.5059105039901454e-09\n", "Overall NLP error.......: 2.5059105039901454e-09 2.5059105039901454e-09\n", "\n", "\n", "Number of objective function evaluations = 7\n", "Number of objective gradient evaluations = 7\n", "Number of equality constraint evaluations = 7\n", "Number of inequality constraint evaluations = 0\n", "Number of equality constraint Jacobian evaluations = 7\n", "Number of inequality constraint Jacobian evaluations = 0\n", "Number of Lagrangian Hessian evaluations = 6\n", "Total CPU secs in IPOPT (w/o function evaluations) = 0.002\n", "Total CPU secs in NLP function evaluations = 0.000\n", "\n", "EXIT: Optimal Solution Found.\n" ] } ], "source": [ "opt1 = pyo.SolverFactory('ipopt')\n", "status1 = opt1.solve(m, tee=True)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.4 Inspect the Solution](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.4-Inspect-the-Solution)", "section": "1.0.3.4 Inspect the Solution" } }, "source": [ "### 1.0.3.4 Inspect the Solution\n", "\n", "Now let's inspect the solution. We'll use the function `value()` to extract the numeric value from the Pyomo variable object." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.4 Inspect the Solution](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.4-Inspect-the-Solution)", "section": "1.0.3.4 Inspect the Solution" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x1 = 0.4999999999666826\n", "x2 = 0.5000000000333175\n", "x3 = -3.499999999966682\n", "\n", "\n" ] } ], "source": [ "## Return the solution\n", "print(\"x1 = \",pyo.value(m.x1))\n", "print(\"x2 = \",pyo.value(m.x2))\n", "print(\"x3 = \",pyo.value(m.x3))\n", "print(\"\\n\")" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.5 Visualize the Solution](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.5-Visualize-the-Solution)", "section": "1.0.3.5 Visualize the Solution" } }, "source": [ "### 1.0.3.5 Visualize the Solution" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.5 Visualize the Solution](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.5-Visualize-the-Solution)", "section": "1.0.3.5 Visualize the Solution" } }, "source": [ "Is our answer correct?\n", "\n", "We can solve this optimization problem with guess and check. If we guess $x_3$, we can then solve the constraints for $x_1$ and $x_2$:\n", "\n", "$$\\begin{bmatrix}\n", "1 & 1 \\\\\n", "1 & 2\n", "\\end{bmatrix} \\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix} = \\begin{bmatrix} 1 \\\\ 5 + x_3 \\end{bmatrix}$$\n", "\n", "We can then evaluate the objective. Let's see the graphical solution to our optimization problem." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.5 Visualize the Solution](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.5-Visualize-the-Solution)", "section": "1.0.3.5 Visualize the Solution" } }, "source": [ "
\n", "Activity: Verify you agree with how to translate the two linear constraints into a linear system of questions.
" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.5 Visualize the Solution](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.5-Visualize-the-Solution)", "section": "1.0.3.5 Visualize the Solution" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEOCAYAAABxdpuaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6dklEQVR4nO3dd3gU5drH8e9NCDU0QaMIElDwSA0EsSBIRJovBw6KghRFUPQIiFhBjoqKCtgRFQsaFSUKNuxYQlEpEqQjCogQmqCAhCIhud8/nk0hpCeb2U3uz3XtldnZ2Z1flrD3zsxTRFUxxhhj8qOM1wGMMcYEHysexhhj8s2KhzHGmHyz4mGMMSbfrHgYY4zJt7JeBygOtWrV0oiIiAI//+DBg1SuXLnoAvlRMGWF4MprWf0nmPIGU1YoXN74+Pg9qnpylg+qaom/RUVFaWHExcUV6vnFKZiyqgZXXsvqP8GUN5iyqhYuL7BUs/lctdNWxhhj8s2KhzHGmHyz4mGMMSbfPL9gLiKvAt2BP1S1qW/dO8DZvk2qA/tUNVJEIoB1wHrfY4tU9aaC7DcpKYmEhASOHDmS67bVqlVj3bp1BdlNsQumrBA4eStUqECdOnUIDQ31OooxQcHz4gHEAFOAN1JXqGqf1GUReQLYn2H7jaoaWdidJiQkUKVKFSIiIhCRHLc9cOAAVapUKewui0UwZYXAyKuq/PnnnyQkJFC/fn1PsxgTLDw/baWq84G/snpM3Kf6VcCMot7vkSNHqFmzZq6Fw5R8IkLNmjXzdBRqjHFEA2BUXd/pqE9ST1tlWN8eeFJVW2fYbg3wC/A38D9VXZDNaw4FhgKEh4dHxcbGHvd4tWrVOOuss/KULzk5mZCQkHz8Rt4JpqwQWHk3bNjA/v37s308MTGRsLCwYkxUcMGUFYIrbzBlTUmBQ4cKnjc6Ojo+9fM3s0A4bZWTqzn+qGMHcIaq/ikiUcCHItJEVf/O/ERVfQl4CaB169baoUOH4x5ft25dnk+XBMKplbwKpqwQWHkrVKhAy5Yts3187ty5ZP47ClTBlBWCK2+wZE1JgVatoF69bbzxxulUq1a0r+/5aavsiEhZ4HLgndR1qvqPqv7pW44HNgKNvElYeA8//DBNmjShefPmREZGsnjx4iJ9/QsvvBCAzZs38/bbb6etX7p0KbfcckuOz506dSpvvOEuQ8XExLB9+/Z873/37t2Ehoby4osv5vu5GcXExDB8+PB8P2/27NlMmDChUPs2Jlh99hmsWAGLF9fEHwdKgXzkcSnws6ompK4QkZOBv1Q1WUQaAA2BTV4FLIyFCxfyySefsGzZMsqXL8+ePXs4evRoke7jhx9+ANKLR79+/QBo3bo1rVtneSSa5qab0huxxcTE0LRpU2rXrp2v/c+cOZPzzz+fGTNmcOONN+YzfeH16NGDHj16FPt+jQkEzzzjfvbqtY2QkDOL/PU9P/IQkRnAQuBsEUkQkSG+h/py4oXy9sBKEVkBzAJuUtUsL7YHuh07dlCrVi3Kly8PQK1atdI+nOPj47n44ouJioqiS5cu7NixA4AOHTpw991306ZNGxo1asSCBe5yz5o1a2jTpg2RkZFccMEF/PrrrwBp5zlHjx7NggULiIyM5KmnnmLu3Ll0796dlJQUIiIi2LdvX1qus846i127djFu3Dgef/xxZs2axdKlS+nfvz+RkZF8+umn9OrVK237r776issvvzzL33HGjBk88cQTJCQksG3btrT1YWFhjB07lhYtWnDJJZewa9cuAD7++GPOO+88WrZsyaWXXpq2PtWBAweoX78+SUlJAPz9999ERESQlJTE5MmTady4Mc2bN6dv377A8UcsM2fOpGnTprRo0YL27dvn95/LmKCydi18/TVUqgSXXbbDL/vwvHio6tWqepqqhqpqHVWd5ls/SFWnZtr2PVVtoqotVLWVqn5cFBlEcr5VrVol122yuuWkc+fObN26lUaNGnHzzTczb948wPU/GTFiBLNmzSI+Pp7BgwczduzYtOcdO3aMJUuW8PTTT/PAAw8A7hTTyJEjWb58OfPmzaNOnTrH7WvChAm0a9eO5cuXM2rUqLT1ZcqUoWfPnnzwwQcALF68mIiICMLDw9O26d27N61bt+att95i+fLlXHbZZaxbt47du3cD8Nprr3Hddded8Ptt3bqVnTt30qZNG6666ireeSft7CMHDx7k/PPPZ8WKFbRt25aXX34ZgIsuuohFixbx008/0bdvXyZNmnTca1apUoUOHTrw6aefAhAbG8sVV1xBaGgoEyZM4KeffmLlypVMnXrcnw0ADz74IF9++SUrVqxg9uzZOf/jGBPkJk92P6+9FqpUOeaXfXhePEqrsLAw4uPjeemllzj55JPp06cPMTExrF+/ntWrV9OpUyciIyMZP348CQlpZ+7SvuVHRUWxefNmAC644AIeeeQRJk6cyJYtW6hYsWKec/Tp0yftgz02NpY+ffrkuL2IMHDgQKZPn86+fftYuHAh3bp1O2G72NhYrrrqKgD69u3LjBnpB5HlypWje/fuAERGRqb9HgkJCXTp0oVmzZrx2GOPsWbNmhNe9/rrr+e1114Dji9czZs3p3///kyfPp2yZU88G9u2bVsGDRrEyy+/THJycm5vizFB66+/wHe5khEj/LcfKx6Aas63v/8+kOs2Wd1yExISQocOHXjggQeYMmUK7733HqpKkyZNWL58OcuXL2fVqlXMmTMn7Tmpp7lCQkI4dsx9o+jXrx+zZ8+mYsWKXH755Xz77bd5/t0vuOACNmzYwO7du/nwww+zPQWV0XXXXcf06dOZMWMGV155ZZYf1jNmzCAmJoaIiAh69OjBihUr0k6nhYaGpvWvyfh7jBgxguHDh7Nq1SpefPHFLPtdtG3bls2bNzNv3jySk5Np2tS17v70008ZNmwY8fHxREVFpb1mqqlTpzJ+/Hi2bt1KZGQkf/75Z57fI2OCySuvwOHD0LkznHOO//ZjxcMj69evT/swBVi+fDn16tXj7LPPZvfu3SxcuBBwp7Gy+gae0aZNm2jQoAG33HIL3bp1Y+XKlcc9XqVKFQ4cOJDlc0WEXr16cdttt3HOOedQs2bNE7bJ/PzatWtTu3Ztxo8fz6BBg7L83Q4ePMi2bdvYvHkzmzdvZsyYMWTua5PZ/v37Of300wF4/fXXs93ummuu4eqrr0476khJSWHr1q1ER0czadIk9u3bR2Ji4nHP2bhxI+eddx4PPvggtWrVYuvWrTlmMSYYHTsGzz3nlkeO9O++rHh4JDExkWuvvTbtIu/atWsZN24c5cqVY9asWdx99920aNGCyMjItFZT2XnnnXdo2rQpkZGR/PLLL1xzzTXHPd68eXPKli1LixYteOqpp054fp8+fZg+fXq2p6wGDRrETTfdRGRkJIcPHwagf//+1K1bl8aNG5+w/YwZM467qA5wxRVXHHfqKivjxo3jyiuvpF27dtSqVSvb7fr378/evXu5+uqrAdfRcMCAATRr1oyWLVsyatQoqlevftxz7rzzTpo1a0bTpk1p3749LVq0yDGLMcHoo49gyxZo2BC6dvXzzrKb6KMk3bKaDGrt2rUnrMvO33//nedtvVZcWYcNG6avvPJKoV+nIHlnzpypAwYMKPS+M8vtbyKYJgEKpqyqwZU3kLO2a+dOmk+enL7OX5NBBXI/DxOgoqKiqFy5Mk888USx73vEiBF8/vnnfPbZZ8W+b2MC2U8/wYIFULUqZHE2uchZ8TD5Fh8f79m+n332Wc/2bUwgS+0UOHgwFMeIP3bNwxhjgtwff8CMGa5/WQFG8ikQKx7GGBPkXnwRjh6F7t3hzKIfiSRLVjyMMSaIHT0Kzz/vlv3dPDcjKx55IJs2wc03uytRZcq4nzffDBs3eh3NGFPKzZwJO3dC06ZwySXFt18rHrn5/HMqX3ih67Z54IDrOn7ggLvfvDl8/nmBXzokJITIyEiaNm3KlVdeyaFDh4oweP7t2rWL7t2706JFCxo3bsxll12W4/abN29O6+Gdk0ceeeS4+6lDxRtjCi91HKtbbsl9TL2iZMUjJxs3Qu/eyKFD4BvJNU1SEhw6BL17F/gIpGLFiixfvpzVq1dTrly5LAf0K0733XcfnTp1YsWKFaxdu7bI5sLIXDxy6/RojMmbRYtgyRI46STo3794923FIydPPHFi0cgsKQmy6LWdX+3atWPDhg3ce++9PJPa5g4YO3YskydPRlW58847adq0Kc2aNUsbzHDu3LlcfPHFXHXVVTRq1Ij777+ft956izZt2tCsWTM2+grb77//TseOHWnevDkdO3Zky5YtJ2TYsWPHcSPyNm/eHCDbfWeUecKm7t27M3fuXEaPHs3hw4eJjIykv++vO3WoeFXlf//7X5a/U4cOHejduzf/+te/6N+/PxoA0yUbE2hSPypuuMENv16ssus9WJJuBe5hXqWK5mkMxKpVc3+tLFSuXFlVVZOSkrRHjx76/PPP62+//aYtW7ZUVdXk5GRt0KCB7tmzR2fNmqWXXnqpHjt2THfu3Kl169bV7du3a1xcnFarVk23b9+uR44c0dNOO03vu+8+VVV9+umndeTIkaqq2r17d42JiVFV1WnTpmnPnj1PyPPFF19otWrVtEOHDjp+/Hjdtm2bqmq2+/7tt9+0SZMmqqr62muv6bBhw9Je6//+7//Seram/p6Zf+9Zs2ZpdHR0lr9T1apVdevWrZqcnKznn3++LliwoEDvcX5YD3PvBFPeQMmakKBatqxqSIjqli3Zb+evHuZ25JGTTIPrFXq7TFK/kbdu3ZozzjiDIUOGEBERQc2aNfnpp5+YM2cOLVu2pGbNmnz33XdcffXVhISEEB4ezsUXX8yPP/4IwLnnnstpp51G+fLlqV+/Pp07dwagWbNmacOdL1y4MG0mwYEDB/Ldd9+dkKdLly5s2rSJG264gZ9//pmWLVuye/fuHPddGN999x29e/fO8nXbtGlDnTp1KFOmzHHDthtjnBdecAMhXn451K1b/Pu3HuY5CQtzF8fzsl0BpF7zyOz6668nJiaGnTt3MnjwYIAcT9ukDtMOboKn1PtlypQ5YWjyVJLNlbWTTjqJfv360a9fP7p37878+fPzdMqobNmypKSkpN3Pajj1zPL6O2Uctt0YA0eOuL4d4C6Ue8GOPHIyYACEhua8TWgoDBxYpLvt1asXX3zxBT/++CNdunQBoH379rzzzjskJyeze/du5s+fT5s2bfL8mhdeeGHakOhvvfUWF1100QnbfPvtt2ktvg4cOMDGjRs544wz8rTviIgIli9fnjY8+pIlS9IeCw0NTZs6NqP27dvz3nvvFfh3Mqa0evtt2LMHWrWCtm29yWBHHjm5/XZ4/fWcL5qHhkKGqV2LQrly5YiOjqZ69eqEhIQArqAsXLiQFi1aICJMmjSJU089lZ9//jlPrzl58mQGDx7MY489xsknn5w2G19G8fHxDB8+PO0o4vrrr+fcc8+ldevWWe4746mktm3bUr9+/bRhz1u1apX22NChQ2nevDmtWrXirbfeSlvfq1cv5s2bV+DfyZjSSDX9QvnIkcXbPDdTEO8vaPv7Vqgh2T/7TFMqVVINDT3+InloqGqlSqqffZa318mH5ORkbdGihf7yyy/5fm4wDR+vGlh57YK5d4Ipr9dZ4+LcR1B4uOqRI3nZPq7A+yKQL5iLyKsi8oeIrM6wbpyIbBOR5b7bZRkeGyMiG0RkvYh08XvAbt04+MMPMHTo8T3Mhw6FlSshi/m7C2Pt2rWcddZZdOzYkYYNGxbpaxtjgl9qp8CbboIMlwaLXSCctooBpgBvZFr/lKo+nnGFiDQG+gJNgNrA1yLSSFWT/RlQGzSAKVPczc8aN27Mpk2b/L4fY0zw2bzZzRYYGuqKh5c8P/JQ1fnAX3ncvCcQq6r/qOpvwAagwFdYNQ+tiEzpYH8LJhhMmQIpKdCnD5x6qrdZAuHIIzvDReQaYClwu6ruBU4HFmXYJsG37gQiMhQYChAeHs7cuXOPezwsLIyEhASqVauWbbPVVMnJyRzIS5PdABBMWSEw8qoq+/fv5+DBgyf8nWSUmJiY4+OBJJiyQnDl9Srr4cMhTJ16AVCWiy6KZ+7cvP2/8VfeQC0eLwAPAer7+QQwGMjqUz7Lr4yq+hLwEkDr1q21Q4cOxz2elJREQkIC27ZtyzXMkSNHqFChQj7ieyeYskLg5K1QoQItWrQgNIem2anDpgSDYMoKwZXXq6wvvAAHD8KFF8KNN0bl+Xn+yhuQxUNVd6Uui8jLwCe+uwlAxr6UdYDtBdlHaGgo9evXz9O2c+fOpWXLlgXZTbELpqwQfHmN8UJKyvGj5wYCz695ZEVETstwtxeQ2hJrNtBXRMqLSH2gIbAk8/ONMaYk+eor+PlnOP10NxxJIPD8yENEZgAdgFoikgDcD3QQkUjcKanNwI0AqrpGRN4F1gLHgGH+bmlljDFeS+0UOGxY7oNeFBfPi4eqXp3F6mk5bP8w8LD/EhljTOBYv97NOVehght6PVAE5GkrY4wxTmr3sv79oVYtb7NkZMXDGGMC1P79EBPjlkeO9DTKCax4GGNMgJo61U0XFB0NzZp5neZ4VjyMMSYAHT4MTz7plu++29ssWbHiYYwxAWjaNPjjD4iKAt/koAHFiocxxgSYo0dh0iS3fM89Hs7ZkQMrHsYYE2CmT4etW6FxY/jPf7xOkzUrHsYYE0CSk2HCBLc8ZoybQigQBWgsY4wpnWbNgl9/hfr1oW9fr9Nkz4qHMcYECFV45BG3fPfdUNbzMUCyZ8XDGGMCxCefuNmta9eGQYO8TpMzKx7GGBMAVOFh36h9d9zh7fzkeWHFwxhjAkBcHCxeDDVrwtChXqfJnRUPY4wJAKlHHbfeCpUrexolT6x4GGOMxxYtgm+/hapVYfhwr9PkjRWPHKxdC336wC+/hHkdxRhTgqUedQwbBtWrexolz6x45GDaNHj3XXj99QivoxhjSqgVK1wrq4oV3SmrYGHFIwd33eX+QX/4oRbx8V6nMcaURI8+6n7ecAOccoq3WfLDikcOwsPh5pvd8rhxnkYxxpRAv/zizm6EhsKdd3qdJn+seOTirrugQoVkPvkEfvzR6zTGmJJkwgTXv+Paa6FOHa/T5I/nxUNEXhWRP0RkdYZ1j4nIzyKyUkQ+EJHqvvURInJYRJb7blP9ne+UU6Bnz20APPCAv/dmjCktfv8d3nzTDXwYiJM95cbz4gHEAF0zrfsKaKqqzYFfgDEZHtuoqpG+203FEbBPn61UqgSffgpLlhTHHo0xJd3jj8OxY65F51lneZ0m/zwvHqo6H/gr07o5qnrMd3cR4OkBXY0aSWltr+3ahzGmsHbtgldeccv33ONtloISVfU6AyISAXyiqk2zeOxj4B1Vne7bbg3uaORv4H+quiCb1xwKDAUIDw+Pio2NLXC+xMREkpNr0Lfv+Rw5EsJzzy2jceO/C/x6/pSYmEhYWPD0SwmmvJbVf4Ipb1FkffHFBsTGnkHbtnsYP3517k8ohMLkjY6OjlfV1lk+qKqe34AIYHUW68cCH5Be5MoDNX3LUcBWoGpurx8VFaWFERcXp6qqo0ergmrXroV6Ob9KzRosgimvZfWfYMpb2Kx//qkaFuY+SxYvLppMOSlMXmCpZvO56vlpq+yIyLVAd6C/75dAVf9R1T99y/HARqBRcWW6/XYIC4MvvoCFC4trr8aYkmTKFEhMhEsvhTZtvE5TcAFZPESkK3A30ENVD2VYf7KIhPiWGwANgU3FlatWLRgxwi3btQ9jTH4lJsIzz7jlsWO9zVJYnhcPEZkBLATOFpEEERkCTAGqAF9lapLbHlgpIiuAWcBNqvpXli/sJ7ffDlWqwJw58MMPxblnY0ywmzoV/voLLrwQLr7Y6zSF4/kkh6p6dRarp2Wz7XvAe/5NlLOaNeGWW9xAZuPGuSJijDG5OXIEnnjCLY8dCyLe5iksz488gtFtt7mjj6++gu+/9zqNMSYYvPYa7NwJkZHQrZvXaQrPikcBnHQSjBzplu+/39ssxpjAl5QEEye65XvuCf6jDrDiUWC33eYmbvnmG1iQZU8TY4xx3n7bDUdy9tlw+eVepykaVjwKqEaN9LH37ejDGJOd5OT0YdfHjIGQEG/zFBUrHoUwahRUq+Ymrp83z+s0xphA9P77sH491KsH/fp5naboWPEohOrV048+rN+HMSazY8fg3nvd8t13u3k7SgorHoV0663u6GPuXHczxphUr73mjjrOPBOGDPE6TdGy4lFI1au7i+fgrn0EwDiTxpgAcOhQ+hmJ8eOhXDlP4xQ5Kx5FYORIV0Tmz7ejD2OMM3kybN8OrVrBVVd5naboWfEoAtWq2dGHMSbdX3+5KWbB9e8oUwI/aUvgr+SNkSNd890FC+Dbb71OY4zx0qOPwv79buTcSy/1Oo1/WPEoIlWrukETwY4+jCnNtmyBZ591y6lHHyWRFY8iNGKEG7rk++9dz3NjTOlz//3wzz/Qty9ERXmdxn+seBShqlXhjjvcsh19GFP6rF4Nb7wBZcvCQw95nca/rHgUseHD3bDtP/zgRt01xpQe99wDKSlw441w1llep/EvKx5FrEoVO/owpjT67jv4+GOoXDm9V3lJZsXDD4YPd1PWLlpkk0UZUxqouuFHwH15DA/3Nk9xsOLhB2FhcOedbnnsWHcYa4wpuWbPdqeqTz45vdVlSWfFw0+GDYPatSE+HqZP9zqNMcZfjh1zQ62DO11VpYq3eYqLFQ8/qVwZHnnELY8ZAwcPepvHGOMfb7wB69ZB/fruQnlpYcXDjwYOdO28t2+Hxx7zOo0xpqgdPgz33eeWS+LghznxvHiIyKsi8oeIrM6w7iQR+UpEfvX9rJHhsTEiskFE1otIF29S502ZMvDUU2550iRISPA2jzGmaD37LGzbBpGRrlNgaeJ58QBigK6Z1o0GvlHVhsA3vvuISGOgL9DE95znRSSgJ3Vs1w5693bfUFLPixpjgt/evenTy5bUwQ9z4vmvq6rzgb8yre4JvO5bfh34T4b1sar6j6r+BmwA2hRHzsKYONEdzk6fDkuWeJ3GGFMUJkyAffugY0fo1MnrNMVPtAC92ESkEe7b/ymAAruB1ar6a4FCiEQAn6hqU9/9fapaPcPje1W1hohMARap6nTf+mnA56o6K4vXHAoMBQgPD4+KjY0tSDQAEhMTCQsLK/DzAV56qQEzZpxBkyb7efbZnxAp1MtlqyiyFqdgymtZ/SeY8iYmJnLoUE0GDDiPpKQyTJ0az9lnH/A6VrYK895GR0fHq2rrLB9U1TzdgHOAZ4DtQLLvluK7pd7fDjwNnJPX1/W9dgSu+KTe35fp8b2+n88BAzKsnwZckdvrR0VFaWHExcUV6vmqqvv3q55yiiqoxsYW+uWyVRRZi1Mw5bWs/hNMeePi4nTwYPd/+aqrvE6Tu8K8t8BSzeZzNdfTViJypojMAlYDQ4AVwAPANcBlwP/5lh/0PXY9sFpEZopIg7xWuEx2ichpvv2fBvzhW58A1M2wXR1cwQp4VaumD5R2993uGogxJvhs3lyJmBg3+OH48V6n8U5ernmsxR0ZDALCVbWbqj6oqm+p6heq+rlv+QFV7QaEA4OBBr7nFsRs4Frf8rXARxnW9xWR8iJSH2gIBM1VhCFDoHlz+P339FZYxpjg8sorDUhJgRtugIYNvU7jnbwUj6tUtbWqvqmquXZ1U9WDqvq6qkYBfXLbXkRmAAuBs0UkQUSGABOATiLyK9DJdx9VXQO8iytKXwDDVDU5D79DQAgJgSefdMuPPgo7d3qbxxiTP99/D99/X4tKldL7d5RWZXPbQFU/ym2bwjxXVa/O5qGO2Wz/MPBwQTN5rWNH+Pe/3eib//sfvPKK14mMMXmhCqNHu+Xbb4dTT/U2j9fy3VRXRK7yR5DS5PHH3fnSV1+F5cu9TmOMyYtPPnHDrlerdjRt2oXSrCD9PN4WkVI0gkvRa9TIDduuCqNG2ZwfxgS65OT0Tr4DBvxO1are5gkEBSkeMbie3WOzelBELhCRBYVKVQrcd5+b73zuXPiowCcGjTHF4c03Yc0aiIiAHj2CooGn3+W7eKjq9cAk4CERSWszJCJni8j7wHdAy6KLWDLVqAEPPOCW77gD/vnH2zzGmKzt359+1PHQQ1CunJ0qgAIOT6KqY4DbgFtEZLqITAVW4fp8vAiU8Nl7i8aNN8K//gUbN8KUKV6nMcZk5d57XcvICy6Afv28ThM4CjO21cvAt0A/4AZgJq5n+c2qao1Q8yA0FJ54wi0/9BDs3u1tHmPM8ZYtg+eec83sp04tfYMf5qQgra1CRWQksBG4BFiGG98qFNhatPFKvm7doHNnd2g8bpzXaYwxqZKT4aab3DTSt9ziOviadAWpo78CT+IGQ+yhbtCsfkAP4DMRCY7RzQKEiOs4WKaM+2azZo3XiYwxAC+/DD/+CKefnn590qQrSPEIwZ2maqGqnwKo6ru44nE+ECcitYouYsnXpIm7/pGS4jofGWO8tWtX+kXyp58uPfOS50dBikdDVX1VVVMyrlTVOcClQH3gh6IIV5o88ABUqwZffgmff+51GmNKtzvvdHN1dO0KV1zhdZrAVJCmukdyeGwx0A6oUJhQpdHJJ7tWHeCOPpKSvM1jTGk1d67r11GhgmsF6a+5d4JdkbcdUNV1QNuift3SYPhwOPNMWLcOXnzR6zTGlD5Hj8LNN7vle+5x/x9N1vIyn0eWAxTmRFW3+p57aUFClVbly8Njj7nl++93cyQbY4rPk0+6L28NG8Jdd3mdJrDl5cjjCxH5VkS6i0hIbhv7mvL2EpF5wGeFj1i6/Oc/0KED/PUXPPig12mMKT02b07/P/f88+7LnMleXopHJHAMNxHTDhF5S0RG+orJhSLSVkT+LSK3ici7wA7gPeCQ77kmH0TcRFEi7nyrjbprTPG45RY3w2ffvnCpnTPJVa7zeQBXAkOB04CbgZ7A1biOgRkJ8DfwPvCCqv5YhDlLlchIGDbMFY8hQ2DxYjeEuzHGPz76yM2xU7Vq+oRtJmd5+Ui6D/hFVd8GFopIb2AL0Bg4GVdEduPmOP8pcxNeUzCPPur+mJctc/N/pE5CY4wpWgcPuqMOcHOSn3aat3mCRV6Kx16gUob77wADVTXGL4kMAGFh8NJL0KWLG7akVy84+2yvUxlT8jz0EGzZAi1bpre0MrnLyzWP1cAQETndd99aPReTzp1h0CA3XPuQIa4HujGm6KxZ4wYnFXHDA4Xk2iTIpMpL8XgAaAVsEZGfcaepLhORbiIS7td0hiefdHMlf/+9awFijCkaqu5I49gxNzxQmzZeJwouuRYPVf0WaAY8AfyJO/LoB3wCbBeR7SLyiYg85GuiG1EUwXyTSy3PcPtbRG4VkXEisi3D+suKYn+BqkaN9KIxerRrTmiMKbw33oD58+GUU+CRR7xOE3zy1IZHVX8B7gIQkRTgTmATEIWbNbAVcBm+FlgisldVCzU4oqqux9fU19e/ZBvwAXAd8JSqPl6Y1w8mvXpB794wa5b7hvTFFzZkgjGF8ddfbgZPcA1SatTwNk8wKkgD0AeBBaq6BPdhDoCInIYrIlEUff+OjsBGVf1dSumn5pQp8O23MGcOvP66uxZijCmYe+6BPXtch9wBA7xOE5xENfDn4xWRV4FlqjpFRMYBg3B9SpYCt6vqCQN5iMhQXP8UwsPDo2JjYwu8/8TERMLCvJ+mZM6ccB599BzCwpKIifmRmjWPnrBNoGTNq2DKa1n9pzjzrl1bheHDW1GmjDJt2lLq1TuUr+eXpvc2Ojo63jdn04lUNaBvQDlgDxDuux+Om1OkDPAw8GpurxEVFaWFERcXV6jnF5WUFNVu3VRB9fLLs94mULLmVTDltaz+U1x5k5JUIyPd/6HRowv2GqXpvQWWajafq8EwI2833FHHLgBV3aWqyeo6I74MlJo2EqnNCcPC4P333TUQY0zePfecG/KnXr30KRBMwQRD8bgamJF6x3dtJVUvXD+UUuOMM2DSJLc8bBj8+ae3eYwJFtu3pxeMZ5+FSpVy3t7kLKCLh4hUAjrhxstKNUlEVonISiAaGOVJOA/deCO0bw9//AG33eZ1GmOCw6hRcOAA9OwJ//6312mCX0AXD1U9pKo1VXV/hnUDVbWZqjZX1R6qusPLjF4oUwZeecXNdPbGGzZtrTG5mTED3n3XHW0884zXaUqGgC4eJnsNG6bPPXDjjfD3397mMSZQ/f47/Pe/bvnpp931DlN4VjyC2KhR0Lo1bN1qo+4ak5XkZBg4EPbvdxOtXX+914lKDiseQaxsWZg2zf184QU31IIxJt3EibBggRtm/eWXbWSGomTFI8g1b+56y4L7VvXPP/ZPagzAkiVw//1u+fXXoVahBkwymdknTQkwdiw0aQK//gqvvRbhdRxjPJeYCP37uxFzb7sNOnXyOlHJY8WjBChXDl591bXCmjmzLkuXep3IGG/deits2OCOzG3EXP+w4lFCtGnj/sOkpAiDB8PRE4e9MqZUeO89dy2wQgV4+20oX97rRCWTFY8S5KGHoHbtw6xaBRMmeJ3GmOKXkAA33OCWH3/cnc41/mHFowSpVAnuuGM9AOPHw7JlHgcyphilpMA118DevXDZZTYfub9Z8ShhWrbcx7BhkJTkJpDae8Jg9caUTE88AXFxbmbAV1+1Zrn+ZsWjBHr8cWjVCn77zU0aFQRTthhTKMuWuVaHAK+9BuHh3uYpDax4lEAVKrjh2qtXh9mzXTExpqQ6dAj69XNH28OHu1NWxv+seJRQ9eu7jlEAY8a4XrbGlES33w7r10PjxunTFRj/s+JRgvXoAXfd5cb36dMHdu3yOpExRWv2bDdBWrlyrlluxYpeJyo9rHiUcA8/DO3awY4d7tA+OdnrRMYUjR07YMgQtzxhArRo4W2e0saKRwlXtizExroWKN9+C+PGeZ3ImMJLSXGNQfbscUOPjBzpdaLSx4pHKVC7tpsMp0wZ1//DJo8ywe7ZZ2HOHKhZE2Ji3N+2KV72lpcSl1ySPnnUgAGwZYu3eYwpqJUr3bU8cDNq1q7tbZ7SyopHKTJmDHTrBn/9BVddZeNfmeBz+LAbLffoURg61E3wZLxhxaMUKVMG3nwT6taFxYvhzju9TmRM/oweDatXQ6NG8OSTXqcp3QK6eIjIZhFZJSLLRWSpb91JIvKViPzq+1nD65zBpGZNmDkTQkNh8mS3bEwweOcd9zdbtqxrllu5steJSreALh4+0aoaqaqtffdHA9+oakPgG999kw/nnefGAQLX1PGXX7zNY0xuFi2Ca691y48/DlFR3uYxwVE8MusJ+PpO8zrwH++iBK/hw+HKK+HAATeA4qFDXicyJmubN0PPnvDPP3DjjXDLLV4nMgCiATxqnoj8BuwFFHhRVV8SkX2qWj3DNntV9YRTVyIyFBgKEB4eHhUbG1vgHImJiYSFhRX4+cUpP1kPHgzhv/+NYuvWSnTtuoO7717v53QnKqnvrdeCKStkn/fgwRBGjGjJb7+F0br1Xzz66CrKlvX2M6ukvLd5ER0dHZ/hrM/xVDVgb0Bt389TgBVAe2Bfpm325vY6UVFRWhhxcXGFen5xym/WlStVK1ZUBdVp0/yTKScl+b31UjBlVc06b1KSateu7m/znHNU9+4t9lhZKgnvbV4BSzWbz9WAPm2lqtt9P/8APgDaALtE5DQA388/vEsY/Jo1gxdecMvDhsGKFd7mMSbVqFHwxRdQqxZ88okbJdoEjoAtHiJSWUSqpC4DnYHVwGzAd+mMa4GPvElYclx7LVx/PRw54q5/7N/vdSJT2k2Z4m7lysGHH0KDBl4nMpkFbPEAwoHvRGQFsAT4VFW/ACYAnUTkV6CT774ppMmTITISNmyAwYNtAinjnc8/Tx+r6tVXoW1bb/OYrJX1OkB2VHUTcMI4mar6J9Cx+BOVbBUruj4fUVHw/vtuKJP77/c6lSltVq920wekpMC997re5CYwBfKRhylmZ53leqCXKeNG333uOa8TmdJk1y7o3t01H+/TBx54wOtEJidWPMxxevSAF190yyNGuNF4jfG3f/4pQ8+e8PvvcP75bh5yEa9TmZxY8TAnuP56mDjRXfe45hobwt34lypMnPgvFi+GevXcBXKbETDwWfEwWbrrLjdw4rFjcMUV8P33XicyJdW4cRAXdwpVqrgmueHhXicyeWHFw2Rr4kTX8urwYXcueuVKrxOZkuatt1zjjDJllHffhaZNvU5k8sqKh8mWiLv+0asX7NsHXbrApk1epzIlxfffuy8nAMOHb6BrV2/zmPyx4mFylDr8dXQ07Nzp5ovescPrVCbYbdrkJnI6etQN0tmr1zavI5l8suJhclWhgruIGRXl/tN36QJ793qdygSrffvcadA9e9zMlk895XUiUxBWPEyeVK3qWl2dfTasWgX//rcN427yLynJTYG8bp27vhEb645uTfCx4mHy7OSTYc4cqFPHna++8kr3YWBMXhw9Cn37wldfwSmnuJZVVat6ncoUlBUPky9nnOEKSM2a8NlnMGiQG0rCmJwcPuyucbz/PlSr5gpHvXpepzKFYcXD5Ns557hTWGFh7mL6yJE2kKLJ3oEDcNll7m+mVi2Ii4Nzz/U6lSksKx6mQM49111EL1fODZ394INeJzKBaO9e10Jv7lw47TSYPx9atvQ6lSkKVjxMgXXs6Ma+Sh1IccoUrxOZQPLHH66Jd+qwIwsWuKNWUzJY8TCFcvnl8NJLbnnECHcay5ht2+Dii93MlI0aucJx5plepzJFyYqHKbQhQ9xQJuBmJfzwQ0/jGI/99hu0awc//+ymOZ4/H+rW9TqVKWpWPEyRuOsud0sdSHHyZK8TGS+sX+8Kx2+/uetic+faQIcllRUPU2QmTHAXzlNSXAusESNcMTGlw8qV0L69O2XVrh18/TWcdJLXqYy/WPEwRUbETR361lvprbB69nRNNU3JtmQJdOjgLpJ37gxffGEdAEs6Kx6myPXrB998k96RsF07SEjwOpXxl/nz4dJLXbPcnj1h9myoVMnrVMbfArZ4iEhdEYkTkXUiskZERvrWjxORbSKy3He7zOus5kQXXQSLFrmWNitWwHnnwbJlXqcyRe3LL6FrV3d0efXVMHMmlC/vdSpTHAK2eADHgNtV9RzgfGCYiDT2PfaUqkb6bp95F9Hk5KyzYOFCdx58+3Z3BPLxx16nMkXlww/dnPeHD7upi998E0JDvU5likvAFg9V3aGqy3zLB4B1wOnepjL5ddJJbiysgQPdKLw9e8Izz9hwJsHu7behd2832OHIka6vT0iI16lMcQrY4pGRiEQALYHFvlXDRWSliLwqIjW8S2byonx5eP111xJLFW691VpiBStVePJJGDAAkpPhnnvcfBwiXiczxU00wL8CikgYMA94WFXfF5FwYA+gwEPAaao6OIvnDQWGAoSHh0fFxsYWOENiYiJhYWEFfn5xCvSs33xzChMn/oukpDKcd96f3HbbEk45paLXsfIk0N/bjPyRde/eUCZO/BeLF9cEYMiQTQwYsKVIXru0v7f+VJi80dHR8araOssHVTVgb0Ao8CVwWzaPRwCrc3udqKgoLYy4uLhCPb84BUPW775TrVlTFVQbNDigW7Z4nShvguG9TVXUWb/6SvXUU92/WY0aqh98UKQvX6rfW38rTF5gqWbzuRqwp61ERIBpwDpVfTLD+tMybNYLWF3c2UzhtG2b3hJr06Ywa4kVwJKSYPRo13dj507X+GHFCjc3hyndArZ4AG2BgcAlmZrlThKRVSKyEogGRnma0hRIakusFi32sWOHa4k1e7bXqUxGmza5JtcTJ7prGg8+CN9+a+NUGSdgZw9W1e+ArC7DWdPcEuKkk+Cxx1YwffrFvPGG+zY7frwbI8vmtfbWjBlw442u/0bduq511UUXeZ3KBJJAPvIwpUBoqBITAw895FryjB0LbdrA0qVeJyudEhPhuuvcKAEHDrhBLlessMJhTmTFw3hOBP73PzdNab168NNPrkf6qFHuw8wUj2XLoFUriImBChXgxRddj/Ea1hjeZMGKhwkYXbvCmjVw++3u/tNPQ+PG1ivd31RdX43zz4dff4WmTd2R39Ch1n/DZM+KhwkolSvD44/Djz9CVBRs3eqGwLjyStixw+t0Jc8ff0D37nDbba5l1c03uxFymzTxOpkJdFY8TEBq1co1533qKVdQZs1y819PnermCzGF9/XX0KKFG/m4Rg344AN47jmoGBx9No3HrHiYgFW2rBvKZO1a9+14/374739ds941a7xOF7x27HDDw1jfDVMYVjxMwDvjDNcHZOZMOPVU+OEHaNnSTTx15IjX6YLHjh2uEUKDBm6iLuu7YQrDiocJCiJuFNd16+Cmm9z5+fHjoXlziIvzOl1g27kzvWg8/bQruL16udZV995ro+GagrHiYYJK9erwwgvw3XeuJdavv8Ill7i+CX/+6XW6wLJzp7sQXr/+8UXjp5/g/ffd9Q5jCsqKhwlKbdu6D8Hx492Q7zEx7kNy2DBYudLrdN7auROef/5MGjRwDQ6OHHHXM1KLRmSk1wlNSWDFwwStcuVcj/RVq6BTJ9cj+vnn3TfqCy90c4gcPux1yuKza5frI9OgAcycWZfDh13RWLbMtaSyomGKkhUPE/QaNnSzFa5c6Y48qlZ1gy4OGgSnn+7O9//8s9cp/WfXLrjjDnfk9eSTrmC2bbsnrWi0bOl1QlMSWfEwJUazZq4V0fbtMG0anHsu7N3rzvefcw5cfLEb8O+ff7xOWjR27kwvGk884YpGz54QHw/jx6+2omH8yoqHKXEqV4bBg11P6fh4N8xG5cowf74b8K9OHTdy74YNXifNn6QkWLDAtZBq0wZq104vGj16uN/1ww9dB0tj/M2KhynRWrVyA/xt3+5aaUVGwp498Nhj7nRXp06u93pSUhHsbONGN75H1apQpoz7efPNbn0BbdjgruP85z9Qs6br0Dd+vBu+JTTUrV+6FD76yIqGKV42a4IpFapWdf1DbrzRffBOnQqxsW6Ijq+/dnOLtGrlTn01b+5+Nm6cj6E6Pv/cdURJSkqvRAcOwCuvuCv3s2ZBt265vsz+/a7T3pw57rZp0/GPn3OO6xneubM7DVe5cv7eB2OKihUPU6qIuFM+bdq4i8tvvumOTNasSS8kqcqUcVPlphYUqEm9em7Y+DIZj9k3bnSF49ChE3eYWkx693ZX9M8887iHk5PdkcOXX7pisWiRW5eqRg249FLo0sUdJZ1xRlG+G8YUnBUPU2pVr+7GeBo+HH7/3TX5XbXKfcavXAm//OJaaf38sxsaBZpx770QFuYKSmpR6fHlE5z+T1KO54CPHU7i80ueYtIZU0hMdAcliYnuSCPjECshIW7ipc6dXcGIirIe4CYwWfEwpZ4IRES427//nb7+yBE3HEpqQZk37y8SEk5i507XFHjhQrfdAKZThpwvmpTVJNpveZMeW6ac8FiDBq5QdO4M0dFQrVrR/W7G+IsVD2OyUaGC6yOR2uR17tyVdOjQgd27jz9KqfJq3qY7rCqJzJsLVaq4o5ewsPRlY4KNFQ9j8unkk914Wpdc4lsxM8ydh8qFVAmjfXv/ZjOmuARtU10R6Soi60Vkg4iM9jqPKcUGDHDtZnMSGgoDBxZPHmOKQVAWDxEJAZ4DugGNgatFpLG3qUypdfvteSseo0YVTx5jikFQFg+gDbBBVTep6lEgFujpcSZTWp15puvHUanSiUUkNNStnzXrhGa6xgQzUVWvM+SbiPQGuqrq9b77A4HzVHV4hm2GAkMBwsPDo2JjYwu8v8TERMKC5KpmMGWF4MqbW9YK27ZRZ+ZMTv3qK0IOHya5YkV2dupEwpVXcuT004sxaXC9rxBceYMpKxQub3R0dLyqts7yQVUNuhtwJfBKhvsDgWez2z4qKkoLIy4urlDPL07BlFU1uPJaVv8JprzBlFW1cHmBpZrN52qwnrZKADLOulwH2O5RFmOMKXWCtXj8CDQUkfoiUg7oC8z2OJMxxpQaQdnPQ1WPichw4EsgBHhVVdd4HMsYY0qNoLxgnl8ishv4vRAvUQvYU0Rx/C2YskJw5bWs/hNMeYMpKxQubz1VPTmrB0pF8SgsEVmq2bU4CDDBlBWCK69l9Z9gyhtMWcF/eYP1mocxxhgPWfEwxhiTb1Y88uYlrwPkQzBlheDKa1n9J5jyBlNW8FNeu+ZhjDEm3+zIwxhjTL5Z8TDGGJNvVjyyISJXisgaEUkRkdaZHhvjm0dkvYh08SpjdkQkUkQWichyEVkqIm28zpQTERnhey/XiMgkr/PkhYjcISIqIrW8zpIdEXlMRH4WkZUi8oGIVPc6U2bBNC+PiNQVkTgRWef7Wx3pdabciEiIiPwkIp8U9Wtb8cjeauByYH7Glb55Q/oCTYCuwPO++UUCySTgAVWNBO7z3Q9IIhKNG06/uao2AR73OFKuRKQu0AnY4nWWXHwFNFXV5sAvwBiP8xwnCOflOQbcrqrnAOcDwwI8L8BIYJ0/XtiKRzZUdZ2qrs/ioZ5ArKr+o6q/ARtw84sEEgWq+parEdiDRv4XmKCq/wCo6h8e58mLp4C7cO9zwFLVOap6zHd3EW4A0UASVPPyqOoOVV3mWz6A+1Au3rH280FE6gD/B7zij9e34pF/pwNbM9xPIPD+gG4FHhORrbhv8gH1jTOTRkA7EVksIvNE5FyvA+VERHoA21R1hddZ8mkw8LnXITIJhv9LWRKRCKAlsNjjKDl5GvclJ8UfLx6UAyMWFRH5Gjg1i4fGqupH2T0ti3XF/g00p+xAR2CUqr4nIlcB04BLizNfRrlkLQvUwJ0GOBd4V0QaqIdtyHPJew/QuXgTZS8vf8MiMhZ3yuWt4syWBwHxfym/RCQMeA+4VVX/9jpPVkSkO/CHqsaLSAd/7KNUFw9VLcgHakDMJZJTdhF5A3euE2Amfjpszatcsv4XeN9XLJaISApuILfdxZUvs+zyikgzoD6wQkTA/dsvE5E2qrqzGCOmye1vWESuBboDHb0syNkIiP9L+SEiobjC8Zaqvu91nhy0BXqIyGVABaCqiExX1QFFtQM7bZV/s4G+IlJeROoDDYElHmfKbDtwsW/5EuBXD7Pk5kNcRkSkEVCOAB2xVFVXqeopqhqhqhG4D79WXhWO3IhIV+BuoIeqHvI6TxaCal4ecd8YpgHrVPVJr/PkRFXHqGod399pX+DboiwcUMqPPHIiIr2AZ4GTgU9FZLmqdlHVNSLyLrAWdypgmKome5k1CzcAz4hIWeAIvrncA9SrwKsisho4ClwbgN+Qg9UUoDzwle9IaZGq3uRtpHRBOC9PW9yU16tEZLlv3T2q+pl3kbxjw5MYY4zJNzttZYwxJt+seBhjjMk3Kx7GGGPyzYqHMcaYfLPiYYwxJt+seBhjjMk3Kx7GGGPyzYqHMcaYfLPiYYyfiUhFEUkQkS0iUj7TY6+ISLKI9PUqnzEFYcXDGD9T1cPA/bhBAG9OXS8ijwJDgBGqGutRPGMKxIYnMaYY+GbNWwGcAjQArsdNKnW/qj7oZTZjCsKKhzHFxDfHwsfAN7iRhKeo6i3epjKmYOy0lTHFRFU/AZbhJut6h/Q5V9KIyPMislVE/haRbSLytG+4cmMCihUPY4qJb1bHSN/dA9kMPT8F+JeqVvVt2wI3J4cxAcWKhzHFQEQ6A28CHwCxwGAROSfzdqq6VlUPZlwFnFU8KY3JOysexviZiJwHvA98D/QH/gekAI9ms/1oETkA/AE0x01KZkxAsQvmxviR7+hiAbAF6KCqf/vWvwDcBFykqt/n8NyBwPOqmlBMkY3JEysexviJiJyBO9o4ClyoqrsyPHYasBH4SVXb5vAaVwFDVfVSf+c1Jj9sDnNj/ERVt+A6Bmb12A6gUh5eJgRoWJS5jCkKds3DmAAhImEicp2IVBenGXAv8KXX2YzJzE5bGRMgRKQy8CEQBZTDXTB/D9cL/ZCH0Yw5gRUPY4wx+WanrYwxxuSbFQ9jjDH5ZsXDGGNMvlnxMMYYk29WPIwxxuSbFQ9jjDH5ZsXDGGNMvv0/K3lxYfhkRtgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def constraints(x3):\n", " ''' Solve the linear constraints\n", " \n", " Args:\n", " x3: Value for the decision variable x3\n", " \n", " Returns:\n", " x1 and x2: Values calculated from the constraints\n", " '''\n", "\n", " # Define the matrices in the above equations\n", " A = np.array([[1, 1],[1, 2]])\n", " b = np.array([1, 5+x3])\n", "\n", " # Solve the linear system of equations\n", " z = np.linalg.solve(A,b)\n", " x1 = z[0]\n", " x2 = z[1]\n", " return x1, x2\n", "\n", "# Define a lambda function to plot the objective\n", "objective = lambda x1, x2, x3: x1**2 + 2*x2**2 - x3\n", "\n", "# Guess many values of x3.\n", "x3_guesses = np.linspace(-10,4,21)\n", "obj = []\n", "for x3 in x3_guesses:\n", " # Solve the constraints to determine x1 and x2\n", " x1, x2 = constraints(x3)\n", " # Calculate the store the objective function value\n", " obj.append(objective(x1,x2,x3))\n", "\n", "# Plot the objective function value versus x3\n", "plt.plot(x3_guesses, obj,color='blue',linewidth=2,label=\"Sensitivity Analysis\")\n", "plt.xlabel(\"$x_3$\",fontsize=18)\n", "plt.ylabel(\"$f(x)$\",fontsize=18)\n", "\n", "# Plot the solution from Pyomo\n", "x3_sln = pyo.value(m.x3)\n", "obj_sln = pyo.value(m.OBJ)\n", "plt.plot(x3_sln, obj_sln,marker='o',color='red',markersize=10,label=\"Pyomo Solution\",linestyle='')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[1.0.3.5 Visualize the Solution](https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html#1.0.3.5-Visualize-the-Solution)", "section": "1.0.3.5 Visualize the Solution" } }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "60268883", "metadata": {}, "source": [ "\n", "| [Contents](toc.html) | [Tag Index](tag_index.html) | [1.1 60 Minutes to Pyomo: An Energy Storage Model Predictive Control Example](https://ndcbe.github.io/CBE60499/01.01-Pyomo-Nuts-and-Bolts.html) >

\"Open

\"Download\"" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" } }, "nbformat": 4, "nbformat_minor": 2 }