LP01_lp-intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dd0HlxDDOJ93"
   },
   "source": [
    "# <center>Linear programming introduction</center>\n",
    "### <center>Alfred Galichon (NYU & Sciences Po)</center>\n",
    "## <center>'math+econ+code' masterclass series</center>\n",
    "#### <center>With python code examples</center>\n",
    "© 2018–2023 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed [here](http://www.math-econ-code.org/team).\n",
    "\n",
    "**If you reuse material from this masterclass, please cite as:**<br>\n",
    "Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kTSWIALlOJ98"
   },
   "source": [
    "## Linear programming: duality\n",
    "\n",
    "### Learning objectives\n",
    "\n",
    "* Linear programming duality\n",
    "\n",
    "* Economic interpretation of the dual\n",
    "\n",
    "* Numerical computation\n",
    "\n",
    "### References\n",
    "\n",
    "* Galichon, *Optimal Transport Methods in Economics*. Appendix B.\n",
    "\n",
    "* Stigler (1945), The cost of subsistence. *Journal of Farm Economics*.\n",
    "\n",
    "* Dantzig (1990), The diet problem. *Interface*.\n",
    "\n",
    "* Complements:\n",
    "\n",
    "    * Gale (1960), *The theory of linear economic models*.\n",
    "\n",
    "    * Vohra (2011), *Mechanism Design: A Linear Programming Approach*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bKgWhZyKOJ98"
   },
   "source": [
    "\n",
    "### The diet problem\n",
    "\n",
    "During World War II, engineers in US Army were wondering how to feed their personnel at minimal cost, leading to what is now called the **optimal diet problem**.\n",
    "\n",
    "* Nutritionists have identified a number of vital nutrients (calories, protein, calcium, iron, etc.) that matter for a person's health, and have determined the minimum daily intake of each nutrient.\n",
    "\n",
    "* For each basic food (pasta, butter, bread, etc), nutritionists have characterized the intake in each of the various nutrients.\n",
    "\n",
    "* Each food has a unit cost, and the problem is to find the optimal diet, i.e. the combination of foods that meets the minimal intake in each of the nutrients and achieves minimal cost.\n",
    "\n",
    "The problem was taken on by G. Stigler, who published a paper about it in 1945, giving a first heuristic solution, exhibiting a diet that costs \\\\$ 39.93 per year in 1939 dollars. Later (in 1947) it was one of the first\n",
    "application of G.B. Dantzig's method (the simplex algorithm), which provided the exact solution (\\$39.67). It then took 120 man-day to perform this operation. However today the computer will perform it for us in a\n",
    "fraction of second.\n",
    "\n",
    "However, don't try this diet at home! Dantzig did so and almost died from it..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DZJ1PauuOJ99"
   },
   "source": [
    "## Motivation\n",
    "\n",
    "### A look at the Data\n",
    "\n",
    "Our dataset was directly taken from Stigler's article. It is a .csv file called `StiglerData1939.txt'. Have a look at the diet problem data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "from tabulate import tabulate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['1. Wheat Flour (Enriched)', '2. Macaroni',\n",
       "       '3. Wheat Cereal (Enriched)', '4. Corn Flakes', '5. Corn Meal',\n",
       "       '6. Hominy Grits', '7. Rice', '8. Rolled Oats',\n",
       "       '9. White Bread (Enriched)', '10. Whole Wheat Bread',\n",
       "       '11. Rye Bread', '12. Pound Cake', '13. Soda Crackers', '14. Milk',\n",
       "       '15. Evaporated Milk (can)', '16. Butters', '17. Oleomargarine',\n",
       "       '18. Eggs', '19. Cheese (Cheddar)', '20. Cream',\n",
       "       '21. Peanut Butter', '22. Mayonnaise', '23. Crisco', '24. Lard',\n",
       "       '25. Sirloin Steak', '26. Round Steak', '27. Rib Roast',\n",
       "       '28. Chuck Roast', '29. Plate', '30. Liver (Beef)',\n",
       "       '31. Leg of Lamb', '32. Lamb Chops', '33. Pork Chops',\n",
       "       '34. Pork Loin Roast', '35. Bacon', '36. Ham-smoked',\n",
       "       '37. Salt Pork', '38. Roasting Chicken ', '39. Veal Cutlets',\n",
       "       '40. Salmon Pink (can)', '41. Apples', '42. Bananas', '43. Lemons',\n",
       "       '44. Oranges', '45. Green Beans', '46. Cabbage', '47. Carrots',\n",
       "       '48. Celery', '49. Lettuce', '50. Onions ', '51. Potatoes',\n",
       "       '52. Spinach ', '53. Sweet Potatoes', '54. Peaches (can)',\n",
       "       '55. Pears (can)', '56. Pineapple (can)', '57. Asparagus (can)',\n",
       "       '58. Green Beans (can)', '59. Pork and Beans (can)',\n",
       "       '60. Corn (can)', '61. Peas (can)', '62. Tomatoes (can)',\n",
       "       '63. Tomato Soup (can)', '64. Peaches Dried', '65. Prunes Dried',\n",
       "       '66. Raisins Dried', '67. Peas Dried', '68. Lima Beans Dried',\n",
       "       '69. Navy Beans Dried', '70. Coffee', '71. Tea', '72. Cocoa',\n",
       "       '73. Chocolate', '74. Sugar', '75. Corn Sirup', '76. Molasses',\n",
       "       '77. Strawberry Preserves '], dtype=object)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/lp_stigler-diet/'\n",
    "filename = 'StiglerData1939.txt'\n",
    "thedata = pd.read_csv(thepath + filename, sep='\\t')\n",
    "thedata = thedata.dropna(how = 'all')\n",
    "commodities = (thedata['Commodity'].values)[:-1]\n",
    "allowance = thedata.iloc[-1, 4:].fillna(0).transpose()\n",
    "commodities\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_stigler_data(verbose=False):\n",
    "    import pandas as pd\n",
    "    thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/lp_stigler-diet/'\n",
    "    filename = 'StiglerData1939.txt'\n",
    "    thedata = pd.read_csv(thepath + filename, sep='\\t')\n",
    "    thedata = thedata.dropna(how = 'all')\n",
    "    commodities = (thedata['Commodity'].values)[:-1]\n",
    "    allowance = thedata.iloc[-1, 4:].fillna(0).transpose()\n",
    "    if verbose:\n",
    "        print('Daily nutrient content:')\n",
    "        print(tabulate(thedata.head()))\n",
    "        print('\\nDaily nutrient requirement:')\n",
    "        print(allowance)\n",
    "    return({'N_i_j':thedata.iloc[:-1, 4:].fillna(0).transpose().values,\n",
    "            'd_i':np.array(allowance),\n",
    "            'c_j':np.ones(len(commodities)),\n",
    "            'nbi':len(allowance),\n",
    "            'nbj':len(commodities),\n",
    "            'names_i': list(thedata.columns)[4:],\n",
    "            'names_j':commodities}) \n",
    "\n",
    "\n",
    "data = load_stigler_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MFgcIMCuOJ-A"
   },
   "source": [
    "Our dataset has the nutritional content of 77 commodities, and 9 nutrients."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IWdBlLrOOJ-C"
   },
   "source": [
    "### The Diet problem\n",
    "\n",
    "Problem setup:\n",
    "\n",
    "* Assume there are nutrients $i \\in \\{1,\\dots,I\\}$ (calories, protein, calcium, iron, etc.) that matter for a person's health, in such way that the minimum daily intake of nutrient $i$ should be $d_{i}$.\n",
    "\n",
    "* Nutrients do not come as standalone elements, but are combined into various foods. One dollar worth of food $j \\in \\{1,\\dots,J\\}$ yields a quantity $N_{ij}$ of nutrient $i \\in \\{1,\\dots,I\\}$.\n",
    "\n",
    "The problem is to find the diet that achieves the minimal intake of each nutrient for the cheapest price. If $q\\in\\mathbb{R}^{J}$ is a vector such that $q_{j}\\geq0$ is the quantity of food $j$ purchased, the quantity of nutrient $i$ ingested is $\\sum_{j=1}^{J}N_{ij}q_{j}$, and the cost of the diet is $\\sum_{j=1}^{J}q_{j}$. Letting $1_J$ be the column vector of ones of size $J$, the optimal diet is therefore given by\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{q\\geq 0} & ~ 1_J^{\\top}q\\\\\n",
    "\\text{s.t.} & ~ Nq \\geq d.\n",
    "\\end{align*}\n",
    "\n",
    "Before we tackle this problem, let's look at the more general instance of  a linear programming problem in standard form."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aVTufYycOJ-C"
   },
   "source": [
    "## A Crash Course on Linear Programming\n",
    "\n",
    "### Linear programming in canonical form\n",
    "\n",
    "Let $c\\in\\mathbb{R}^{n}$, $d\\in\\mathbb{R}^{m}$, $A$ be a $m\\times n$ matrix, and consider the following problem\n",
    "\n",
    "\\begin{align}\n",
    "V_{P} = \\max_{x\\in\\mathbb{R}_{+}^{n}}  & ~ c^{\\top} x \\\\\n",
    "\\text{s.t.} & ~ Ax \\leq d.\n",
    "\\end{align}\n",
    "\n",
    "This problem is a *linear programming problem*, as the objective function $x \\mapsto c^{\\top}x$ is linear, and as the constraints, namely $x \\in \\mathbb{R}_{+}^{n}$ and $Ax \\leq d$, are also linear (or more accurately, affine). This problem is called the *primal program*, for reasons to be explained soon. The set of $x$'s that satisfy the constraints are called *feasible solutions*; the set of solutions of the primal problem are called *optimal solutions*.\n",
    "\n",
    "#### Remarks\n",
    "\n",
    "* The previous diet problem can be reformulated into this problem — why?\n",
    "\n",
    "* A problem doesn't necessarily have a feasible solution (e.g. if $A=0$ and $d_j<0$ for some $j$), in which case (by convention) $V_{P}=-\\infty$.\n",
    "\n",
    "* If there are feasible solutions with $c^\\top x$ arbitrarily large, the primal problem is unbounded and $V_{P}=+\\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3HlLYyZYOJ-D"
   },
   "source": [
    "\n",
    "### Duality\n",
    "\n",
    "There is a powerful tool called *duality* which provides much insight into the analysis of the primal problem. The idea is to rewrite the problem as\n",
    "\n",
    "\\begin{align*}\n",
    "V_{P} = \\max_{x\\in\\mathbb{R}_{+}^{n}} \\left\\{ c^{\\top}x + L_{P}\\left(d-Ax\\right) \\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "where $L_{P}\\left(z\\right)$ is a penalty function whose value is zero if the constraint is met, that is if $z\\geq0$, and $-\\infty$ if it is not, namely if $z_j<0$ for at least one $j$. The simplest choice of such a penalty function is given by $L_{P}\\left(  z\\right)  =\\min_{y\\in\\mathbb{R}^{m}_{+}}\\left\\{  z^{\\top}y\\right\\}$. One has\n",
    "\n",
    "\\begin{align*}\n",
    "V_{P} = \\max_{x\\in\\mathbb{R}_{+}^{n}} \\min_{y\\in\\mathbb{R}^{m}_{+}} \\left\\{ c^{\\top}x + (d-Ax)^{\\top}y \\right\\}.\n",
    "\\end{align*}\n",
    "\n",
    "However, the minimax inequality $\\max_{x}\\min_{y}\\leq\\min_{y}\\max_{x}$ always holds, thus\n",
    "\n",
    "\\begin{align*}\n",
    "V_{P}\n",
    "&\\leq \\min_{y\\in\\mathbb{R}^{m}_{+}} \\max_{x\\in\\mathbb{R}_{+}^{n}} \\left\\{\n",
    "c^{\\top}x + (d-Ax)^{\\top}y \\right\\}\n",
    "= \\min_{y\\in\\mathbb{R}^{m}_{+}} \\max_{x\\in\\mathbb{R}_{+}^{n}} \\left\\{ x^{\\top}(c-A^{\\top}y) + d^{\\top}y \\right\\} \\\\\n",
    "&\\leq \\min_{y\\in\\mathbb{R}^{m}_{+}} \\left\\{ d^{\\top}y + L_{D} \\left(c-A^{\\top}y\\right) \\right\\} =:V_{D}\n",
    "\\end{align*}\n",
    "\n",
    "where $L_{D}\\left(z\\right) = \\max_{x\\in\\mathbb{R}_{+}^{n}}\\left\\{x^{\\top}z\\right\\}$ is equal to $0$ if $z\\in\\mathbb{R}_{-}^{n}$, and to $+\\infty$ if not. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qxdT3DNAOJ-D"
   },
   "source": [
    "Therefore, the value $V_{D}$ is expressed by the *dual program*\n",
    "\n",
    "\\begin{align}\n",
    "V_{D} = \\min_{y\\in\\mathbb{R}^{m}_{+}}  & ~ d^{\\top}y, \\\\\n",
    "\\text{s.t.} & ~ A^{\\top}y \\geq c\n",
    "\\end{align}\n",
    "\n",
    "and the weak duality inequality $V_{P}\\leq V_{D}$ holds. \n",
    "\n",
    "It turns out that as soon as either the primal or dual program has an optimal solution, then both\n",
    "programs have an optimal solution and the values of the two programs coincide,\n",
    "so the weak duality becomes an equality $V_{P}=V_{D}$ called *strong duality*.\n",
    "Further, if $x^{\\ast}\\in\\mathbb{R}_{+}^{n}$ is an optimal primal solution, and\n",
    "$y^{\\ast}\\in\\mathbb{R}^{m}_{+}$ is an optimal dual solution, then complementary\n",
    "slackness holds, that is $x_{i}^{\\ast}>0$ implies $\\left(  A^{\\top}y^{\\ast\n",
    "}\\right)  _{i}=c_{i}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gkHfQIegOJ-D"
   },
   "source": [
    "\n",
    "### Duality theorem\n",
    "We summarize these results into the following statement.\n",
    "\n",
    "---\n",
    "**Theorem.** In the setting described above:\n",
    "\n",
    "1. The weak duality inequality holds:\n",
    "\n",
    "\\begin{align}\n",
    "V_{P}\\leq V_{D}.\n",
    "\\end{align}\n",
    "\n",
    "2. As soon as the primal or the dual program has an optimal solution, then both programs have an optimal solution, and strong duality holds:\n",
    "\n",
    "\\begin{align}\n",
    "V_{P}=V_{D}.\n",
    "\\end{align}\n",
    "\n",
    "3. If $x^{\\ast}\\in\\mathbb{R}_{+}^{n}$ is an optimal primal solution, and $y^{\\ast}\\in\\mathbb{R}^{m}_{+}$ is an optimal dual solution, then complementary slackness holds:\n",
    "\n",
    "\\begin{align}\n",
    "x_{i}^{\\ast}>0 &\\implies (A^{\\top}y^{\\ast})_{i} = c_{i}, \\\\\n",
    "y_{j}^* > 0 &\\implies (Ax^*)_j = d_j.\n",
    "\\end{align}\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "46dTAOLDOJ-E"
   },
   "source": [
    "\n",
    "# Back to the diet problem\n",
    "\n",
    "Recall the optimal diet problem\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{q\\geq0} & ~ c^{\\top}q \\\\\n",
    "\\text{s.t.} & ~  Nq \\geq d\n",
    "\\end{align*}\n",
    "\n",
    "which has minimax formulation $\\min_{q\\geq0} \\max_{\\pi\\geq0} c^{\\top}q + d^{\\top}\\pi - q^{\\top}N^{\\top}\\pi$, so the dual is\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{\\pi\\geq0} & ~ d^{\\top}\\pi \\\\\n",
    "\\text{s.t.} ~ & N^{\\top}\\pi \\leq c.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XMd_IEfwOJ-E"
   },
   "source": [
    "## Interpretation of duality\n",
    "\n",
    "Imagine that there is a new firm called Nutrient Shoppe, who sells raw nutrients. Let $\\pi_{i}$ be the price of nutrient $i$. The cost of the diet is $d^{\\top}\\pi$. Consumers purchase raw nutrients and can generate\n",
    "\"synthetic foods\". The cost of the synthetic version of food $j$ is $\\sum_{i=1}^{m} N_{ij}\\pi_{i} = (N^{\\top}\\pi)_{j}$. The constraint thus means that each \"synthetic food\" is more affordable than its natural counterpart.\n",
    "\n",
    "The duality means that it is possible to price the nutrients so that the\n",
    "synthetic foods are cheaper than the natural ones, in such a way that the\n",
    "price of the synthetic diet equals the price of the natural diet.\n",
    "\n",
    "\n",
    "Complementary slackness yields:\n",
    "\n",
    "* $q_{j}>0$ implies $(N^{\\top}\\pi)_{j} = c_{j}$; that\n",
    "is, if natural food $j$ is actually purchased, then the prices of its\n",
    "synthetic and natural versions coincide.\n",
    "\n",
    "* $\\pi_{i}>0$ implies $\\left(  Nq\\right)  _{i}=d_{i}$; that is, if\n",
    "nutrient $i$ has a positive price, then the natural diet has the\n",
    "\"just right\" amount.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "S8X3uyGGOJ-E"
   },
   "source": [
    "## Numerically solving the diet problem\n",
    "\n",
    "To solve the primal problem we need to construct the objects $N$ and $d$. $c$ is simply a vector of ones, the size of the number of commodities. $N$ is a matrix of amounts of nutrients in each commodity. $d$ is the required daily allowance of each nutrient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RwM-wDLSOJ-F"
   },
   "source": [
    "### Finding a solver\n",
    "\n",
    "There are many solvers available. Some are open-souce, some are commercially available, even though the latter often come for free for an academic use. `scipy` has a linear programming solver called `scipy.optimize.linprog`, which we don't recommend except for small problems.\n",
    "\n",
    "Actually, with a problem of the (tiny) size that we are discussing in this introduction, the choice of the solver does not matter: all of them will return a solution whithin a fraction of seconds. But later on, we shall deal with problems with hundreds of thousands of variables or constraints, and the choice of the right linear programming solver will become crucial.\n",
    "We will need *large scale* solvers, which deal with sparse constraints (often the case as we shall see). \n",
    "\n",
    "See a benchmark of large scale solvers at:\n",
    "\n",
    "http://plato.asu.edu/ftp/lpsimp.html\n",
    "\n",
    "As we can see, there is a factor 100 between the average running time of the best performer (COPT, a commercial solver) and the worse performer of the list (GLPK, an open-source solver). [Gurobi](http://www.gurobi.com/), also a commercial solver with a free academic license, is among the best performers, and widely available in Python, R, C++ and many other languages. We shall use it as our linear programming solver of choice for the entirety of this course. \n",
    "\n",
    "To access Gurobi from a docker container, you need a token server license. These are available free of charge for academic users.\n",
    "\n",
    "If you do not have access to a full Gurobi license, or if you are running this notebook from Colab, you may still install Gurobi using `pip` with limited functionality -- namely, up to 2,000 variables and 2,000 constraints. Remove the commenting dash # in the cell below and install using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "aERxvLNAOJ-F",
    "outputId": "4209b31b-ba1d-438b-8cc6-bbd62e4c6f28"
   },
   "outputs": [],
   "source": [
    "# !python -m pip install -i https://pypi.gurobi.com gurobipy  ## only if Gurobi not here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uGkrZy_nOJ-F"
   },
   "source": [
    "Assuming we've installed the Gurobi solver and solved the license issue, let's first load up the Gurobi library and get some help. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "lnfYyvvVOJ-G"
   },
   "outputs": [],
   "source": [
    "import gurobipy as grb\n",
    "\n",
    "#grb.Model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eazveEm1OJ-G"
   },
   "source": [
    "Define the model and call the solver by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "DgHCRMp7OJ-G",
    "outputId": "afb1c3bb-42e6-4698-a394-e9d995902ca7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter Username\n",
      "Academic license - for non-commercial use only - expires 2023-12-23\n",
      "Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (win64)\n",
      "\n",
      "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, 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 9 rows, 77 columns and 570 nonzeros\n",
      "Model fingerprint: 0x7ae6b743\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e-01, 5e+03]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [8e-01, 8e+01]\n",
      "Presolve removed 0 rows and 47 columns\n",
      "Presolve time: 0.02s\n",
      "Presolved: 9 rows, 30 columns, 240 nonzeros\n",
      "\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0    0.0000000e+00   1.384688e+01   0.000000e+00      0s\n",
      "       5    1.0866228e-01   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 5 iterations and 0.04 seconds (0.00 work units)\n",
      "Optimal objective  1.086622782e-01\n"
     ]
    }
   ],
   "source": [
    "m = grb.Model()\n",
    "varq_j = m.addMVar(data['nbj'] )\n",
    "m.setObjective(varq_j.sum(), grb.GRB.MINIMIZE)\n",
    "c_j = m.addConstr(data['N_i_j'] @ varq_j >= data['d_i'] )\n",
    "m.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H1P3XFHfOJ-G"
   },
   "source": [
    "We are after the optimal solutions `X`, the dual solution `pi` and the value function `objval`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "Y0NhNJDlOJ-G"
   },
   "outputs": [],
   "source": [
    "q_j,π_i = varq_j.X,c_j.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iY61K1sCOJ-H"
   },
   "source": [
    "We can print the optimal diet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "xxdxyN62OJ-H",
    "outputId": "3e2e3c8d-1477-4bfb-823b-a473b3ed33b7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "***Optimal solution***\n",
      "-------------------------  -----------\n",
      "1. Wheat Flour (Enriched)   0.0295191\n",
      "30. Liver (Beef)            0.00189256\n",
      "46. Cabbage                 0.0112144\n",
      "52. Spinach                 0.00500766\n",
      "69. Navy Beans Dried        0.0610286\n",
      "Total cost (optimal):      39.6617\n",
      "-------------------------  -----------\n"
     ]
    }
   ],
   "source": [
    "def print_optimal_diet(q_j):\n",
    "    print('***Optimal solution***')\n",
    "    total,thelist = 0.0, []\n",
    "    for j,commodity in enumerate(data['names_j']):\n",
    "        if q_j[j] > 0:\n",
    "            total += q_j[j] * 365\n",
    "            thelist.append([data['names_j'][j],q_j[j]])\n",
    "    thelist.append(['Total cost (optimal):', total])\n",
    "    print(tabulate(thelist))\n",
    "\n",
    "print_optimal_diet(q_j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MHZK_OUEOJ-H"
   },
   "source": [
    "As promised, we achieve the minimum cost bundle at \\$39.67 per year in 1939 dollars. We can compare this to Stigler's solution which was:\n",
    "\n",
    "|Food| Annual Quantities| Annual Cost|\n",
    "| ---------- | ------------------ | ------------ |\n",
    "| Wheat Flour | \t370 lb.|   \\$13.33 |\n",
    "| Evaporated Milk | \t57 cans |\t  \\$3.84 |\n",
    "|Cabbage| \t111 lb. \t  |\\$4.11|\n",
    "|Spinach| \t23 lb. \t  |\\$1.85|\n",
    "|Dried Navy Beans| \t285 lb. |\t\\$16.80|\n",
    "|Total Annual Cost| \t&nbsp; \t| \\$39.93 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kwdLkCwyOJ-I"
   },
   "source": [
    "**Exercise**. Recover the solution with `scipy`'s linear programming solver, namely `scipy.optimize.linprog`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "D642CY2KOJ-I",
    "outputId": "14fe0a54-3660-45c0-f4ea-82fe400e2aae"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "***Optimal solution***\n",
      "-------------------------  -----------\n",
      "1. Wheat Flour (Enriched)   0.0295191\n",
      "30. Liver (Beef)            0.00189256\n",
      "46. Cabbage                 0.0112144\n",
      "52. Spinach                 0.00500766\n",
      "69. Navy Beans Dried        0.0610286\n",
      "Total cost (optimal):      39.6617\n",
      "-------------------------  -----------\n"
     ]
    }
   ],
   "source": [
    "# answer:\n",
    "from scipy.optimize import linprog as lp\n",
    "\n",
    "res = lp(data['c_j'], A_ub=-data['N_i_j'], b_ub=-data['d_i'])\n",
    "qbis_j = res.x\n",
    "\n",
    "print_optimal_diet(qbis_j)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}