{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuous Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This code cell installs packages on Colab\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()\n",
    "    helper.install_glpk()\n",
    "    helper.download_data(['student_diet.csv'])\n",
    "    helper.download_figures(['pack1.png','pack2.png','pack3.png'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import pyomo.environ as pyo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Programs: Student Diet Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reference: https://docs.mosek.com/modeling-cookbook/linear.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You want to save money eating while remaining healthy. A healthy diet requires at least P=6 units of protein, C=15 units of carbohydrates, F=5 units of fats and V=7 units of vitamins. Due to compounding factors (blizzard during Lent), our campus only has these options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>P</th>\n",
       "      <th>C</th>\n",
       "      <th>F</th>\n",
       "      <th>V</th>\n",
       "      <th>price</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>takeaway</th>\n",
       "      <td>3.0</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>vegtables</th>\n",
       "      <td>1.0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bread</th>\n",
       "      <td>0.5</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             P  C  F  V  price\n",
       "takeaway   3.0  3  2  1      5\n",
       "vegtables  1.0  2  0  4      1\n",
       "bread      0.5  4  1  0      2"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load data from file, use the first column (0, recall Python starts counting at 0) as the index\n",
    "food_options = pd.read_csv('https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/student_diet.csv',index_col=0)\n",
    "\n",
    "# Print up the the first 10 rows of data\n",
    "food_options.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's build a Python dictionary to store the nutrient requirements. (I strongly recommend not touching Python until we write the model on paper. I am including this here to avoid scrolling between the problem description and this cell.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment and fill in with all of the data\n",
    "# nutrient_requirements = {'P':6, 'C':15 }\n",
    "\n",
    "# Add your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Propose an Optimization Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sets**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Parameters**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Variables**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Constraints**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Degree of Freedom Analysis**\n",
    "\n",
    "We will later learn more about how to factor inequality constraints into degree of freedom analysis. For now, please count the number of equality and inequality constraints separately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve in Pyomo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our optimization model written on paper, we can proceed to solve in Pyomo. Before we start, let's review a few pandas tricks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['P', 'C', 'F', 'V', 'price']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Extract the column names, convert to a list\n",
    "food_options.columns.to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['P', 'C', 'F', 'V']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Same as above, but drop the last entry\n",
    "nutrients = food_options.columns.to_list()[0:4]\n",
    "nutrients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['takeaway', 'vegtables', 'bread']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Extract the index names, convert to a list\n",
    "foods = food_options.index.to_list()\n",
    "foods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('takeaway', 'P'): 3.0,\n",
       " ('takeaway', 'C'): 3.0,\n",
       " ('takeaway', 'F'): 2.0,\n",
       " ('takeaway', 'V'): 1.0,\n",
       " ('vegtables', 'P'): 1.0,\n",
       " ('vegtables', 'C'): 2.0,\n",
       " ('vegtables', 'F'): 0.0,\n",
       " ('vegtables', 'V'): 4.0,\n",
       " ('bread', 'P'): 0.5,\n",
       " ('bread', 'C'): 4.0,\n",
       " ('bread', 'F'): 1.0,\n",
       " ('bread', 'V'): 0.0}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a dictionary with keys such as ('takeaway', 'P')\n",
    "# Do not include 'price'\n",
    "food_info = food_options[nutrients].stack().to_dict()\n",
    "food_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'takeaway': 5, 'vegtables': 1, 'bread': 2}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create dictionary of only prices\n",
    "price = food_options['price'].to_dict()\n",
    "price"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's build our Pyomo model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 Set Declarations\n",
      "    FOOD : Size=1, Index=None, Ordered=Insertion\n",
      "        Key  : Dimen : Domain : Size : Members\n",
      "        None :     1 :    Any :    3 : {'takeaway', 'vegtables', 'bread'}\n",
      "    NUTRIENTS : Size=1, Index=None, Ordered=Insertion\n",
      "        Key  : Dimen : Domain : Size : Members\n",
      "        None :     1 :    Any :    4 : {'P', 'C', 'F', 'V'}\n",
      "    food_info_index : Size=1, Index=None, Ordered=True\n",
      "        Key  : Dimen : Domain         : Size : Members\n",
      "        None :     2 : FOOD*NUTRIENTS :   12 : {('takeaway', 'P'), ('takeaway', 'C'), ('takeaway', 'F'), ('takeaway', 'V'), ('vegtables', 'P'), ('vegtables', 'C'), ('vegtables', 'F'), ('vegtables', 'V'), ('bread', 'P'), ('bread', 'C'), ('bread', 'F'), ('bread', 'V')}\n",
      "\n",
      "3 Param Declarations\n",
      "    food_info : Size=12, Index=food_info_index, Domain=Any, Default=None, Mutable=False\n",
      "        Key                : Value\n",
      "            ('bread', 'C') :   4.0\n",
      "            ('bread', 'F') :   1.0\n",
      "            ('bread', 'P') :   0.5\n",
      "            ('bread', 'V') :   0.0\n",
      "         ('takeaway', 'C') :   3.0\n",
      "         ('takeaway', 'F') :   2.0\n",
      "         ('takeaway', 'P') :   3.0\n",
      "         ('takeaway', 'V') :   1.0\n",
      "        ('vegtables', 'C') :   2.0\n",
      "        ('vegtables', 'F') :   0.0\n",
      "        ('vegtables', 'P') :   1.0\n",
      "        ('vegtables', 'V') :   4.0\n",
      "    needs : Size=4, Index=NUTRIENTS, Domain=Any, Default=None, Mutable=False\n",
      "        Key : Value\n",
      "          C :    15\n",
      "          F :     5\n",
      "          P :     6\n",
      "          V :     7\n",
      "    price : Size=3, Index=FOOD, Domain=Any, Default=None, Mutable=False\n",
      "        Key       : Value\n",
      "            bread :     2\n",
      "         takeaway :     5\n",
      "        vegtables :     1\n",
      "\n",
      "1 Var Declarations\n",
      "    food_eaten : Size=3, Index=FOOD\n",
      "        Key       : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "            bread :     0 :   1.0 :  None : False : False : NonNegativeReals\n",
      "         takeaway :     0 :   1.0 :  None : False : False : NonNegativeReals\n",
      "        vegtables :     0 :   1.0 :  None : False : False : NonNegativeReals\n",
      "\n",
      "1 Objective Declarations\n",
      "    cost : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Sense    : Expression\n",
      "        None :   True : minimize : 5*food_eaten[takeaway] + food_eaten[vegtables] + 2*food_eaten[bread]\n",
      "\n",
      "1 Constraint Declarations\n",
      "    diet_min : Size=4, Index=NUTRIENTS, Active=True\n",
      "        Key : Lower : Body                                                                         : Upper : Active\n",
      "          C :  15.0 : 3.0*food_eaten[takeaway] + 2.0*food_eaten[vegtables] + 4.0*food_eaten[bread] :  +Inf :   True\n",
      "          F :   5.0 :                                 2.0*food_eaten[takeaway] + food_eaten[bread] :  +Inf :   True\n",
      "          P :   6.0 :     3.0*food_eaten[takeaway] + food_eaten[vegtables] + 0.5*food_eaten[bread] :  +Inf :   True\n",
      "          V :   7.0 :                             food_eaten[takeaway] + 4.0*food_eaten[vegtables] :  +Inf :   True\n",
      "\n",
      "9 Declarations: FOOD NUTRIENTS needs food_info_index food_info price food_eaten diet_min cost\n"
     ]
    }
   ],
   "source": [
    "# Add your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition note\"> \n",
    "<p class=\"title\"><b>Activity</b></p>\n",
    " Check the Pyomo model. Specifically, are the input (parameter) data correct? Do the equations match our model written on paper?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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",
      "\n",
      "This is Ipopt version 3.13.2, running with linear solver ma27.\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:        0\n",
      "Number of nonzeros in inequality constraint Jacobian.:       10\n",
      "Number of nonzeros in Lagrangian Hessian.............:        0\n",
      "\n",
      "Total number of variables............................:        3\n",
      "                     variables with only lower bounds:        3\n",
      "                variables with lower and upper bounds:        0\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:        0\n",
      "Total number of inequality constraints...............:        4\n",
      "        inequality constraints with only lower bounds:        4\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  8.0000000e+00 6.00e+00 1.10e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0\n",
      "   1  8.6444860e+00 5.02e+00 9.73e-01  -1.0 1.16e+00    -  2.18e-01 1.54e-01h  1\n",
      "   2  1.3016656e+01 0.00e+00 5.33e-01  -1.0 8.85e-01    -  6.53e-01 1.00e+00h  1\n",
      "   3  1.2884986e+01 0.00e+00 6.91e-02  -1.7 1.53e-01    -  7.53e-01 8.71e-01f  1\n",
      "   4  1.2512801e+01 0.00e+00 1.19e-01  -2.5 4.08e+00    -  1.17e-01 6.52e-01f  1\n",
      "   5  1.2513485e+01 0.00e+00 2.83e-08  -2.5 5.33e-02    -  1.00e+00 1.00e+00f  1\n",
      "   6  1.2500398e+01 0.00e+00 1.50e-09  -3.8 4.46e-02    -  1.00e+00 1.00e+00f  1\n",
      "   7  1.2500005e+01 0.00e+00 1.84e-11  -5.7 5.47e-04    -  1.00e+00 1.00e+00f  1\n",
      "   8  1.2500000e+01 0.00e+00 2.54e-14  -8.6 1.25e-05    -  1.00e+00 1.00e+00f  1\n",
      "\n",
      "Number of Iterations....: 8\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:   1.2499999882508366e+01    1.2499999882508366e+01\n",
      "Dual infeasibility......:   2.5375692660596042e-14    2.5375692660596042e-14\n",
      "Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Complementarity.........:   2.5136445446423303e-09    2.5136445446423303e-09\n",
      "Overall NLP error.......:   2.5136445446423303e-09    2.5136445446423303e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 9\n",
      "Number of objective gradient evaluations             = 9\n",
      "Number of equality constraint evaluations            = 0\n",
      "Number of inequality constraint evaluations          = 9\n",
      "Number of equality constraint Jacobian evaluations   = 0\n",
      "Number of inequality constraint Jacobian evaluations = 9\n",
      "Number of Lagrangian Hessian evaluations             = 8\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": [
    "# Specify the solver\n",
    "solver = pyo.SolverFactory('ipopt')\n",
    "\n",
    "# Solve\n",
    "results = solver.solve(m, tee=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition note\"> \n",
    "<p class=\"title\"><b>Activity</b></p>\n",
    " Does your degree of freedom analysis match Ipopt?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analyze Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's extract the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Units of takeaway eaten = 0.9999999904584187\n",
      "Units of vegtables eaten = 1.4999999892477958\n",
      "Units of bread eaten = 2.9999999704842386\n"
     ]
    }
   ],
   "source": [
    "# Add your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: After we discuss optimization theory, add discussion of shadow prices and multipliers here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear Programs: Circle Packing Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the smallest rectangle you can use to enclose three given circles? Reference: Example 4.4 in Biegler (2010).\n",
    "\n",
    "![picture](https://ndcbe.github.io/optimization/_images/pack1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Propose an Optimization Model\n",
    "\n",
    "The following optimization model is given in Biegler (2010):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![picture](https://ndcbe.github.io/optimization/_images/pack2.png)\n",
    "![picture](https://ndcbe.github.io/optimization/_images/pack3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition note\"> \n",
    "<p class=\"title\"><b>Activity</b></p>\n",
    " Identify the sets, parameters, variables, and constraints.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sets**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Parameters**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Variables**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Constraints**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition note\"> \n",
    "<p class=\"title\"><b>Activity</b></p>\n",
    " Perform degree of freedom analysis.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Degree of Freedom Analysis**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implement in Pyomo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will define functions to create and intialize the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "\n",
    "def create_circle_model(circle_radii):\n",
    "    ''' Create circle optimization model in Pyomo\n",
    "    \n",
    "    Arguments:\n",
    "        circle_radii: dictionary with keys=circle name and value=radius (float)\n",
    "        \n",
    "    Returns:\n",
    "        model: Pyomo model\n",
    "    '''\n",
    "\n",
    "    # Number of circles to consider\n",
    "    n = len(circle_radii)\n",
    "\n",
    "    # Create a concrete Pyomo model.\n",
    "    model = pyo.ConcreteModel()\n",
    "\n",
    "    # Initialize index for circles\n",
    "    model.CIRCLES = pyo.Set(initialize = circle_radii.keys())\n",
    "    \n",
    "    # Create parameter\n",
    "    model.R = pyo.Param(model.CIRCLES, domain=pyo.PositiveReals, initialize=circle_radii)\n",
    "\n",
    "    # Create variables for box\n",
    "    model.a = pyo.Var(domain=pyo.PositiveReals)\n",
    "    model.b = pyo.Var(domain=pyo.PositiveReals)\n",
    "\n",
    "    # Set objective\n",
    "    model.obj = pyo.Objective(expr=2*(model.a + model.b), sense = pyo.minimize)\n",
    "\n",
    "    # Create variables for circle centers\n",
    "    model.x = pyo.Var(model.CIRCLES, domain=pyo.PositiveReals)\n",
    "    model.y = pyo.Var(model.CIRCLES, domain=pyo.PositiveReals)\n",
    "\n",
    "    # \"In the box\" constraints\n",
    "    def left_x(m,c):\n",
    "        return m.x[c] >= model.R[c]\n",
    "    model.left_x_con = pyo.Constraint(model.CIRCLES, rule=left_x)\n",
    "\n",
    "    def left_y(m,c):\n",
    "        return m.y[c] >= model.R[c]\n",
    "    model.left_y_con = pyo.Constraint(model.CIRCLES, rule=left_y)\n",
    "\n",
    "    def right_x(m,c):\n",
    "        return m.x[c] <= m.b - model.R[c]\n",
    "    model.right_x_con = pyo.Constraint(model.CIRCLES, rule=right_x)\n",
    "\n",
    "    def right_y(m,c):\n",
    "        return m.y[c] <= m.a - model.R[c]\n",
    "    model.right_y_con = pyo.Constraint(model.CIRCLES, rule=right_y)\n",
    "\n",
    "    # No overlap constraints\n",
    "    def no_overlap(m,c1,c2):\n",
    "        if c1 < c2:\n",
    "            return (m.x[c1] - m.x[c2])**2 + (m.y[c1] - m.y[c2])**2 >= (model.R[c1] + model.R[c2])**2\n",
    "        else:\n",
    "            return pyo.Constraint.Skip\n",
    "    model.no_overlap_con = pyo.Constraint(model.CIRCLES, model.CIRCLES, rule=no_overlap)\n",
    "    \n",
    "    return model\n",
    "\n",
    "def initialize_circle_model(model, a_init=25, b_init=25):\n",
    "    ''' Initialize the x and y coordinates using uniform distribution\n",
    "    \n",
    "    Arguments:\n",
    "        a_init: initial value for a (default=25)\n",
    "        b_init: initial value for b (default=25)\n",
    "        \n",
    "    Returns:\n",
    "        Nothing. But per Pyomo scoping rules, the input argument `model`\n",
    "        can be modified in this function.\n",
    "    \n",
    "    '''\n",
    "    # Initialize \n",
    "    model.a = 25\n",
    "    model.b = 25\n",
    "\n",
    "    for i in model.CIRCLES:\n",
    "        # Adding circle radii ensures the remains in the >0, >0 quadrant\n",
    "        model.x[i] = random.uniform(0,10) + model.R[i]\n",
    "        model.y[i] = random.uniform(0,10) + model.R[i]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will create a dictionary containing the circle names and radii values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'A': 10.0, 'B': 5.0, 'C': 3.0}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create dictionary with circle data\n",
    "circle_data = {'A':10.0, 'B':5.0, 'C':3.0}\n",
    "circle_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['A', 'B', 'C'])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Access the keys\n",
    "circle_data.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 Set Declarations\n",
      "    CIRCLES : Size=1, Index=None, Ordered=Insertion\n",
      "        Key  : Dimen : Domain : Size : Members\n",
      "        None :     1 :    Any :    3 : {'A', 'B', 'C'}\n",
      "    no_overlap_con_index : Size=1, Index=None, Ordered=True\n",
      "        Key  : Dimen : Domain          : Size : Members\n",
      "        None :     2 : CIRCLES*CIRCLES :    9 : {('A', 'A'), ('A', 'B'), ('A', 'C'), ('B', 'A'), ('B', 'B'), ('B', 'C'), ('C', 'A'), ('C', 'B'), ('C', 'C')}\n",
      "\n",
      "1 Param Declarations\n",
      "    R : Size=3, Index=CIRCLES, Domain=PositiveReals, Default=None, Mutable=False\n",
      "        Key : Value\n",
      "          A :  10.0\n",
      "          B :   5.0\n",
      "          C :   3.0\n",
      "\n",
      "4 Var Declarations\n",
      "    a : Size=1, Index=None\n",
      "        Key  : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "        None :     0 :  None :  None : False :  True : PositiveReals\n",
      "    b : Size=1, Index=None\n",
      "        Key  : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "        None :     0 :  None :  None : False :  True : PositiveReals\n",
      "    x : Size=3, Index=CIRCLES\n",
      "        Key : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "          A :     0 :  None :  None : False :  True : PositiveReals\n",
      "          B :     0 :  None :  None : False :  True : PositiveReals\n",
      "          C :     0 :  None :  None : False :  True : PositiveReals\n",
      "    y : Size=3, Index=CIRCLES\n",
      "        Key : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "          A :     0 :  None :  None : False :  True : PositiveReals\n",
      "          B :     0 :  None :  None : False :  True : PositiveReals\n",
      "          C :     0 :  None :  None : False :  True : PositiveReals\n",
      "\n",
      "1 Objective Declarations\n",
      "    obj : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Sense    : Expression\n",
      "        None :   True : minimize : 2*(a + b)\n",
      "\n",
      "5 Constraint Declarations\n",
      "    left_x_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body : Upper : Active\n",
      "          A :  10.0 : x[A] :  +Inf :   True\n",
      "          B :   5.0 : x[B] :  +Inf :   True\n",
      "          C :   3.0 : x[C] :  +Inf :   True\n",
      "    left_y_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body : Upper : Active\n",
      "          A :  10.0 : y[A] :  +Inf :   True\n",
      "          B :   5.0 : y[B] :  +Inf :   True\n",
      "          C :   3.0 : y[C] :  +Inf :   True\n",
      "    no_overlap_con : Size=3, Index=no_overlap_con_index, Active=True\n",
      "        Key        : Lower : Body                                : Upper : Active\n",
      "        ('A', 'B') : 225.0 : (x[A] - x[B])**2 + (y[A] - y[B])**2 :  +Inf :   True\n",
      "        ('A', 'C') : 169.0 : (x[A] - x[C])**2 + (y[A] - y[C])**2 :  +Inf :   True\n",
      "        ('B', 'C') :  64.0 : (x[B] - x[C])**2 + (y[B] - y[C])**2 :  +Inf :   True\n",
      "    right_x_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body              : Upper : Active\n",
      "          A :  -Inf : x[A] - (b - 10.0) :   0.0 :   True\n",
      "          B :  -Inf :  x[B] - (b - 5.0) :   0.0 :   True\n",
      "          C :  -Inf :  x[C] - (b - 3.0) :   0.0 :   True\n",
      "    right_y_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body              : Upper : Active\n",
      "          A :  -Inf : y[A] - (a - 10.0) :   0.0 :   True\n",
      "          B :  -Inf :  y[B] - (a - 5.0) :   0.0 :   True\n",
      "          C :  -Inf :  y[C] - (a - 3.0) :   0.0 :   True\n",
      "\n",
      "13 Declarations: CIRCLES R a b obj x y left_x_con left_y_con right_x_con right_y_con no_overlap_con_index no_overlap_con\n"
     ]
    }
   ],
   "source": [
    "# Create model\n",
    "model = create_circle_model(circle_data)\n",
    "model.pprint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's initialize the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 Set Declarations\n",
      "    CIRCLES : Size=1, Index=None, Ordered=Insertion\n",
      "        Key  : Dimen : Domain : Size : Members\n",
      "        None :     1 :    Any :    3 : {'A', 'B', 'C'}\n",
      "    no_overlap_con_index : Size=1, Index=None, Ordered=True\n",
      "        Key  : Dimen : Domain          : Size : Members\n",
      "        None :     2 : CIRCLES*CIRCLES :    9 : {('A', 'A'), ('A', 'B'), ('A', 'C'), ('B', 'A'), ('B', 'B'), ('B', 'C'), ('C', 'A'), ('C', 'B'), ('C', 'C')}\n",
      "\n",
      "1 Param Declarations\n",
      "    R : Size=3, Index=CIRCLES, Domain=PositiveReals, Default=None, Mutable=False\n",
      "        Key : Value\n",
      "          A :  10.0\n",
      "          B :   5.0\n",
      "          C :   3.0\n",
      "\n",
      "4 Var Declarations\n",
      "    a : Size=1, Index=None\n",
      "        Key  : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "        None :     0 :    25 :  None : False : False : PositiveReals\n",
      "    b : Size=1, Index=None\n",
      "        Key  : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "        None :     0 :    25 :  None : False : False : PositiveReals\n",
      "    x : Size=3, Index=CIRCLES\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          A :     0 :  19.37245589877164 :  None : False : False : PositiveReals\n",
      "          B :     0 : 14.175030366901597 :  None : False : False : PositiveReals\n",
      "          C :     0 : 3.9363588537058454 :  None : False : False : PositiveReals\n",
      "    y : Size=3, Index=CIRCLES\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          A :     0 : 14.452708848019624 :  None : False : False : PositiveReals\n",
      "          B :     0 :  7.543623830819569 :  None : False : False : PositiveReals\n",
      "          C :     0 :  6.234252703186074 :  None : False : False : PositiveReals\n",
      "\n",
      "1 Objective Declarations\n",
      "    obj : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Sense    : Expression\n",
      "        None :   True : minimize : 2*(a + b)\n",
      "\n",
      "5 Constraint Declarations\n",
      "    left_x_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body : Upper : Active\n",
      "          A :  10.0 : x[A] :  +Inf :   True\n",
      "          B :   5.0 : x[B] :  +Inf :   True\n",
      "          C :   3.0 : x[C] :  +Inf :   True\n",
      "    left_y_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body : Upper : Active\n",
      "          A :  10.0 : y[A] :  +Inf :   True\n",
      "          B :   5.0 : y[B] :  +Inf :   True\n",
      "          C :   3.0 : y[C] :  +Inf :   True\n",
      "    no_overlap_con : Size=3, Index=no_overlap_con_index, Active=True\n",
      "        Key        : Lower : Body                                : Upper : Active\n",
      "        ('A', 'B') : 225.0 : (x[A] - x[B])**2 + (y[A] - y[B])**2 :  +Inf :   True\n",
      "        ('A', 'C') : 169.0 : (x[A] - x[C])**2 + (y[A] - y[C])**2 :  +Inf :   True\n",
      "        ('B', 'C') :  64.0 : (x[B] - x[C])**2 + (y[B] - y[C])**2 :  +Inf :   True\n",
      "    right_x_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body              : Upper : Active\n",
      "          A :  -Inf : x[A] - (b - 10.0) :   0.0 :   True\n",
      "          B :  -Inf :  x[B] - (b - 5.0) :   0.0 :   True\n",
      "          C :  -Inf :  x[C] - (b - 3.0) :   0.0 :   True\n",
      "    right_y_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body              : Upper : Active\n",
      "          A :  -Inf : y[A] - (a - 10.0) :   0.0 :   True\n",
      "          B :  -Inf :  y[B] - (a - 5.0) :   0.0 :   True\n",
      "          C :  -Inf :  y[C] - (a - 3.0) :   0.0 :   True\n",
      "\n",
      "13 Declarations: CIRCLES R a b obj x y left_x_con left_y_con right_x_con right_y_con no_overlap_con_index no_overlap_con\n"
     ]
    }
   ],
   "source": [
    "# Initialize model\n",
    "initialize_circle_model(model)\n",
    "model.pprint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition note\"> \n",
    "    <p class=\"title\"><b>Activity</b></p>\n",
    " Compare the initial values for <tt>x</tt> and <tt>y</tt> with and without initialization. What is the default initial value in Pyomo?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize Initial Point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll define a function to plot the solution (or initial point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAFlCAYAAADyArMXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAApx0lEQVR4nO3dWYxk130e8O9/19p7756FMxwNRQuinVi2CCGAgkCBacOhH2Q/OIgeHAUZgX6wAhvwQwy/WEhgwAhEOy+BEVojWAFsxwIsxwJMJBYFI0qAQDCHoCRKQ4kccjhbr9Pdtd/95KFqRsNh7333+/2AQXXfqe46VX3vV6fOKkopEBFR8WhZF4CIiE6GAU5EVFAMcCKigmKAExEVFAOciKigGOBERAVlpPlgi4uL6tKlS2k+JBFR4V27dm1LKbX0+PFUA/zSpUt49dVX03xIIqLCE5H39jrOJhQiooJigBMRFRQDnIiooBjgREQFxQAnIiooBjgRUUExwImICurQABeRCyLy9yJyXUS+LyK/OT3+BRG5KyKvT/89n3xxiYjogaNM5AkA/LZS6jURaQO4JiLfmP7fHymlvphc8YiIaD+HBrhSahXA6vTrvohcB3A+6YIREdHBjjWVXkQuAfgZAN8G8EkAnxeRfw3gVUxq6Tt7/MwLAF4AAF3XISKnLTMRUWE9+eSTuHnzZiy/S466J6aItAD8bwC/r5T6moisANgCoAD8RwBnlVL/9pDfobgHJxFVmYjguDkoIteUUs8+fvxIo1BExATwVwD+TCn1NQBQSq0rpUKlVATgTwB84lglIiKiUznKKBQBcBXAdaXUHz5y/Owjd/sVAG/EXzwiItrPUdrAPwng1wB8T0Renx77XQCfEZGPYdKEchPArydQPiIi2seR28BjeTC2gRNRxaXeBk5ERPmT6o48RJR/UaTghRH8MIIfKvhhBC+IfnwsUAiiCArApCKpoBQefi8CCCa3gDz83tS16T+BZWiwHnxvTI/pGocZHxMDnKiCgjDC0A0x9AKMvBBe8CCwJ6Ednx//LsePDrynCGDq8jDoLUND0zLQtHU0LAO6xnB/HAOcqOS8IMLICzBwJ2E9cAO4h4RpFpQCvEDBC0IA4fSoC2AS7jVTR2sa5k3bQNPSYejVbgVmgBOViB9G6DsBhm6AoTe59YLiDxxQChh7IcZeCMB7eLxmapMwnwZ6u2ZWqqbOACcquJEXYGfkY2foYeAGqNJAL8eP4Pge7g8moa4J0KmbmGtYmG2YqJl6xiVMFgOcqGCiSKE79rEz8rA79nPZHJKVSAG7Ix+7Ix8A0LD0SZg3TbRto3SdpAxwogJwgxC7o0lo98YBwqhC1exTGHkhRt4Yd3fHMHXBbMPEbMPCbN0sRfs5A5wop/wwwmbfxdbAxdAND/8BOpAfKmz2PWz2PYgA7ZqBpZaNhZZd2HZzBjhRznTHPjZ6DraHHljRToZSQG8coDcO8N72CIstGysdGw2rWJFYrNISlZQfRtgauFjvudORFpSWIFRY6zpY6zpo1wwsd2wsNm1oBaiVM8CJMtRzJrXt+wPWtvOg7wToOwFu6Q9q5TXUrfyOZGGAE6UsCCNsDTys9xyMWNvOJT9UWO06WO066NQNrHRqmG9YuauVM8CJUhKE0cNQ4CiS4pi0lQ9gGYLzsw0st/PTvMIAJ0pYGCmsdsdY7ToIYl1nhNLkBQrvbg1xrzvGE3N1LLXszMeVM8CJEhJFCut9B/d2x6WYzk4Trh/hxsYQ93YdPDFXx2LLzqwsDHCimCmlsNl3cXtnDC/gLMmyGnsh3lof4N7uGBfmGphrWqmXgQFOFBOlFLYGHu7sjA5dOpXKY+iGeHOtj3bNwIW5BmYaZmqPzQAnisH20MPt7RFHlVRY3wnwg9UeZuomLszX0a4lH+QMcKJTcIMQ724NsTP0sy4K5UR37KN718eZmRouzjcSnabPACc6ofWeg1vbI44soT2tdR3sjDxcXmxitpFM+zgDnOiYHD/Ejc0BeuMg66JQzrl+hOurfSy1bTy50IAZ8wqIDHCiI1JqMjvv9vaI097pWDb7LrpjD5cWmrH+XgY40RGMvADvbA7Rd1jrppPxAoUfrQ9gLl6EF0SwjNPXxhngRAeIIoW7u5MNAaq0VRklR6vP4Dt3dvHkQgPL7dqpfhcDnGgfQzfA2xsDDg2k2AWhwo2NIbb6Hp5absI2TrbiYfH3FCJKwGbfxRt3uwxvSlR37OONu130nJMNQ2WAEz1CKYWbW0O8vTFgRyWlwgsUfnCvh7Wuc+yfZRMK0ZQfRvjRep/DAyl1SgHvbg0x9AJ8aKF55OVqGeBEAAZugB+u9bn4FGVqo+di5Ib4iTOtI7WLswmFKm+j7+D7d7sMb8qFgRvgjbtddMeHt4szwKmylJos0H9jY8j2bsoVL1C4vtrDand84P3YhEKV5AWT9m5OzKG8Ugq4uTXC0N1/JBQDnCpn6AZ4k+3dVBCbfXff/2OAU6X0HR9vrvW5giCVAgOcKqM78vHD9T53hKfSYIBTJWwPPby13mdnJZUKA5xKb2vg4u2NARejotJhgFOpbfZd3NhkeFM5cRw4lRbDm8qOAU6ltDVgeFP5McCpdNjmTVXBAKdS2Rl6DG+qDAY4lcbIC/AWw5sqhAFOpeAFEd5c4yQdqhYGOBVeFCn8aL0P1+faJlQtDHAqvHe2BlxVkCqJAU6FdmdnhM2+l3UxiDLBAKfC2hq4uL198IL3RGXGAKdCGrgBbmwMsi4GUaYY4FQ4bhDih2s9rixIlccAp0IJIzXdPZ7pTcQAp0J57/7wwD0Ciark0AAXkQsi8vcicl1Evi8ivzk9Pi8i3xCRt6a3c8kXl6psd+Rhvbf//oBEVXOUGngA4LeVUh8F8E8A/IaIPAPgdwB8Uyn1NIBvTr8nSkQQRrixOcy6GES5cmiAK6VWlVKvTb/uA7gO4DyATwP4yvRuXwHwywmVkQg37w+5izzRY47VBi4ilwD8DIBvA1hRSq0Ck5AHsBx76YgA3B+4nKxDtIcjb6kmIi0AfwXgt5RSPRE56s+9AOCFkxWPqs4LIry7xaYTor0cqQYuIiYm4f1nSqmvTQ+vi8jZ6f+fBbCx188qpV5SSj2rlHo2jgJTtbyzNYAfcsgg0V6OMgpFAFwFcF0p9YeP/NfXAXx2+vVnAfxN/MWjKtvoOdgZ+lkXgyi3jtKE8kkAvwbgeyLy+vTY7wL4AwBfFZErAG4B+NVESkiV5Pghbt4fZV0Molw7NMCVUv8XwH4N3j8Xb3GIJm5sDrg5A9EhOBOTcmej76A35vreRIdhgFOuRJHiErFER8QAp1y51x1zwg7RETHAKTf8MMJq18m6GESFwQCn3LizM0bAMd9ER8YAp1xw/BDrPda+iY6DAU65cGt7BMXKN9GxMMApc33Hx/0BF6siOi4GOGXuPc64JDoRBjhlanvooe9w0g7RSTDAKTNKKdzaZu2b6KQY4JSZ7aGHsccNiolOigFOmVnjsEGiU2GAUyaGbsAFq4hOiQFOmWDtm+j0GOCUuiCMOO6bKAYMcErdRt/lZg1EMWCAU6qUUmw+IYoJA5xStTvy4fpc75soDgxwShXX+yaKDwOcUjP2QnTHftbFICoNBjilhm3fRPFigFMqokhha+BmXQyiUmGAUyp6js/t0ohixgCnVGwPOXGHKG4McErFzogBThQ3Bjglru/48AI2nxDFjQFOidsZcuggURIY4JS4bTafECWCAU6JGnshd90hSggDnBLF2jdRchjglKgdDh8kSgwDnBLjBRH6DrdNI0oKA5wSw7HfRMligFNielx5kChRDHBKzMBl8wlRkhjglIggjOBw5x2iRDHAKRFDl2O/iZLGAKdEDDw2nxAljQFOiRiy/ZsocQxwSgQ7MImSxwCn2HlBBJcdmESJY4BT7Nh8QpQOBjjFjs0nROlggFPshhyBQpQKBjjFjk0oROlggFOsgjDi/pdEKWGAU6z8kOFNlBYGOMXKCzh8kCgtDHCKlRcywInSwgCnWDHAidLDAKdY+WxCIUoNA5xixRo4UXoY4BQrdmISpefQABeRL4vIhoi88cixL4jIXRF5ffrv+WSLSUXBGjhReo5SA/9TAL+4x/E/Ukp9bPrv5XiLRUXFNnCi9Bwa4EqpbwHYTqEsVHB+GCHiPB6i1JymDfzzIvLdaRPLXGwlosLy2XxClKqTBvgfA3gKwMcArAJ4cb87isgLIvKqiLx6wseiguA0eqJ0nSjAlVLrSqlQKRUB+BMAnzjgvi8ppZ5VSj170kJSMSjFACdK04kCXETOPvLtrwB4Y7/7UnUwv4nSZRx2BxH5CwCfArAoIncA/B6AT4nIxwAoADcB/HpyRaSiYH5X1yuvvIKrV7+E9fV1rKys4MqVz+G5557Lulild2iAK6U+s8fhqwmUhQqOTSjV9Morr+DFF78Ix3EBAGtr63jxxS8CAEM8YZyJSbFhfFfT1atfehjeDziOi6tXv5RRiaqDAU6xYQW8mtbX1491nOLDAKfYiGRdAsrCysrKsY5TfBjgFBvmdzVdufI51Gr2+47VajauXPlcRiWqjkM7MYmOSlgFr6QHHZUchZI+BjjFhvFdXc899xwDOwNsQqHYsAJOlC4GOMWGTShE6WKAU2xMnQFOlCYGOMXG1Hk6EaWJVxzFxtQ1aKyEE6WGAU6xMg2eUkRp4dVGsbLYjEKUGl5tFCuLNXCi1PBqo1ixBk6UHl5tFCu2gROlh1cbxYo1cKL08GqjWDHAidLDq41ixU5MovTwaqNYcTo9UXoY4BQrQ9dgGQxxojQwwCl2TZvLzBOlgQFOsWtaDHCiNDDAKXYt1sCJUsEAp9ixCYUoHbzSKHaWocE2Nbh+lHVRYhFFCkGkEEQRglDBDyNESiFSH7yvJoCuaTA0gaELDE2DoQs07lZECWCAUyJatgHX97IuxrFESsHxQzh+NL0N4YURor2S+pg0TWDpGuqmDtvUUTM11EydwU6nwgCnRDRtA/cH+Q7wKFIYeAEGToCxH8INIkCdPqz3eywnmrwpPCQC29DQMHW0agaatsFAp2NhgFMiWjkdieIFEQZugL7jY+SFUAkF9pEoBdcP4fohdkYeRARN20DLNtCuGdyijg6Vz6uMCq9p61kX4aEoUtgd+9gdee+vAeeMUgoDx8fA8bHWBeqWgdmGiZmaCY171dEeGOCUCEPXUDM1OBl2ZDp+iJ2hh+7YR5RlTfuExl6AsRdgXXMwUzcx17BQM/PzxkjZY4BTYlq2ASeDjsze2Mf9oYexF6T+2EmIIoWdoYedoYe6ZWCxZaFdM7MuFuUAA5wSM1M3sZViR+bACbDRd3LdTHJaYy/A7e0AdcvAUtvmpKmK41+fEjPbsAAME3+coRtgo++WpsZ9FGMvwK37ARq2geW2jUZOO40pWfyrU2IsQ0O7ZqDvJBOsXhBhtTvG0K1OcD9u5Aa46QZo10ycmalx5ErFMMApUXNNK/YAV0phe+hhs+8WsnMyCX3Hx9ANsNypYb5pZV0cSgnfrilR8414w8TxQ9y8P8J6z2F4PyZSCmvdMW7eH8INytsPQD/GAKdE1S0ddev0Q9+UUtjsu3h3a1iptu6TGLkB3tkcYmvgZl0UShgDnBJ32lp4EEW4tT3CZt/JduZkgSilsNFzcGt7hDCGtVwonxjglLi55snHLDt+iHc3h5XuqDyNgePj3S02qZQVA5wS166ZJ9qtvjv2cXNrCD8sx7K0WfGCyZtgz/GzLgrFjAFOqZhrHK8Wvt5zcHdnxI7KmERK4c72CBs9J+uiUIwY4JSKow5tU0rh7s4Y99kBl4itgYvV3TH7EkqCAU6pmKmbMPSDV9SLlMKdnTG643yvI150OyMPdxnipcAAp1SICBZb9r7//yC8+2ynTUVv7OPODkO86BjglJozndqexx+E94Dhnaq+47MmXnAMcEpN3dIxu0dn5uquw/DOSG/sY7XLjs2iYoBTqh6vhW/2XbZ5Z2x35LHTuKAY4JSq2YYJ25ycdr2xj80+a395sN532f9QQAxwSpWI4EynhrEX4u7uOOvi0APT4Ztl3gyjjBjglLqZusHOsxyKlMLt7RECznwtDAY4pe7m/REaOdq1nn7MDyPcY6dmYTDAKVWr3TF64wAL3HQgtwaOj50RO5aLgAFOqRl7IW7dHwEAbENnLTzH1rsOvIBNKXl3aICLyJdFZENE3njk2LyIfENE3preziVbTCo6pRRubA7w6NLUce/WQ/GJlMK9LjuZ8+4oNfA/BfCLjx37HQDfVEo9DeCb0++J9nV3d/yBvTE7dROWyVp4Xo3cgOPDc+7QAFdKfQvA9mOHPw3gK9OvvwLgl+MtFpWJ44e4u7NXbU6w0t5/fRTK3mbf5aiUHDtpG/iKUmoVAKa3y/vdUUReEJFXReTVEz4WFdydnRH229WrXTPZFp5j0XQvUsqnxDsxlVIvKaWeVUo9m/RjUf4M3ACb/YNHNKy0917kivJhZ+xzS7acOmmAr4vIWQCY3m7EVyQqkwejTg5Stwx06iffN5MSphTWe6yF59FJA/zrAD47/fqzAP4mnuJQmeyOPHTHR1tfY7ljQw7e74EyNHB8biydQ0cZRvgXAP4fgI+IyB0RuQLgDwD8vIi8BeDnp98Tvc+t7cNr3w9Yuo45DivMtQ22heeOcdgdlFKf2ee/fi7mslCJdMc+hu7x2k2X2ja6Yx/hfj2elKmxF2Dshahb7HTOC87EpESsn2D3c13TsHDAtmuUvW1Osc8VBjjFzg1CbA9PdqEvNC0YOk/LvOqNfQQRx4XnBa8Uit1Gz8VJV4oVESx3WAvPK6UUdkfc+CEvGOAnEEWKa1nvI4oUNk65y85s3ULTPrR7hjKyM/R4/ucEr5I9+GGEoRtg6IWTWzdAEClEkXrfjEIRQBOBrgkalo6mZaBp62jaBmoVXeNjd+zDC05/cZ+dreHdzSE7NHPIDyMMvRAtvslmjn8BTNbq2Bq4GLohBm5w5GU0lQJCpRBGCl4Qve+jpaFPQr1lG5hvWmjXqjFR5aRt34+zdB0rMzXc23MNFcpa3/EZ4DlQ2b+AUgo7Ix/rPQfdsX/iNtv9BKFCbxygNw5wb9dB09ax3K5hsVXeTrpJ+2h8oxRm6xb644Cb7eZQ3wlwdibrUlDlAtwNQmz0XGz03VQXrB+6Id51h7i1PcJCy8JKp1a6GkzPCeCH8b4Tnp2pYeQFbErJmSCMOCY8B8qVIAdwgxDv3R9he+jFXts+jjBSkzeQnouWbeDiQgMzJVkHZCem5pNHGbqGszM13GFTSu70HZ8BnrFyfpZ/zFrXwXdud3F/kG14P27gBvjBvR7e2RyUYs3lpCZ5dOoWZhrleJMrkz7XRslcqWvgYy/Ejc3BB3aCyZv1noudkY/Li03MFXSzX8cP4frJvQmd6dQwdMNSvNGVheuHCKIIhlaJemAulfKVV0rh7u4Y372zm/vwfsALIry51sfbG334BQyppFeq0zUN52a5bnjeOF7xztUyKV2Au0GIN+72cOv+/rvA5Nlm38N3bu+iW7DZbsdduOokWrZZ2E8oZTX2udFDlkoV4CMvwBt3exgUvG3ODxWur/WwVaANZdN6zVc6tcpOksojhwGeqdIEeN/x8YN7vVSHBiZJKeCt9QHWuqeblp6WoZdOgGsiuDBfL+1Y+qJhDTxbpbgK+o6P66v92Mcg58G7W0Pc2833EDrHDxGk+Nqbuo4L8w3u4JMDQRixYzlDhQ/wgRvgzbV+qSd6vHd/dKL1tdOSxcfouqnj/Gw99celD3JL8qm3iAod4GMvxJurvVRrf1l5Z3OIzZxuaZVVs1WnbmGJS89mLihx5SnvChvgUaTw1kY5m032887mAGMvf22OXoYfoZdaNU7yyRibULJT2AC/szNOZehankQKuLE5yN1azFl3HJ+bqXNKd4ZYA89OIQO87/i41813x15S+k6Auznr1MyyBg5MdvG5MNeAaRTydC68Ik48K4vCnfFhpPD2xiBXa5qkbfLpIz9j3f0YNnA4LUPXcHG+AU3j0JS0sQaencKthXJrewQnwTU3ikAp4O2NAf7R+ZlcBFZeNrm1DR1PzNZxa2cEpJAp1669hpdf/lts72xjfm4ezz//S/j4x382+QfOmTKPAMu7QtXAuyO/MBNbkjbywtwssZqny7dVM/HEbB1I+H3t2rXX8NWv/iW2t7cBBWxvb+OrX/1LXLv2WrIPnENV/jSctUIF+O2dUdZFyJXV7jgX7Y95u4A7dWsyRjzBEH/55b+F571/+VzP8/Dyy3+b3IPmVs5OgAopTIAP3KAwKwumJVLARk7HhmdtJuEQ397ZPtZxoiQUJsDZdLK3ta6Tu2GFeZFkiM/PzR/reLll3w9TVYUIcD+McL9AK/OlyQui2HaCP6k8r0mSVIg///wvwbLev7StZVl4/vlfiveBCiDPf/+yK8QolPWeU8i1vdOy1nOw0MpuSrme8yt4pm5BF8HtnVFs7fUPRptwFMpkhUjKRu4DXCmF9R5r3wfpjQOMvAANK5s/p6lrAPI9K7ZVM3FxoYnb2yNEMdUGPv7xn61kYD/OyMFQ1qrKfRPKzsjPfKp2EWTZR2AVZAZk0zLw5HwDOgMnVobO1zMrub/yuuNibS2WlV6GI3SsAm2uULcMXFpsckOIGHFT4+zk/pUfcOjgkYy9MLMx4UWpgT9gGzqeWmqiYXMBrDiYrIEnaqm9f/9Wrq+8MFKpbdVVBlm92RXxAtY1DZcWmpjnJsmnxhp4MkSADy028eHl1r73yfUrP3CD3M3yy7OsNnO2C7vJsODMTB3n5+ocCncKpsEXL26WIXjmXAdnZmoH3i/Xo1CKvrt82rKaqdowdYjkb0r9Uc3ULViGjjs7I/jsMD8WTaRQfSBF0K4ZeHqlBds4vGKU61d+zOaTYxn72bxemiZoZjSEMS51U8flRbaLH5dt6hB+fInNcsfGM2c7RwpvIOcB7uVgneki8UOV2bT6ZgmCj+3ix1cvbPNZvmgCXF5q4qml1rGWiM53gOdgpb0iUSq716xlF7sG/mNsFz8ObmV3epYh+Oi5DlY6B7d37yXXV10elkp93CuvvIKrV7+E9fV1rKys4MqVz+G5557LulgP+aFCFlnaLE2AT8zULdiGjnu7Yzh+vmeZZok18NOZbZi4vNQ8cpPJ43J91eVtp49XXnkFL774RTjOZGr/2to6XnzxiwCQmxDP6jVrWDp0TXL3NzuNmqnj8lITmwMXW323sJ20SdE1rXBzAPLC1AUXFxpYbh+/1v0ovvrHcPXqlx6G9wOO4+Lq1S9lVKI9ZBQyIoLZhpnNgydKsNSq4fJSi80Fj2nVcl3/y61o3MU/fmL21OENMMCPZX19/VjHM5Fhu+1co7ydf7ah40OLTZyZqeViH9I8aDPAj8UyBB8504a/dSu2Ty65DvC8LVO5srJyrONZyDJb5hpmyTv+BPNNG5eXmqVr8z8uEUGr4ENH07TUtvHTT8zGPsIp1wGetynaV658DrXa+9clqNVsXLnyuYxK9EFmhpMqDF1Dp1bGZpT3s3QdTy40cW62XtmVDVu2wU8iR2CbGp4528GHl1uJLKCW67dQy9Dg+PkZifKgozLPo1CynhU337Qqs4LkbMNC09ax1nXRd6rxnB9g88nBRICVTg0XE16+WNKc+CEi6jiP99Z6H1uDbLcLKxJTFzx7Kds9Gb0gwmu3dio3YqPv+Njou3ArMORQRPATK+3Kfvo4zEzdxMWFxr5zI0Tk2BPuROSaUurZx4/n+m20xjGmx5KH18syNCw0rcq98bZrJto1A92xj82+W+pNSGbqJsN7D+2agQvzDczU02tGzHWA82Pa8eRlNuTKTK1yAT4hmKlbmKmb2B1NauRBDiejnRaXGni/pq3jwlwDcxm8Lvm44veRl0Aqiry84XVqJpq2jqFb/uaEvQlmG5Mg3xl52Bq4CMJytCnVLSMXn/TyoG7puDBXz3RD8Xxc8fswdA11S8fYq2oQHE+eJlac6dRwY3OYdTEyJTIZdjjbsLA99HB/4BZ+pipr35ORJU/M1bHUsjNfiTE/V/w+WrbBAD8Cy9BOvJ5CEhZbNm5tj+CXpOZ5GpoIFls25hrmJMiHHqICBvlkmGjuIyMxljEJ7uV29sH9QO7/Gp26gc2+e/gdK26mnq8/paYJzs7Wcev+KOui5IauaVhq16ZDLQPsjrxCLZS1mIMaZxY6dQMrnRrmG1buxr7n66rfw0LTxnv6qDRtiElZimFdhbid7dSw1nVKPSLjJHRNw3zTwnzTwtgLsD3y0Rt7uR56aRka5kq51s3eTF2w1Lax3K7leg2cUwW4iNwE0AcQAgj2Gqd4WromWGrZWO06cf/q0mhYeqpDl45K0wQX5uu4sVHttvCD1C0D5y0DZzo2umMf2yMfXg5r5cvtWiVq3+3apLa90MxfbXsvcdTA/7lSaiuG37OvMzM1BvgBDtv4NEtLLRtrXafCI1KOZlIrtzHftDH0Js0rvbGfi1p53TLQyWEFIS6GPqkkLndsNAq2vkshSlszdcw1TewMqzVd+SgMfdJBllcigovzDVxf7WddlMJoWgaaloEznQi7Ix+9sY9xhrXylU5+z6+TEpnUtpfaNhabdiFq23s5bYArAH8nIgrAf1VKvfT4HUTkBQAvnPJxcKZTY4DvYall535W3GzD4hvwCeiahoWWjYWWjSCMMHAD9J0AQy9IbRRLp24Wrla6H1OfrFk/27AwWzcTWVwqbadaC0VEziml7onIMoBvAPh3SqlvHXD/Y62F8rjXb+9ySOEjRICPXZgtxMQKL4jw3Tu7HFYYA6UUhl6AgROg7wbwE+okNnQNl5eaMLTiBl3T1jHXsDDbMNGyjVy04+dmLRSl1L3p7YaI/DWATwDYN8BP69xsjR1ij1hoWoUIb2AyiuHSYhNvrQ+yLkrhiQhatomWbeIMADcI0XcC9J1pU0tM75FnZ2qFC29NgJmG+TC08zQ3IgknDnARaQLQlFL96de/AOA/xFayPSy3a9geevwojsnuHpcWm1kX41gWWzZ2hl5F10lJjm3osFs6Fls2wijC2Asx9kM4foSxH55oPZaZhoV2AdZ2r1s6mpaOhm2gZRlo16q1TvlpauArAP56+pHEAPDnSqn/GUupDnB5sYXvOvwofnmxlenmDSd1abGJnuPDC6r990uKrmlo1TS0HgnfMIoeBrrjh3D88MCx+aau4UwnXyObRIC6qaNpG2ja01vLyH3/T9JyvR74frYGbqU/ii93bDy11Mq6GCe2O/Lw5lo/F0PkqiqMIjhBBNcPMfYjhFEEP1QIowjnZxuZbBknMulotHQdpiGwdA1N20DD0tG0ylOzzk0beFaq/FHcNjVcWihW08njZhsWLs438B6n2WdG1zQ0LQ3Nx0aYXF5qYrltwwsnge4HEbwwghdE8B8cCyMEkYJSCgqYvhErKDVpfhdgujeqQOTB9wJDE1iGBlPXpkGtPfL95FgeOhmLpJABDlT3o/hTS61SfGw8N1vH2A+x0eM6N3lxdqaGlWnTiW3osA0A5RsCXirFa0SdMnUNH15ql3wX9Pd7Yq6eyynzJ3V5sYlOzhbhqqrZhoknFxpZF4OOqbABDkyGCz293KpEiK90bFyYL9cF9mBvxZpZ6NOw8BqWPr2OKnAhlUzhr5yF1qRDr8zn3lLbxuUCd1oexNQ1fPRsB5ZR4j9gjtmmho+caZdiVmIVleKvNgm4ZilDfDLipNidloepmTqeOTvDEE+ZZWh45mynMJPB6INKEeDAZJLPR1bapejge+D8bH366aI8z2k/dYshnibb1PCT5xjeRVeaAAeAuaaFj55tw9SLHQIiwKXFBi5WrFOpbun4yXMzsNkmnqjJ68zwLoPSXSntmomfvjCLxVYxN19t2jp+6vwMzs7Usy5KJmrmJFwaOd4Fpciato5nznZKv0ZIVRRyJuZR7Qw9vLM1LMSWXpoA5+fqOD9br0STyWHCSOHG5gD3KzhZKymLLQuXSzKPoMjinIlZ6gAHgCCMcGt7hPUcTxhp1ww8tdTK9d57WbmzM8KdnTGn3Z+CCHBxvoFzs9X8VJc3DPAT6I59vLs1zNV64oYuuDDXwEqnmrt9H9XO0MPbmwNubH0Cpi54ermNmQptSJx3DPATUkqhO/ax3nOxM8puF/CWbWClM9lphR9nj8bxQ/xovc+9NY+hXTPw4eUWOytzhgEeAzeYrMOx0XdTaSPXNcFCy8JKp4ZWBiu9lYFSCve6Du5sj5DSjmKFpGuCJ+bqODtTjZ3ki4YBHiOlFHZHPrYGLvpuANePL8xNXdCwDCy0LCyyth2bsRfixuYAfSfIuii506lP+lNY684vBniC/DDCyA0x9AIM3QBDLzxSu7llyMNF5h+sYcyLKFlrXQe3tkcIWR2HoQuenG9gOWcbMdAHMcBTFkaTNZCVAiKlECkFEYEmgCYCXZNC7o5TBl4Q4e7uGBs9p5LNKpoAK50azs3WYRk8B4ug8hs6pE3XBLrG2nQeWYaGDy02cXamhjs7Y2wN3EoMORQBlts2zs/VOSmnwhjgVAo1U8eHl1s4P1vHnZ0R7g+zG2WUJJHJhJwn5hpsoiMGOJVL3dLx9EobF/wQm30XG32nFLs2WYaG5baN5Y7NGjc9xACnUqqZOi7MN/DEXB3bQw/rPRfdsZ91sY5FBJipm1jp1DDXMDkkkD6AAU6lJiJYaE0mTTl+iO2hh52Rh74T5LKJRWQyAWeuYWG+abGZhA7EAKfKqJk6zs3WcW62jiCMsDv2sTvysDvy4Wc4Td/UBbMNC7MNE7N1k7vj0JExwKmSDF3DYsvGYsuGUgqOH2HgBhh5wfQ2TGTtFUOX6VwB/eG8gZqpsXmEToQBTpUnIqhb+nQ1SPvhcccP4fghvDCCHyr4QQQvjOAFEYJIQSkFBUCpSdOHTH+XoQlsQ4OpazANDaYusHQNNZOTuyheDHCifTBwKe/Y2EZEVFAMcCKigmKAExEVFAOciKigGOBERAXFACciKigGOBFRQTHAiYgKigFORFRQDHAiooJigBMRFRQDnIiooBjgREQFxQAnIiooBjgRUUExwImICooBTkRUUAxwIqKCYoATERUUA5yIqKAY4EREBcUAJyIqKAY4EVFBMcCJiAqKAU5EVFAMcCKigmKAExEV1KkCXER+UUR+KCJvi8jvxFUoIiI63IkDXER0AP8FwL8A8AyAz4jIM3EVjIiIDnaaGvgnALytlHpHKeUB+O8APh1PsYiI6DCnCfDzAG4/8v2d6TEiIkqBcYqflT2OqQ/cSeQFAC9Mv3VF5I1TPGbRLQLYyroQGav6a8DnX+3nDwCLInLc1+DJvQ6eJsDvALjwyPdPALj3+J2UUi8BeAkARORVpdSzp3jMQqv68wf4GvD5V/v5A/G+BqdpQvkHAE+LyIdExALwrwB8PY5CERHR4U5cA1dKBSLyeQD/C4AO4MtKqe/HVjIiIjrQaZpQoJR6GcDLx/iRl07zeCVQ9ecP8DXg86fYXgNR6gP9jkREVACcSk9EVFCpBDin3AMiclNEvicir4vIq1mXJ2ki8mUR2Xh02KiIzIvIN0TkrentXJZlTNo+r8EXROTu9Dx4XUSez7KMSRKRCyLy9yJyXUS+LyK/OT1eifPggOcf2zmQeBPKdMr9jwD8PCZDD/8BwGeUUj9I9IFzRkRuAnhWKVWJMbAi8s8ADAD8N6XUT02P/ScA20qpP5i+kc8ppf59luVM0j6vwRcADJRSX8yybGkQkbMAziqlXhORNoBrAH4ZwL9BBc6DA57/v0RM50AaNXBOua8gpdS3AGw/dvjTAL4y/formJzMpbXPa1AZSqlVpdRr06/7AK5jMlu7EufBAc8/NmkEOKfcTygAfyci16azU6toRSm1CkxObgDLGZcnK58Xke9Om1hK2XzwOBG5BOBnAHwbFTwPHnv+QEznQBoBfqQp9xXwSaXUz2KyeuNvTD9eU/X8MYCnAHwMwCqAFzMtTQpEpAXgrwD8llKql3V50rbH84/tHEgjwI805b7slFL3prcbAP4ak6alqlmftgs+aB/cyLg8qVNKrSulQqVUBOBPUPLzQERMTMLrz5RSX5sersx5sNfzj/McSCPAKz/lXkSa004MiEgTwC8AqOKiXl8H8Nnp158F8DcZliUTD4Jr6ldQ4vNARATAVQDXlVJ/+Mh/VeI82O/5x3kOpDKRZzpM5j/jx1Pufz/xB80REbmMSa0bmMx+/fOyvwYi8hcAPoXJ6nPrAH4PwP8A8FUAFwHcAvCrSqnSdvLt8xp8CpOPzgrATQC//qA9uGxE5J8C+D8Avgcgmh7+XUzagUt/Hhzw/D+DmM4BzsQkIioozsQkIiooBjgRUUExwImICooBTkRUUAxwIqKCYoATERUUA5yIqKAY4EREBfX/AaWotS6ml5KmAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot initial point\n",
    "\n",
    "def plot_circles(m):\n",
    "    ''' Plot circles using data in Pyomo model\n",
    "    \n",
    "    Arguments:\n",
    "        m: Pyomo concrete model\n",
    "    \n",
    "    Returns:\n",
    "        Nothing (but makes a figure)\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    # Create figure\n",
    "    fig, ax = plt.subplots(1,figsize=(6,6))\n",
    "    \n",
    "    # Adjust axes\n",
    "    l = max(m.a.value,m.b.value) + 1\n",
    "    ax.set_xlim(0,l)\n",
    "    ax.set_ylim(0,l)\n",
    "    \n",
    "    # Draw box\n",
    "    art = mpatches.Rectangle((0,0), width=m.b.value, height=m.a.value,fill=False)\n",
    "    ax.add_patch(art)\n",
    "\n",
    "    # Draw circles and mark center\n",
    "    for i in m.CIRCLES:\n",
    "        art2 = mpatches.Circle( (m.x[i].value,m.y[i].value), radius=m.R[i],fill=True,alpha=0.25)\n",
    "        ax.add_patch(art2)\n",
    "        \n",
    "        plt.scatter(m.x[i].value,m.y[i].value,color='black')\n",
    "    \n",
    "    # Show plot\n",
    "    plt.show()\n",
    "    \n",
    "plot_circles(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve and Inspect the Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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",
      "\n",
      "This is Ipopt version 3.13.2, running with linear solver ma27.\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:        0\n",
      "Number of nonzeros in inequality constraint Jacobian.:       30\n",
      "Number of nonzeros in Lagrangian Hessian.............:       12\n",
      "\n",
      "Total number of variables............................:        8\n",
      "                     variables with only lower bounds:        8\n",
      "                variables with lower and upper bounds:        0\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:        0\n",
      "Total number of inequality constraints...............:       15\n",
      "        inequality constraints with only lower bounds:        9\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        6\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0  1.0000000e+02 1.50e+02 1.01e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0\n",
      "   1  1.0313665e+02 8.56e+01 1.10e+00  -1.0 1.11e+01    -  4.84e-01 3.16e-01h  1\n",
      "   2  1.0122024e+02 4.44e+01 1.04e+00  -1.0 4.29e+01    -  2.73e-01 1.46e-01h  1\n",
      "   3  9.9356199e+01 3.24e-01 4.47e-01  -1.0 1.00e+01    -  9.36e-01 8.45e-01h  1\n",
      "   4  9.8949272e+01 0.00e+00 1.23e-01  -1.0 1.61e+01    -  8.63e-01 1.00e+00h  1\n",
      "   5  9.8907634e+01 0.00e+00 3.14e-02  -1.0 5.93e+01    -  1.00e+00 1.00e+00h  1\n",
      "   6  9.8897118e+01 0.00e+00 1.01e-02  -1.0 3.46e+01    -  1.00e+00 1.00e+00h  1\n",
      "   7  9.8899831e+01 0.00e+00 3.28e-03  -1.0 3.59e+01    -  1.00e+00 1.00e+00h  1\n",
      "   8  9.8899882e+01 0.00e+00 1.40e-04  -1.0 1.53e+01    -  1.00e+00 1.00e+00h  1\n",
      "   9  9.8394322e+01 0.00e+00 2.11e-04  -1.7 1.24e+00    -  1.00e+00 1.00e+00h  1\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "  10  9.8284887e+01 0.00e+00 4.44e-05  -3.8 2.68e-01    -  1.00e+00 9.98e-01f  1\n",
      "  11  9.8284281e+01 0.00e+00 6.66e-09  -5.7 7.01e-03    -  1.00e+00 1.00e+00h  1\n",
      "  12  9.8284270e+01 0.00e+00 3.24e-12  -8.6 5.49e-04    -  1.00e+00 1.00e+00h  1\n",
      "\n",
      "Number of Iterations....: 12\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:   9.8284270438748507e+01    9.8284270438748507e+01\n",
      "Dual infeasibility......:   3.2404583559775637e-12    3.2404583559775637e-12\n",
      "Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Complementarity.........:   2.5221946198982285e-09    2.5221946198982285e-09\n",
      "Overall NLP error.......:   2.5221946198982285e-09    2.5221946198982285e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 13\n",
      "Number of objective gradient evaluations             = 13\n",
      "Number of equality constraint evaluations            = 0\n",
      "Number of inequality constraint evaluations          = 13\n",
      "Number of equality constraint Jacobian evaluations   = 0\n",
      "Number of inequality constraint Jacobian evaluations = 13\n",
      "Number of Lagrangian Hessian evaluations             = 12\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      0.003\n",
      "Total CPU secs in NLP function evaluations           =      0.000\n",
      "\n",
      "EXIT: Optimal Solution Found.\n"
     ]
    }
   ],
   "source": [
    "# Specify the solver\n",
    "solver = pyo.SolverFactory('ipopt')\n",
    "\n",
    "# Solve the model\n",
    "results = solver.solve(model, tee = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we can inspect the solution. Because Pyomo is a Python extension, we can use Pyoth (for loops, etc.) to programmatically inspect the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Name\tValue\n",
      "a \t 29.14213541618474\n",
      "b \t 19.99999980318951\n",
      "x[A] \t 9.99999990125201\n",
      "x[B] \t 14.999999849644157\n",
      "x[C] \t 4.723219051632863\n",
      "y[A] \t 19.1421355149319\n",
      "y[B] \t 4.999999951252102\n",
      "y[C] \t 5.4079499719072\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAFoCAYAAACotWuNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAArJklEQVR4nO3daYwk530e8OdfVV1999zX3uTyEkXJVLBQAisIZJkxFCOAbAMytB8MBliB/mABMuAPFvzFcoAAQiDayYfACM0VzAC2HAGyLEEw4mgFG4qCQDZXoEVSy0urJbm7c199V3VV//Ohe5bL5cxO90x3nc8PWs1Ms3fq7a2up996T1FVEBFRMhhhF4CIiEaHoU5ElCAMdSKiBGGoExElCEOdiChBGOpERAlyaKiLSE5E/lFE/llEXhWRP+o/Pi0i3xORN/tfp8ZfXCIiuh85bJy6iAiAoqrWRSQD4IcAvgjgNwBsqepXRORLAKZU9ffHXmIiIjrQoTV17an3f8z0/yiAzwB4of/4CwB+bRwFJCKiwVmDPElETABXATwE4L+p6o9EZEFVlwFAVZdFZP6w3zM7O6vnzp07TnmJiFLn6tWrG6o6N8hzBwp1VfUBPCkikwC+JSJPDFoYEXkGwDMAcObMGbz44ouD/lUiIgIgIm8P+tyhRr+o6g6AfwDwaQCrIrLUP+ASgLUD/s5zqnpBVS/MzQ30QUNEREc0yOiXuX4NHSKSB/AUgNcAfAfA0/2nPQ3g22MqIxERDWiQ5pclAC/029UNAN9Q1e+KyP8D8A0RuQTgHQCfHWM5iYhoAIeGuqr+BMDH9nl8E8Avj6NQRER0NJxRSkSUIAx1IqIEYagTESUIQ52IKEEGmnw0Ki+//DJ6S8lQZIgBMTOAaUHMDMTsD3ISAJAPnK/eWkHaWygCXajvA34H2vWgXgfQbuAvge7v7NmzuHHjRtjFoIAEGuqu64IbXQdHVdF0fTRdH47no+MrOn4XrteF63fR8brojvh0GAJkLAO2aSBjGrAtAxlTYFsGiraFgm3ygz1g/PdOl0BDncan21U0Oz4ajtf/46PpeiMP7UPLoYDT6cLp7F9jNwQo2BYKWROlbC/ki7YFw2DwEI0CQz2mGo6HWttD3fHQdD00XR9xuAnqKlB3euVegwMAEAHyGRPFrIVi1kQ5l0Epy7cm0VHwyomJblex2+pgu+liu9mB6yWn7VoVd5qJ1mu9x2xLMFmwMVWwMZHPwGRNnmggDPUIa3f8O0G+2+wE3pQSJtdTrFUdrFUdGAJU8hlMFWxMFjLIZcywi0cUWQz1iKk7HrbqLrabLpquH3ZxIqGrwE6zg51mBwBQsE1MFWxMFTMo5zIhl44oWhjqEeD5XWzUXazV2mg4DPLD9JpqWri100LBNjFfyWKulIVlctoF0aF7lI70YCLKIY3vqTseVqttbNZd+GlqWxkDQ4CZUhYLlSxr7/cQEQ4ljjkRuaqqFwZ5LmvqAfO7is26g9Wqg7rjhV2cxOgqsF5zsF5zUMyamC/nMFuyWXun1GGoB6TpeljZbWODtfKxazg+fu408M5WEzMlG4uVHIocIkkpwXf6mLVcHze3m9iou2EXJXX87nsjaKaLNk5P51Gw+ZanZOM7fEzaHR83t1vYqDuxmBSUdFuN3oii2ZKNU1MFDoukxGKoj5jrdXFrp4W1ajtV48rjQBVYr7nYqLuYK2dxaiqPrMVwp2RhqI9Ix+/i9k4LK7sM86hTBdaqDjZqDhYqOZyYzMO22KFKycBQPya/q7i908LybpsdoDHTVWB5t421moPFSg4nJnMcLUOxx1A/hp2mi+sbjQNXJKR48LuKWzstrNcdPDhbxFTRDrtIREfGUD8Cz+/ixmYT6zUn7KLQCLleF6+t1DBbsnFutogMa+0UQwz1IW3WHdzYbMD12NSSVBt1F7utDs7OFDFXzoZdHKKhMNQH5Hpd3NhsYJPjzVOh4yveWqtjs+HggdkiR8lQbDDUB7BWa+PtzSY8n7XztNludFBr7+LMdAELlVzYxSE6FEP9PhzPx/X1xp0lXymdPF9xfb2BjbqD83MlTlyiSGNP0AGq7Q5eubXLQKc7qi0Pr9zaxS7fExRhDPV9rOy28dPbVXaG0gd0fMW1lSpu77TCLgrRvtj8cpduV/HzzQbWqhyqSAdTBd7ebKLheHhwrsT9UylSGOp9jufjjZU61zingW3UXbQ6u3hkocx2dooMNr8A2G312s8Z6DSshuOznZ0iJfWhvrLbxrVltp/T0e21s99iOztFQGqbX1QV1zfYfk6joQq8s9lE0/Hw0HwJImxnp3Cksqbe7SreXKsz0GnkNuouXlupocsVOykkqQt1v6t4fbXG6f40NjvNDq6tVLkUM4UiVaHu+V1cW65yQhGNXbXl4ae3q+j4XJaZgpWaUO8Feg21Nke4UDDqDoOdgpeKUN8LdA5ZpKA1XZ/BToFKfKh7fm/jAwY6haXp+ri2zGCnYCQ61PcCnU0uFLaGw2CnYCQ21FUVb6zWGegUGQ3Hx+sc7khjlthQv77RwG6Lo1woWmptD9c36mEXgxIskaG+vNvixCKKrPWai5vbzbCLQQmVuFDfbrh4e5MXDEXbu1stbNZZ8aDRS1SoN10Pb67VoWyypBh4a41LPdPoHRrqInJaRP5eRK6JyKsi8sX+418WkVsi8lL/z6+Ov7gHc73eSBdOzaa46Crw+koVjueHXRRKkEFWafQA/J6q/lhEygCuisj3+v/tT1T1q+Mr3mC6XcUbqzU4HQ4Xo3hxPcXrKzV8+MQEd1CikTi0pq6qy6r64/73NQDXAJwcd8GGcX2DQxcpvhqOj7fWOCKGRmOoNnUROQfgYwB+1H/oCyLyExH5mohMjbpwg7i908J6jSsuUrxtNVy8u8UOfjq+gUNdREoAvgngd1W1CuBPAZwH8CSAZQDPHvD3nhGRF0XkxeMX9/2arscLgRLj1k6LHad0bKIDDBURkQyA7wL4O1X9433++zkA31XVJw75PTrI8QahqnjlVpUXASVKwTbxkZMTMEbYvi4iGNV1R+EQkauqemGQ5w4y+kUAXAZw7e5AF5Glu5726wBeGbagx3Fzm7UaSp6m6+NdTkyiYxhk9MsnAPwWgJdF5KX+Y38A4KKIPAlAAdwA8NtjKN++6o7HTX4psZZ325gs2JjIZ8IuCsXQQM0vIzvYCJpful3FT27touVybC8lVzZj4BdOTY5kmCObX+JvpM0vUfPOVpOBTonndLq4sdkIuxgUQ7EK9d1WB8u77bCLQRSItaqD7QaH69JwYhPqnt/Fz9Y5QYPS5fpGnRtr0FBiE+o3t1tcBoBSx/UU73AuBg0hFqHe7vhYrbLZhdJpveawH4kGFotQf3erCS6+SGmlCry9xU5TGkzkQ73ueNios7OI0m270eH2jDSQyIf62xzWRQQAeIc7etEAIh3q2w0X1RaXAiAC9u5auQUe3V9kQ12Vvf5E93p3q4kuO5joPiIb6us1B032+BO9T7vTxWqNI8HoYJEMdb+rXKmO6AC3tlvwOCGJDhDJUF+ptuF6vMUk2k/HVy6XQQeKXKirKlb4hiW6r9Vqm23rtK/IhfpWw4Xr8daS6H46vmKjwZEw9EGRC/UVLgdANJDVXYY6fVCkQr3heByXTjSguuOh1uYsU3q/SIU6a+lEw+FCd3SvyIS653exyTVeiIayUWcfFL1fZEJ9rebAZ28+0VBUWVun94tEqKsqm16Ijmit1ubG0nRHJEJ9p9nhrkZER+R6ik3uZUp9kQj1tRqHZhEdx1qV1xD1hB7qflex02Qtg+g4qu0O14MhABEI9d1Wh1vVER2TKrDd5Jh1ikCob7EtkGgktnnHSwg51FXZ9EI0KjvNDhf5onBDvdr20PH5JiQaBb+r3Jyawg31bTa9EI3UFu98Uy/UUOcbkGi0dpouJyKlXGih3nA8TjgiGjHXU9QcrnSaZqGFOke9EI0HmzXTLbRQZ4cO0XjscLx6qoUS6qqKBm8Ricai1fG54mmKhRLqDdfnLFKiMVHt7YpE6RROqPMNRzRWvMbSK5RQZy2CaLwY6unFmjpRArHilF6Bh3q3q2i6ftCHJUqVdqfLpXhTKvBQb7geOOGNaPwaDitPaRR4qPO2kCgYNYfj1dMo+Jo6Q50oEKypp1MIoc43GlEQeFecToGHuuOx84YoCB2/yxUbUyjwUOf0ZaJgqAIuR8CkzqGhLiKnReTvReSaiLwqIl/sPz4tIt8TkTf7X6cOPZrICIpMRINyeWecOoPU1D0Av6eqHwLwrwD8jog8DuBLAL6vqg8D+H7/50Mw1ImCxO0i0+fQUFfVZVX9cf/7GoBrAE4C+AyAF/pPewHArx16NNbUiQLFmnr6DNWmLiLnAHwMwI8ALKjqMtALfgDzA/yGoQtIREfXYZt66liDPlFESgC+CeB3VbUqA9a6ReQZAM8AAAzzCEUkoqPiaLP0GaimLiIZ9AL9L1T1r/sPr4rIUv+/LwFY2+/vqupzqnpBVS8IQ50oUKypp88go18EwGUA11T1j+/6T98B8HT/+6cBfPvww7H5hShIbFNPn0GaXz4B4LcAvCwiL/Uf+wMAXwHwDRG5BOAdAJ899Dcx04kCxZp6+hwa6qr6Qxwcx7883OGY6kRB4pDG9AllkwwiChIrU2kScKjzzUUUOM4PSRXW1Om+rly5gosXP4dPfeqXcPHi53DlypWwi0RE9zHwOHVKnytXruDZZ7+KdtsBAKysrOLZZ78KAHjqqafCLBoNgzX1VGFNnQ50+fLzdwJ9T7vt4PLl50MqER0NQz1NGOp0oNXV1aEep6jiCJg0YajTgRYWFoZ6nCKKG2WkCkOdDnTp0ueRy2Xf91gul8WlS58PqUREdJiAO0pZY4iTvc7Qy5efx+rqKhYWFnDp0ufZSRo3rKmnigS5h6Fh5/X/XrsZ2PGICPjFh2a5V2nMichVVb0wyHPZ/EKUYKbBkS9pE2yos7ZAFKiMyVBPm4Br6gx1oiDZFm/G0ybQM852PaJg2SZDPW3Y/EKUYKyppw+bX4gSLMOaeuqwpk6UYKyppw9r6kQJxpp6+rCmTpRgWdbUUyfwM86lnYmCw5p6+gR+xvkmIwqGZQpnlKZQ4AlbzJpBH5IolQo2r7U0Cj7Ube6gRxSEUpbXWhoFHup8oxEFo8hrLZVCaH7hG40oCKxApVPgoW5bBidEEI1ZxhTkMmxTT6NQ0rWcYw2CaJx4R5xeoYQ633BE48Wml/QKJdRLHAFDNFasOKVXSDV1tvURjROvsfQKJdQt00Auw85SonGwLUHWYqinVWjJWslnwjo0UaJVcry20iy0UJ8u2GEdmijRpoq8ttIstFCfyGe42BDRiBkCTLHClGqhhbphCCYLvE0kGqUKK0upF2pvJWsURKM1zaaX1As51DPcNINoRIRNL4SQQ90yDfbUE41IKWtxXSUKN9QB3i4SjQpHvRAQgVCfKrKmTjQKHCZMQARCPWuZXHyI6Jjytok8t68jRCDUAWC2zBoG0XHMlngNUU8kQn2ulOXYWqIjMgSYL+fCLgZFRCRC3TIN1jSIjmimZHPUC90RmXfC4gRrGkRHsVDhtUPviUyoF2wLlTw7TImGUcpaKHOuB93l0FAXka+JyJqIvHLXY18WkVsi8lL/z6+OojCLrHEQDYV3uHSvQWrqfw7g0/s8/ieq+mT/z9+OojDTRRtZbp5BNBDbEsxwwhHd49AEVdUfANgKoCwQEbYPEg1ovpyDwVFjdI/jVIu/ICI/6TfPTB30JBF5RkReFJEXB/ml8+Us+D4luj8RYL6SDbsYFEFHDfU/BXAewJMAlgE8e9ATVfU5Vb2gqhcG+cUZ08BMiW9WovuZLtrch5T2daRQV9VVVfVVtQvgzwB8fJSFOjWV55K8RAcQ6V0jRPs5UqiLyNJdP/46gFcOeu5R5DIm29aJDjBbyqJgc/gv7e/Qd4aIfB3AJwHMishNAH8I4JMi8iQABXADwG+PumCnpvJYrznwuzrqX00UW4YAp6dZS6eDiWpwoSkiOszxbm438e5Wa4wlIoqXk5N5nJkpDPV3RARBXuc0eiJyddB+yUgPCj8xkeeaFkR9GVNwYpLNknR/kU5MwxCcZocQEQDg5FQelhnpS5YiIPLvkLlyFgUu/k8pl8sYXEaDBhL5UBcRnJkerg2RKGnOTBcgHOdLA4h8qAO9DXUn8lyJjtKpnLM4IY8GFotQB4BzswUuH0CpIwKcHXK0C6VbbEK9YFs4xWYYSpmTk3mul05DiU2oA8CJiRzKOc6ko3QoZk0uB0BDi1Woiwgemi9xk2pKPEOAh+ZL7BylocUq1IHeujAcDUNJd2q6wPVd6EhiF+pAbwsvjoahpCrnLJzgNnV0RLEMdQA4P1+EZfLWlJLFNITNLnQssQ31rGVyqBclztmZAnIZzqCmo4ttqAO9PRqnufEuJcRkIcN9BOjYYh3qAPDgXBG5TOxfBqVcNmPg/Fwp7GJQAsQ+DTOmgccWK2xfp9gyDcGjC2UuM00jkYh3Ud428ch8mfuaUuyIAA/Pl1DMcvgijUYiQh0AJgoZPDBbDLsYREM5M13AFPuFaIQSE+oAsFDJYYnjeykm5itZnJjkMgA0WokKdaA3JGyywIlJFG2VvIUHeWdJY5C4UBcRPDxf4m5JFFm5jIFHFsqcYERjkbhQBwDLNPDoYhkZjoihiLFMwWOLFWS41yiNSWLfWbmMiUcXy1zRkSLDEOCRhTLyvIukMUpsqANAOZfBY0sMdgqfIcBjixUuREdjl+hQB4AKg51CdifQ2YFPAUh8qAMMdgqPIcCji2UGOgUmFaEOvBfsXE6AgmIavU7RyQInF1FwUhPqQC/YP7RU4agYGjvLFDy2xBo6BS9VoQ4ApayFD5+Y4OJJNDYZU/D4UgWVHAOdgpfKZMvbJj58osIle2nkbMvAh09McIEuCk1qUy2XMfHEyQlU8rz4aDTKOQsfOTnBcegUqtSGOtBbi/3xpQoXAaNjW6hk8eETFTbrUehSX00VEZybLaKYtXB9vY6uhl0iihNDgAdmi5jnNnQUEakP9T1z5SwKtonXV2twOt2wi0MxYFsGHlkoocwOUYoQ3ivepZjttYlyKjcdZq/9nIFOUcNQv0fGNPChpTJOTPJ2mva3OJFj+zlFFptf9iEiODvTa2f/+UYDns+GdurNED03W8B8mR/4FF0M9fuYLWVRzvWCfbvRCbs4FKKJfAYPzhWRy3C4IkUbQ/0QWcvEY4sVrNccvL3ZQIe19lSxTMHZ6QJHt1BsMNQHNFfOYrKQwY2NBjbqbtjFoQBMF208MFtk2znFCkN9CBnTwMMLZcyUXPx8owHX49DHJLItwbmZImZK2bCLQjQ0hvoRTBdtVHIW3t5qYq3qhF0cGqG5so2zM0XuIUqxxVA/Iss0cH6uhNliFjc2G2i6fthFomPI2ybOzRS49jnFHkP9mCYKGXw0P4GNuoub2020ORs1VrIZA6em8pgrZSHCdfYp/hjqIyAimCtnMVuysV5z8O52i+3tEWdbgpOTBcyXszC4zSElyKGhLiJfA/DvAayp6hP9x6YB/E8A5wDcAPCbqro9vmLGg4hgvpLDbCmL1Vobt3dacD0OgYySjCk4MZnHYiXHMKdEGqQ36M8BfPqex74E4Puq+jCA7/d/pj7DECxN5PHk6SmcmSlwX9QIsEzBqak8njw9iROTeQY6JZaoHl6TFJFzAL57V039dQCfVNVlEVkC8A+q+ugAv0cHOV7SeH4Xy7ttrNXarLkHLGMK5ss5LE3mUjuiRUSQxusuSUTkqqpeGOS5R21TX1DVZQDoB/v8EX9PKlimgdPTBZyaymOr4WK16mC3xWUHxqmSt7BQyWG6YN+plftdRcfvwvW76Hhd+KpA738Q9P7PEEHGNGCbBjKmwErpBwHF19g7SkXkGQDPjPs4cSAimCllMVPKot3xsVZ1sFZrc+mBEcmYvQ7rybwNhaLueHhrvY6m68P1uvCPsAOKaQgypiBvmyjaFopZC8WsiazFNWAomtj8ErJuV7HVdLFabaPa8sIuTiwVbBO5jAkBUHO8QEYeZUxBJZ/BZD6DyYId6aUE2PwSf0E0v3wHwNMAvtL/+u0j/p7UMwzBbCmL2VIWLdfHWq2NrYbL8e6HEOk1mZimoNXxA5/81fEVm3UXm3UXQAOlrIXJQgYzJRsFmyOFKTyH1tRF5OsAPglgFsAqgD8E8DcAvgHgDIB3AHxWVbcOPRhr6gNruT62my62my5qbQ9p/2cTAYq2CTEEHa8b6Q+9cs7CfCWL2WI0xsCzph5/w9TUB2p+GRWG+tF4fhc7rQ52mi52mp3UtMFnTMFkIYOCbcHp+NhqurEaPWSZgrlSFosTuVDXYWeoxx9DPcFUex2AO80OdlsdNF3/SB2AUWQagoJtopLLYLKYQc4ycHunjdVqG3F+iSK9pZtPTuZDCXeGevwx1FOm5fpouB4ajoeG46PpepGvzVtmL8BLWQsF20IpayGXMSAid8b1L++2E/OBBQCGAPOVHE5O5gPtWGWoxx9DndDudx42HA8NtzcipON3A2++yJgC2zJgWwaKtoWCbaKYtQ6ssa5W23hnq5nofWFNQ3BishfuQSwixlCPP4Y6HUhVe5NvfEWnH/RO/2vH703O6c3J0f5X3AkEEYFgb+SJQKS3cUjG7E/YsfqTdqzeY7ZpDBxa7Y6Pn63XUzWss5g18eBcCaXseEfLMNTjj6FOsaGqWN5t492tZqzbzY9KBDgxkcepqfGtR8NQj78gxqkTHZvj+XhztY5aOz2183upArd2Wthqunh0oYy8zZmqdDzRnQZHiVZtd/DKrd1UB/rdWq6PV27vYqvBTc3peBjqFLiV3TZ+ersaqzHnQfB8xesrNby71Qy7KBRjbH6hwKgqrm80uFn3IW5ut9BwPTw8X4YZgRmpFC+sqVMgul3FG6t1BvqAthsdXFuuwvOjuxwCRRNDncau21W8tlJje/GQam0PP12uosNgpyEw1Ifkel20Oz5abm/mZsv14Xg+L7wD7AU6NwU5mobj4xqDnYbANvUDqCpaHR91x0PT6X89ZJ2VjCl3prwXsmZ/6nt6h6ipKl5fZaAfV8Px8dpyDY+fqLCNnQ7FUL/L3oYVa1UHtXZn6MkwHV+x2+q8L8QsUzCRz2ChnMNEITPiEkfbjc0mdpoM9FGoOx7eWqvjkYVSIEsLUHwx1IGxbi3n3bWZQi5jYL6Sw3w5m/hNkFd221jZbYddjETZarh4d6uFMzOFsItCEZbqUN9tdnBrpxVY80C708U7m03c3Gpiumjj5FQ+kbvk7DRd3NhshF2MRLq100LONjBfzoVdFIqo5CXKADy/ixubTazXwhle11Vgo+5is+Hi5GQeJyfHt+5H0NodH2+u1VO/U9M4XV9vIJ8xUc6lqzmPBpPsNoB9bNYd/PPNndAC/W6qvYkmL9/aRa0d/7ZnVcVba/VEL5sbBarAW2v1RK01T6OTmlB3vS7eWK3hjdV65KanN10fr96u4sZGI9YX6u3dNtdyCUi708XbbOKifaQi1GvtDn5yc6e/83s0qQLLu228fGsX7Y4fdnGG1nA83OSaJYFarTrYaUb3PU3hSHyo7zRdXFuuRX57tz2tfq296canxtvtKn62Xk/leuhh+9l6g0sJ0PskOtQ36w5eW6nFrknD9bp49XYVdScewb5SbaPhxO/uIglcr4ub262wi0ERkthQ3264sR6F4fmK15ajX2P3/C5u7zBUwrRabceyyY7GI5Ghvtvq4I3VWmwDfU/HV1xbrkb6gr2104pN01ZSdRVcg53uSFyod/wu3lqrJaZ91/V6wwSjuMdku+Nz1mhEbNTd2DTX0XglLtRvbDQiN2TxuGptD7ci2MRxc7uVmA/PJOAQRwISFuobdQcbER62eBy3tltoRKgm5npdbNbDn8BF76m2vEi9RygciQl11+vixkZyayrd/izCbkSqxqvVNmvpEbRSZXNY2iUm1K9v1BPfYdd0/UgMX1NVrNUYHlG0UXO4oUbKJSLUd1sdbDfiv3bKIG7vtuB64V60mw03cf0WSdFVYC0C6xpReBIR6mkagaHaa/oIU9jHp/vj+Um32Id6u+NjO2XrX6zV2qENcXS9Lhftijin0+XwxhSLfaivVZ3YTzIalutpaKN8dppu6v6942i7ka6KDr0n1qHe7aa3wy6sW+ytlN0VxVXa7l7pPbEO9Y2Gk/gRLweptb3Ab7H9rmKXG0nHQsPxI728BI1PrEM97QET9Frau60Ox6bHCGvr6RTrUE97Z1DQy93W2UEaKzxf6RTbUO/4XbQ76Z5kEfSHWto/ROOG5yudYhvqXOOiN7wwyIlIjYiv7U7v1+50uStSCsU21FkL6Qnqw63d8eGltFM6zrgjVfrENtT5Zu0J6sONH6LxVOfdVerENtQdL7xQv3LlCi5e/Bw+9alfwsWLn8OVK1dCK0tQ/w5OyOvN0NE4HNaYOlbYBTiqsDaTvnLlCp599qtot3uLJq2srOLZZ78KAHjqqacCL09Q/wxhLyJGR5PWeRxpFtuaeljjpS9ffv5OoO9ptx1cvvx8KOXpBjRnn8u5xhM/jNMntqEOhJPqq6urQz0+bkGtw8JwiCeXH8apc6zmFxG5AaAGwAfgqeqFURRqwGMjjGBfWFjAysoHA3xhYSHwsgCAIRLIcRgO8cQ7rPQZRU39l1T1ySADHQguzO516dLnkctl3/dYLpfFpUufD6U8RkD/DGH1YYxblDq9x0EVgMT4hpyGFtuOUiuoNLvHXmfo5cvPY3V1FQsLC7h06fOhdJICgBHQv0MSl9uNWqf32IRUAaJwyHE2WxCRnwPYRq8d5L+r6nOHPF9HtbnD9fU6VqvctuuB2SIWJ3JjP86Prm8mbjGvixc/t29T2uLiAr7+9b8KoUTj8YuPLEL9dC9+F3cicnXQ1pDj1tQ/oaq3RWQewPdE5DVV/cE9hXkGwDPHPM4HlLIWVsFQL2bNsIsQW1Hr9B6fhH0a030dq7FNVW/3v64B+BaAj+/znOdU9cKo29yL2di2HI2MCFC0g/l3kATewh/UuR1Wp/fYJLHtjA505FAXkaKIlPe+B/ArAF4ZVcEOU7DNwDoJo6poW4G1qZsJ7GuLWqf32DDUU+U41bwFAN/q1+AsAH+pqv9rJKUagIigmLVSvQlykE0vGdOAG+LSDOMQtU7vcbDMcIb+UniO1VE69MFG2FEKADc2GljeTecepQBwfq6I+cr4O0kB4NpyFTsp32kqjgq2iSfPTCHI65xGb5iO0ljfVM+Ws4c/KaFMQzBdtAM7nm3F+q2SWpkktpvRfcX6jJeyFsq5dHaYzpZsWAFesDbDIZb4YZw+sT/jQYzRjqKgXzfDIZ74YZw+sT/jM0UbtpWuYTCVvIVCQEMZ9+RtjoePowLnMaRO7ENdRDBfTldtfWkiH/gxS7bF2eYxVOJ8jtSJfagDwEIlBzMlg9ZzGQNThUzgxzUMQYG19VixTEEuw3OWNokIddsycHamEHYxxk4EOD9fCm12J2fxxktQs40pWhIR6kCvtj4ZQg02SEsTOVRy4b1G3srHS1pHhqVdYkIdAB6cK/Zn0CVPwTZxeircu5GJfLI/NJOmwvOVSokK9axl4txMMexijNxes0tQ67wcJJcx2a4eExlTUGFNPZUSFeoAMFfOYqYU3EzLIJyczEem6SPIWax0dJMFO5Era9LhEhfqAHB+rpSY9sS5so1TU8EPYTzIFEM9Fvjhm16JDHXTEDy2WI5M7faoZko2zs+FN9plP6WsxdmlEWcI+z/SLLFXp2UaeGypHNsa+2zJxsMhDl+8n/kUL6QWBzOlbGrmbdAHJTbUgd4KdR9aqsRuqOPiRA4PRTTQAWC+kuXs0ghL63pI1JPoUAfea4qJwxvdEODsTAEPzBYjG+hAb5TRDNtsI6mcs2Lf7EjHk/hQB3rrwzwwW8TjJyrIZaL5kss5Cx89NYkTk9HpFL2fhRh8SKbRQkCbplB0RTPhxmQin8EvnJrEiclcZJoPTKP3gfPEyYlYrYRYyWUC3U6PDmdbgtmEDeel4aUq1IHewlRnZ3ohGnYoTRYy+OipiVg0De0n7Bmu9H4nJwuRbrajYKS28a2UtfCRkxPYbXWwWnWw3XQD2XR9bxu6hUoW5RDXcRmFqaKNSt5CtZXezb+jIm+bWKhwVBKlONSBXlv7ZMHGZMGG4/lYqzpYqzlwve7Ij1WwTSxUcoFvQzduZ2eKePnmbtjFSL3TU3nW0glAykP9blnLxOnpAk5N5bHT7KDa7qDh+Gi4Hjx/+Cq8bQmKWQtF28JEIRPq6orjVMpamC3Z2Ki7YRcltco5CzMl1tKph6F+DxHBVNF+33T4dsdHw/HQdH04ng+/C3RVoQoYBmCIwJDeglelbG+ruTTNujw9XcB2swO/G0D7Fb2P9IfBEu1hqA8glzGRy5iYCbsgEZXLmDgzXcDPNxphFyV1liZyse+bodFKT3WSxmpxIsf1RgIWhTX2KXoY6jQy5+eTu0lJ1IgAD0VgjX2KHoY6jUzWMtm+G5BTU3nuGUv7YqjTSM2XcxwvPWZTxQxOxmQ5CQoeQ51G7oHZItvXx6SYNfHwfJlj0ulADHUaORHBIwulWK1lEwe2JXhkocy10um+GOo0FpZp4LHFMjLsOB0JQ4BHFsrIZfhBSffHUKexyWVMPLpY5oiYY9oLdI5Hp0Ew1GmsyrkMHmOwH9leoHPDbxoUQ53GrpzL4ENLFTbFDMkQ4NFFBjoNh6FOgShlLTx+opKqNXGOwzIFjy1VMFlgoNNweIVRYAp2bw37co6TZu4nb5t44sQEh4XSkTDUKVC2ZeDxpQonKB1gumjjIzHb2pCihVUmCpxhCB6cK6GYtXBjowGu2NtzaiqP09NcZoGOh6FOoVmo5FDMWvjZWh1N1w+7OKHJZgycny1hosDmFjo+hjqFqpS18NFTE7i53cLtnVaqau0ivQ+2M9MFzhKlkWGoU+hEBKenC5gp2bi+3kCtnfyNrPO2iQfniond5pDCw1CnyCjYFj58ooL1moObOy04ndFvAB62jCk4MZnHYiXHtdBpLBjqFCkigvlKDrOlLNZqDm7tNOF68W+TsUzB0kQOSxN5NrXQWDHUKZIMQ7A4kcNcOYuVahsru224Xvxq7hlTsFDJYWkiB8vkCGIaP4Y6RZppCE5O5nFiIoftZger1TZ2Wx1oxCvv5ZyFhUoOM0WbzSwUqGOFuoh8GsB/BWACeF5VvzKSUhHdQ0QwXbQxXbTR7vhYrznYqDtoR6jd3bYMzBRtzFeyKNisL1E4RI9Y5RERE8AbAP4tgJsA/gnARVX96X3+jh71eET7aXd8bDddbDc6qLaDrcGL9IZkThYymCrYkd0zVETA6y7eROSqql4Y5LnHeRd+HMBbqnq9f9C/AvAZAAeGOtGo5TImlibyWJrIw/O7qDseGq6PhuOh4XgjrcnbloFS1kLBNlHMWijnLGTYTk4Rc5xQPwng3bt+vgngXx6vOERHZ5kGJgs2Ju+aae/5XTQ7PjpeF67fRcfT3le/C78/00m1V+sGAEMEtiXImAZsy0DG7P0p2CYDnGLhOKG+X+/PB+7xROQZAM/0f3RE5JVjHDPKZgFshF2IMeLri69ZEUnqawMSfu7Qe21nB/0Lxwn1mwBO3/XzKQC3732Sqj4H4DkAEJEXB20XipskvzaAry/OkvzagGS/vqO8tuPcT/4TgIdF5AERsQF8DsB3jvH7iIjomI5cU1dVT0S+AODv0BvS+DVVfXVkJSMioqEdawyWqv4tgL8d4q88d5zjRVySXxvA1xdnSX5tQLJf39Cv7cjj1ImIKHo4RouIKEECCXUR+bSIvC4ib4nIl4I4ZpBE5IaIvCwiL4nIi2GX57hE5Gsisnb38FMRmRaR74nIm/2vU2GW8agOeG1fFpFb/fP3koj8aphlPCoROS0ify8i10TkVRH5Yv/xpJy7g15fUs5fTkT+UUT+uf/6/qj/+FDnb+zNL0dZTiBuROQGgAuqmoixsiLybwDUAfwPVX2i/9h/BrClql/pfzBPqervh1nOozjgtX0ZQF1Vvxpm2Y5LRJYALKnqj0WkDOAqgF8D8B+QjHN30Ov7TSTj/AmAoqrWRSQD4IcAvgjgNzDE+Quipn5nOQFVdQHsLSdAEaWqPwCwdc/DnwHwQv/7F9C7mGLngNeWCKq6rKo/7n9fA3ANvZnfSTl3B72+RNCeev/HTP+PYsjzF0So77ecQGJORJ8C+N8icrU/gzaJFlR1GehdXADmQy7PqH1BRH7Sb56JZfPE3UTkHICPAfgREnju7nl9QELOn4iYIvISgDUA31PVoc9fEKE+0HICMfcJVf0XAP4dgN/p3+JTfPwpgPMAngSwDODZUEtzTCJSAvBNAL+rqtWwyzNq+7y+xJw/VfVV9Un0Zuh/XESeGPZ3BBHqAy0nEGeqerv/dQ3At9Brckqa1X6b5l7b5lrI5RkZVV3tX0xdAH+GGJ+/flvsNwH8har+df/hxJy7/V5fks7fHlXdAfAPAD6NIc9fEKGe6OUERKTY77SBiBQB/AqAJC5a9h0AT/e/fxrAt0Msy0jtXTB9v46Ynr9+R9tlANdU9Y/v+k+JOHcHvb4Enb85EZnsf58H8BSA1zDk+Qtk8lF/iNF/wXvLCfynsR80ICLyIHq1c6A3Q/cv4/76ROTrAD6J3gpxqwD+EMDfAPgGgDMA3gHwWVWNXYfjAa/tk+jduiuAGwB+e68NM05E5F8D+D8AXgawt5D8H6DX7pyEc3fQ67uIZJy/j6LXEWqiV+H+hqr+RxGZwRDnjzNKiYgShDNKiYgShKFORJQgDHUiogRhqBMRJQhDnYgoQRjqREQJwlAnIkoQhjoRUYL8f8luUKIfKnfCAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Print variable values\n",
    "print(\"Name\\tValue\")\n",
    "for c in model.component_data_objects(pyo.Var):\n",
    "    print(c.name,\"\\t\", pyo.value(c))\n",
    "\n",
    "# Plot solution\n",
    "plot_circles(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "left_x_con[A] \t 10.0 \t 9.99999990125201 \t None\n",
      "left_x_con[B] \t 5.0 \t 14.999999849644157 \t None\n",
      "left_x_con[C] \t 3.0 \t 4.723219051632863 \t None\n",
      "left_y_con[A] \t 10.0 \t 19.1421355149319 \t None\n",
      "left_y_con[B] \t 5.0 \t 4.999999951252102 \t None\n",
      "left_y_con[C] \t 3.0 \t 5.4079499719072 \t None\n",
      "right_x_con[A] \t None \t 9.806250034216646e-08 \t 0.0\n",
      "right_x_con[B] \t None \t 4.6454646351890005e-08 \t 0.0\n",
      "right_x_con[C] \t None \t -12.276780751556647 \t 0.0\n",
      "right_y_con[A] \t None \t 9.87471615587765e-08 \t 0.0\n",
      "right_y_con[B] \t None \t -19.142135464932636 \t 0.0\n",
      "right_y_con[C] \t None \t -20.73418544427754 \t 0.0\n",
      "no_overlap_con[A,B] \t 225.0 \t 224.99999778541843 \t None\n",
      "no_overlap_con[A,C] \t 169.0 \t 216.47226866513608 \t None\n",
      "no_overlap_con[B,C] \t 64.0 \t 105.77864678972617 \t None\n"
     ]
    }
   ],
   "source": [
    "# Print constraints\n",
    "for c in model.component_data_objects(pyo.Constraint):\n",
    "    print(c.name,\"\\t\", pyo.value(c.lower),\"\\t\", pyo.value(c.body),\"\\t\", pyo.value(c.upper))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reinitialize and Resolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition note\"> \n",
    "    <p class=\"title\"><b>Activity</b></p>\n",
    " Reinitialize the model, plot the initial point, resolve, and plot the solution. Is there more than one solution?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 Set Declarations\n",
      "    CIRCLES : Size=1, Index=None, Ordered=Insertion\n",
      "        Key  : Dimen : Domain : Size : Members\n",
      "        None :     1 :    Any :    3 : {'A', 'B', 'C'}\n",
      "    no_overlap_con_index : Size=1, Index=None, Ordered=True\n",
      "        Key  : Dimen : Domain          : Size : Members\n",
      "        None :     2 : CIRCLES*CIRCLES :    9 : {('A', 'A'), ('A', 'B'), ('A', 'C'), ('B', 'A'), ('B', 'B'), ('B', 'C'), ('C', 'A'), ('C', 'B'), ('C', 'C')}\n",
      "\n",
      "1 Param Declarations\n",
      "    R : Size=3, Index=CIRCLES, Domain=PositiveReals, Default=None, Mutable=False\n",
      "        Key : Value\n",
      "          A :  10.0\n",
      "          B :   5.0\n",
      "          C :   3.0\n",
      "\n",
      "4 Var Declarations\n",
      "    a : Size=1, Index=None\n",
      "        Key  : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "        None :     0 :    25 :  None : False : False : PositiveReals\n",
      "    b : Size=1, Index=None\n",
      "        Key  : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "        None :     0 :    25 :  None : False : False : PositiveReals\n",
      "    x : Size=3, Index=CIRCLES\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          A :     0 : 19.584576155980667 :  None : False : False : PositiveReals\n",
      "          B :     0 :  7.314989643761216 :  None : False : False : PositiveReals\n",
      "          C :     0 :  6.471478399524947 :  None : False : False : PositiveReals\n",
      "    y : Size=3, Index=CIRCLES\n",
      "        Key : Lower : Value              : Upper : Fixed : Stale : Domain\n",
      "          A :     0 : 14.902286214475343 :  None : False : False : PositiveReals\n",
      "          B :     0 : 10.489228593009017 :  None : False : False : PositiveReals\n",
      "          C :     0 :   9.00880553584265 :  None : False : False : PositiveReals\n",
      "\n",
      "1 Objective Declarations\n",
      "    obj : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Sense    : Expression\n",
      "        None :   True : minimize : 2*(a + b)\n",
      "\n",
      "5 Constraint Declarations\n",
      "    left_x_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body : Upper : Active\n",
      "          A :  10.0 : x[A] :  +Inf :   True\n",
      "          B :   5.0 : x[B] :  +Inf :   True\n",
      "          C :   3.0 : x[C] :  +Inf :   True\n",
      "    left_y_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body : Upper : Active\n",
      "          A :  10.0 : y[A] :  +Inf :   True\n",
      "          B :   5.0 : y[B] :  +Inf :   True\n",
      "          C :   3.0 : y[C] :  +Inf :   True\n",
      "    no_overlap_con : Size=3, Index=no_overlap_con_index, Active=True\n",
      "        Key        : Lower : Body                                : Upper : Active\n",
      "        ('A', 'B') : 225.0 : (x[A] - x[B])**2 + (y[A] - y[B])**2 :  +Inf :   True\n",
      "        ('A', 'C') : 169.0 : (x[A] - x[C])**2 + (y[A] - y[C])**2 :  +Inf :   True\n",
      "        ('B', 'C') :  64.0 : (x[B] - x[C])**2 + (y[B] - y[C])**2 :  +Inf :   True\n",
      "    right_x_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body              : Upper : Active\n",
      "          A :  -Inf : x[A] - (b - 10.0) :   0.0 :   True\n",
      "          B :  -Inf :  x[B] - (b - 5.0) :   0.0 :   True\n",
      "          C :  -Inf :  x[C] - (b - 3.0) :   0.0 :   True\n",
      "    right_y_con : Size=3, Index=CIRCLES, Active=True\n",
      "        Key : Lower : Body              : Upper : Active\n",
      "          A :  -Inf : y[A] - (a - 10.0) :   0.0 :   True\n",
      "          B :  -Inf :  y[B] - (a - 5.0) :   0.0 :   True\n",
      "          C :  -Inf :  y[C] - (a - 3.0) :   0.0 :   True\n",
      "\n",
      "13 Declarations: CIRCLES R a b obj x y left_x_con left_y_con right_x_con right_y_con no_overlap_con_index no_overlap_con\n"
     ]
    }
   ],
   "source": [
    "# Initialize and print the model\n",
    "\n",
    "# Add your solution here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAFlCAYAAADyArMXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAn40lEQVR4nO3da4xk6Vkf8P9zzqlT175OT/dc9jJe48sam2B7hSKBIhAGJUskwwdM/CE4ypLlA44g4gOWv8RShGQhG/KBCGVhLYwEBEvYYMEqwUYoJhJC7Gwc1uvF9u5s7+7s9P1S93N/8qGqZ3tm+lLVdS516vx/0qi7T1d3vVVz6t9vved531dUFURElD9G1g0gIqKLYYATEeUUA5yIKKcY4EREOcUAJyLKKQY4EVFOWWne2crKit64cSPNuyQiyr2bN2/uqurl+4+nGuA3btzA888/n+ZdEhHlnoi8ftJxDqEQEeUUA5yIKKcY4EREOcUAJyLKKQY4EVFOMcCJiHKKAU5ElFPnBriIPCwifyMiL4vISyLyy8PjnxGRt0Tkm8N/TybfXCIiOjLKRJ4AwK+q6gsiMgfgpoh8bfi931LVzyXXPCIiOs25Aa6qGwA2hp+3ReRlANeTbhgREZ1trKn0InIDwAcB/D2AHwbwSRH5eQDPY9BLPzjhZ54G8DQAmKYJEZm0zUREufXoo49ifX09lt8lo+6JKSINAP8bwK+r6pdFZA3ALgAF8F8AXFXVf3/O71DuwUmUX34YDf4FCjcM4YeKKBq8plUBhQ4/vk0AiAACGXwUwDYNlEwDJctAyRTYplGYzp2IYNwcFJGbqvrE/cdH6oGLSAnAnwL4Q1X9MgCo6tax7/8ugL8Yq0VENFX6XgjHD+GFEbxgGNShwgsieMPgTqr/JQKUTBmE+vCfbRqwhwFfLpmolUwYRjFCflTnBrgM/iw+C+BlVf3NY8evDsfHAeBnAHwrmSYSUZxUFX0/RMcN0HOHH70QYZTdu2NVwAsUXhACCE+8jQhQs03UbAuNsoVa2UTdtmAWONTPHUIRkR8B8LcAXgQQDQ9/GsDHAfwgBu+W1gH84rFAP+13cQiFKEVRpOj5IXpucDeou26ADLM6ViJApWSiUR4Ee71soW6bsMzpneIS5xDKyGPgcWCAEyWv4wY46Hpo9n103CCxYY9pVikZWKzZWKqVMF8pTdXQS+pj4EQ0vcJI0ez7OOh5OOz58ILo/B+acY4fYbPpYLPpwDQEC9USlmolLNZs2Nb09s7HxQAnyiHHD3HYG4R2q+/PzJBIEsJIsd/1sN/1AHTRKFtYrJWwVLfRKOc7AvPdeqIC6bgB9jseDnoeet7JF/rofJ3h9YDbB33YlmCxZmO5ZmOxVspdKSMDnGiKhZFit+Niq+Wg6zK04+YFiu2Wi+2WC9sSrM5VcHmujErJzLppI2GAE02hrhtgq+Vgt+NlWt5XJF6guH3Qx1uHfSzWSlibq0x9r5wBTjQlwkix13Gx1XLRcYOsm1NYqsBB18dB14dtGVidK2N1voyyNX29cgY4UcZ6XoCtlovdjosgZG97mnhBdLdXvlSzsTZfxkJ1enrlDHCijBx0Pbx12EfbYW972qnibiVLuWTg6kIFa3OVzOvLGeBEKWv2fby532Nw55TrR1jf7eHOoYOHlqpYnStn1iNngBOlpO34eHO/j2bfz7opFAMviHBrp4s7h308tFTDSsNOPcgZ4EQJ67oB3jzo4aDL4J5Fjh/hle0O7hyaeGipikuNcmr3zQAnSkjfC3H7oIfdjpd1UygFPS/Ed7c6aBw6eHi5isWanfh9MsCJYub4IW4f9LHbcQu5kFTRddwAL2+0MVex8MilGuYrpcTuiwFOFJMoUrx12Medwz7XJiG0nQAvvdXCct3GjZVaInXkDHCiGLQdH6/udNHnGiV0n/2uh5bj49HlGlbnK7H+bgY40QTCSPHmfg+bLYfDJXSqIFS8utPFbseDWPGNjTPAiS6o2fPx6m4Hrs/1t2k0zb4P+8q7sNHs48p8ZeKyQwY40ZiCMML6Xg87bTfrplAeiYH13R72Oh4eu1xHzb54DDPAicaw3/Xw2m4HXsDxEppM2wnw4u0mri1W8dBS9UK9cQY40Qi8IML6Xhd7rOmmGEUK3D7oY7/r4Z2rjbF3CJqdzeGIEtJyfLz41iHDmxLT80J8660mNpr9sX6OPXCiM2w2HazvdVlhQolTBdZ3e+i6AR5baYy00iEDnOgEUaS4tdvlhUpK3U7bQ89r4t1rc+du7cYhFKL7uEGIl+60GN6Uma47GFJp9s5eAI09cKJjmn0fr2y3WWVCmfNDxcubLTyyXDv1NgxwoqGNZh+v7/U43k1TQxV4fa936vcZ4FR4HO+mvGKAU6G5QYjvbna4CzzlEgOcCqvvhfj2RgtewLVMKJ8Y4FRIXTfAP222eLGSco0BToXTcnx8Z7ONIGR4U74xwKlQmj0f39lqI+SWOTQDGOBUGIc9D9/ZbHO7M5oZnIlJhcDwplnEAKeZx/CmWcUAp5nG8KZZxgCnmdXs+wxvmmkMcJpJfS/Ed7cY3jTbGOA0c/wwwj9ttljnTTOPAU4zRVXx3a02HJ/T42n2McBpptza7aLV58JUVAwMcJoZdw772G5xSVgqDgY4zYT9roc39k9f+J5oFjHAKfe6boBXtjvcSYcKhwFOueYFEf5pk4tTUTExwCm3omhQccINGaioGOCUW6/tddF2WHFCxcUAp1w66HqsOKHCY4BT7vhhhFu7naybQZQ5Bjjlzmu7Xe5lSQQGOOXMTtvFXsfLuhlEU4EBTrnhBiHW97pZN4NoajDAKTdu7XS5wiDRMecGuIg8LCJ/IyIvi8hLIvLLw+PLIvI1Efne8ONS8s2lotpqOTjs+Vk3g2iqjNIDDwD8qqo+DuCfA/glEXkfgE8B+GtVfReAvx5+TRQ7xw/x+h7XOSG637kBrqobqvrC8PM2gJcBXAfwUQBfHN7siwB+OqE2UoGpKl7Z7nCqPNEJxhoDF5EbAD4I4O8BrKnqBjAIeQCrsbeOCm+z5XC2JdEprFFvKCINAH8K4FdUtSUio/7c0wCevljzqMiCMMLtg37WzSCaWiP1wEWkhEF4/6Gqfnl4eEtErg6/fxXA9kk/q6rPqOoTqvpEHA2m4njrsM+qE6IzjFKFIgCeBfCyqv7msW99FcAnhp9/AsCfx988KirHD7HZdLJuBtFUG2UI5YcB/FsAL4rIN4fHPg3gswC+JCJPAXgDwM8m0kIqpNsHPfC6JdHZzg1wVf0/AE4b8P7xeJtDNNhhZ6fN6fJE5+FMTJo6rPkmGg0DnKbKYc9Ds88Zl0SjYIDT1FBV9r6JxsAAp6mx03HR88Ksm0GUGwxwmgpRpHhzn5N2iMbBAKepsNV2uLs80ZgY4JQ5VeWkHaILYIBT5g57PhyfvW+icTHAKXObLfa+iS6CAU6Z6nshd9ohuiAGOGVqi71vogtjgFNmwkix03GzbgZRbjHAKTM7bZfrfRNNgAFOmeHFS6LJMMApE82ejz6nzRNNhAFOmdhocdo80aQY4JQ6x2fpIFEcGOCUuoOeB+W1S6KJMcApdftdbpdGFAcGOKXKDyO0nSDrZhDNBAY4pYrDJ0TxYYBTqg66vHhJFBcGOKUmjBSHPY5/E8WFAU6pafZ9RBw+IYoNA5xSw+oTongxwCkVqhw+IYobA5xS0XIC+Fx5kChWDHBKxQGHT4hixwCnVBz2WT5IFDcGOCUuCCMuHUuUAAY4Ja7rMryJksAAp8R1PK59QpQEBjglrusywImSwACnxHUY4ESJYIBTovwwgutHWTeDaCYxwClRHD4hSg4DnBLF4ROi5DDAKVEsISRKDgOcEtVxOQOTKCkMcEqMG4TwAi5gRZQUBjglhtPniZLFAKfEeAHLB4mSxACnxLgMcKJEMcApMX7IACdKEgOcEuMxwIkSxQCnxPisQCFKFAOcEuOFrEIhShIDnBKhqtzEmChhDHBKhBdGUOY3UaIY4JQI9r6JkscAp0RwEg9R8qysG0CTiyKFF0bwwwh+qPDDCEGkUNW7wxgigIjAEKBkGiiZBmzTgG0ZMA2JvU2sASdKHgM8R8JI0fUCdN0AXTdEzwvgBhGCCYcrTENQMgU120K9bKJuW6iXLdjWxd+gsQdOlDwG+BTzwwiHPR/NvoeOG8Lxw0QuDIaRIowUju9hv/v2cdsS1MsW5islLNVsVG1z5N8Z8QomUeLODXAR+QKAfw1gW1XfPzz2GQD/AcDO8GafVtXnkmpkkfS8AAc9HwddDx03yLSSwwsUXuDjoOvj9b0eKiUDSzUbSzUb81ULIqcPvTC/iZI3Sg/89wH8NoA/uO/4b6nq52JvUQF5QYTttoPttjvVGwA7foSNpoONpgPLFFyq21idr6BRfvA0Yn4TJe/cAFfVb4jIjRTaUjiHPQ/bbRf7XS93PdYgVGy1XGy1XDTKFtbmy7jUKN+9IKp5e0AUm69//et49tnfw9bWFtbW1vDUU7+Aj3zkI1k3ayZNMgb+SRH5eQDPA/hVVT2IqU0zLYoUW20Hm00HzhT3tsfRcQN0dgK8vt/DSqOMa4sV9sAL6utf/zo+//nPwXFcAMDm5hY+//nBG3WGePwuWmbwOwDeCeAHAWwA+PxpNxSRp0XkeRF5/oL3NRNUFVstB//3zUOs7/ZmJryPC0LFZtPBN984xO39PgKWEhbOs8/+3t3wPuI4Lp599vcyatFsu1APXFW3jj4Xkd8F8Bdn3PYZAM8Mb1u4jpmqYqfj4vZBf6rHt+MUKbDTcdHu+1iu2/cMrdBs29raGus4TeZCPXARuXrsy58B8K14mjNbmj0f/+92E69udwsT3kdEBqWEux0X39tuY6/jcly8ANbW1sY6TpM5N8BF5I8B/B2A94jIbRF5CsBviMiLIvKPAH4MwH9KuJ25EoQRXt3p4NsbrcJu7Hu8vx1Fg+Gj9b0eHL+Yz0dRPPXUL6BSKd9zrFIp46mnfiGjFs22UapQPn7C4WcTaMtM2O96eG23W/iZiCcNmPS9AK/tdrHSsLHSKJ9ZR075dHShklUo6eBMzJj4YYT13S52O17WTZkKp2WzqmKn7aLlBLi2UB1rdiflw0c+8hEGdkq4GmEMOm6AF99qMryPkRP74G9z/RDre10c9PicEV0Ue+AT2mm7uLXTQcTrc/eQEboGqoqNwz4cL8TaQgUGh1SIxsIAvyBVxfpeD5tNJ+umTKWSOfqbu4OeByeI8PBSFdYYP0dUdHy1XIAfRvj2RovhfQZrzLrvvhfg1m4XPS9IqEVEs4cBPiY3CPHtOy20+gyas4zTAz8ShBHe2Ouh4/K5JRoFA3wMjh/ipTst9Apa2z2OcXvgRyJVvLnfQ8vxY24R0exhgI/oKLyLNqPyokzDOLWU8DyqitsHfbT6DHGiszDAR+D4Ib690Sr85JxxTXRBUhW3D/vsiROdgQF+Dj+M8PIGe94XcdFhlLuGPfEux8SJTsQAP0MUKb6z2Z7JpV/TcJELmQ8Yhjjf/RA9iAF+hlu7XbQd9v4uauIe+FAYRXhjv4eQs6WI7sEAP8Vbh33stN3zb0ininNSjheEuH3Q45K0RMcwwE+w3/Xwxl4v62bkXsmMd2p81w04eYroGAb4fdwgxKs7naybMRNiGQO/z0HPQ5PlhUQAGOAPeHW7iyDk2/Q4VErmyQuDT2iz6XC/TSIwwO+x2XTYu4uRIQLbin+97zCKcIdDKUQM8COOH+KNfY57x61aSuYU6zg+1xKnwmOAYzB1+5XtDsvUElAtJbfjzlbTYX04FRoDHMB222W9d0KSDPBIB5slExVV4QM8jBS3Dzh0kpSkLmQeaTs+1xCnwip8gN857MMLOHSSFBFBJYELmcdttTjhioqp0AHuBRE2WM2QuErCO8/3vYBLz1IhFTrAbx9wfY00JDkOfmS77XCaPRVOYQPc8UNsc62TVFRSCHAviHDQYy+ciqWwAb7ZdMAOWzqqpYvvzjOO/S7rwqlYChngYaTY6bD3nR5JqRceckNkKpRCBvhux+V6Jymbq5RSuZ8D9sKpQAoZ4FySNH3zFSuV+2m7AWdnUmEULsCbfR89L8y6GYVjWybsFIZRoMo1UqgwChfg3GUnO/PldHrhzb7PkkIqhEIFuKrikL2zzMylNIwShBH6Pt9l0exL5xU1JVr9AH6OLl66QQjHj+D4IdwgQhgpVBUKQBUQGay5LQKUDAOVkoFyyUS1ZMKMaUPhOFVtC5YpqVxAbjsBanahTm8qoEKd4ftT3PsOwggdN4DjD3qPrh8iGnMYoNl/+/OSaaBSMlEpmajaJuq2CUmjGPscc5VSKpUibcfH2nwl8fshylKhAnwaL2513ACHPQ8tJ0CcM4v8MIIfRmg7g9mJJdPAYs3GUq0U627x45qrWKkEuBcM3rmkUX9OlJXCBHjXDeD601FeFkQRmj0fBz0fXpDOWK0fRthpO9jpuJgrW1iq2WikNCZ9XN22YBiCKIU1aNpOwACnmVaYAJ+GvS69IMJO20XLybBKQhVtx0fb8WFbBpZqNpbrdmrDKyKCubKVyv/HYJ3wcuL3Q5SVwgR4N8Mp1qqK/a6H7bY7VeVtXhBhqzXYyPnaYjW13mqjkk6AsxKFZl1hygizWiPD8UOs7/Ww1Zre5U4dP8St3S62W87YF04vYq5SgpFClUwUKWdl0kwrRIAHYQQn5fHvSBU7bQe3drvo52HLL1Xsdlzc2ukmvkWZIYLFajpro7AXTrOsEAHeddN9ETt+iNd2u4NZn1Pa6z6NFwzeMQyW202u7ct1O7HffZzDAKcZVogA76TYA+66Adb3unDzHByq2O+6eGO/l1i1iG2ZqVTBsAdOs6wQAd5PafGqVt9PNPTS1nUDvL7fQxAlM/y0XEu+F84xcJplhQhwP0z+RXzY83D7sD+1Fyovqu8FeH0vmRBvVCzYVrKnYDBcfoBoFhUiwJPuhbX6Pu40ndyNd4/K9UO8vpfEBtCCpaR74aoIZuQdEdH9ChHgSfbA246P24f9mQ3vI64fJjI8tFhLvqSQuy/RrJr5AI8iTWwFQi+I8FYBwvtI3wuw0Yp3NyPTMLCQcEmhn9AYPlHWZj7AvQR733ea/Zm5YDmqZs9Dy4l3FmXSJYXsgdOsmvmp9PGP2w7sdVz0pnwH9Js3X8Bzz/0l9g/2sby0jCef/Cl8+MMfmvj3bjYd1GwTlhHP3/+yZaJWNtFLqF4/qXOAKGsz3wNPYmq4G4TYnvKt2W7efAFf+tKfYH9/H1Bgf38fX/rSn+DmzRcm/t1BGGHjMN6hlJV6kotOMcBpNs18gMf90lVV3Dmc3nVNjjz33F/C8+5dd9vzPDz33F/G8vvbjh/r9nSNSglVO5nFtKb8v4rowmY+wOO22/FysbbJ/sH+WMcvYrPlxFrhc4U76BCNZeYDPM4CtSCKsNuZ7qGTI8tLy2Mdv4go0sF6LzGp2hbmE6hImYKd5IgSMfsBHuOr97CX4UYMY3ryyZ+Cbd9b3WHbNp588qdivZ9m3491lubqfDnev7oAEviFRFPh3AAXkS+IyLaIfOvYsWUR+ZqIfG/4cSnZZl5cXHNEjjZlyIsPf/hD+NjHfg7Ly8uAAMvLy/jYx34uliqU41QVh734ygpt04x9jZQUlh4nysQoZYS/D+C3AfzBsWOfAvDXqvpZEfnU8Otfi795kyvFtIFv2wkQpLCmSpw+/OEPxR7YJznoergU47Zsl+fKOOz7sdXYm0xwmlHnppuqfgPA/Ve+Pgrgi8PPvwjgp+NtVnxKphFLD2x/Cne0nxZ+GKEdY028aRi41IivFx7XH3GiaXPRM3tNVTcAYPhx9bQbisjTIvK8iDx/wfuaWGnCFe8cP5z6STtZO4h5eOlSvQwrpuC1TPbAaTYl3jVR1WdU9QlVfSLp+zqNPWEQxB1Os6jrBnCD+GZSGiJYnYtnck8pphmjRNPmomf2lohcBYDhx+34mhS/8oQ98DiHB2ZZx4n3eVqslVAuTTa5xxBJZQNloixcNNm+CuATw88/AeDP42lOMibZNMAPo9xdvMxKL/adj2TiyT2TDp8RTbNRygj/GMDfAXiPiNwWkacAfBbAT4jI9wD8xPDrqVWzL75mV/yhNLuS2H+yXrYmWq2wYiUzPZ9oGpybbKr68VO+9eMxtyUxjfLFAzyt/TRnQRBG8IIo9m3SVucr6LjBhXZWqpTYA6fZVYizu2qbF65E4K7m40ni+TJEcH2xeqEJlUktkEU0DQoR4ABQv+AwisMAH0tSz1fVtsZfclYElQkvghJNs8IE+EWGUYIoys3aJ9Miyd1vLs+VxwrksmXA4EpWNMOKE+CVCwQ4t+IaW5wLW91PRHBtsTry6oI19r5pxs38lmpHFqolGAKMs7zGNAZ4UtukxSWpDaSPVEomLs9VsD3C5soX+aNNlCeF6YGbhmChNt5a00lsxzaJJLdJi0saz9lKwz734qQhgvoE1UdEeVCYAAcw9jKl0xXfyW+Tlh+C60vVM2dY1ssWx79p5hUqwBdrdq53Z0ljm7S8sE0Ta2fM0pzj8AkVQKEC3LaMsapRpi3s09gmbVJpPmVLNfvkWZoimKvEvzUb0bQpVIADGGudaWvKFkFKa5u0SaS9ecKVhcoDFyvrtslNHKgQChfglxvlkV/c1pQtQ5rWNmmTSH/zhMEsTftYyeAka6cQ5UnhBgot08BKw8ZW6/zd1KdxI4C0tkm7qCzetZiGgUeWq3htpwtDZKK1b4jyZLq6mCm5sjDaEqWGSGy7whRFVsu32qaJh5druNSIb29OomlXyHSq2RYWqqNd5KpyNt9Ysny+5iolfPDRxczunyhthQxwYPReOFezG4NIpgG+0rBxbaE2WLmQqAAKG+BLtRJqI4TzKLehgYplZLZ9mQhwbRjcj1yq8UImFUJhA1xE8Oil2rm3q5TM6SsIn1JZvltZm6/cs1Lh9602UC/zjy/NtsIGODCYmXneWLiR8bBAnkyydd0kLFPw0NK9wyamIXjvlXkOgdFMK3SAA4O32+cZ9YJnkRmGYC6j8r2rC5UT689ty8DjV+cY4jSzCh/gjbKFy3Nnj5cuVkuZje3mxWLVzuQ5si0D1xZOv2hZtkw8fnWOe2PSTOJZDeChpdqZszMNQ9gLP4tIZhcNH1munfuHo2yZeN+1eYY4zRye0RhcqDzvgiarGk7XKFux70Q/iuW6jctzo+2TyRCnWcSzeWhtvoLFMzZ8KFsmNwg4xbjrrMehZAresVIf62fKlonvv7bA0lCaGQzwYx67XD9z/RP2wh9kW2YmW5c9drlxoV6/bRn4/mvzXC+cZgID/JiyZZ7Zq5urlFBmSeE9svijdnnulHXAR2SZBh6/Oo/5KkOc8o0Bfp+VRvnMqpRrC1VO7Bmqla3UA7xcMnDj0nhDJycxDcHjV+axVOfFacovBvgJHltpnLokadU2scKhFBgiZ5bvJcEyBe+9MhfbCpGGIXjP2hyuLY62Lg7RtGGAn8AwBO+5MnfqGOvlufI907aLaG2hkmrliQjwrtVG7LM9B0sq1PGutQZ38aHcYYCfwrYMvPfK3IkvahHBtcVqYdedblRKWEq58uTRSzUsJnifK40yvv/aPMosM6Qpc1a/gmfrGeplC9+32jjxe5WSidURa5BniWkYuDbiUrxxWZsv42oKwzX1soUPXF/gpC2aGuWSgfdfXzj1+wzwcyzXbdxYOXmSz3LdRq1gteFXFyqp7lK0VC+NXe89iZI5WD+Fa4pT1haqJXzg+sKZ808Y4CO4ulA9caamiODhpVphxsOvLFQxn2LvdKlewrtX51IfqhIRPHKpxnFxysy1xQoevzp37ibhDPARXVs8OcRNY/Bin/X68LX5Sqolg0fhneUiYiuNMt5/ndPvKT2GDNayf/RSfaSOC8/MMVxbrOKxy/UHysAtw8Cjy7Mb4pfnKrjUSG+8f6Vh4z1r2Yb3kZo9GBdnvTglrTIc7x51fR+AAT62tfkK3nm58cCVYcscTDCpZrSpQSJEcGWhOtYJNanV+TK+b7UxVRU+lmngvVfm8c7Vs5daILoIkcG1pR94aHHs9ZZmKG3SM6gDN/DdrTa8QO8eN43BNm1vHfTRdvwMWzg5EcH1xfTGvEWAG5fqI282nYXVuQoWqzbW97rY63hZN4dmQM028djlOuYqF3udsQd+QXOVEj5wffGBRZEMGWzvdWWhCmOKepHjqNoW3rFSTy28bUvw+NX5qQ7vI7Zl4N1rc3j3WgO2lc//X8qeCPDQUhUfuL5w4fAG2AOfiG0ZeN/VeazvdbHVcu8el+EGB42yhY1mH103yLCVozNEcHmujOW6ndoQRqNs4d1XGihb+bp+cKlRxny1hNf3ethpu+f/ANGQej38wEMLscwqFlU9/1YxERFN8/7StN1ysL7XQxg9+PgOex42Ww6iE743LeplC1cXqqlOj1+dL+Mdl+pTcbFyEoc9D7d2u3D9KOum0BQzBHh4uYbrSzWMm4MiclNVn3jgOAM8Po4f4tZOF83+g+PfQRhho+lM3di4YQjW5iupTo0vlww8tlJPdGp82sJI8cZ+D5tNJ+um0BSar1p45+UGKiUTIsIAn2bbLQev7/cQhA8+1p4X4KDno9X3x/5PjJNtmViqlbBQK8Ey0ut1X1mo4JHls/cgzbOW42N9t4uuG2bdFJoCtiV4aKmGtfm3r+8wwHPADUKs7/aw3z25WiGMFM2+j4OeB9dP58UuIpirWFiq2alvD1cdXm2fn+CCTZ7sdVy8edBH32OQF5FlDha8uzJfeaCzwgDPkcOehzf3++iccSGz5wU46Ppou34i4+TlkonFavq9bWDQA7m2WMXaXCX3Y93jUlXsdFzcPuhzfLwgTENwdaFy5ppBDPAc2hu+kHvn9Mi8IILjh+j7IRw/hONHCKMRX/wiKFsGKpaJSslApWSiUjIzGa44qwdSNKqK7fbg/98LGOSzyJDB8OC1xeq565cwwHNKVbHb8XD7oAdnjB6ZH0Zw/QiRKiJVKADVwUkjIhABLENQKZmZ156P0gMpqihSbLYc3Dnswz/h+gjljwiwOlfG9aXqyKWwDPCcU1Uc9nxstR0c9nzMwlNSL5tYnatgpWEzuM8RRoqNZh8bTefEC900/UQGi509tFQdezVSBvgMcfwQO20X223nnmn5eWDIYELL2nx5otlkRRWEETZbDrbbLsfIc8I0BCsNG1cXqqjaF5t8xgCfQUe98v2eh8OeP7VjpaYhWKiWsFQrYbnO3nYcVAcVSVstFwc9bybekc2aRtnC6nwZK43yxNd0GOAF0HEDHHQHYX5WBUsayiUDSzUbS7US5iulwlWTpMkNQmy3XGy33an9I14UpiG41LCxNl9BI8ayWwZ4wfhhhK4boOuFg49uMNZF0HGUTEG9bKFuW6iXTdTLVmF2HJomqoqDno/tGbpOkhdJX89hgBOCMELXC+EFEfxw8M8LInhhBD/U4ZosClVAAQgGF14MEZRMA7ZloGQaKJkCe/h11TZzt6hUEeT5OklepHk9hwFOVECqirYb4LA7mMF73pwCOpttCRZrNharJSzW7NTmK8QZ4FxOlignRATzlcF1iEcu1eD4IQ57gzBv9X1M8WKXU6NRtrBYK2FpuNxz3uX/ERAVVKVk4sqCiSsLlXvW1jnseRxqGTpeNbVYs1NdLjkNEwW4iKwDaAMIAQQndfGJKHmmMdhEZLk+WKL3qIqp4wboeUFhAt00BDXbvNvTnvWqqTh64D+mqrsx/B4iikmjbN0zROAGIXpuOAz0wce8lylaptytlqrZg8dbKRlTtSF20jiEQlQAZWtQYbRUf3sTDT+MBqHuBei5AToJlqdOyrbkbkjXbJa3Hpk0wBXAX4mIAvjvqvrM/TcQkacBPD3h/RBRzEqmgYWagYXa22Vzqnq3FNUflqW+Xaqqg4XVgghBGE180VRkMO9gUM46+Fc+Vt5asgzYw+NFX9HyNBOVEYrINVW9IyKrAL4G4D+q6jfOuD3LCIlmxNH8gzB6e4VMKKDD+QfAIKQFAsjR54O5CNZw/kGRhjuOTE0ZoareGX7cFpGvAPghAKcGOBHNjqNeM2Xnws++iNRFZO7ocwA/CeBbcTWMiIjONkkPfA3AV4ZvgSwAf6Sq/zOWVhER0bkuHOCqegvAP4uxLURENAYOYBER5RQDnIgopxjgREQ5xQAnIsopBjgRUU4xwImIcooBTkSUUwxwIqKcYoATEeUUA5yIKKcY4EREOcUAJyLKKQY4EVFOMcCJiHKKAU5ElFMMcCKinGKAExHlFAOciCinGOBERDnFACciyikGOBFRTjHAiYhyigFORJRTDHAiopxigBMR5RQDnIgopxjgREQ5xQAnIsopBjgRUU4xwImIcooBTkSUUwxwIqKcYoATEeUUA5yIKKcY4EREOcUAJyLKKQY4EVFOMcCJiHKKAU5ElFMMcCKinGKAExHlFAOciCinGOBERDnFACciyikGOBFRTjHAiYhyigFORJRTDHAiopxigBMR5RQDnIgopxjgREQ5xQAnIsopBjgRUU5NFOAi8i9F5Dsi8oqIfCquRhER0fkuHOAiYgL4bwD+FYD3Afi4iLwvroYREdHZJumB/xCAV1T1lqp6AP4HgI/G0ywiIjrPJAF+HcCbx76+PTxGREQpsCb4WTnhmD5wI5GnATw9/NIVkW9NcJ95twJgN+tGZKzozwEff7EfPwCsiMi4z8GjJx2cJMBvA3j42NcPAbhz/41U9RkAzwCAiDyvqk9McJ+5VvTHD/A54OMv9uMH4n0OJhlC+QcA7xKRd4iIDeDfAPhqHI0iIqLzXbgHrqqBiHwSwP8CYAL4gqq+FFvLiIjoTJMMoUBVnwPw3Bg/8swk9zcDiv74AT4HfPwU23Mgqg9cdyQiohzgVHoiopxKJcA55R4QkXUReVFEvikiz2fdnqSJyBdEZPt42aiILIvI10Tke8OPS1m2MWmnPAefEZG3hufBN0XkySzbmCQReVhE/kZEXhaRl0Tkl4fHC3EenPH4YzsHEh9CGU65/y6An8Cg9PAfAHxcVb+d6B1PGRFZB/CEqhaiBlZE/gWADoA/UNX3D4/9BoB9Vf3s8A/5kqr+WpbtTNIpz8FnAHRU9XNZti0NInIVwFVVfUFE5gDcBPDTAP4dCnAenPH4P4aYzoE0euCccl9AqvoNAPv3Hf4ogC8OP/8iBifzzDrlOSgMVd1Q1ReGn7cBvIzBbO1CnAdnPP7YpBHgnHI/oAD+SkRuDmenFtGaqm4Ag5MbwGrG7cnKJ0XkH4dDLDM5fHA/EbkB4IMA/h4FPA/ue/xATOdAGgE+0pT7AvhhVf0QBqs3/tLw7TUVz+8AeCeAHwSwAeDzmbYmBSLSAPCnAH5FVVtZtydtJzz+2M6BNAJ8pCn3s05V7ww/bgP4CgZDS0WzNRwXPBof3M64PalT1S1VDVU1AvC7mPHzQERKGITXH6rql4eHC3MenPT44zwH0gjwwk+5F5H68CIGRKQO4CcBFHFRr68C+MTw808A+PMM25KJo+Aa+hnM8HkgIgLgWQAvq+pvHvtWIc6D0x5/nOdAKhN5hmUy/xVvT7n/9cTvdIqIyGMY9LqBwezXP5r150BE/hjAj2Kw+twWgP8M4M8AfAnAIwDeAPCzqjqzF/lOeQ5+FIO3zgpgHcAvHo0HzxoR+REAfwvgRQDR8PCnMRgHnvnz4IzH/3HEdA5wJiYRUU5xJiYRUU4xwImIcooBTkSUUwxwIqKcYoATEeUUA5yIKKcY4EREOcUAJyLKqf8PSSo5ZZRO+GcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot initial point\n",
    "# Add your solution here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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",
      "\n",
      "This is Ipopt version 3.13.2, running with linear solver ma27.\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:        0\n",
      "Number of nonzeros in inequality constraint Jacobian.:       30\n",
      "Number of nonzeros in Lagrangian Hessian.............:       12\n",
      "\n",
      "Total number of variables............................:        8\n",
      "                     variables with only lower bounds:        8\n",
      "                variables with lower and upper bounds:        0\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:        0\n",
      "Total number of inequality constraints...............:       15\n",
      "        inequality constraints with only lower bounds:        9\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        6\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0  1.0000000e+02 6.11e+01 1.02e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0\n",
      "   1  1.0059888e+02 1.59e+01 2.23e+01  -1.0 1.80e+01    -  1.67e-01 2.64e-01h  1\n",
      "   2  1.0278855e+02 0.00e+00 2.06e+00  -1.0 1.89e+01    -  1.42e-01 1.00e+00h  1\n",
      "   3  1.0025727e+02 0.00e+00 4.27e-01  -1.0 7.37e+00    -  9.58e-01 6.22e-01f  1\n",
      "   4  9.9403733e+01 0.00e+00 1.44e-01  -1.0 2.60e+01    -  7.40e-01 1.00e+00h  1\n",
      "   5  9.8843064e+01 0.00e+00 6.36e-03  -1.0 8.38e+00    -  1.00e+00 1.00e+00h  1\n",
      "   6  9.8390750e+01 0.00e+00 3.92e-04  -1.7 2.22e+00    -  1.00e+00 1.00e+00h  1\n",
      "   7  9.8284949e+01 0.00e+00 5.50e-05  -3.8 2.33e-01    -  1.00e+00 9.98e-01h  1\n",
      "   8  9.8284281e+01 0.00e+00 1.32e-07  -5.7 1.40e-01    -  1.00e+00 1.00e+00h  1\n",
      "   9  9.8284281e+01 0.00e+00 4.01e-10  -5.7 1.26e+01    -  1.00e+00 1.00e+00h  1\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "  10  9.8284281e+01 0.00e+00 3.54e-11  -5.7 4.54e-01    -  1.00e+00 1.00e+00h  1\n",
      "  11  9.8284281e+01 0.00e+00 1.84e-11  -5.7 2.74e-03    -  1.00e+00 1.00e+00h  1\n",
      "  12  9.8284270e+01 0.00e+00 1.43e-13  -8.6 2.69e-05    -  1.00e+00 1.00e+00h  1\n",
      "\n",
      "Number of Iterations....: 12\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:   9.8284270438748536e+01    9.8284270438748536e+01\n",
      "Dual infeasibility......:   1.4274501609950705e-13    1.4274501609950705e-13\n",
      "Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Complementarity.........:   2.5071318587575989e-09    2.5071318587575989e-09\n",
      "Overall NLP error.......:   2.5071318587575989e-09    2.5071318587575989e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 13\n",
      "Number of objective gradient evaluations             = 13\n",
      "Number of equality constraint evaluations            = 0\n",
      "Number of inequality constraint evaluations          = 13\n",
      "Number of equality constraint Jacobian evaluations   = 0\n",
      "Number of inequality constraint Jacobian evaluations = 13\n",
      "Number of Lagrangian Hessian evaluations             = 12\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      0.003\n",
      "Total CPU secs in NLP function evaluations           =      0.000\n",
      "\n",
      "EXIT: Optimal Solution Found.\n"
     ]
    }
   ],
   "source": [
    "# Solve the model\n",
    "# Add your solution here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAFoCAYAAACotWuNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqoklEQVR4nO3da4hk6X0e8Od/LnWvrr739Nx3ZiVLluyswqIEbILsbIxiArINMtoPYQOzrD9YIIM/WPiL5UBAhB07gQTDekd4DY4TgexIJCKxRtjIgiBrV6z2ol15b3Pv+63ut3P++VDVuz2z3dNVXXVubz0/aLqnprvrPX3qPPWe9yqqCiIiMoMVdQGIiGh8GOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAY5NtRFJCMi/yAiPxaR10XkD/uPz4rId0Tkrf7nmeCLS0REDyPHjVMXEQGQV9WqiLgAvg/gSwB+A8C2qn5VRL4MYEZVfy/wEhMR0ZGOralrT7X/T7f/oQA+B+CF/uMvAPi1IApIRESDcwb5JhGxAbwE4FEA/1VVfyAiS6q6AgCquiIii8f9nvn5eb148eIo5SUimjgvvfTSpqouDPK9A4W6qnoAHhORaQB/LSKfHLQwIvIMgGcA4Pz583jxxRcH/VEiIgIgIjcH/d6hRr+o6i6AvwPwWQBrIrLcf8JlAOtH/Mxzqvq4qj6+sDDQGw0REZ3QIKNfFvo1dIhIFsATAN4E8C0AT/W/7SkA3wyojERENKBBml+WAbzQb1e3AHxdVf+XiPw/AF8XkSsAbgH4fIDlJCKiARwb6qr6CoBPHfL4FoB/GUShiIjoZDijlIjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIAx1IiKDMNSJiAzCUCciMghDnYjIIMeGuoicE5G/FZE3ROR1EflS//GviMhdEXm5//GrwReXiIgexhnge7oAfldVfyQiRQAvich3+v/3x6r6bHDFIyKiYRwb6qq6AmCl/3VFRN4AcCboghER0fCGalMXkYsAPgXgB/2Hvigir4jI10RkZtyFIyKi4Qwc6iJSAPANAL+jqmUAfwLgMoDH0KvJXz3i554RkRdF5MWNjY3RS0xEREcaKNRFxEUv0P9CVf8KAFR1TVU9VfUB/CmATx/2s6r6nKo+rqqPLywsjKvcRER0iEFGvwiAawDeUNU/OvD48oFv+3UAr42/eERENIxBRr/8AoB/C+BVEXm5/9jvA3hSRB4DoABuAPitAMpHRERDGGT0y/cByCH/9e3xF4eIiEbBGaVERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAZhqBMRGYShTkRkEIY6EZFBGOpERAY5NtRF5JyI/K2IvCEir4vIl/qPz4rId0Tkrf7nmeCLS0REDzNITb0L4HdV9eMA/jmA3xaRnwXwZQDfVdWPAPhu/99ERBShY0NdVVdU9Uf9rysA3gBwBsDnALzQ/7YXAPxaQGUkIqIBDdWmLiIXAXwKwA8ALKnqCtALfgCLYy8dERENxRn0G0WkAOAbAH5HVcsiMujPPQPgGQCwbRuD/lzsiAWxXYjtAP3PYjuAWMD+Man2P/tQr9v/6AD7n9WPrvyUHJbzwevLdiCWC7H2X2fywevtQ687BaBQ3+u/5rpQ/4PXIMXLhQsXcOPGjbH/3oFCXURc9AL9L1T1r/oPr4nIsqquiMgygPXDflZVnwPwXP/3qO6/AGOs1fVQa3motbqotbuotTy0u6MHcsqxUEg7yKVs5NMO8mkbacceQ4kpSRptD7V2F422h7bno+P56HT1/a+DuEREANcWuLb1/kfasZB2LeRTvddkYitcCRXU3/vYUJfeM18D8Iaq/tGB//oWgKcAfLX/+ZuBlDAEvq8oNzvYqXewU2+j1QmmRt3u+tjutrFd++CxjGthOpfCTM7FVMaFZfHCMoWqotHxUG11UW/1P7c9eH74FRtVoN1VtLseAO9D/28JkEv1Khq9CoeDnGvz9ZhAclzNWUR+EcDfA3gVwH7a/T567epfB3AewC0An1fV7WN+V2xq6r6v2Kq1sV1rY6/RieRCe5BtCUpZF3OFFGZzKV5QCdPqethrdFBt9sK71uoiBi+rExMBcikbuZSDQtrBdM5FxuWd5biICAbNQxF5SVUfH+h7wwzZOIR6o+1hrdzEZrWFjhffK861BQvFNJamMryQYqzS7GC3f4dXa324BmyabMrGTM7FdC6FqYzDJpsRMNRHtF1rY2WvgXIjeR1GpayL5VIGM/lU1EWZeJ6v2K23sVPvYK/RRrsb34pB0BxbMJ3tBfx0zoVrc4L6MBjqJ7Rbb+P2dgPVVvLC/EHFjINzMzmUcm7URZko7a6PrVoLO7UOKs1OoptUgiKC95to5gtp3l0OgKE+pHKzg1tbdVSayQ/zB01lHZyfzaGYYbgHaa/ewVqlie1aO5ARKSYrZV0sTaUxm0+xieYIDPUBtbs+bm7VsFltB/o8cbBQTOHCXJ63vWPU8XysV1pYLzfRDGgU1CRJOYLFYgYLRdbeH8RQH8BGpYWbW7VYd4COm2sLLs7nMV9IR12URNtrdLBe7tXK2bwyfiLAdM7FYjGDmZzL2jsY6g/V6np4b7OGnVpn7L87KWbzKVycz3Ey0xA8X7FeaWKt3EKjbf7IlbhIORYW+yO7Us7k3mUy1I+wW2/jrfUquhNUOz+Kaws+slhkR+oxfF+xWm7i3m5jou7q4sa2BMulDJZLGTgT2ITIUD/EnZ067uw02Il1gAhwfjaH09PZqIsSO6qK9UoLd3YaY1n2gcbDsffDPQt7gibcMdQP8HzFOxtVbE1AZ+hJzRdSuLRQmKiL5Ciqio1qL8yDWgKCRufagjMzWSwVMxMxm5qh3tfqenhzpYI620CPlU/b+NipqYlut9yqtnB7p8E28wRJORbOzmSxWEwb3aHKUAfQ7Hj4yUqZta0hZFwLH1+emrjhZLv1Nm5t1ydi6r6pMq6FszM5LBTNHNk18aHeaHv4ycreRE/LPqmUY+ETpycj2Dteb57CRoVNc6aYzrl4ZD5v3Os3qFBPxH35fg2dgX4y7a6P1++V0eyYXWvdrLbw49u7DHTD7NY7eOXOHlb3mgOH4CSLfajvBxJHK4ym3fXxk5UyOp55f8dW18NPVyt4a63KIYqG8nzFe5s1vH6vzP6RY8Q61H1f8dPVCgN9TFodHz9drcA3aMrkermJV+7sYbvG2vkkqDS7eOXOLu7s1FlrP0KsQ/3tjaoRqyvGSaXZxbub1aiLMbJmx8NP7pXxzkaNE88mjK/A7e0GXr27x3w4RGxD/fZ2nePQA7JRaePOTj3qYpzY6l6vdr7XmNxlIQiotTy8dncPt7ZYaz8olqG+XWvjzk4j6mIY7fZ2A7v1ZL1per7irbUK3tusxWL7QYqeKnB3t2Fsf9FJxC7U210f724kv3kgCd7ZqCbmQmh2PLx+b28illSm4ZUbXTbH9MUu1N/bnKylc6PU7ipubNaiLsax9uodvHZ3jxOJ6KFaHR+v393DRqUVdVEiFatQX+/vMkPh2ay2sVmN70Vwd7eBN1bLfKOngfgKvL1exY3N2sS2s8cm1FtdDze3ktt5l2Q3NmuxGza6337e6wSLujSUNCt7zYltZ49NqN/ernNoWkQ6nuJ2jEbDsP2cxqHc6OKVO5PXzh6LUK+1uryAI7ZRicfuP+Um289pfNrdXjv7VoybGMctFqF+k7fYkVMFbm5H22m6W2/jzZUK289prHwF3lqvYr3cjLoooYg81HfrbU4iiYmdWieyc7FVbeGnqxWOP6dAqALvbNSwsmf+/JfIQ/32tvl/5CS5vR1+2/p6pYm31qtgnlPQbmzWI3mNhynSUC83OxPXiRF3lWY31HOyXm7infUam98oNHd2Grhl8Ei7SEN9bW8y2riSZjWk87JebuKdjfhPfiLz3N1tGFtjjyzUW10PW5xoFEtb1Vbg43sZ6BS1OztmBntkob5ebvGWO6Z8BdYCHCmwXmni3QQsT0DmMzHYowv1CV+fIe6COj87tTbe3WAbOsXHnZ1GoJWYsEUS6pVmJ3bT0ul+rY6P2pg7TOvtLt5arzLQKXbe26xhr27G0OpIQn2nZsYfz3TjXFyt4/l4k+PQKaZUgX9cr8RiVvWoIgn17YRtzjCpdsZ0nvb3mm11eHdG8dX1FG+uJn8RsNBDvdH2jHg3nAS1lodmZ/Rz9e5mFZUm5yNQ/DU7Pv5xrZLoZXtDD/XdBmvpSTLqsgF3durYqPCcU3KUG91Ej84KPdTH3flGwRpldulWtcVlICiR1sst3NtN5ms39FCvcknVRKmesNmk2uri7XXuNUvJdWu7nsid2CJpU6fkaHS8oUesdD0fP12tcIEuSjTtb403jn6lMIUb6hL5opA0JFWg1h6utn5jK37b4xGdhOcr3l6vJqrjlKFOx6oP0WS2VW2xY5SMUml2cS9Biw+GHOoS6tPReLS6g4V6u+vjvQSPGiA6yp3tOupD3rFGJdRQF4Z6Ig06GePdzSq3oiMj+f329SQ0w4TcHsJQT6LWAO3j65Uml38go9VaHu7sxH+YY8jNL6E+G43JcbXvZsfDTYN3kiHad3e3gUoz3pWXY0NdRL4mIusi8tqBx74iIndF5OX+x68O9nTmpfr169fx5JNfwC//8i/hySe/gOvXr0ddpLE7rvnlnY0qumx2oQmwv4G1H+PxuoPU1P8MwGcPefyPVfWx/se3B3s6s0L9+vXruHr1WayurkEVWF1dw9WrzxoX7A8bp76y10C5kYwOJKJxaLQ93IzxxhrHhrqqfg/AdghlSZxr155Hs3n/ZhLNZgvXrj0fUYmCcVTfULvrcxkAmkire81QN2gfxiht6l8UkVf6zTMzR32TiDwjIi+KyIvqx/OPcFJra2tDPZ5kh/X6396pc310mli3YtqPdNJQ/xMAlwE8BmAFwNWjvlFVn1PVx1X1cbGcEz5dPC0tLQ31eJI9OBy10fawwS0JaYLtNTrYjeHeECcKdVVdU1VPVX0Afwrg0+MtVjJcufI0Mpn0fY9lMmlcufJ0RCUKxmHTC25uc59Roptb9diNXT9R1VlEllV1pf/PXwfw2sO+/wPxOvhRPfHEEwB6betra2tYWlrClStPv/+4KWzr/lTfa3Q4Jp0IQL3tYaPawmIxE3VR3ndsqIvIXwL4DIB5EbkD4A8AfEZEHkMvpW8A+K3Bns6sUAd6wW5aiD/Ite+/oYtrWyJRFG5vNzCfT8Oy4jG679hQV9UnD3n42omeLWa3KTSY1IFQ36y2YtvrTxSFdtfHSrmJM9PZqIsCIKKNpylZUk6vBuL7itsxHp9LFJV7u43YbFgdaqjHrUOBBpOybQDAWqWJZiceL1yiOOl6Gpt1YcKtqTPUEynlWPB9xd2YvGiJ4mit3Bx4meoghRzqrOUlUS5tY7Pa4rK6RA+h2tuwOmoMdXooEaCQcrBaTs7OL0RRWSs3I1/sK/SO0mzKDvspaQS5lI1qu4vaEFvaEU2qjqfYrEVbWw891AtphnqS5NMOVhO0PyNR1Nb2Ji7U3bCfkkbg2oLtWvzWtyCKq2qrG+lGGqGH+nSOoZ4kna5y0BLRkKK8uw091DOujRzb1RMhl7KxE8NV6IjibqvWRnuAvX2DEMmM0tl8KoqnpWGJchgj0Qmo9kbCRCGSUJ9hqCdCu8NAJzqp9Uozkln0kYR6Ie0g5XDZmXhTdLmrEdGJtbuKnXr4HaaRJetiMX38N1FkXL7pEo0sipFjkV25S1OZQ3fUoehZAtjgySEa1W69HXoTTGShnnIszLFtPZaKWQfNiHruiUzS8RTlZrj7D0R6j32qFJ8toOgDGYdDTonGJezNqSMN9WLGRTFzom1SKSDFjMM104nGKOx29ch7w87N5KIuAh1wZjqDcoRTnIlM0+z4qLfDa4KJPNRLOZdLB8TEbD6Frs+9TIjGLczaeuShDgAX5nIcCRMxEeD8bI7LAhAFYKcW3t1vLEI9l3IwX+C49SgtFtNIOxZ2I5gsQWS6aqsb2lZ3sQh1oFdLdGxW16Pg2oJzszlUWl14nEVKFIi9kCpMsQn1lGPh4lw+6mJMpEfm83BtC9VWuONpiSZJWNdXbEIdABaKacwVOCEpTAvFFOb6TV81hjpRYCYy1IFerTHlsBkmDA/eHbGmThScetsLZVPq2IW6a1u4NF+IuhjGEwEeXSjAsXsvgY7no8VJR0SBUQVqIYxXj12oA7311s/NZqMuhtHOz+ZQOjA/gE0vRMGrtYIfARPLUAeAszM5LBTZvh6Exak0Tk/f/6ZZCXnRIaJJFEYTZ2xDHQAuzRe4NsyYTWUdXJr/8CijMG4LiSZdGHfEsQ51yxJ8dKmItBvrYiZGxrXw0aUi5JDpu2x+IQpeo+MFPhck9mmZciz87PIUt78bUdq18PHlKbj2h/+Ora6HdpeTjoiCFkZnaSKSMuPa+MRpBvtJpd3eG2PGPXyd9HoInTdE1BP0XXFiUnI/2DNsihlKNrX/dzt644sWdzkiCk3QQ4cTlZC9YC8hn+bOPIMopB184vQU0sfsZNRmqBOFpuMx1O+Tcix84nQJ81xO4KEWiml84vThbegPagf8IiOiDwR9Z5zI8YK2JfjIUhH5dAO3tuvc1OEAEeDiXH6o/V9ZUycKT9A19USG+r7T01nkUw7eWq+g4zHZU47g0cUiStnhdpIK+kVGRB8IuhKVuOaXB5VyLn7ubAmz+clujpkrpPBzZ6aHDnSAzS9EYfIV6AZ4zSW6pr4v7dj4mVNFbFZbuLlVm6gx1ylH8Mh84cRvar6v6PIuhyhUQVakjAj1ffOFNEpZFze3atiomL/X5kIxjYtzufdXWjwJ1tKJwtcJsOJpVKgDvaV7H10sYmmqg9vbDew1zNtzczrn4txsDoX06KePY9SJwtfygpvwZ1yo7ytmXPzsaRd7jQ5ub9eNWIWwmHFwfi6Hqczw7eZHYScpUfiC7Cw1NtT3lbIuSmdK2Km1sbLXTFzNXaR3DMulDKZz4+8MZns6UfiCXNTL+FDfN5NPYSafQrPjYb3cwnqlGethkClHsFjMYKGYfugU/1Ep4vs3IDJVkAs1HhvqIvI1AP8GwLqqfrL/2CyA/wHgIoAbAH5TVXeCK+b4ZFwb5+dyODuTxXa9je1aG3uNTixqrI4tmM66mM2nMJtPHbpE7rhx4hZR+DTAC2+QmvqfAfgvAP78wGNfBvBdVf2qiHy5/+/fG3/xgmNZgvlCGvOFNFQV5WYXu/U2duodNNrhrVqYTdmYybmYzqUwlXFCCfKDmOnJdP36dVy79jzW1tawtLSEK1eexhNPPBF1sWhAQV53x4a6qn5PRC4+8PDnAHym//ULAP4OCQv1g0Sk1/aedXFhrteJUW93UW11UW97qLa6Y1lZLeNayKed3kfKRj7tDLQ2S5CCrDFQMK5fv46rV59Fs9kCAKyuruHq1WcBgMGeEEFedidtU19S1RUAUNUVEVkcY5kil3IspJzUfR2Tvq9oez7ano9O10fHU3Q8H57fa5VWVYgIBL21aVzbgmsLXMdCyu59WFa4tfBBMNOT59q1598P9H3NZgvXrj3PUE+MBHeUisgzAJ4J+nmCZlmCjGUH2mkZhZBbe2gM1tbWhnqc4ii4C++k9/5rIrIMAP3P60d9o6o+p6qPq+rjJ3wuClDYbfg0uqWlpaEep/gJ8rI7aah/C8BT/a+fAvDN8RSHwsZIT54rV55GJpO+77FMJo0rV56OqEQ0rCCvu0GGNP4lep2i8yJyB8AfAPgqgK+LyBUAtwB8PsAyUoBYUU+e/XZzjn5JriDvkCXM0Q8iohxtES+re028t1mLuhhEE2W5lMEjC4WBR5+JyEuDNmEnfj11Go0dwxE5RKYL8rpjqE+4VMTj5IkmUdoJ7rrjFT3hUgG+uIjocEFOOuQVPeFcm80vRGELsjLFUJ9wjm2xXZ0oZKypU6BYWycKjwhr6hQwtqsThSfoRfx4NRNHwBCFKMiRLwBDncCaOlGYWFOnwEW9pjvRJAm6EsWrmQK/HSSiDzDUKXD59MTsP04UuXwq2D0ZGOqEjGvD4bBGolAEXYliqBMAIJ9ibZ0oaGnXYkcphaPAJhiiwIVxnTHUCQCQT5u19ypRHDHUKTTsLCUKXhjXGUOdAPQ6S7kGDFGwgh75AjDU6YBChrV1oqBkUzacECb68Soekecrau0uaq0u6m0Pnq/wVbG/9aAlAkt6szZzaRv5lINcyg5049mTyqcc7NQ6UReDyEiFkPqtGOpDanY8bNfaqLW6qLU9NDseht1L2xIgl3KQT9soZBzM5lKhvIMfhyNgiIITVr8Vr+IBqCp26h2slZvYrY9ek/UVqLa6qLa6WCu3cMOqY66QwtJUJtJgLWYciGDoNykiOl4p64byPAz1h2h2PGxUWlivtNDu+oE9j+cr1sstrJdbKKQdLE2lMVdIh74jkWNbmMq42GuwCYZonDKuhVxIE/wY6ofwfcXtnTpW9pqh11qrrS6qG13c2q7jwlweC8V0qM8/m08x1InGbDafCu25om/IjZm9Rgc/vrOLe7vhB/pBHU/x9noVb66W0ep6oT3vTD6cW0SiSTITYqizpt7X9Xzc2q5jrdyKuij32al1UG7s4cJcDktTmcCfL+3YKKQdVFvdwJ+LaBKkHEExxL4yhjp6tfO316uBtpuPwvMV727UsFlt4dHFAtJOsEOjZvIuQ51oTKZzqVCHME9888tmtYU3V8qxDfSDyo0uXr9XRqMdbHNMmO1/RKabzYV7PU10qK/uNfH2ehV+gobwtTo+Xr+3F2hNOpdykHEn+qVBNBa2JaENZdw3sVfuWrmJ9zZriRyT3fEUb6yUUQsw2FlbJxrddM6FFfLQ5IkM9fVKE+9u1KIuxki6/WCvt4MJ9jB764lMNRNy0wswgaFebnYSH+j7Op7izdUKut74+wOmMi6bYIhGYFuCmVz4Q4Qn6qr1fMU769VENrkcpdXxcWOrHsjvDmMIJZGpForpSNZ0mqhQv7lVQ7MT/1Euw9qotLBda4/99y4Ww1+qgMgUpyKqFE1MqO/W27GbWDRO721W0RlzM4xjW5gvsG2daFilrItsCBtiHGYiQr3r+XjHkHb0o7S7ivc2x3+Mp0psgiEa1nKE181EhPqt7XoiJheNaqvaHnszTC7lYCrLicdEg0q7FqYj6CDdZ3yodzwfGxVzm10edG+3MfbfuVzKjv13Epnq1FQm0p3NjA/19UorUTNGR1Vpdsc+KWkm5yLN4Y1Ex7ItwWLIy2U/yOgrVVWxuteMuhihWy2P95hFhMMbiQYwX4h+a0qjQ3271p6ItvQHbVZaYx8Js1hMg6MbiR4uDgMLjA71cddYk8LXXrPTOLm2xbZ1ooeYK6RC27LuYYwN9WbHQ7kxuWuCbwbQOXx6OgPXZnWd6EEiwPnZXNTFAGBwqAe5gmESNDoevDH3EDu2hTMzrK0TPWhpKoOMG81kowcZHOrh7esZR6oIZM31pWKGC30RHeDYgrMxquyMdHWKyA0ReVVEXhaRF8dVqHHgdmzB3K1YluBcTG4zieJguZSBG/GIl4PG0ar/S6q6OYbfM1a1gNYZT5KgmqDmC2ms7Db5xkkTL+VYOB2zAQTxeXsZo0bbQ9eboBlHRwgydM/PsbZOdG4mG/rORscZNdQVwN+IyEsi8sw4CjQOQe0GlDTNjj/2ztJ9payLmXx061sQRS2XsrEQ8ezRw4za/PILqnpPRBYBfEdE3lTV7x38hn7Yhxr4QQXZoK5fv45r157H2toalpaWcOXK03jiiSciKUvX92FbwfTKn5/NYbe+Z9SmI0SDOj+bi3SNl6OMVFNX1Xv9z+sA/hrApw/5nudU9XFVfXyU5xpGlJl+/fp1XL36LFZX16AKrK6u4erVZ3H9+vVIyhNk4OZSTmQbARBFaSbvxnYf3xOHuojkRaS4/zWAXwHw2rgKNgo/wqrjtWvPo9m8f+JPs9nCtWvPR1KeoP8W52dzkW0GQBQF1xZcmi9EXYwjjVJTXwLwfRH5MYB/APC/VfX/jKdYo4nyjmhtbW2ox4MmCPaPYVmCRxcLkf7NicL0yHweKSe+Y0xO3Kauqu8C+CdjLMvYWBEmzNLSElZXPxzgS0tLEZQmnDe4QtrBmeks7uyMfy13ojhZKKYwV4hf5+hB8X27GUGUtcYrV55GJnP/Sc9k0rhy5elIyhPWxtFnZ7IopKNfzIgoKCnHwsW5fNTFOJaRV2Haia6Nd3+USxxGv9iWhDbTTURweTGPV+/sTdSmJDQ5Li/kI18rfRCiIXYqioiG8Xyer/jhje2JH2o3lXXwidOlUJ/z3m4DN7fqoT4nUdCWptK4tDDezlERwaB5KCIvDTqCMP5vOydgW4JsTFZMi1IUzSHLpQyKGSNvAGlCZVwLFxLQ7LLPyFAHgDzbdyP5G4j0RsOE1ZZPFCQR4HLCXs/Ghjo77aL7G2RcG5cXklOzITrK+dkcpjLJWg7D2FDPpye7+cW1JdJF++cKaZybjdfqdUTDWJxK4/R08l7DxoZ6Ie0g5STnlmncpnPR1y7OzuSwUIznVGqih5nKOrg0n8y7TWNDXUSwWJzcdUmWYrImy6X5AjtOKVEyroWPLhVjuVjXIIwNdaB3+5TQ8zKSQtpBMSbtgJYl+JlTRaS5BR4lgGMLPnZqKlY7GQ0ruSUfQNqxMRvTldSCtFSK1zRm17bwsVNFOPYEvsNSYogAH10sJn6BOqNDHQBOleLRDBEW1xbM5+MV6kBvmd6PcOEvirFH5vMoxaAvalTGh/pUxp2okTCLxUzsttfaN51L4QK3waMYOlXKxKYfalTGhzoAXJhNZi/2sFKOYHk63i/M5VIWp2NeRposc4UULhpU2ZiIUC/l3Ilohrk0X0hEB8+FuTyWJ+B8UPzNFVL9ZsF43t2eRPwTYEzOz+aQMXgExuJUOrbbax3m4jyDnaI1mzcv0IEJCnXb4B160m4y1nl+0MX5/ETcQVH8zOZT+OiSeYEOTFCoA0Ax4+J0KXnTfo9zeSFZCw4d9Mh8nm3sFKq5grmBDkxYqAO9HXqmsubMcDw7k0Upm+xhWBfm8jg7Y96bLcXP4lTayCaXgyYu1C2rN2PMhKnrp6czODdrRq/9udkcLs6bcSwUT8ulDC4vmB3owASGOtBrX//YqWKil+c9VcokauH+QSyXsri8kDey34OidXYmi4sJXaBrWBMZ6gDg2BY+vlxMZFPM2ZksHjH0Bbo4lcHHT03B5ZICNAaWAI8uFoy5ox3ExIY60A/2U1OYKyRjKKBIr2PR9BdoKefik2dKEzUTmMYv7Vr45JkSForxWzYjSEZuPH0S6+Umbm7X0fXiWb5cysalhXxsVl8Mg+8r3t2sYqPSjroolDClrIuPLMV7Ml5QG08nr+0hIItTGZRyLm5s1rFdi0+IiABnprM4M52N7ZouQbEswaOLReTTDdzcqiOm9QGKmdPTGZyfzRnfIXoU1tQPsVlt4eZWDe1utGUtpB1cXswjl+J7716jg7fWKujE9E6KomdbgksLecwXktHcElRNnaF+hI7n495uAxuVVuhBkk3ZWC5lsFhMT2xt4zDNjoe31qqotrpRF4ViJu1a+JmlIvIJGtHGUI+I7yu2am2slZuoNIMLE0t6M90WpzKJ2708TL6veG+rhvVyK+qiUEzM5F1cXoh3+/lhGOoxUG93sV5uYavWRrvrj/z7RICsa2OhmMZCMZ24F2WUduttvLtZQ6sz+nmgZHJtwYW5fGJHtzDUY6bd9VFrdVFrd1Fve6i2ug8NGJHeCJZcykEh7SCXtpFPOYldsyUOPF9xa7uO1b1m1EWhkPXWQM8j5SS3IsRQTwDPV3R9H6qArwqBQASwROBYMnGjV8JSbnbw7kYNjbYXdVEoYClH8Mh8wYi9hzmkMQFsS2BbnDATtqmMi58/U8KdnQbu7TU49NFQi1NpXJjNwWEz5UMx1MkIliU4P5fDXCGFdzaqqLVYazdF2rVweb5gxKbQYWCok1HyaQc/d6aElb0m7u42YjtDmI5nSW/hurMzOfY9DYGhTsYREZyezmKxmMbKXhMre014PsM9KUSAxWIaZ2aySDtszhwWQ52M5dgWzs3mcKqUwd2dBtbKTTDb40sEmC+kcXYmi4zLMD8phjoZz7Wt3kbX071wX6+02JkaM3OFFM7N5JBNMcxHxVCniZF2bFxaKOD0dBZ3durYrLYZ7hGbybs4N5NL1PT+uONfkiZOxrXx6GIRp6e7uLfbxFa1xWaZEIkAM7kUlqe5JEYQOPmIJl7X87FRbWGt3OIEpgClXQuL/SUx2AHKyUdEgXFsC8ulLJZLWew1OtioNLFVbbP2PgYiwHTOxVIxg+mcy1VHQ8BQJzqglHVRyrq4MOdjo9LCWrmJJhcNG1rK6dXKF6dYKw8bQ53oEK5t4fR0Fqens9ird7BRbWGv0Y5845Q4c2zBdNbFXCGNGdbKI8NQJzpGKee+P0W90uxgt97BTr3NpQjQ29BlJudiOpfCVMZhkMcAQ51oCMWMi2LGxbnZHFpdD3v1DnbqHew1OhMxa1Wkt4DaTN7FTC7FSUIxNFKoi8hnAfxnADaA51X1q2MpFVECpB0bi1M2Fqcy8H1FuV+Lr7Z6a+ybEPKWALmUg3zaRinbq5FzHZZ4O/GQRhGxAfwjgH8F4A6AHwJ4UlV/8pCf4ZBGmgiqikbHQ63l3beZSpwXGLMtQS5lI592kO9/zqVsNqkEJI5DGj8N4G1Vfbf/pP8dwOcAHBnqRJNCRJBLOcilnPu2W2t2+iHf8tDoeOh4Ptqej07XD2UIpUivEzhlW3AdQcbph3jaRtZlgJtglFA/A+D2gX/fAfDPRisOkdkyro2Ma2Ou8OH/63h+76OraB8I+47no+srVAGFvr+0gWovpIHe5/2dtmxLesHtWHBt6QX4+x/C4DbcKKF+2CvjQ3UNEXkGwDP9f7ZE5LURnjPO5gFsRl2IAPH4ksvkYwMSfHwDvMHuH9uFQX/nKKF+B8C5A/8+C+Deg9+kqs8BeA4AROTFQduFksbkYwN4fElm8rEBZh/fSY5tlM3+fgjgIyLyiIikAHwBwLdG+H1ERDSiE9fUVbUrIl8E8H/RG9L4NVV9fWwlIyKioY00Tl1Vvw3g20P8yHOjPF/MmXxsAI8vyUw+NsDs4xv62EJdepeIiII1Sps6ERHFTCihLiKfFZGfisjbIvLlMJ4zTCJyQ0ReFZGXReTFqMszKhH5moisHxx+KiKzIvIdEXmr/3kmyjKe1BHH9hURuds/fy+LyK9GWcaTEpFzIvK3IvKGiLwuIl/qP27KuTvq+Ew5fxkR+QcR+XH/+P6w//hQ5y/w5peTLCeQNCJyA8DjqprIsbIPEpF/AaAK4M9V9ZP9x/4jgG1V/Wr/jXlGVX8vynKexBHH9hUAVVV9NsqyjUpElgEsq+qPRKQI4CUAvwbg38GMc3fU8f0mzDh/AiCvqlURcQF8H8CXAPwGhjh/YdTU319OQFXbAPaXE6CYUtXvAdh+4OHPAXih//UL6F1MiXPEsRlBVVdU9Uf9rysA3kBv5rcp5+6o4zOC9lT7/3T7H4ohz18YoX7YcgLGnIg+BfA3IvJSfwatiZZUdQXoXVwAFiMuz7h9UURe6TfPJLJ54iARuQjgUwB+AAPP3QPHBxhy/kTEFpGXAawD+I6qDn3+wgj1gZYTSLhfUNV/CuBfA/jt/i0+JcefALgM4DEAKwCuRlqaEYlIAcA3APyOqpajLs+4HXJ8xpw/VfVU9TH0Zuh/WkQ+OezvCCPUB1pOIMlU9V7/8zqAv0avyck0a/02zf22zfWIyzM2qrrWv5h8AH+KBJ+/flvsNwD8har+Vf9hY87dYcdn0vnbp6q7AP4OwGcx5PkLI9SNXk5ARPL9ThuISB7ArwAwcdGybwF4qv/1UwC+GWFZxmr/gun7dST0/PU72q4BeENV/+jAfxlx7o46PoPO34KITPe/zgJ4AsCbGPL8hTL5qD/E6D/hg+UE/kPgTxoSEbmEXu0c6M3Q/W9JPz4R+UsAn0Fvhbg1AH8A4H8C+DqA8wBuAfi8qiauw/GIY/sMerfuCuAGgN/ab8NMEhH5RQB/D+BVAH7/4d9Hr93ZhHN31PE9CTPO38+j1xFqo1fh/rqq/nsRmcMQ548zSomIDMIZpUREBmGoExEZhKFORGQQhjoRkUEY6kREBmGoExEZhKFORGQQhjoRkUH+PyJafAf3BSA1AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot solution\n",
    "# Add your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take Away Messages\n",
    "* Linear programs are convex. We will learn this means all local optima are global optima.\n",
    "* Nonlinear programs may be nonconvex. For nonconvex problems, there often existings many local optima that are not also global optima.\n",
    "* We will learn how to mathematically define convexity and analyze this property.\n",
    "* Initialization is really important in optimization problems with nonlinear objectives or constraints!\n",
    "* There are specialize solves for linear programs, quadratic programs, and convex programs. In this class, we will focus on more general algorithms for (non)convex nonlinear programs including the algorithms used by the `ipopt` solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}