diet.py

Can linear programming save money on the food budget of the US Army without damaging the nutritional health of members of the armed forces?

This example solves a simple variation of the well-known diet problem that was posed by George Stigler and George Dantzig: how to choose foods that satisfy nutritional requirements while minimizing costs or maximizing satiety.

Stigler solved his model “by hand” because technology at the time did not yet support more sophisticated methods. However, in 1947, Jack Laderman, of the US National Bureau of Standards, applied the simplex method (an algorithm that was recently proposed by George Dantzig) to Stigler’s model. Laderman and his team of nine linear programmers, working on desk calculators, showed that Stigler’s heuristic approximation was very close to optimal (only 24 cents per year over the optimum found by the simplex method) and thus demonstrated the practicality of the simplex method on large-scale, real-world problems.

The problem that is solved in this example is to minimize the cost of a diet that satisfies certain nutritional constraints.

  1# --------------------------------------------------------------------------
  2# Source file provided under Apache License, Version 2.0, January 2004,
  3# http://www.apache.org/licenses/
  4# (c) Copyright IBM Corp. 2015, 2018
  5# --------------------------------------------------------------------------
  6
  7# The goal of the diet problem is to select a set of foods that satisfies
  8# a set of daily nutritional requirements at minimal cost.
  9# Source of data: http://www.neos-guide.org/content/diet-problem-solver
 10
 11from collections import namedtuple
 12
 13from docplex.mp.model import Model
 14from docplex.util.environment import get_environment
 15# ----------------------------------------------------------------------------
 16# Initialize the problem data
 17# ----------------------------------------------------------------------------
 18
 19FOODS = [
 20    ("Roasted Chicken", 0.84, 0, 10),
 21    ("Spaghetti W/ Sauce", 0.78, 0, 10),
 22    ("Tomato,Red,Ripe,Raw", 0.27, 0, 10),
 23    ("Apple,Raw,W/Skin", .24, 0, 10),
 24    ("Grapes", 0.32, 0, 10),
 25    ("Chocolate Chip Cookies", 0.03, 0, 10),
 26    ("Lowfat Milk", 0.23, 0, 10),
 27    ("Raisin Brn", 0.34, 0, 10),
 28    ("Hotdog", 0.31, 0, 10)
 29]
 30
 31NUTRIENTS = [
 32    ("Calories", 2000, 2500),
 33    ("Calcium", 800, 1600),
 34    ("Iron", 10, 30),
 35    ("Vit_A", 5000, 50000),
 36    ("Dietary_Fiber", 25, 100),
 37    ("Carbohydrates", 0, 300),
 38    ("Protein", 50, 100)
 39]
 40
 41FOOD_NUTRIENTS = [
 42    ("Roasted Chicken", 277.4, 21.9, 1.8, 77.4, 0, 0, 42.2),
 43    ("Spaghetti W/ Sauce", 358.2, 80.2, 2.3, 3055.2, 11.6, 58.3, 8.2),
 44    ("Tomato,Red,Ripe,Raw", 25.8, 6.2, 0.6, 766.3, 1.4, 5.7, 1),
 45    ("Apple,Raw,W/Skin", 81.4, 9.7, 0.2, 73.1, 3.7, 21, 0.3),
 46    ("Grapes", 15.1, 3.4, 0.1, 24, 0.2, 4.1, 0.2),
 47    ("Chocolate Chip Cookies", 78.1, 6.2, 0.4, 101.8, 0, 9.3, 0.9),
 48    ("Lowfat Milk", 121.2, 296.7, 0.1, 500.2, 0, 11.7, 8.1),
 49    ("Raisin Brn", 115.1, 12.9, 16.8, 1250.2, 4, 27.9, 4),
 50    ("Hotdog", 242.1, 23.5, 2.3, 0, 0, 18, 10.4)
 51]
 52
 53Food = namedtuple("Food", ["name", "unit_cost", "qmin", "qmax"])
 54Nutrient = namedtuple("Nutrient", ["name", "qmin", "qmax"])
 55
 56
 57# ----------------------------------------------------------------------------
 58# Build the model
 59# ----------------------------------------------------------------------------
 60
 61
 62def build_diet_model(mdl, **kwargs):
 63    ints = kwargs.pop('ints', False)
 64
 65    # Create tuples with named fields for foods and nutrients
 66    foods = [Food(*f) for f in FOODS]
 67    nutrients = [Nutrient(*row) for row in NUTRIENTS]
 68
 69    food_nutrients = {(fn[0], nutrients[n].name):
 70                          fn[1 + n] for fn in FOOD_NUTRIENTS for n in range(len(NUTRIENTS))}
 71
 72    # Decision variables, limited to be >= Food.qmin and <= Food.qmax
 73    ftype = mdl.integer_vartype if ints else mdl.continuous_vartype
 74    qty = mdl.var_dict(foods, ftype, lb=lambda f: f.qmin, ub=lambda f: f.qmax, name=lambda f: "q_%s" % f.name)
 75
 76    # Limit range of nutrients, and mark them as KPIs
 77    for n in nutrients:
 78        amount = mdl.sum(qty[f] * food_nutrients[f.name, n.name] for f in foods)
 79        mdl.add_range(n.qmin, amount, n.qmax)
 80        mdl.add_kpi(amount, publish_name="Total %s" % n.name)
 81
 82    # Minimize cost
 83    total_cost = mdl.sum(qty[f] * f.unit_cost for f in foods)
 84    mdl.add_kpi(total_cost, 'Total cost')
 85
 86    # add a functional KPI , taking a model and a solution as argument
 87    # this KPI counts the number of foods used.
 88    def nb_products(mdl_, s_):
 89        qvs = mdl_.find_matching_vars(pattern="q_")
 90        return sum(1 for qv in qvs if s_[qv] >= 1e-5)
 91
 92    mdl.add_kpi(nb_products, 'Nb foods')
 93    mdl.minimize(total_cost)
 94
 95    return mdl
 96
 97
 98# ----------------------------------------------------------------------------
 99# Solve the model and display the result
100# ----------------------------------------------------------------------------
101
102if __name__ == '__main__':
103    with Model(name="diet", log_output=True, float_precision=6) as mdl:
104        build_diet_model(mdl, ints=True)
105        mdl.print_information()
106
107        s = mdl.solve()
108        if s:
109            qty_vars = mdl.find_matching_vars(pattern="q_")
110            for fv in qty_vars:
111                food_name = fv.name[2:]
112                print("Buy {0:<25} = {1:9.6g}".format(food_name, fv.solution_value))
113
114            mdl.report_kpis()
115            # Save the CPLEX solution as "solution.json" program output
116            with get_environment().get_output_stream("solution.json") as fp:
117                mdl.solution.export(fp, "json")
118        else:
119            print("* model has no solution")