ampl-lecture/diet_case_study.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vrsB0mVy8Iyu"
      },
      "source": [
        "# Diet and Other Input Models: Minimizing Costs\n",
        "[![diet_case_study.ipynb](https://img.shields.io/badge/github-%23121011.svg?logo=github)](https://github.com/ampl/colab.ampl.com/blob/master/ampl-lecture/diet_case_study.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/ampl/colab.ampl.com/blob/master/ampl-lecture/diet_case_study.ipynb) [![Kaggle](https://kaggle.com/static/images/open-in-kaggle.svg)](https://kaggle.com/kernels/welcome?src=https://github.com/ampl/colab.ampl.com/blob/master/ampl-lecture/diet_case_study.ipynb) [![Gradient](https://assets.paperspace.io/img/gradient-badge.svg)](https://console.paperspace.com/github/ampl/colab.ampl.com/blob/master/ampl-lecture/diet_case_study.ipynb) [![Open In SageMaker Studio Lab](https://studiolab.sagemaker.aws/studiolab.svg)](https://studiolab.sagemaker.aws/import/github/ampl/colab.ampl.com/blob/master/ampl-lecture/diet_case_study.ipynb) [![Hits](https://h.ampl.com/https://github.com/ampl/colab.ampl.com/blob/master/ampl-lecture/diet_case_study.ipynb)](https://colab.ampl.com)\n",
        "\n",
        "Description: Diet case study, Chapter 2 from the AMPL book adapted to Python\n",
        "\n",
        "Tags: amplpy, ampl-lecture\n",
        "\n",
        "Notebook author: Marcos Dominguez Velad <<marcos@ampl.com>>\n",
        "\n",
        "Model author: N/A\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "wJe2Z4eu8Iyw",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "b4f0970f-1e14-45cf-d779-3f5503a21895"
      },
      "outputs": [],
      "source": [
        "# Install dependencies\n",
        "%pip install -q amplpy numpy pandas"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "XLCYKvPM8Iyx",
        "outputId": "1c01ea21-92b3-46fb-fe88-fb4b82392257"
      },
      "outputs": [],
      "source": [
        "# Google Colab & Kaggle integration\n",
        "from amplpy import AMPL, ampl_notebook\n",
        "\n",
        "ampl = ampl_notebook(\n",
        "    modules=[\"highs\"],  # modules to install\n",
        "    license_uuid=\"default\",  # license to use\n",
        ")  # instantiate AMPL object and register magics"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Nv9DV-Cj8Iyx"
      },
      "source": [
        "This notebook provides the implementation of the production problem described in the book *AMPL: A Modeling Language for Mathematical Programming* by Robert Fourer, David M. Gay, and Brian W. Kernighan. This implementation is Python-based using amplpy rather than pure Ampl style.\n",
        "\n",
        "## Diet problem\n",
        "As an intuitive example of a cost-minimizing model, we use the well-known \"diet problem\", which finds a mix of foods that satisfies requirements on the amounts of various vitamins. We will construct a small, explicit linear program, and then show how a general model can be formulated for all linear programs of that kind.\n",
        "\n",
        "After formulating the diet model, we will discuss a few changes that might make it more realistic. The full power of this model, however, derives from its applicability to many situations that have nothing to do with diets. A general model derived from the diet formulation can be also applied to blending, economics, and scheduling.\n",
        "\n",
        "### Solving a diet problem instance\n",
        "\n",
        "Consider the problem of choosing prepared foods to meet certain nutritional requirements. Suppose that precooked dinners of the following kinds are available for the following prices per package:\n",
        "\n",
        "|      | Name              | Price $ |\n",
        "| ---- | ----------------- | ------- |\n",
        "| BEEF | beef              | 3.19    |\n",
        "| CHK  | chicken           | 2.59    |\n",
        "| FISH | fish              | 2.29    |\n",
        "| HAM  | ham               | 2.89    |\n",
        "| MCH  | macaroni & cheese | 1.89    |\n",
        "| MTL  | meat loaf         | 1.99    |\n",
        "| SPG  | spaghetti         | 1.99    |\n",
        "| TUR  | turkey            | 2.49    |\n",
        "\n",
        "These dinners provide the following percentages, per package, of the minimum daily requirements for vitamins A, C, B1 and B2:\n",
        "\n",
        "\n",
        "|      | A   | C   | B1  | B2  |\n",
        "| ---- | --- | --- | --- | --- |\n",
        "| BEEF | 60% | 20% | 10% | 15% |\n",
        "| CHK  | 8   | 0   | 20  | 20  |\n",
        "| FISH | 8   | 10  | 15  | 10  |\n",
        "| HAM  | 40  | 40  | 35  | 10  |\n",
        "| MCH  | 15  | 35  | 15  | 15  |\n",
        "| MTL  | 70  | 30  | 15  | 15  |\n",
        "| SPG  | 25  | 50  | 25  | 15  |\n",
        "| TUR  | 60  | 20  | 15  | 10  |\n",
        "\n",
        "The problem is to find the cheapest combination of packages that will meet a week's requirements - that is, at least 700% of the daily requirement for each nutrient.\n",
        "\n",
        "Let us write $X_{beef}$ for the number of packages of beef dinner to be purchased, $X_{chk}$ for the number of packages of chicken dinner, and so forth. Then the total cost of the diet will be:\n",
        "```ampl\n",
        "total cost =\n",
        "      3.19 Xbeef + 2.59 Xchk + 2.29 Xfish + 2.89 Xham +\n",
        "      1.89 Xmch + 1.99 Xmtl + 1.99 Xspg + 2.49 Xtur\n",
        "```\n",
        "\n",
        "The total percentage of the vitamin A requirement is given by a similar formula, except that X BEEF , X CHK , and so forth are multiplied by the percentage per package instead of the cost per package:\n",
        "```ampl\n",
        "total percentage of vitamin A daily requirement met =\n",
        "      60 Xbeef + 8 Xchk + 8 Xfish + 40 Xham +\n",
        "      15 Xmch + 70 Xmtl + 25 Xspg + 60 Xtur\n",
        "```\n",
        "\n",
        "This amount needs to be greater than or equal to 700 percent. There is a similar formula for each of the other vitamins, and each of these also needs to be $\\geq$ 700.\n",
        "\n",
        "Putting these all together, we have the following linear program:\n",
        "\n",
        "```ampl\n",
        "Minimize\n",
        "    3.19 Xbeef + 2.59 Xchk + 2.29 Xfish + 2.89 Xham +\n",
        "    1.89 Xmch + 1.99 Xmtl + 1.99 Xspg + 2.49 Xtur\n",
        "\n",
        "Subject to\n",
        "     60 Xbeef + 8 Xchk + 8 Xfish + 40 Xham +\n",
        "    15 Xmch + 70 Xmtl + 25 Xspg + 60 Xtur >= 700\n",
        "\n",
        "      20 Xbeef + 0 Xchk + 10 Xfish + 40 Xham +\n",
        "      35 Xmch + 30 Xmtl + 50 Xspg + 20 Xtur >= 700\n",
        "\n",
        "      10 Xbeef + 20 Xchk + 15 Xfish + 35 Xham +\n",
        "      15 Xmch + 15 Xmtl + 25 Xspg + 15 Xtur >= 700\n",
        "\n",
        "      15 Xbeef + 20 Xchk + 10 Xfish + 10 Xham +\n",
        "      15 Xmch + 15 Xmtl + 15 Xspg + 10 Xtur >= 700\n",
        "\n",
        "      Xbeef >= 0, Xchk >= 0, Xfish >= 0, Xham >= 0,\n",
        "      Xmch >= 0, Xmtl >= 0, Xspg >= 0, Xtur >= 0\n",
        "```\n",
        "\n",
        "At the end we have added the common-sense requirement that no fewer than zero packages of a food can be purchased. And now, we can transcribe the model to an AMPL statement of the explicit diet LP:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "NghxO-oM8Iyy",
        "outputId": "9aaae3cc-bb62-4851-bb89-1a4fda43bec5"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Writing diet0.mod\n"
          ]
        }
      ],
      "source": [
        "%%writefile diet0.mod\n",
        "# Variables\n",
        "var Xbeef >= 0; var Xchk >= 0; var Xfish >= 0;\n",
        "var Xham >= 0; var Xmch >= 0; var Xmtl >= 0;\n",
        "var Xspg >= 0; var Xtur >= 0;\n",
        "\n",
        "# Objective function (minimizing)\n",
        "minimize cost:\n",
        "\t3.19*Xbeef + 2.59*Xchk + 2.29*Xfish + 2.89*Xham +\n",
        "\t1.89*Xmch + 1.99*Xmtl + 1.99*Xspg + 2.49*Xtur;\n",
        "\n",
        "# Constraints\n",
        "subject to A:\n",
        "\t60*Xbeef + 8*Xchk + 8*Xfish + 40*Xham +\n",
        "\t15*Xmch + 70*Xmtl + 25*Xspg + 60*Xtur >= 700;\n",
        "subject to C:\n",
        "\t20*Xbeef + 0*Xchk + 10*Xfish + 40*Xham +\n",
        "\t35*Xmch + 30*Xmtl + 50*Xspg + 20*Xtur >= 700;\n",
        "subject to B1:\n",
        "\t10*Xbeef + 20*Xchk + 15*Xfish + 35*Xham +\n",
        "\t15*Xmch + 15*Xmtl + 25*Xspg + 15*Xtur >= 700;\n",
        "subject to B2:\n",
        "\t15*Xbeef + 20*Xchk + 10*Xfish + 10*Xham +\n",
        "\t15*Xmch + 15*Xmtl + 15*Xspg + 10*Xtur >= 700;"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "laQz-8lv8Iyz"
      },
      "source": [
        "A few AMPL commands then suffice to read the file, send the LP to a solver, and retrieve the results (using HiGHS solver)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "rGWGL6Oa8Iyz",
        "outputId": "8fe71fbb-f167-413d-d432-13efb94bd838"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "HiGHS 1.10.0: \b\b\b\b\b\b\b\b\b\b\b\b\b\bHiGHS 1.10.0: optimal solution; objective 88.2\n",
            "2 simplex iterations\n",
            "0 barrier iterations\n",
            "Xbeef = 0\n",
            "Xchk = 0\n",
            "Xfish = 0\n",
            "Xham = 0\n",
            "Xmch = 46.6667\n",
            "Xmtl = 0\n",
            "Xspg = 0\n",
            "Xtur = 0\n",
            "\n"
          ]
        }
      ],
      "source": [
        "ampl.read('diet0.mod')\n",
        "ampl.solve(solver='highs')\n",
        "ampl.display('Xbeef,Xchk,Xfish,Xham,Xmch,Xmtl,Xspg,Xtur')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "BVOAhrwRtG76"
      },
      "outputs": [],
      "source": [
        "assert ampl.solve_result == \"solved\", ampl.solve_result"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BWWFP-eN8Iy0"
      },
      "source": [
        "The optimal solution is found quickly, but it is hardly what we might have hoped for. The cost is minimized by a monotonous diet of 46 and 2/3 packages of macaroni and cheese! You can check that this neatly provides $15\\% \\times 46\\frac{2}{3} = 700\\%$ of the requirement for vitamins A, B1 and B2, and a lot more vitamin C than necessary; the cost is only $\\$1.89 × 46\\frac{2}{3} = \\$88.20.$\n",
        "\n",
        "You might guess that a better solution would be generated by requiring the amount of each vitamin to equal 700% exactly. Such a requirement can easily be imposed by changing each >= to = in the AMPL constraints. If you go ahead and solve the changed LP, you will find that the diet does indeed become more varied: approximately 19.5 packages of chicken, 16.3 of macaroni and cheese, and 4.3 of meat loaf. But since equalities are more restrictive than inequalities, the cost goes up to $89.99.\n",
        "\n",
        "## An AMPL model for the diet problem\n",
        "\n",
        "Clearly we will have to consider more extensive modifications to our linear program in order to produce a diet that is even remotely acceptable. We will probably want to change the sets of food and nutrients, as well as the nature of the constraints and bounds. As in the production example of the previous chapter, this will be much easier to do if we rely on a general model that can be coupled with a variety of specific data files.\n",
        "\n",
        "This model deals with two things: nutrients and foods. Thus we begin an AMPL model by declaring sets of each:\n",
        "```ampl\n",
        "set NUTR;\n",
        "set FOOD;\n",
        "```\n",
        "Next we need to specify the numbers required by the model. Certainly a positive cost should be given for each food:\n",
        "```ampl\n",
        "param cost {FOOD} > 0;\n",
        "```\n",
        "We also specify that for each food there are lower and upper limits on the number of packages in the diet:\n",
        "```ampl\n",
        "param f_min {FOOD} >= 0;\n",
        "param f_max {j in FOOD} >= f_min[j];\n",
        "```\n",
        "Notice that we need a dummy index `j` to run over `FOOD` in the declaration of `f_max`, in order to say that the maximum for each food must be greater than or equal to the corresponding minimum.\n",
        "\n",
        "To make this model somewhat more general than our examples so far, we also specify similar lower and upper limits on the amount of each nutrient in the diet:\n",
        "\n",
        "```ampl\n",
        "param n_min {NUTR} >= 0;\n",
        "param n_max {i in NUTR} >= n_min[i];\n",
        "```\n",
        "Finally, for each combination of a nutrient and a food, we need a number that represents the amount of the nutrient in one package of the food. Such a \"product\" of two sets is written by listing them both:\n",
        "```ampl\n",
        "param amt {NUTR,FOOD} >= 0;\n",
        "```\n",
        "References to this parameter require two indices. For example, `amt[i,j]` is the amount of nutrient `i` in a package of food `j`.\n",
        "\n",
        "The decision variables for this model are the numbers of packages to buy of the different foods:\n",
        "```ampl\n",
        "var Buy {j in FOOD} >= f_min[j], <= f_max[j];\n",
        "```\n",
        "The number of packages of some food `j` to be bought will be called `Buy[j]`; in any acceptable solution it will have to lie between `f_min[j]` and `f_max[j]`.\n",
        "\n",
        "The total cost of buying a food `j` is the cost per package, `cost[j]`, times the number of packages, `Buy[j]`. The objective to be minimized is the sum of this product over all foods `j`:\n",
        "```ampl\n",
        "minimize Total_Cost: sum {j in FOOD} cost[j] * Buy[j];\n",
        "```\n",
        "Similarly, the amount of a nutrient `i` supplied by a food `j` is the nutrient per package, `amt[i,j]`, times the number of packages `Buy[j]`. The total amount of nutrient `i` supplied is the sum of this product over all foods `j`:\n",
        "```ampl\n",
        "sum {j in FOOD} amt[i,j] * Buy[j];\n",
        "```\n",
        "To complete the model, we need only specify that each such sum must lie between the appropriate bounds. Our constraint declaration begins\n",
        "```ampl\n",
        "subject to Diet {i in NUTR}:\n",
        "```\n",
        "to say that a constraint named `Diet[i]` must be imposed for each member `i` of `NUTR`. The rest of the declaration gives the algebraic statement of the constraint for nutrient `i`: the variables must satisfy\n",
        "```ampl\n",
        "n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i]\n",
        "```\n",
        "A \"double inequality\" like this is interpreted in the obvious way: the value of the sum in the middle must lie between `n_min[i]` and `n_max[i]`. We can write all together into a model file `diet.mod`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "61jys-Xc8Iy0",
        "outputId": "d1c5c362-26c0-4ac5-c0e4-803ad820e105"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Writing diet.mod\n"
          ]
        }
      ],
      "source": [
        "%%writefile diet.mod\n",
        "set NUTR;\n",
        "set FOOD;\n",
        "\n",
        "param cost {FOOD} > 0;\n",
        "param f_min {FOOD} >= 0;\n",
        "param f_max {j in FOOD} >= f_min[j];\n",
        "\n",
        "param n_min {NUTR} >= 0;\n",
        "param n_max {i in NUTR} >= n_min[i];\n",
        "\n",
        "param amt {NUTR,FOOD} >= 0;\n",
        "\n",
        "var Buy {j in FOOD} >= f_min[j], <= f_max[j];\n",
        "\n",
        "minimize Total_Cost:  sum {j in FOOD} cost[j] * Buy[j];\n",
        "\n",
        "subject to Diet {i in NUTR}:\n",
        "   n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i];"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qx0zHP1qtG78"
      },
      "source": [
        "We can use the `read()` function from Amplpy to read a model file by specifying the name (or the path) of the model. This must happen before reading the concrete data for the problem."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "id": "0tK9cvPMtG79"
      },
      "outputs": [],
      "source": [
        "# Reset Ampl to delete the previous model. Another Ampl object could be created as well\n",
        "ampl.reset()\n",
        "ampl.read('diet.mod')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CT-nHc-k8Iy1"
      },
      "source": [
        "## Data\n",
        "\n",
        "By specifying appropriate data, we can solve any of the linear programs that correspond to the above model. In this example we are assuming you are familiar to basic data structures in Python like lists, dictionaries, Pandas DataFrames, and Numpy Matrices.\n",
        "\n",
        "Let's begin by using the data from the previous example.\n",
        "\n",
        "The sets of `FOOD` and `NUTR` are defined by Python lists (`foods` and `nutrs`). The parameters related to foods are `cost`, `f_min` and `f_max`, which are in fact the columns of the dataframe `food_df`. There is another column in this dataframe to specify the food, which be the index for each row in the dataframe. Something equivalent happens with the parameters related to the nutrients, all the data is gathered in `nutr_df`. Notice that the values of `f_min` and `n_min` are as given originally, while `f_max` and `n_max` are set, for the time being, to large values that won't affect the optimal solution. The data for `amt` is given by another dataframe, generated through a numpy matrix to show a more compact representation.\n",
        "\n",
        "This is not the only way to load data into the Ampl model, and these will be covered in other resources."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "id": "vBx37ZFc8Iy1"
      },
      "outputs": [],
      "source": [
        "# Defining the data in Python data structures\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "\n",
        "# Indexing sets (as Python lists)\n",
        "foods = ['BEEF', 'CHK', 'FISH', 'HAM', 'MCH', 'MTL', 'SPG', 'TUR']\n",
        "nutrs = ['A', 'B1', 'B2', 'C']\n",
        "\n",
        "# Pandas Dataframe with 4 columns: FOOD (indexing set in the model), cost, f_min, f_max\n",
        "# Notice that since FOOD is the indexing set in the model, we set it as the index in the dataframe\n",
        "food_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"BEEF\", 3.19, 0, 100),\n",
        "        (\"CHK\", 2.59, 0, 100),\n",
        "        (\"FISH\", 2.29, 0, 100),\n",
        "        (\"HAM\", 2.89, 0, 100),\n",
        "        (\"MCH\", 1.89, 0, 100),\n",
        "        (\"MTL\", 1.99, 0, 100),\n",
        "        (\"SPG\", 1.99, 0, 100),\n",
        "        (\"TUR\", 2.49, 0, 100),\n",
        "    ],\n",
        "    columns=[\"FOOD\", \"cost\", \"f_min\", \"f_max\"],\n",
        ").set_index(\"FOOD\")\n",
        "\n",
        "# Pandas Dataframe with 3 columns: NUTR (indexing set in the model), n_min, n_max\n",
        "nutr_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"A\", 700, 10000),\n",
        "        (\"C\", 700, 10000),\n",
        "        (\"B1\", 700, 10000),\n",
        "        (\"B2\", 700, 10000),\n",
        "    ],\n",
        "    columns=[\"NUTR\", \"n_min\", \"n_max\"],\n",
        ").set_index(\"NUTR\")\n",
        "\n",
        "# amt_df is a Pandas Dataframe generated from a numpy matrix, each row corresponds\n",
        "# to a food, and each column to a nutrient.\n",
        "amt_df = pd.DataFrame(\n",
        "    np.array(\n",
        "        [\n",
        "            [60, 8, 8, 40, 15, 70, 25, 60],\n",
        "            [20, 0, 10, 40, 35, 30, 50, 20],\n",
        "            [10, 20, 15, 35, 15, 15, 25, 15],\n",
        "            [15, 20, 10, 10, 15, 15, 15, 10],\n",
        "        ]\n",
        "    ),\n",
        "    columns=food_df.index.to_list(),\n",
        "    index=nutr_df.index.to_list(),\n",
        ")\n",
        "\n",
        "# Load the data into Ampl:\n",
        "\n",
        "ampl.set[\"FOOD\"] = foods\n",
        "ampl.set[\"NUTR\"] = nutrs\n",
        "\n",
        "# Load data related to foods:\n",
        "ampl.set_data(food_df, \"FOOD\")\n",
        "# Or alternatively:\n",
        "#ampl.param['cost'] = df_food['cost']\n",
        "#ampl.param['f_min'] = df_food['f_min']\n",
        "#ampl.param['f_max'] = df_food['f_max']\n",
        "\n",
        "# Load data related to nutrients:\n",
        "ampl.set_data(nutr_df, \"NUTR\")\n",
        "# Or alternatively:\n",
        "#ampl.param['n_min'] = df_nutr['n_min']\n",
        "#ampl.param['n_max'] = df_nutr['n_max']\n",
        "\n",
        "ampl.param[\"amt\"] = amt_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YYUt_whZ8Iy1"
      },
      "source": [
        "Once we have read the model and loaded the data in AMPL, it is used to solve the resulting linear program and show the result for the `Buy` variables."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0EuIFJQy8Iy1",
        "outputId": "04dc944d-16fc-4409-bd78-ab78ab5c8d68"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "HiGHS 1.10.0: \b\b\b\b\b\b\b\b\b\b\b\b\b\bHiGHS 1.10.0: optimal solution; objective 88.2\n",
            "2 simplex iterations\n",
            "0 barrier iterations\n",
            "{'BEEF': 0, 'CHK': 0, 'FISH': 0, 'HAM': 0, 'MCH': 46.66666666666666, 'MTL': 0, 'SPG': 0, 'TUR': 0}\n"
          ]
        }
      ],
      "source": [
        "ampl.solve(solver='highs')\n",
        "buy = ampl.var['Buy'].to_dict()\n",
        "print(buy)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qgiJ-C9w8Iy2"
      },
      "source": [
        "Naturally, the result is the same as before. This time we have saved the `Buy` data into a Python dictionary later displayed.\n",
        "\n",
        "Now suppose that we want to make the following enhancements. To promote variety, the weekly diet must contain between 2 and 10 packages of each food. The amount of sodium and calories in each package is also given; total sodium must not exceed 40,000 mg, and total calories must be between 16,000 and 24,000. All of these changes can be made through a few modifications to the data. Putting this new data in the Python data structures and run Ampl again:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "UgySW6oZ8Iy2",
        "outputId": "2f17cd16-e0c7-488a-99f2-560022281af2"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "HiGHS 1.10.0: \b\b\b\b\b\b\b\b\b\b\b\b\b\bHiGHS 1.10.0: infeasible problem\n",
            "4 simplex iterations\n",
            "0 barrier iterations\n",
            "\n",
            "suffix dunbdd OUT;\n",
            "Solve status: infeasible\n"
          ]
        }
      ],
      "source": [
        "# New data (diet2)\n",
        "foods = ['BEEF', 'CHK', 'FISH', 'HAM', 'MCH', 'MTL', 'SPG', 'TUR']\n",
        "nutrs = ['A', 'B1', 'B2', 'C', 'NA', 'CAL']\n",
        "\n",
        "food_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"BEEF\", 3.19, 2, 10),\n",
        "        (\"CHK\", 2.59, 2, 10),\n",
        "        (\"FISH\", 2.29, 2, 10),\n",
        "        (\"HAM\", 2.89, 2, 10),\n",
        "        (\"MCH\", 1.89, 2, 10),\n",
        "        (\"MTL\", 1.99, 2, 10),\n",
        "        (\"SPG\", 1.99, 2, 10),\n",
        "        (\"TUR\", 2.49, 2, 10),\n",
        "    ],\n",
        "    columns=[\"FOOD\", \"cost\", \"f_min\", \"f_max\"],\n",
        ").set_index(\"FOOD\")\n",
        "\n",
        "# Create a pandas.DataFrame with data for n_min, n_max\n",
        "nutr_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"A\", 700, 20000),\n",
        "        (\"C\", 700, 20000),\n",
        "        (\"B1\", 700, 20000),\n",
        "        (\"B2\", 700, 20000),\n",
        "        (\"NA\", 0, 40000),\n",
        "        (\"CAL\", 16000, 24000),\n",
        "    ],\n",
        "    columns=[\"NUTR\", \"n_min\", \"n_max\"],\n",
        ").set_index(\"NUTR\")\n",
        "\n",
        "amt_df = pd.DataFrame(\n",
        "    np.array(\n",
        "        [\n",
        "            [60, 8, 8, 40, 15, 70, 25, 60],\n",
        "            [20, 0, 10, 40, 35, 30, 50, 20],\n",
        "            [10, 20, 15, 35, 15, 15, 25, 15],\n",
        "            [15, 20, 10, 10, 15, 15, 15, 10],\n",
        "            [928, 2180, 945, 278, 1182, 896, 1329, 1397],\n",
        "            [295, 770, 440, 430, 315, 400, 370, 450],\n",
        "        ]\n",
        "    ),\n",
        "    columns=food_df.index.to_list(),\n",
        "    index=nutr_df.index.to_list(),\n",
        ")\n",
        "\n",
        "# Load data into Ampl\n",
        "ampl.set[\"FOOD\"] = foods\n",
        "ampl.set[\"NUTR\"] = nutrs\n",
        "ampl.set_data(food_df, \"FOOD\")\n",
        "ampl.set_data(nutr_df, \"NUTR\")\n",
        "ampl.param[\"amt\"] = amt_df\n",
        "\n",
        "# Run Ampl\n",
        "ampl.solve(solver=\"highs\")\n",
        "print(\"Solve status: \" + ampl.solve_result)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EnnKpEAq8Iy2"
      },
      "source": [
        "The message infeasible problem tells us that we have constrained the diet too tightly; there is no way that all of the restrictions can be satisfied. This can be checked by examining the attribute `ampl.solve_result`.\n",
        "\n",
        "AMPL lets us examine a variety of values produced by a solver as it attempts to find a solution. We can use marginal (or dual) values to investigate the sensitivity of an optimum solution to changes in the constraints. Here there is no optimum, but the solver does return the last solution that it found while attempting to satisfy the constraints. We can look for the source of the infeasibility by displaying some values associated with this solution:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SnbLAGXk8Iy3",
        "outputId": "f67a9e21-5dd0-4aef-c4f7-64beadcb7c03"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            ":   Diet.lb  Diet.body Diet.ub    :=\n",
            "A       700     700      20000\n",
            "B1      700     700      20000\n",
            "B2      700     700      20000\n",
            "C       700    1633.33   20000\n",
            "CAL   16000   14700      24000\n",
            "NA        0   55160      40000\n",
            ";\n",
            "\n"
          ]
        }
      ],
      "source": [
        "ampl.display(\"Diet.lb, Diet.body, Diet.ub\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CvQLW2Kn8Iy3"
      },
      "source": [
        "For each nutrient, `Diet.body` is the sum of the terms `amt[i,j] * Buy[j]` in the constraint `Diet[i]`. The `Diet.lb` and `Diet.ub` values are the \"lower bounds\" and \"upper bounds\" on the sum in `Diet[i]` - in this case, just the values `n_min[i]` and `n_max[i]`. We can see that the diet returned by the solver supply the minimum vitamin B2, while the amount of sodium (NA) has overpassed its upper bound.\n",
        "\n",
        "At this point, there are two obvious choices: we could require less B2 or we could allow more sodium. If we try the latter, and relax the sodium limit to 50,000 mg, a feasible solution becomes possible:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "yHrX1ctN8Iy3",
        "outputId": "96c635fe-d96e-476c-9467-361d20a229d0"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "HiGHS 1.10.0: \b\b\b\b\b\b\b\b\b\b\b\b\b\bHiGHS 1.10.0: optimal solution; objective 117.8989859\n",
            "1 simplex iterations\n",
            "0 barrier iterations\n",
            "Buy [*] :=\n",
            "BEEF   5.22693\n",
            " CHK   2\n",
            "FISH   2\n",
            " HAM  10\n",
            " MCH  10\n",
            " MTL  10\n",
            " SPG   9.43973\n",
            " TUR   2\n",
            ";\n",
            "\n",
            "Solve status: solved\n"
          ]
        }
      ],
      "source": [
        "ampl.param[\"n_max\"][\"NA\"] = 50000;\n",
        "ampl.solve(solver=\"highs\")\n",
        "ampl.display(\"Buy\")\n",
        "print(\"Solve status: \" + ampl.solve_result)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U-E4sGiZ8Iy3"
      },
      "source": [
        "This is at least a start toward a palatable diet, although we have to spend \\$117.89, compared to \\$88.20 for the original, less restricted case. Clearly it would be easy, now that the model is set up, to try many other possibilities. Section 11.3 of the AMPL book describes ways to quickly change the data and re-solve.\n",
        "\n",
        "One still disappointing aspect of the solution is the need to buy 5.22693 packages of beef, and 9.43973 of spaghetti. How can we find the best possible solution in terms of whole packages? You might think that we could simply round the optimal values to whole numbers - or integers, as they're often called in the context of optimization - but it is not so easy to do so in a feasible way. Using AMPL to modify the reported solution , we can observe that rounding up to 6 packages of beef and 10 of spaghetti, for example, will violate the sodium limit:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Cy2xYRSJ8Iy3",
        "outputId": "79820019-450b-4337-f814-5c3fd7038ef0"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            ":   Diet.lb Diet.body Diet.ub    :=\n",
            "A       700     2012    20000\n",
            "B1      700     1060    20000\n",
            "B2      700      720    20000\n",
            "C       700     1730    20000\n",
            "CAL   16000    20240    24000\n",
            "NA        0    51462    50000\n",
            ";\n",
            "\n"
          ]
        }
      ],
      "source": [
        "ampl.var[\"Buy\"][\"BEEF\"] = 6\n",
        "ampl.var[\"Buy\"][\"SPG\"] = 10\n",
        "\n",
        "ampl.display(\"Diet.lb, Diet.body, Diet.ub\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gBN3ZXRN8Iy3"
      },
      "source": [
        "You can similarly check that rounding the solution down to 5 of beef and 9 of spaghetti will provide insufficient vitamin B2. Rounding one up and the other down doesn't work either. With enough experimenting you can find a nearby all-integer solution that does satisfy the constraints, but still you will have no guarantee that it is the least-cost allinteger solution.\n",
        "\n",
        "AMPL does provide for putting the integrality restriction directly into the declaration of the variables:\n",
        "```ampl\n",
        "var Buy {j in FOOD} integer >= f_min[j], <= f_max[j];\n",
        "```\n",
        "This will only help, however, if you use a solver that can deal with problems whose variables must be integers. The HiGHS solver can handle these so called integer programs. If we add *integer* to the declaration of variable `Buy` as above, save the resulting model in the file `dieti.mod`, and add the higher sodium limit to the data (diet2a), then we can re-solve as follows:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FUyT_G1b8Iy3",
        "outputId": "a9d499fe-4430-40c9-fe2a-82b77d3a0faf"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Writing dieti.mod\n"
          ]
        }
      ],
      "source": [
        "%%writefile dieti.mod\n",
        "set NUTR;\n",
        "set FOOD;\n",
        "\n",
        "param cost {FOOD} > 0;\n",
        "param f_min {FOOD} >= 0;\n",
        "param f_max {j in FOOD} >= f_min[j];\n",
        "\n",
        "param n_min {NUTR} >= 0;\n",
        "param n_max {i in NUTR} >= n_min[i];\n",
        "\n",
        "param amt {NUTR,FOOD} >= 0;\n",
        "\n",
        "var Buy {j in FOOD} integer >= f_min[j], <= f_max[j];\n",
        "\n",
        "minimize Total_Cost:  sum {j in FOOD} cost[j] * Buy[j];\n",
        "\n",
        "subject to Diet {i in NUTR}:\n",
        "   n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i];"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "R76alwsQ8Iy4",
        "outputId": "59a51138-9203-4e43-b776-49873cc8136c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "HiGHS 1.10.0: \b\b\b\b\b\b\b\b\b\b\b\b\b\bHiGHS 1.10.0: optimal solution; objective 119.3\n",
            "9 simplex iterations\n",
            "1 branching nodes\n",
            "Buy [*] :=\n",
            "BEEF   9\n",
            " CHK   2\n",
            "FISH   2\n",
            " HAM   8\n",
            " MCH  10\n",
            " MTL  10\n",
            " SPG   7\n",
            " TUR   2\n",
            ";\n",
            "\n"
          ]
        }
      ],
      "source": [
        "ampl.reset()\n",
        "ampl.read(\"dieti.mod\")\n",
        "\n",
        "# New data (diet2a)\n",
        "foods = ['BEEF', 'CHK', 'FISH', 'HAM', 'MCH', 'MTL', 'SPG', 'TUR']\n",
        "nutrs = ['A', 'B1', 'B2', 'C', 'NA', 'CAL']\n",
        "\n",
        "food_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"BEEF\", 3.19, 2, 10),\n",
        "        (\"CHK\", 2.59, 2, 10),\n",
        "        (\"FISH\", 2.29, 2, 10),\n",
        "        (\"HAM\", 2.89, 2, 10),\n",
        "        (\"MCH\", 1.89, 2, 10),\n",
        "        (\"MTL\", 1.99, 2, 10),\n",
        "        (\"SPG\", 1.99, 2, 10),\n",
        "        (\"TUR\", 2.49, 2, 10),\n",
        "    ],\n",
        "    columns=[\"FOOD\", \"cost\", \"f_min\", \"f_max\"],\n",
        ").set_index(\"FOOD\")\n",
        "\n",
        "# Create a pandas.DataFrame with data for n_min, n_max\n",
        "nutr_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"A\", 700, 20000),\n",
        "        (\"C\", 700, 20000),\n",
        "        (\"B1\", 700, 20000),\n",
        "        (\"B2\", 700, 20000),\n",
        "        (\"NA\", 0, 50000),\n",
        "        (\"CAL\", 16000, 24000),\n",
        "    ],\n",
        "    columns=[\"NUTR\", \"n_min\", \"n_max\"],\n",
        ").set_index(\"NUTR\")\n",
        "\n",
        "amt_df = pd.DataFrame(\n",
        "    np.array(\n",
        "        [\n",
        "            [60, 8, 8, 40, 15, 70, 25, 60],\n",
        "            [20, 0, 10, 40, 35, 30, 50, 20],\n",
        "            [10, 20, 15, 35, 15, 15, 25, 15],\n",
        "            [15, 20, 10, 10, 15, 15, 15, 10],\n",
        "            [928, 2180, 945, 278, 1182, 896, 1329, 1397],\n",
        "            [295, 770, 440, 430, 315, 400, 370, 450],\n",
        "        ]\n",
        "    ),\n",
        "    columns=food_df.index.to_list(),\n",
        "    index=nutr_df.index.to_list(),\n",
        ")\n",
        "\n",
        "# Load data into Ampl\n",
        "ampl.set[\"FOOD\"] = foods\n",
        "ampl.set[\"NUTR\"] = nutrs\n",
        "ampl.set_data(food_df, \"FOOD\")\n",
        "ampl.set_data(nutr_df, \"NUTR\")\n",
        "ampl.param[\"amt\"] = amt_df\n",
        "\n",
        "# Run Ampl\n",
        "ampl.solve(solver=\"highs\")\n",
        "ampl.display(\"Buy\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9ZV55KCI8Iy4"
      },
      "source": [
        "Since integrality is an added constraint, it is no surprise that the best integer solution costs about \\$1.40 more than the best \"continuous\" one. But the difference between the diets is unexpected; the amounts of 3 foods change, each by two or more packages. In general, integrality and other \"discrete\" restrictions make solutions for a model much harder to find.\n",
        "\n",
        "## Generalizations to blending, economics and scheduling\n",
        "\n",
        "Your personal experience probably suggests that diet models are not widely used by people to choose their dinners. These models would be much better suited to situations in which packaging and personal preferences don't play such a prominent role — for example , the blending of animal feed or perhaps food for college dining halls.\n",
        "\n",
        "The diet model is a convenient, intuitive example of a linear programming formulation that appears in many contexts. Suppose that we rewrite the model in a more general way:\n",
        "\n",
        "```ampl\n",
        "set INPUT;            # inputs\n",
        "set OUTPUT;           # outputs\n",
        "param cost {INPUT} > 0;\n",
        "param in_min {INPUT} >= 0;\n",
        "param in_max {j in INPUT} >= in_min[j];\n",
        "param out_min {OUTPUT} >= 0;\n",
        "param out_max {i in OUTPUT} >= out_min[i];\n",
        "param io {OUTPUT,INPUT} >= 0;\n",
        "var X {j in INPUT} >= in_min[j], <= in_max[j];\n",
        "minimize Total_Cost: sum {j in INPUT} cost[j] * X[j];\n",
        "subject to Outputs {i in OUTPUT}:\n",
        "\tout_min[i] <= sum {j in INPUT} io[i,j] * X[j] <= out_max[i];\n",
        "```\n",
        "\n",
        "The objects that were called *foods* and *nutrients* in the diet model are now referred to more generically as *inputs* and *outputs*. For each input `j`, we must decide to use a quantity `X[j]` that lies between `in_min[j]` and `in_max[j]`; as a result we incur a cost equal to `cost[j]·X[j]`, and we create `io[i,j]·X[j]` units of each output `i`. Our goal is to find the least-cost combination of inputs that yields, for each output `i`, an amount between `out_min[i]` and `out_max[i]`.\n",
        "\n",
        "In one common class of applications for this model, the inputs are raw materials to be mixed together. The outputs are qualities of the resulting **blend**. The raw materials could be the components of an animal feed, but they could equally well be the crude oil derivatives that are blended to make gasoline, or the different kinds of coal that are mixed as input to a coke oven. The qualities can be amounts of something (sodium or calories for animal feed), or more complex measures (vapor pressure or octane rating for gasoline), or even physical properties such as weight and volume.\n",
        "\n",
        "In another well-known application, the inputs are production activities of some sector of an **economy**, and the outputs are various products. The `in_min` and `in_max` parameters are limits on the levels of the activities, while `out_min` and `out_max` are regulated by demands. Thus the goal is to find levels of the activities that meet demand at the lowest cost. This interpretation is related to the concept of an economic equilibrium.\n",
        "\n",
        "In still another, quite different application, the inputs are **work schedules**, and the outputs correspond to hours worked on certain days of a month. For a particular work schedule `j`, `io[i,j]` is the number of hours that a person following schedule `j` will work on day `i` (zero if none), `cost[j]` is the monthly salary for a person following schedule `j`, and `X[j]` is the number of workers assigned that schedule. Under this interpretation, the objective becomes the total cost of the monthly payroll, while the constraints say that for each day `i`, the total number of workers assigned to work that day must lie between the limits `out_min[i]` and `out_max[i]`. The same approach can be used in a variety of other scheduling contexts, where the hours, days or months are replaced by other periods of time.\n",
        "\n",
        "Although linear programming can be very useful in applications like these, we need to keep in mind the assumptions that underlie the LP model. We have already mentioned the \"continuity\" assumption whereby `X[j]` is allowed to take on any value between `in_min[j]` and `in_max[j]`. This may be a lot more reasonable for blending than for scheduling.\n",
        "\n",
        "As another example, in writing the objective as\n",
        "```ampl\n",
        "sum {j in INPUT} cost[j] * X[j]\n",
        "```\n",
        "we are assuming \"linearity of costs\", that is, that the cost of an input is proportional to the amount of the input used, and that the total cost is the sum of the inputs' individual costs.\n",
        "\n",
        "In writing the constraints as\n",
        "```ampl\n",
        "out_min[i] <= sum {j in INPUT} io[i,j] * X[j] <= out_max[i]\n",
        "```\n",
        "we are also assuming that the yield of an output `i` from a particular input is proportional to the amount of the input used, and that the total yield of an output `i` is the sum of the yields from the individual inputs. This \"linearity of yield\" assumption poses no problem when the inputs are schedules, and the outputs are hours worked. But in the blending example, linearity is a physical assumption about the nature of the raw materials and the qualities, which may or may not hold. In early applications to refineries, for example, it was recognized that the addition of lead as an input had a nonlinear effect on the quality known as octane rating in the resulting blend.\n",
        "\n",
        "AMPL makes it easy to express discrete or nonlinear models, but any departure from continuity or linearity is likely to make an optimal solution much harder to obtain. At the least, it takes a more powerful solver to optimize the resulting mathematical programs.\n",
        "\n",
        "## Bibliography\n",
        "\n",
        "* George B. Dantzig, \"The Diet Problem.\" Interfaces 20, 4 (1990) pp. 43-47. An entertaining account of the origins of the diet problem.\n",
        "\n",
        "* Robert Fourer, David M. Gay, and Brian W. Kernighan, \"AMPL: A Modeling Language for Mathematical Programming (2nd edition).\" Cengage Learning (2002).\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "* Susan Garner Garille and Saul I. Gass, \"Stigler's Diet Problem Revisited.\" Operations Research 49, 1 (2001) pp. 1-13. A review of the diet problem's origins and its influence over the years on linear programming and on nutritionists.\n",
        "\n",
        "* Said S. Hilal and Warren Erikson, \"Matching Supplies to Save Lives: Linear Programming the Production of Heart Valves.\" Interfaces 11, 6 (1981) pp. 48-56. A less appetizing equivalent of the diet problem, involving the choice of pig heart suppliers."
      ]
    }
  ],
  "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.12"
    },
    "vscode": {
      "interpreter": {
        "hash": "ee4198aac13dae9b6fccd443041680d99c2643c4956ef7c80ff6dcc0057cb523"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}