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
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
# --------------------------------------------------------------------------
# Source file provided under Apache License, Version 2.0, January 2004,
# http://www.apache.org/licenses/
# (c) Copyright IBM Corp. 2015, 2018
# --------------------------------------------------------------------------

# The goal of the diet problem is to select a set of foods that satisfies
# a set of daily nutritional requirements at minimal cost.
# Source of data: http://www.neos-guide.org/content/diet-problem-solver

from collections import namedtuple

from docplex.mp.model import Model
from docplex.util.environment import get_environment
# ----------------------------------------------------------------------------
# Initialize the problem data
# ----------------------------------------------------------------------------

FOODS = [
    ("Roasted Chicken", 0.84, 0, 10),
    ("Spaghetti W/ Sauce", 0.78, 0, 10),
    ("Tomato,Red,Ripe,Raw", 0.27, 0, 10),
    ("Apple,Raw,W/Skin", .24, 0, 10),
    ("Grapes", 0.32, 0, 10),
    ("Chocolate Chip Cookies", 0.03, 0, 10),
    ("Lowfat Milk", 0.23, 0, 10),
    ("Raisin Brn", 0.34, 0, 10),
    ("Hotdog", 0.31, 0, 10)
]

NUTRIENTS = [
    ("Calories", 2000, 2500),
    ("Calcium", 800, 1600),
    ("Iron", 10, 30),
    ("Vit_A", 5000, 50000),
    ("Dietary_Fiber", 25, 100),
    ("Carbohydrates", 0, 300),
    ("Protein", 50, 100)
]

FOOD_NUTRIENTS = [
    ("Roasted Chicken", 277.4, 21.9, 1.8, 77.4, 0, 0, 42.2),
    ("Spaghetti W/ Sauce", 358.2, 80.2, 2.3, 3055.2, 11.6, 58.3, 8.2),
    ("Tomato,Red,Ripe,Raw", 25.8, 6.2, 0.6, 766.3, 1.4, 5.7, 1),
    ("Apple,Raw,W/Skin", 81.4, 9.7, 0.2, 73.1, 3.7, 21, 0.3),
    ("Grapes", 15.1, 3.4, 0.1, 24, 0.2, 4.1, 0.2),
    ("Chocolate Chip Cookies", 78.1, 6.2, 0.4, 101.8, 0, 9.3, 0.9),
    ("Lowfat Milk", 121.2, 296.7, 0.1, 500.2, 0, 11.7, 8.1),
    ("Raisin Brn", 115.1, 12.9, 16.8, 1250.2, 4, 27.9, 4),
    ("Hotdog", 242.1, 23.5, 2.3, 0, 0, 18, 10.4)
]

Food = namedtuple("Food", ["name", "unit_cost", "qmin", "qmax"])
Nutrient = namedtuple("Nutrient", ["name", "qmin", "qmax"])


# ----------------------------------------------------------------------------
# Build the model
# ----------------------------------------------------------------------------


def build_diet_model(mdl, **kwargs):
    ints = kwargs.pop('ints', False)

    # Create tuples with named fields for foods and nutrients
    foods = [Food(*f) for f in FOODS]
    nutrients = [Nutrient(*row) for row in NUTRIENTS]

    food_nutrients = {(fn[0], nutrients[n].name):
                          fn[1 + n] for fn in FOOD_NUTRIENTS for n in range(len(NUTRIENTS))}

    # Decision variables, limited to be >= Food.qmin and <= Food.qmax
    ftype = mdl.integer_vartype if ints else mdl.continuous_vartype
    qty = mdl.var_dict(foods, ftype, lb=lambda f: f.qmin, ub=lambda f: f.qmax, name=lambda f: "q_%s" % f.name)

    # Limit range of nutrients, and mark them as KPIs
    for n in nutrients:
        amount = mdl.sum(qty[f] * food_nutrients[f.name, n.name] for f in foods)
        mdl.add_range(n.qmin, amount, n.qmax)
        mdl.add_kpi(amount, publish_name="Total %s" % n.name)

    # Minimize cost
    total_cost = mdl.sum(qty[f] * f.unit_cost for f in foods)
    mdl.add_kpi(total_cost, 'Total cost')

    # add a functional KPI , taking a model and a solution as argument
    # this KPI counts the number of foods used.
    def nb_products(mdl_, s_):
        qvs = mdl_.find_matching_vars(pattern="q_")
        return sum(1 for qv in qvs if s_[qv] >= 1e-5)

    mdl.add_kpi(nb_products, 'Nb foods')
    mdl.minimize(total_cost)

    return mdl


# ----------------------------------------------------------------------------
# Solve the model and display the result
# ----------------------------------------------------------------------------

if __name__ == '__main__':
    with Model(name="diet", log_output=True, float_precision=6) as mdl:
        build_diet_model(mdl, ints=True)
        mdl.print_information()

        s = mdl.solve()
        if s:
            qty_vars = mdl.find_matching_vars(pattern="q_")
            for fv in qty_vars:
                food_name = fv.name[2:]
                print("Buy {0:<25} = {1:9.6g}".format(food_name, fv.solution_value))

            mdl.report_kpis()
            # Save the CPLEX solution as "solution.json" program output
            with get_environment().get_output_stream("solution.json") as fp:
                mdl.solution.export(fp, "json")
        else:
            print("* model has no solution")