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")