optimization101/Modeling_Session_1/completed_modeling1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "89eeb5ad",
   "metadata": {},
   "source": [
    "# Translating a Decision Problem to an Optimization Model\n",
    "\n",
    "In the first video we discussed a few key concepts that are necessary for mathematical optimization:\n",
    "- parameters \n",
    "- decision variables\n",
    "- constraints\n",
    "- objective function\n",
    "\n",
    "In this first modeling example we will see how these are used to formulate a decision problem as an optimization model and code the formulation using `gurobipy`. For more information on all of the commands in the Python API check out our [documentation](https://www.gurobi.com/documentation/10.0/refman/py_python_api_details.html).\n",
    "\n",
    "## The Decision Problem\n",
    "We make widgets. Have a set of production facilities that produce boxes of widgets. There is also a set of distribution locations that will then distribute the widgets for sale. Each distribution center has a forecasted demand and each production facility has a min and max number of widgets it can make during this period. We need to ensure that each distribution facility receives enough widgets to satisfy demand from production and we want to do this at minimal cost. The minimum production is 75% of the production facilities max value. \n",
    "\n",
    "## Sets and Define Model\n",
    "Our sets are:\n",
    "- $P = \\{\\texttt{'Baltimore', 'Cleveland', 'Little Rock', 'Birmingham', 'Charleston'}\\} \\quad\\quad\\quad\\quad\\quad\\quad\\quad\\space\\space \\texttt{production}$\n",
    "- $D = \\{\\texttt{'Columbia', 'Indianapolis', 'Lexington', 'Nashville', 'Richmond', 'St. Louis'}\\} \\quad\\quad\\quad \\texttt{distribution}$\n",
    "\n",
    "To index each set, we'll use the lowercase letter of each set. Letters used for sets and indices are up to you. Typically, capital letters are for sets and corresponding lowercase will be the index. Single letters are used mainly for conciseness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "28e16593-5c56-43d6-b2fb-d05835d38ea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy\n",
    "\n",
    "# First, import packages\n",
    "import pandas as pd\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# Sets P and D, respectively\n",
    "# When we code sets we can be more descriptive in the name\n",
    "production = ['Baltimore','Cleveland','Little Rock','Birmingham','Charleston']\n",
    "distribution = ['Columbia','Indianapolis','Lexington','Nashville','Richmond','St. Louis']\n",
    "\n",
    "# Define a gurobipy model for the decision problem\n",
    "m = gp.Model('widgets')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "728aaa2d-b631-4ce3-9bd2-c619732a291c",
   "metadata": {},
   "source": [
    "## Parameters\n",
    "\n",
    "Parameters of a math optimization problem are values treated as constants in the model and are associated with the decision variables. For this decision problem these values are the limits of each production facility, the demand for each distribution center, and the pairwise costs between production and distribution locations. \n",
    "\n",
    "- $m_p$ is the max production in location $p$, $\\forall p \\in P \\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\space\\space \\texttt{max}\\_\\texttt{prod[p]}$\n",
    "- $n_d$ is the number of customers for a distribution center $d$, $\\forall d \\in D \\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad \\texttt{n}\\_\\texttt{demand[d]}$\n",
    "- $c_{p,d}$ is the cost to ship a widget between location $p$ and location $d$, $\\forall p \\in P, d \\in D \\quad\\quad\\quad \\texttt{cost[p,d]}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "924fa2e7-7f88-4f69-a5e1-74feed55adf5",
   "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>distribution</th>\n",
       "      <th>Columbia</th>\n",
       "      <th>Indianapolis</th>\n",
       "      <th>Lexington</th>\n",
       "      <th>Nashville</th>\n",
       "      <th>Richmond</th>\n",
       "      <th>St. Louis</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Baltimore</th>\n",
       "      <td>4.50</td>\n",
       "      <td>5.09</td>\n",
       "      <td>4.33</td>\n",
       "      <td>5.96</td>\n",
       "      <td>1.96</td>\n",
       "      <td>7.30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <td>3.33</td>\n",
       "      <td>4.33</td>\n",
       "      <td>3.38</td>\n",
       "      <td>1.53</td>\n",
       "      <td>5.95</td>\n",
       "      <td>4.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Charleston</th>\n",
       "      <td>3.02</td>\n",
       "      <td>2.61</td>\n",
       "      <td>1.61</td>\n",
       "      <td>4.44</td>\n",
       "      <td>2.36</td>\n",
       "      <td>4.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Cleveland</th>\n",
       "      <td>2.43</td>\n",
       "      <td>2.37</td>\n",
       "      <td>2.54</td>\n",
       "      <td>4.13</td>\n",
       "      <td>3.20</td>\n",
       "      <td>4.88</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <td>6.42</td>\n",
       "      <td>4.83</td>\n",
       "      <td>3.39</td>\n",
       "      <td>4.40</td>\n",
       "      <td>7.44</td>\n",
       "      <td>2.92</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "distribution  Columbia  Indianapolis  Lexington  Nashville  Richmond  \\\n",
       "production                                                             \n",
       "Baltimore         4.50          5.09       4.33       5.96      1.96   \n",
       "Birmingham        3.33          4.33       3.38       1.53      5.95   \n",
       "Charleston        3.02          2.61       1.61       4.44      2.36   \n",
       "Cleveland         2.43          2.37       2.54       4.13      3.20   \n",
       "Little Rock       6.42          4.83       3.39       4.40      7.44   \n",
       "\n",
       "distribution  St. Louis  \n",
       "production               \n",
       "Baltimore          7.30  \n",
       "Birmingham         4.01  \n",
       "Charleston         4.60  \n",
       "Cleveland          4.88  \n",
       "Little Rock        2.92  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Use .squeeze(\"columns\") to make the costs a series\n",
    "path = 'https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/Modeling_Session_1/'\n",
    "transp_cost = pd.read_csv(path + 'cost.csv', index_col=[0,1]).squeeze(\"columns\")\n",
    "# transp_cost = pd.read_csv('cost.csv', index_col=[0,1]).squeeze(\"columns\")\n",
    "# Pivot to view the costs a bit easier\n",
    "transp_cost.reset_index().pivot(index='production', columns='distribution', values='cost')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cf7664cc-4547-4b72-948f-a02e138c45d9",
   "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>max_production</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Baltimore</th>\n",
       "      <td>180</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Cleveland</th>\n",
       "      <td>200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <td>140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <td>80</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Charleston</th>\n",
       "      <td>180</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             max_production\n",
       "Baltimore               180\n",
       "Cleveland               200\n",
       "Little Rock             140\n",
       "Birmingham               80\n",
       "Charleston              180"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_prod = pd.Series([180,200,140,80,180], index = production, name = \"max_production\")\n",
    "n_demand = pd.Series([89,95,121,101,116,181], index = distribution, name = \"demand\") \n",
    "max_prod.to_frame()\n",
    "#n_demand.to_frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd659abc-ab18-49d4-b42e-43e43f19cccf",
   "metadata": {},
   "source": [
    "We also have the requirement that each production facility needs to produce at 75% of this maximum output. We'll denote this value by $a$ in the formulation \"frac\" for the fraction of maximum production required. Initially we set $a = 0.75$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dd77683c-47d1-4ee2-af3d-e9ce2003c209",
   "metadata": {},
   "outputs": [],
   "source": [
    "frac = 0.75"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cca98af8-c5c2-4883-8c48-f104ff798010",
   "metadata": {},
   "source": [
    "## Decision Variables\n",
    "This is what the optimization solver determines, which are the actions you have control over. As a reminder, they come in three main flavors:\n",
    "- `Continuous`: Price of a product\n",
    "- `Integer`: The number of food trucks to use for an event\n",
    "- `Binary`: Yes/no decision to include a certain stock in a portfolio\n",
    "\n",
    "Decision variables (and parameters) are indexed using elements of sets that we define for the problem. In this example, let's start with a set of cities that produce our widget, which we call set $P$ for the formulation but can define as 'production' in the code. And a set of cities that distribute the widget $D$ and 'distribution' similarly. The decision here is to determine the number of boxes to send from each production facility to each distribution location. \n",
    "\n",
    "Let $x_{p,d}$ be the number of widgets that are produced at facility $i$ and shipped to location $j$.\n",
    "\n",
    "### Add Variables in gurobipy\n",
    "`gurobipy` let's you add decision variables primarily with two (similar) commands:\n",
    "- [addVar()](https://www.gurobi.com/documentation/10.0/refman/py_model_addvar.html) adds a single variable\n",
    "- [addVars()](https://www.gurobi.com/documentation/10.0/refman/py_model_addvar.html) adds a group of variables by sets/indices\n",
    "\n",
    "When using `addVars` you have to provide the indices of the variables you want to add, which for us are the production and distribution locations. There are other arguments we can use and will cover a couple of them later on.  \n",
    "\n",
    "### Our Decision Variables\n",
    "As is often the case in writing code, there are several ways to get to the same point. Below we can see three different ways to create the decision variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ed5c511b-6da0-4edd-a647-7cf91c67aeda",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('Baltimore', 'Columbia'): <gurobi.Var Baltimore_to_Columbia>,\n",
       " ('Baltimore', 'Indianapolis'): <gurobi.Var Baltimore_to_Indianapolis>,\n",
       " ('Baltimore', 'Lexington'): <gurobi.Var Baltimore_to_Lexington>,\n",
       " ('Baltimore', 'Nashville'): <gurobi.Var Baltimore_to_Nashville>,\n",
       " ('Baltimore', 'Richmond'): <gurobi.Var Baltimore_to_Richmond>,\n",
       " ('Baltimore', 'St. Louis'): <gurobi.Var Baltimore_to_St. Louis>,\n",
       " ('Cleveland', 'Columbia'): <gurobi.Var Cleveland_to_Columbia>,\n",
       " ('Cleveland', 'Indianapolis'): <gurobi.Var Cleveland_to_Indianapolis>,\n",
       " ('Cleveland', 'Lexington'): <gurobi.Var Cleveland_to_Lexington>,\n",
       " ('Cleveland', 'Nashville'): <gurobi.Var Cleveland_to_Nashville>,\n",
       " ('Cleveland', 'Richmond'): <gurobi.Var Cleveland_to_Richmond>,\n",
       " ('Cleveland', 'St. Louis'): <gurobi.Var Cleveland_to_St. Louis>,\n",
       " ('Little Rock', 'Columbia'): <gurobi.Var Little Rock_to_Columbia>,\n",
       " ('Little Rock', 'Indianapolis'): <gurobi.Var Little Rock_to_Indianapolis>,\n",
       " ('Little Rock', 'Lexington'): <gurobi.Var Little Rock_to_Lexington>,\n",
       " ('Little Rock', 'Nashville'): <gurobi.Var Little Rock_to_Nashville>,\n",
       " ('Little Rock', 'Richmond'): <gurobi.Var Little Rock_to_Richmond>,\n",
       " ('Little Rock', 'St. Louis'): <gurobi.Var Little Rock_to_St. Louis>,\n",
       " ('Birmingham', 'Columbia'): <gurobi.Var Birmingham_to_Columbia>,\n",
       " ('Birmingham', 'Indianapolis'): <gurobi.Var Birmingham_to_Indianapolis>,\n",
       " ('Birmingham', 'Lexington'): <gurobi.Var Birmingham_to_Lexington>,\n",
       " ('Birmingham', 'Nashville'): <gurobi.Var Birmingham_to_Nashville>,\n",
       " ('Birmingham', 'Richmond'): <gurobi.Var Birmingham_to_Richmond>,\n",
       " ('Birmingham', 'St. Louis'): <gurobi.Var Birmingham_to_St. Louis>,\n",
       " ('Charleston', 'Columbia'): <gurobi.Var Charleston_to_Columbia>,\n",
       " ('Charleston', 'Indianapolis'): <gurobi.Var Charleston_to_Indianapolis>,\n",
       " ('Charleston', 'Lexington'): <gurobi.Var Charleston_to_Lexington>,\n",
       " ('Charleston', 'Nashville'): <gurobi.Var Charleston_to_Nashville>,\n",
       " ('Charleston', 'Richmond'): <gurobi.Var Charleston_to_Richmond>,\n",
       " ('Charleston', 'St. Louis'): <gurobi.Var Charleston_to_St. Louis>}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# loop through each p and d combination to create a decision variable\n",
    "m = gp.Model('widgets')\n",
    "x = {}\n",
    "for p in production:\n",
    "    for d in distribution:\n",
    "        x[p,d] = m.addVar(name = p+\"_to_\"+d)\n",
    "m.update()\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b66187a6-2b42-4d9b-8a96-6969da5150fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('Baltimore', 'Columbia'): <gurobi.Var prod_ship[Baltimore,Columbia]>,\n",
       " ('Baltimore', 'Indianapolis'): <gurobi.Var prod_ship[Baltimore,Indianapolis]>,\n",
       " ('Baltimore', 'Lexington'): <gurobi.Var prod_ship[Baltimore,Lexington]>,\n",
       " ('Baltimore', 'Nashville'): <gurobi.Var prod_ship[Baltimore,Nashville]>,\n",
       " ('Baltimore', 'Richmond'): <gurobi.Var prod_ship[Baltimore,Richmond]>,\n",
       " ('Baltimore', 'St. Louis'): <gurobi.Var prod_ship[Baltimore,St. Louis]>,\n",
       " ('Cleveland', 'Columbia'): <gurobi.Var prod_ship[Cleveland,Columbia]>,\n",
       " ('Cleveland', 'Indianapolis'): <gurobi.Var prod_ship[Cleveland,Indianapolis]>,\n",
       " ('Cleveland', 'Lexington'): <gurobi.Var prod_ship[Cleveland,Lexington]>,\n",
       " ('Cleveland', 'Nashville'): <gurobi.Var prod_ship[Cleveland,Nashville]>,\n",
       " ('Cleveland', 'Richmond'): <gurobi.Var prod_ship[Cleveland,Richmond]>,\n",
       " ('Cleveland', 'St. Louis'): <gurobi.Var prod_ship[Cleveland,St. Louis]>,\n",
       " ('Little Rock', 'Columbia'): <gurobi.Var prod_ship[Little Rock,Columbia]>,\n",
       " ('Little Rock',\n",
       "  'Indianapolis'): <gurobi.Var prod_ship[Little Rock,Indianapolis]>,\n",
       " ('Little Rock', 'Lexington'): <gurobi.Var prod_ship[Little Rock,Lexington]>,\n",
       " ('Little Rock', 'Nashville'): <gurobi.Var prod_ship[Little Rock,Nashville]>,\n",
       " ('Little Rock', 'Richmond'): <gurobi.Var prod_ship[Little Rock,Richmond]>,\n",
       " ('Little Rock', 'St. Louis'): <gurobi.Var prod_ship[Little Rock,St. Louis]>,\n",
       " ('Birmingham', 'Columbia'): <gurobi.Var prod_ship[Birmingham,Columbia]>,\n",
       " ('Birmingham',\n",
       "  'Indianapolis'): <gurobi.Var prod_ship[Birmingham,Indianapolis]>,\n",
       " ('Birmingham', 'Lexington'): <gurobi.Var prod_ship[Birmingham,Lexington]>,\n",
       " ('Birmingham', 'Nashville'): <gurobi.Var prod_ship[Birmingham,Nashville]>,\n",
       " ('Birmingham', 'Richmond'): <gurobi.Var prod_ship[Birmingham,Richmond]>,\n",
       " ('Birmingham', 'St. Louis'): <gurobi.Var prod_ship[Birmingham,St. Louis]>,\n",
       " ('Charleston', 'Columbia'): <gurobi.Var prod_ship[Charleston,Columbia]>,\n",
       " ('Charleston',\n",
       "  'Indianapolis'): <gurobi.Var prod_ship[Charleston,Indianapolis]>,\n",
       " ('Charleston', 'Lexington'): <gurobi.Var prod_ship[Charleston,Lexington]>,\n",
       " ('Charleston', 'Nashville'): <gurobi.Var prod_ship[Charleston,Nashville]>,\n",
       " ('Charleston', 'Richmond'): <gurobi.Var prod_ship[Charleston,Richmond]>,\n",
       " ('Charleston', 'St. Louis'): <gurobi.Var prod_ship[Charleston,St. Louis]>}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Provide each set for the indices \n",
    "m = gp.Model('widgets')\n",
    "x = m.addVars(production, distribution, name = 'prod_ship')\n",
    "m.update()\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4420b1b9-844b-4dc1-90f0-955938f5885c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('Baltimore', 'Columbia'): <gurobi.Var prod_ship[Baltimore,Columbia]>,\n",
       " ('Baltimore', 'Indianapolis'): <gurobi.Var prod_ship[Baltimore,Indianapolis]>,\n",
       " ('Baltimore', 'Lexington'): <gurobi.Var prod_ship[Baltimore,Lexington]>,\n",
       " ('Baltimore', 'Nashville'): <gurobi.Var prod_ship[Baltimore,Nashville]>,\n",
       " ('Baltimore', 'Richmond'): <gurobi.Var prod_ship[Baltimore,Richmond]>,\n",
       " ('Baltimore', 'St. Louis'): <gurobi.Var prod_ship[Baltimore,St. Louis]>,\n",
       " ('Cleveland', 'Columbia'): <gurobi.Var prod_ship[Cleveland,Columbia]>,\n",
       " ('Cleveland', 'Indianapolis'): <gurobi.Var prod_ship[Cleveland,Indianapolis]>,\n",
       " ('Cleveland', 'Lexington'): <gurobi.Var prod_ship[Cleveland,Lexington]>,\n",
       " ('Cleveland', 'Nashville'): <gurobi.Var prod_ship[Cleveland,Nashville]>,\n",
       " ('Cleveland', 'Richmond'): <gurobi.Var prod_ship[Cleveland,Richmond]>,\n",
       " ('Cleveland', 'St. Louis'): <gurobi.Var prod_ship[Cleveland,St. Louis]>,\n",
       " ('Little Rock', 'Columbia'): <gurobi.Var prod_ship[Little Rock,Columbia]>,\n",
       " ('Little Rock',\n",
       "  'Indianapolis'): <gurobi.Var prod_ship[Little Rock,Indianapolis]>,\n",
       " ('Little Rock', 'Lexington'): <gurobi.Var prod_ship[Little Rock,Lexington]>,\n",
       " ('Little Rock', 'Nashville'): <gurobi.Var prod_ship[Little Rock,Nashville]>,\n",
       " ('Little Rock', 'Richmond'): <gurobi.Var prod_ship[Little Rock,Richmond]>,\n",
       " ('Little Rock', 'St. Louis'): <gurobi.Var prod_ship[Little Rock,St. Louis]>,\n",
       " ('Birmingham', 'Columbia'): <gurobi.Var prod_ship[Birmingham,Columbia]>,\n",
       " ('Birmingham',\n",
       "  'Indianapolis'): <gurobi.Var prod_ship[Birmingham,Indianapolis]>,\n",
       " ('Birmingham', 'Lexington'): <gurobi.Var prod_ship[Birmingham,Lexington]>,\n",
       " ('Birmingham', 'Nashville'): <gurobi.Var prod_ship[Birmingham,Nashville]>,\n",
       " ('Birmingham', 'Richmond'): <gurobi.Var prod_ship[Birmingham,Richmond]>,\n",
       " ('Birmingham', 'St. Louis'): <gurobi.Var prod_ship[Birmingham,St. Louis]>,\n",
       " ('Charleston', 'Columbia'): <gurobi.Var prod_ship[Charleston,Columbia]>,\n",
       " ('Charleston',\n",
       "  'Indianapolis'): <gurobi.Var prod_ship[Charleston,Indianapolis]>,\n",
       " ('Charleston', 'Lexington'): <gurobi.Var prod_ship[Charleston,Lexington]>,\n",
       " ('Charleston', 'Nashville'): <gurobi.Var prod_ship[Charleston,Nashville]>,\n",
       " ('Charleston', 'Richmond'): <gurobi.Var prod_ship[Charleston,Richmond]>,\n",
       " ('Charleston', 'St. Louis'): <gurobi.Var prod_ship[Charleston,St. Louis]>}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The index of the tranporation costs have each combination of prodiction and distribution location\n",
    "m = gp.Model('widgets')\n",
    "x = m.addVars(transp_cost.index, name = 'prod_ship')\n",
    "m.update()\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94a42856-5458-4501-959e-c0265af517c2",
   "metadata": {},
   "source": [
    "The command `m.update()` updates the model to include any changes that have been made, like adding variables. It doesn't need to be run in every cell but if you see *Awaiting Model Update* in the output of a cell, then this should prevent that from happening."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dde6c86",
   "metadata": {},
   "source": [
    "## Constraints\n",
    "We outlined production and demand constraints at the beginning of this example; now we formulate and code them. Note that it doesn't matter the order in which constraints (and/or decision variables) are added to the model.\n",
    "\n",
    "### Add Constraints in gurobipy\n",
    "Adding constraints to a model is similar to adding variables:\n",
    "- [addConstr()](https://www.gurobi.com/documentation/10.0/refman/py_model_addconstr.html) adds a single constraint\n",
    "- [addConstrs()](https://www.gurobi.com/documentation/10.0/refman/py_model_addconstrs.htmll) adds a group of constraints using a Python `generator` expression\n",
    " \n",
    "### Our Constraints\n",
    "To start, we'll formulate the demand constraints for each distribution location first and add them to the model.\n",
    "\n",
    "\\begin{align*} \n",
    "\\sum_{p}x_{p,d} \\ge n_d, \\quad \\forall d \\in D \\quad\\quad \\texttt{meet}\\_\\texttt{demand[d]}\\\\ \n",
    "\\end{align*}\n",
    "\n",
    "This will be the first time we use [gp.quicksum()](https://www.gurobi.com/documentation/10.0/refman/py_quicksum.html). There are other ways to sum expressions in gurobipy and while this method isn't the most concise to code, it is easy to compare it to the summation in the formulation to see how it works. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5fee7987-ff5e-4b00-9e28-c93121c44ec3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Columbia': <gurobi.Constr meet_demand[Columbia]>,\n",
       " 'Indianapolis': <gurobi.Constr meet_demand[Indianapolis]>,\n",
       " 'Lexington': <gurobi.Constr meet_demand[Lexington]>,\n",
       " 'Nashville': <gurobi.Constr meet_demand[Nashville]>,\n",
       " 'Richmond': <gurobi.Constr meet_demand[Richmond]>,\n",
       " 'St. Louis': <gurobi.Constr meet_demand[St. Louis]>}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meet_demand = m.addConstrs((gp.quicksum(x[p,d] for p in production) >= n_demand[d] for d in distribution), name = 'meet_demand')\n",
    "#meet_demand = m.addConstrs(x.sum('*', j) >= demand[j] for j in distribution)\n",
    "m.update()\n",
    "meet_demand"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecca1f83-c111-47f1-9a34-6eeb2a1f63b6",
   "metadata": {},
   "source": [
    "Next we have the maximum number of widgets each production facility can make. We also have that each facility must make at least 75% of its max production. \n",
    "\n",
    "$$\n",
    "\\begin{align*} \n",
    "\\sum_{d}x_{p,d} &\\le m_p, &\\forall p \\in P \\quad\\quad &\\texttt{can}\\_\\texttt{produce[p]}\\\\ \n",
    "\\sum_{d}x_{p,d} &\\ge a*m_p,&\\forall p \\in P \\quad\\quad &\\texttt{must}\\_\\texttt{produce[p]}\\\\ \n",
    "\\end{align*}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d70dbb23-db8f-40fb-af49-1682071c6c62",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Baltimore': <gurobi.Constr can_produce[Baltimore]>,\n",
       " 'Cleveland': <gurobi.Constr can_produce[Cleveland]>,\n",
       " 'Little Rock': <gurobi.Constr can_produce[Little Rock]>,\n",
       " 'Birmingham': <gurobi.Constr can_produce[Birmingham]>,\n",
       " 'Charleston': <gurobi.Constr can_produce[Charleston]>}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "can_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')# (x.sum(i, '*')\n",
    "must_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "m.update()\n",
    "can_produce"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5be8505-9359-4004-8ac0-b54502db661e",
   "metadata": {},
   "source": [
    "## Objective Function\n",
    "We were told to **reduce** the transportation costs and we'll use this to determine our objective function as minimizing the total cost to ship widgets from production to distribution locations. \n",
    "\n",
    "### Setting the Objective in gurobipy\n",
    "This is done using [setObjective()](https://www.gurobi.com/documentation/10.0/refman/py_model_setobjective.html). The second argument (in this case `GRB.MINIMIZE`) is called the model's *sense*. For a maximization problem we would use `GRB.MAXIMIZE`. \n",
    "\n",
    "### Our Objective Function\n",
    "\\begin{align*} \n",
    "{\\rm minimize} \\space \\sum_{p,d}c_{p,d}x_{p,d}, \\quad \\forall p \\in P, d \\in D\\\\ \n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2756333c-6a92-4a3a-ba4d-9ec706d900f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.setObjective(gp.quicksum(transp_cost[i,j]*x[i,j] for i in production for j in distribution), GRB.MINIMIZE) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b77c7b0-ddd1-48de-8d41-71a2e0fd77f3",
   "metadata": {},
   "source": [
    "## Find, Extract, and Analyze the Solution\n",
    "Before running the optimization, it is a good idea to write an `lp` file. This is a text file that prints out the variables, constraints, and object like we would see in the *formulation*, just without the summation symbols and using the names we designated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f2f93cde-16ae-4e1d-a9ab-835761587e7f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: variable name \"prod_ship[Baltimore,St. Louis]\" has a space\n",
      "Warning: constraint name \"meet_demand[St. Louis]\" has a space\n",
      "Warning: to let Gurobi read it back, use rlp format\n"
     ]
    }
   ],
   "source": [
    "m.write('widget_shipment.lp')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12763405-b27f-4a6d-897f-db30df211caa",
   "metadata": {},
   "source": [
    "### Run the Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "119fc40e-d1fc-4d13-b227-90255fd806ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 10.0.2 build v10.0.2rc0 (win64)\n",
      "\n",
      "CPU model: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "\n",
      "Optimize a model with 16 rows, 30 columns and 90 nonzeros\n",
      "Model fingerprint: 0x20186c14\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [2e+00, 7e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [6e+01, 2e+02]\n",
      "Presolve removed 5 rows and 0 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 11 rows, 35 columns, 65 nonzeros\n",
      "\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0    0.0000000e+00   1.610000e+02   0.000000e+00      0s\n",
      "      15    1.7048900e+03   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 15 iterations and 0.01 seconds (0.00 work units)\n",
      "Optimal objective  1.704890000e+03\n"
     ]
    }
   ],
   "source": [
    "m.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ee05c74-fe8d-445f-9ffe-4a25d123a17b",
   "metadata": {},
   "source": [
    "### Extract the Solution\n",
    "There are many ways to get the values of decision variables out of gurobipy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "69d92eda-0937-4aff-8e60-ff012d766fa0",
   "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></th>\n",
       "      <th>cost</th>\n",
       "      <th>shipment</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Baltimore</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>5.96</td>\n",
       "      <td>19.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>116.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.37</td>\n",
       "      <td>95.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nashville</th>\n",
       "      <td>4.13</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>1.53</td>\n",
       "      <td>80.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Charleston</th>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.60</td>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost  shipment\n",
       "production  distribution                \n",
       "Baltimore   Nashville     5.96      19.0\n",
       "            Richmond      1.96     116.0\n",
       "Cleveland   Columbia      2.43      89.0\n",
       "            Indianapolis  2.37      95.0\n",
       "            Nashville     4.13       2.0\n",
       "Little Rock St. Louis     2.92     140.0\n",
       "Birmingham  Nashville     1.53      80.0\n",
       "Charleston  Lexington     1.61     121.0\n",
       "            St. Louis     4.60      41.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_values = pd.Series(m.getAttr('X', x), name = \"shipment\", index = transp_cost.index)\n",
    "sol = pd.concat([transp_cost, x_values], axis=1)\n",
    "#sol \n",
    "sol[sol.shipment > 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09af17eb-34c3-44f0-a65b-5835fd313f24",
   "metadata": {},
   "source": [
    "Here are a couple of other ways to get the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4143dc50-9073-4131-896f-ec7261160c0c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'prod_ship[Baltimore,Columbia]': 0.0,\n",
       " 'prod_ship[Baltimore,Indianapolis]': 0.0,\n",
       " 'prod_ship[Baltimore,Lexington]': 0.0,\n",
       " 'prod_ship[Baltimore,Nashville]': 19.0,\n",
       " 'prod_ship[Baltimore,Richmond]': 116.0,\n",
       " 'prod_ship[Baltimore,St. Louis]': 0.0,\n",
       " 'prod_ship[Cleveland,Columbia]': 89.0,\n",
       " 'prod_ship[Cleveland,Indianapolis]': 95.0,\n",
       " 'prod_ship[Cleveland,Lexington]': 0.0,\n",
       " 'prod_ship[Cleveland,Nashville]': 2.0,\n",
       " 'prod_ship[Cleveland,Richmond]': 0.0,\n",
       " 'prod_ship[Cleveland,St. Louis]': 0.0,\n",
       " 'prod_ship[Little Rock,Columbia]': 0.0,\n",
       " 'prod_ship[Little Rock,Indianapolis]': 0.0,\n",
       " 'prod_ship[Little Rock,Lexington]': 0.0,\n",
       " 'prod_ship[Little Rock,Nashville]': 0.0,\n",
       " 'prod_ship[Little Rock,Richmond]': 0.0,\n",
       " 'prod_ship[Little Rock,St. Louis]': 140.0,\n",
       " 'prod_ship[Birmingham,Columbia]': 0.0,\n",
       " 'prod_ship[Birmingham,Indianapolis]': 0.0,\n",
       " 'prod_ship[Birmingham,Lexington]': 0.0,\n",
       " 'prod_ship[Birmingham,Nashville]': 80.0,\n",
       " 'prod_ship[Birmingham,Richmond]': 0.0,\n",
       " 'prod_ship[Birmingham,St. Louis]': 0.0,\n",
       " 'prod_ship[Charleston,Columbia]': 0.0,\n",
       " 'prod_ship[Charleston,Indianapolis]': 0.0,\n",
       " 'prod_ship[Charleston,Lexington]': 121.0,\n",
       " 'prod_ship[Charleston,Nashville]': 0.0,\n",
       " 'prod_ship[Charleston,Richmond]': 0.0,\n",
       " 'prod_ship[Charleston,St. Louis]': 41.0}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# You can get the name and value of all the decision variables:\n",
    "all_vars = {v.varName: v.x for v in m.getVars()} \n",
    "all_vars"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80e2541e-d147-477c-af49-965e6aa835d0",
   "metadata": {},
   "source": [
    "Or you can only iterate over a specific variable and only return values that are of interest to you. Remember, x is a dict in python. So, iterate over it, the same way you iterate over any dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8b2613a8-2239-4e81-af59-86420dbafe91",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('Baltimore', 'Nashville'): 19.0,\n",
       " ('Baltimore', 'Richmond'): 116.0,\n",
       " ('Cleveland', 'Columbia'): 89.0,\n",
       " ('Cleveland', 'Indianapolis'): 95.0,\n",
       " ('Cleveland', 'Nashville'): 2.0,\n",
       " ('Little Rock', 'St. Louis'): 140.0,\n",
       " ('Birmingham', 'Nashville'): 80.0,\n",
       " ('Charleston', 'Lexington'): 121.0,\n",
       " ('Charleston', 'St. Louis'): 41.0}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xvals = {k: v.x for k,v in x.items() if v.x > 0} \n",
    "xvals "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7572f7c5-41ac-46fa-a2ff-a1a242cd7b98",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(703.0, 703)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_values.sum(), n_demand.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bde35449-0b63-4bb4-9926-d20013b1b95b",
   "metadata": {},
   "source": [
    "### Solution Analysis\n",
    "While determining the optimal transportation of widgets was our goal, we may want to dig a little deeper into the solution. For example we can aggregate the total production by facility to see which locations (if any) did not produce their maximum capacity of widgets and which (if any) production facilities are at the lower bound of their production. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "53a42bc7-1e45-4245-84b9-f648d7e74fa5",
   "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>Remaining</th>\n",
       "      <th>Utilization</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Baltimore</th>\n",
       "      <td>45.0</td>\n",
       "      <td>0.75</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Charleston</th>\n",
       "      <td>18.0</td>\n",
       "      <td>0.90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Cleveland</th>\n",
       "      <td>14.0</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Remaining  Utilization\n",
       "Baltimore         45.0         0.75\n",
       "Birmingham         0.0         1.00\n",
       "Charleston        18.0         0.90\n",
       "Cleveland         14.0         0.93\n",
       "Little Rock        0.0         1.00"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sum the shipment amount by production facility\n",
    "ship_out = sol.groupby('production')['shipment'].sum()\n",
    "pd.DataFrame({'Remaining':max_prod - ship_out, 'Utilization':ship_out/max_prod})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "105d6054-67e4-4237-b00a-705b9d34a153",
   "metadata": {},
   "source": [
    "In mathematical optimization, when the left-hand and right-hand sides of an inequality constraint are equal, we say the constraint is `binding`. When this *doesn't happen* then there is `slack` or `surplus` in that constraint. We can get this value by calling the `Slack` attribute of a constraint. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "30213e7c-bc28-4d31-be7c-b57925ee68b1",
   "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>Remaining</th>\n",
       "      <th>Utilization</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Baltimore</th>\n",
       "      <td>45.0</td>\n",
       "      <td>0.75</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Cleveland</th>\n",
       "      <td>14.0</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Charleston</th>\n",
       "      <td>18.0</td>\n",
       "      <td>0.90</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Remaining  Utilization\n",
       "Baltimore         45.0         0.75\n",
       "Cleveland         14.0         0.93\n",
       "Little Rock        0.0         1.00\n",
       "Birmingham         0.0         1.00\n",
       "Charleston        18.0         0.90"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({'Remaining':[can_produce[p].Slack for p in production], \n",
    "              'Utilization':[1-can_produce[p].Slack/max_prod[p] for p in production]}, \n",
    "             index = production)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f2d53b5-54bc-48c7-97e6-26c1c0f329b0",
   "metadata": {},
   "source": [
    "# Using Binary Variables\n",
    "As we described in the first session and also at the top of this notebook, binary variables are used to choose alternatives in mathematical optimization. They can be interpreted as a yes/no decision or an on/off switch. \n",
    "\n",
    "In the original problem Birmingham's production is much lower than the rest of the facilities. Suppose we have the option to expand that facilities max capacity by either 25 or 50 widgets, but there is a cost of \\\\$50 and \\\\$75, respectively, to choose one of these options and we can choose at most one. We'll use a binary decision variable for each option named $xprod$. \n",
    "\n",
    "Let $xprod_0 = 1$ if we choose the first option and expand production capacity by 25 and $0$ otherwise.\n",
    "Let $xprod_1 = 1$ if we choose the second option and expand production capacity by 50 and $0$ otherwise.\n",
    "\n",
    "While it's fairly common to use single lowercase letters as decision variables, it is not necessary and you'll see variables defined as above (where they are more descriptive) quite often. We will formulate a new model that contains the same decision variables and demand constraints as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "dd9a333d-ae1d-4b29-a69c-760a7414ef59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We use m2 for the second model\n",
    "# These parts are the same as above outside of the new model name\n",
    "m2 = gp.Model('widgets2')\n",
    "x = m2.addVars(production, distribution, obj = transp_cost, name = 'prod_ship')\n",
    "meet_demand = m2.addConstrs((gp.quicksum(x[p,d] for p in production) >= n_demand[d] for d in distribution), name = 'meet_demand')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bc726e6-7982-46bd-9568-925351b22306",
   "metadata": {},
   "source": [
    "In the cell above we did use new argument of the `addVars()` function: `obj`. This will set the coefficient of the added decision variables in the objective function and is equivalent to what we did earlier by attaching the transportation costs between each production and distribution location to the appropriate decision variable.  \n",
    "\n",
    "Next, we'll add the same constraints for production limits as before for each production facility other than Birmingham. The formulation is basically the same other than the set the constraints hold for. \n",
    "$$\n",
    "\\begin{align*} \n",
    "\\sum_{d}x_{p,d} &\\le m_p, &\\forall p \\in P -\\{\\texttt{Birmingham}\\} \\\\ \n",
    "\\sum_{d}x_{p,d} &\\ge a*m_p,&\\forall p \\in P -\\{\\texttt{Birmingham}\\} \\\\ \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "In gurobipy, this is done by adding a condition in the `generator` expression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "38e30563-b8c8-454b-9b61-e6ef6fa5fe83",
   "metadata": {},
   "outputs": [],
   "source": [
    "can_produce = m2.addConstrs((gp.quicksum(x[p,d] for d in distribution) <= max_prod[p] for p in production if p != 'Birmingham'), name = 'can_produce')\n",
    "must_produce =  m2.addConstrs((gp.quicksum(x[p,d] for d in distribution) >= frac*max_prod[p] for p in production if p != 'Birmingham'), name = 'must_produce')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed539f12-4ce2-4de6-820f-076e45ac0a73",
   "metadata": {},
   "source": [
    "Now, add the new binary variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "5eda68b8-6f7c-4245-8553-42f41fff5a7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "xprod = m2.addVars(range(2), vtype = GRB.BINARY, obj = [50, 75], name = 'expand_Birmingham_prod')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf53d755-2ff2-46e3-a38a-2c8fb1434c64",
   "metadata": {},
   "source": [
    "Let's breakdown each of the arguments in the cell above -- there are a few new things there.\n",
    "1. `range(n)` is used to add $n$ decision variables. In this case we add 2 variables.\n",
    "2. We need to declare this as a binary variable using `vtype`.\n",
    "3. We again use the `obj` capability to immediately set the objective function coefficients for these variables.\n",
    "\n",
    "The objective and new binary variables look like this in the formulation: \n",
    "\n",
    "\\begin{align*} \n",
    "{\\rm minimize} \\space &\\sum_{p,d}c_{p,d}x_{p,d} + 50*xprod_0 + 75*xprod_1, \\quad &\\forall p \\in P, d \\in D\\\\ \n",
    "&xprod_i \\in \\{0,1\\}, &{\\rm for} \\space i \\in \\{0,1\\}\n",
    "\\end{align*}\n",
    "\n",
    "Next we have the production constraints that are specific to the Birmingham facility.\n",
    "\n",
    "$$\n",
    "\\begin{align*} \n",
    "\\sum_{d}x_{p,d} &\\le m_p + 25*xprod_0 + 50*xprod_1, & p = \\texttt{Birmingham} \\\\ \n",
    "\\sum_{d}x_{p,d} &\\ge a*(m_p+ 25*xprod_0 + 50*xprod_1),& p = \\texttt{Birmingham} \\\\ \n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "487ee3df-6cf0-421a-95dc-54440248d81f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Birmingham_max = m2.addConstr(gp.quicksum(x['Birmingham',d] for d in distribution) <= max_prod['Birmingham'] + 25*xprod[0] + 50*xprod[1], name = 'Birmingham_add_max')\n",
    "Birmingham_min = m2.addConstr(gp.quicksum(x['Birmingham',d] for d in distribution) >= frac*(max_prod['Birmingham'] + 25*xprod[0] + 50*xprod[1]), name = 'Birmingham_add_min')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ed63a6b-ef71-4eb6-b237-637096518f0a",
   "metadata": {},
   "source": [
    "It was stated above that we can select at most one of the expansion options which means we cannot allow both $xprod_0$ and $xprod_1$ to equal one. To model this we add a constraint limiting the sum of these two binary variables to at most one.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{i}xprod_i \\le 1\n",
    "\\end{align*}\n",
    "$$\n",
    "The corresponding constraint in gurobipy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "a48175c5-6b22-402f-a207-9f590b7dc5cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "Birmingham_lim = m2.addConstr(gp.quicksum(xprod[i] for i in range(2)) <= 1, name = 'expansion_choice')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14e96d35-1d54-4c65-be2c-0b4c6d26d964",
   "metadata": {},
   "source": [
    "Now we can run this optimization model and see if this potential expansion will help us reduce overall costs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "fb182d4c-c533-4cc1-a559-05c07ef18975",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 10.0.2 build v10.0.2rc0 (win64)\n",
      "\n",
      "CPU model: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "\n",
      "Optimize a model with 17 rows, 32 columns and 96 nonzeros\n",
      "Model fingerprint: 0x80ef8db6\n",
      "Variable types: 30 continuous, 2 integer (2 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 5e+01]\n",
      "  Objective range  [2e+00, 8e+01]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 2e+02]\n",
      "Presolve time: 0.00s\n",
      "Presolved: 17 rows, 32 columns, 96 nonzeros\n",
      "Variable types: 30 continuous, 2 integer (2 binary)\n",
      "Found heuristic solution: objective 1704.8900000\n",
      "\n",
      "Root relaxation: objective 1.684260e+03, 14 iterations, 0.00 seconds (0.00 work units)\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 1684.26000    0    1 1704.89000 1684.26000  1.21%     -    0s\n",
      "H    0     0                    1700.4000000 1684.26000  0.95%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Implied bound: 1\n",
      "  MIR: 1\n",
      "\n",
      "Explored 1 nodes (14 simplex iterations) in 0.02 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 2: 1700.4 1704.89 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 1.700400000000e+03, best bound 1.700400000000e+03, gap 0.0000%\n"
     ]
    }
   ],
   "source": [
    "m2.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "34d94ec0-6ab4-4b68-8181-0c432595642e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The original model had a total cost of 1704.89\n",
      "The new formualtion has a total cost of 1700.4\n"
     ]
    }
   ],
   "source": [
    "obj1 = m.getObjective()\n",
    "obj2 = m2.getObjective()\n",
    "print(f\"The original model had a total cost of {round(obj1.getValue(),2)}\")\n",
    "print(f\"The new formualtion has a total cost of {round(obj2.getValue(),2)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6813d599-d8c1-49c3-8f09-f495c3b1af16",
   "metadata": {},
   "source": [
    "What does the change in objective function value tell us? \n",
    "\n",
    "Let's look at values of our binary variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "3177169c-67b9-4db3-804d-8b17b32fe70d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    1.0\n",
       "1    0.0\n",
       "dtype: float64"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.Series(m2.getAttr('X', xprod))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3410274c-55fd-4724-9f76-92cfc8178647",
   "metadata": {},
   "source": [
    "The model selected the first expansion option since $xprod_0 = 1$, which was increasing production by 25 widgets in Birmingham. We can see the rest of the solution, which will include the increase in Birmingham's production capacity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "05f918a8-ed81-4210-a341-c2d401f1276c",
   "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></th>\n",
       "      <th>cost</th>\n",
       "      <th>shipment</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Baltimore</th>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>135.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.37</td>\n",
       "      <td>95.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>1.53</td>\n",
       "      <td>101.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.01</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Charleston</th>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.60</td>\n",
       "      <td>37.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost  shipment\n",
       "production  distribution                \n",
       "Baltimore   Richmond      1.96     135.0\n",
       "Cleveland   Columbia      2.43      89.0\n",
       "            Indianapolis  2.37      95.0\n",
       "Little Rock St. Louis     2.92     140.0\n",
       "Birmingham  Nashville     1.53     101.0\n",
       "            St. Louis     4.01       4.0\n",
       "Charleston  Lexington     1.61     121.0\n",
       "            St. Louis     4.60      37.0"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x2_values = pd.Series(m2.getAttr('X', x), name = \"shipment\", index = transp_cost.index)\n",
    "sol2 = pd.concat([transp_cost, x2_values], axis=1)\n",
    "sol2[sol2.shipment > 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "3f750129",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Freeing default Gurobi environment\n"
     ]
    }
   ],
   "source": [
    "m.dispose()\n",
    "m2.dispose()\n",
    "gp.disposeDefaultEnv()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90a9b8d1-fe39-410a-a551-4ddb165ffd5a",
   "metadata": {},
   "source": [
    "### Homework! (not really, but something to look into)\n",
    "Analyze how the optimal solution changes between the two models. You'll notice something weird. \n",
    "- What is it that's odd?\n",
    "- Why do you think it happened?\n",
    "- How would you address it from formulation perspective and a business perspective?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "vscode": {
   "interpreter": {
    "hash": "70d79b410022258e996d14a47d5687c6c473c4d5d0c04ea47c92533eee9f501f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}