{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#
Block 1: Linear programming: the diet problem
\n", "###
Alfred Galichon (NYU & Sciences Po)
\n", "##
'math+econ+code' masterclass on optimal transport and economic applications
\n", "
© 2018-2021 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274, and contributions by Jules Baudet, Pauline Corblet, Gregory Dannay, and James Nesbit are acknowledged.
\n", "\n", "####
With python code examples
" ] }, { "cell_type": "markdown", "metadata": {}, "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": {}, "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 = combination of foods that meet 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": {}, "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": 5, "metadata": {}, "outputs": [], "source": [ "import scipy.sparse as sp\n", "import pandas as pd\n", "import numpy as np\n", "\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')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our dataset has the nutritional content of $77$ commodities, and in the final row, the daily minimum requirement of each of these nutrients." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CommodityUnitPrice Aug.15 1939(cents)Edible Weight per $1.00 (grams)Calories (1000)Protein(grams)Calcium(grams)Iron(mg.)Vitamin A(1000 I.U)Thiamine(mg.)Riboflavin(mg.)Niacin(mg.)Asorbic Acid (mg.)
01. Wheat Flour (Enriched)10 lb.36.012600.044.71411.02.0365.0NaN55.433.3441.0NaN
12. Macaroni1 lb.14.13217.011.6418.00.754.0NaN3.21.968.0NaN
23. Wheat Cereal (Enriched)28 oz.24.23280.011.8377.014.4175.0NaN14.48.8114.0NaN
34. Corn Flakes8 oz.7.13194.011.4252.00.156.0NaN13.52.368.0NaN
45. Corn Meal1 lb.4.69861.036.0897.01.799.030.917.47.9106.0NaN
\n", "
" ], "text/plain": [ " Commodity Unit Price Aug.15 1939(cents) \\\n", "0 1. Wheat Flour (Enriched) 10 lb. 36.0 \n", "1 2. Macaroni 1 lb. 14.1 \n", "2 3. Wheat Cereal (Enriched) 28 oz. 24.2 \n", "3 4. Corn Flakes 8 oz. 7.1 \n", "4 5. Corn Meal 1 lb. 4.6 \n", "\n", " Edible Weight per $1.00 (grams) Calories (1000) Protein(grams) \\\n", "0 12600.0 44.7 1411.0 \n", "1 3217.0 11.6 418.0 \n", "2 3280.0 11.8 377.0 \n", "3 3194.0 11.4 252.0 \n", "4 9861.0 36.0 897.0 \n", "\n", " Calcium(grams) Iron(mg.) Vitamin A(1000 I.U) Thiamine(mg.) \\\n", "0 2.0 365.0 NaN 55.4 \n", "1 0.7 54.0 NaN 3.2 \n", "2 14.4 175.0 NaN 14.4 \n", "3 0.1 56.0 NaN 13.5 \n", "4 1.7 99.0 30.9 17.4 \n", "\n", " Riboflavin(mg.) Niacin(mg.) Asorbic Acid (mg.) \n", "0 33.3 441.0 NaN \n", "1 1.9 68.0 NaN \n", "2 8.8 114.0 NaN \n", "3 2.3 68.0 NaN \n", "4 7.9 106.0 NaN " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "thedata.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Diet problem\n", "\n", "Problem setup:\n", "\n", "* Assume there are nutrients $i\\in\\left\\{ 1,...,I\\right\\} $ (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\\left\\{ 1,...,J\\right\\}$ yields a quantity $N_{ij}$ of nutrient $i\\in\\left\\{1,...,I\\right\\}$.\n", "\n", "The problem is to find the diet that achieves the minimal intake of each nutrient at a cheapest price. If $q\\in\\mathbb{R}^{J}$ is a vector such that $q_{j}\\geq0$ is the quantity of food $i$ 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", "s.t.~ & Nq\\geq d.\\nonumber\n", "\\end{align*}\n", "\n", "Before we tackle this problem, let's look at the more general instance of linear programming problem in standard form." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A Crash Course on Linear Programming\n", "\n", "### Linear programming in standard 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", "s.t.~Ax & \\leq d\n", "\\end{align}\n", "\n", "This problem is a *linear programming problem*, as the objective function, namely $x\\rightarrow c^{\\top}x$ is linear, and as the constraint, namely $x\\in\\mathbb{R}_{+}^{n}$ and $Ax=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 constraint 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 reformulate into this problem - why?\n", "\n", "* A problem does not necessarly have a feasible solution (e.g. if $A=0$ and $d\\neq0$), in which case (by convention) $V_{P}=-\\infty$.\n", "\n", "* The whole space may be solution (e.g. if $A=0$ and $d=0$), in which case $V_{P}=+\\infty$." ] }, { "cell_type": "markdown", "metadata": {}, "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(\n", "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=0$, and $-\\infty$ if it is not, namely if $z\\neq0$. The simplest choice of such 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+\\left( d-Ax\\right) ^{\\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} & \\leq\\min_{y\\in\\mathbb{R}^{m}}\\max_{x\\in\\mathbb{R}_{+}^{n}}\\left\\{\n", "c^{\\top}x+\\left( d-Ax\\right) ^{\\top}y\\right\\} =\\min_{y\\in\\mathbb{R}^{m}\n", "}\\max_{x\\in\\mathbb{R}_{+}^{n}}\\left\\{ x^{\\top}\\left( c-A^{\\top}y\\right)\n", "+d^{\\top}y\\right\\} \\\\\n", "& \\leq\\min_{y\\in\\mathbb{R}^{m}}\\left\\{ d^{\\top}y+L_{D}\\left( c-A^{\\top\n", "}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": {}, "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", "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": {}, "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 have 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\\text{ implies }\\left( A^{\\top}y^{\\ast}\\right) _{i}=c_{i}.\n", "\\end{align}\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "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", "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", "s.t.~ & N^{\\top}\\pi\\leq c\n", "\\end{align*}\n", "\n", "Complementary slackness yields:\n", "\n", "* $q_{j}>0$ implies $\\left( N^{\\intercal}\\pi\\right) _{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." ] }, { "cell_type": "markdown", "metadata": {}, "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$. Consumer 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}=\\left(N^{\\intercal}\\pi\\right)_{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" ] }, { "cell_type": "markdown", "metadata": {}, "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 amount of nutrients in each commodity. $d$ is the required daily allowance of each nutrient." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Calories (1000) 3.0\n", "Protein(grams) 70.0\n", "Calcium(grams) 0.8\n", "Iron(mg.) 12.0\n", "Vitamin A(1000 I.U) 5.0\n", "Thiamine(mg.) 1.8\n", "Riboflavin(mg.) 2.7\n", "Niacin(mg.) 18.0\n", "Asorbic Acid (mg.) 75.0\n", "Name: 77, dtype: float64" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "commodity = (thedata['Commodity'].to_numpy())[:-1]\n", "intake = thedata.iloc[:-1, 4:].fillna(0).transpose().to_numpy()\n", "allowance = thedata.iloc[-1, 4:].fillna(0).transpose()\n", "allowance" ] }, { "cell_type": "markdown", "metadata": {}, "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 problem 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. The system administrator of your university can help you setup one. For the purpose of this course, the Gurobi team has graciously provided us with educational licenses so that you can start working immediately. Please keep in mind that these licenses can only be used for the purpose of this course. They should not be communicated to anyone else. They expire January 31, 2021.\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": "markdown", "metadata": {}, "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": 22, "metadata": {}, "outputs": [], "source": [ "import gurobipy as grb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the model and call the solver by:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Academic license - for non-commercial use only - expires 2022-11-14\n", "Using license file c:\\gurobi911\\gurobi.lic\n", "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)\n", "Thread count: 6 physical cores, 12 logical processors, using up to 12 threads\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.01s\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.03 seconds\n", "Optimal objective 1.086622782e-01\n" ] } ], "source": [ "m = grb.Model('optimalDietMatrix')\n", "x = m.addMVar(shape=commodity.shape)\n", "m.setObjective(x.sum(), grb.GRB.MINIMIZE)\n", "c = m.addConstr(sp.csr_matrix(intake) @ x >= np.array(allowance), name=\"c\")\n", "m.optimize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are after the optimal solutions `x`, the dual solution `pi` and the value function `objval`" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "x = m.getAttr('X')\n", "pi = m.getAttr('pi')" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "***Optimal solution***\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01
01. Wheat Flour (Enriched)10.774458
130. Liver (Beef)0.690783
246. Cabbage4.093269
352. Spinach1.827796
469. Navy Beans Dried22.275426
5Total cost (optimal):39.661732
\n", "
" ], "text/plain": [ " 0 1\n", "0 1. Wheat Flour (Enriched) 10.774458\n", "1 30. Liver (Beef) 0.690783\n", "2 46. Cabbage 4.093269\n", "3 52. Spinach 1.827796\n", "4 69. Navy Beans Dried 22.275426\n", "5 Total cost (optimal): 39.661732" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print('***Optimal solution***')\n", "total = 0\n", "thelist = []\n", "for i,thecom in enumerate(commodity):\n", " if x[i] > 1e-8:\n", " total += x[i] * 365\n", " thelist.append([thecom,x[i]*365])\n", "\n", "thelist.append(['Total cost (optimal):', total])\n", "pd.DataFrame(thelist)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As promised, we achieve the minimum cost bundle at $\\$39.67$ per year in $1939$ dollars. We can compare this to Stigler's solutions 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  \t| \\$39.93 |" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**. Recover the solution with `scipy`'s linear programming solver, namely `scipy.optimize.linprog`." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0 1\n", "0 1. Wheat Flour (Enriched) 0.029519\n", "1 30. Liver (Beef) 0.001893\n", "2 46. Cabbage 0.011214\n", "3 52. Spinach 0.005008\n", "4 69. Navy Beans Dried 0.061029\n", "Total cost:39.661731030395636\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\antoi\\Anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:87: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray\n", " return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n" ] } ], "source": [ "# answer:\n", "\n", "from scipy.optimize import linprog as lp\n", "\n", "res = lp(np.ones(len(intake[0])), A_ub=-intake, b_ub=-allowance)\n", "x = res.get('x')\n", "result = res.get('fun')\n", "total = 0\n", "diet = []\n", "for i, com in enumerate(commodity):\n", " if x[i] >1e-8:\n", " total = total + x[i]*365\n", " diet.append([com, x[i]])\n", "print(pd.DataFrame(diet))\n", "print(\"Total cost:\" + str(total))" ] } ], "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.7.4" } }, "nbformat": 4, "nbformat_minor": 2 }