The diet problem

This tutorial solves the classic "diet problem", also known as the Stigler diet.

Required packages

using JuMP
import DataFrames
import GLPK

Formulation

Suppose we wish to cook a nutritionally balanced meal by choosing the quantity of each food $f$ to eat from a set of foods $F$ in our kitchen.

Each food $f$ has a cost, $c_f$, as well as a macronutrient profile $a_{m,f}$ for each macronutrient $m \in M$.

Because we care about a nutritionally balanced meal, we set some minimum and maximum limits for each nutrient, which we denote $l_m$ and $u_m$ respectively.

Furthermore, because we are optimizers, we seek the minimum cost solution.

With a little effort, we can formulate our dinner problem as the following linear program:

\[\begin{aligned} \min & \sum\limits_{f \in F} c_f x_f \\ \text{s.t.}\ \ & l_m \le \sum\limits_{f \in F} a_{m,f} x_f \le u_m, && \forall m \in M \\ & x_f \ge 0, && \forall f \in F \end{aligned}\]

In the rest of this tutorial, we will create and solve this problem in JuMP, and learn what we should cook for dinner.

Data

First, we need some data for the problem:

foods = DataFrames.DataFrame(
    [
        "hamburger" 2.49 410 24 26 730
        "chicken" 2.89 420 32 10 1190
        "hot dog" 1.50 560 20 32 1800
        "fries" 1.89 380 4 19 270
        "macaroni" 2.09 320 12 10 930
        "pizza" 1.99 320 15 12 820
        "salad" 2.49 320 31 12 1230
        "milk" 0.89 100 8 2.5 125
        "ice cream" 1.59 330 8 10 180
    ],
    ["name", "cost", "calories", "protein", "fat", "sodium"],
)

9 rows × 6 columns

namecostcaloriesproteinfatsodium
AnyAnyAnyAnyAnyAny
1hamburger2.494102426730
2chicken2.8942032101190
3hot dog1.556020321800
4fries1.89380419270
5macaroni2.093201210930
6pizza1.993201512820
7salad2.4932031121230
8milk0.8910082.5125
9ice cream1.59330810180

Here, $F$ is foods.name and $c_f$ is foods.cost. (We're also playing a bit loose the term "macronutrient" by including calories and sodium.)

Tip

Although we hard-coded the data here, you could also read it in from a file. See Getting started with data and plotting for details.

We also need our minimum and maximum limits:

limits = DataFrames.DataFrame(
    [
        "calories" 1800 2200
        "protein" 91 Inf
        "fat" 0 65
        "sodium" 0 1779
    ],
    ["name", "min", "max"],
)

4 rows × 3 columns

nameminmax
AnyAnyAny
1calories18002200
2protein91Inf
3fat065
4sodium01779

JuMP formulation

Now we're ready to convert our mathematical formulation into a JuMP model.

First, create a new JuMP model. Since we have a linear program, we'll use GLPK as our optimizer:

model = Model(GLPK.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: GLPK

Next, we create a set of decision variables x, indexed over the foods in the data DataFrame. Each x has a lower bound of 0.

@variable(model, x[foods.name] >= 0)
1-dimensional DenseAxisArray{VariableRef,1,...} with index sets:
    Dimension 1, Any["hamburger", "chicken", "hot dog", "fries", "macaroni", "pizza", "salad", "milk", "ice cream"]
And data, a 9-element Vector{VariableRef}:
 x[hamburger]
 x[chicken]
 x[hot dog]
 x[fries]
 x[macaroni]
 x[pizza]
 x[salad]
 x[milk]
 x[ice cream]

Our objective is to minimze the total cost of purchasing food. We can write that as a sum over the rows in data.

@objective(
    model,
    Min,
    sum(food["cost"] * x[food["name"]] for food in eachrow(foods)),
)

\[ 2.49 x_{hamburger} + 2.89 x_{chicken} + 1.5 x_{hot dog} + 1.89 x_{fries} + 2.09 x_{macaroni} + 1.99 x_{pizza} + 2.49 x_{salad} + 0.89 x_{milk} + 1.59 x_{ice cream} \]

For the next component, we need to add a constraint that our total intake of each component is within the limits contained in the limits DataFrame. To make this more readable, we introduce a JuMP @expression

for limit in eachrow(limits)
    intake = @expression(
        model,
        sum(food[limit["name"]] * x[food["name"]] for food in eachrow(foods)),
    )
    @constraint(model, limit.min <= intake <= limit.max)
end

What does our model look like?

print(model)
Min 2.49 x[hamburger] + 2.89 x[chicken] + 1.5 x[hot dog] + 1.89 x[fries] + 2.09 x[macaroni] + 1.99 x[pizza] + 2.49 x[salad] + 0.89 x[milk] + 1.59 x[ice cream]
Subject to
 410 x[hamburger] + 420 x[chicken] + 560 x[hot dog] + 380 x[fries] + 320 x[macaroni] + 320 x[pizza] + 320 x[salad] + 100 x[milk] + 330 x[ice cream] ∈ [1800.0, 2200.0]
 24 x[hamburger] + 32 x[chicken] + 20 x[hot dog] + 4 x[fries] + 12 x[macaroni] + 15 x[pizza] + 31 x[salad] + 8 x[milk] + 8 x[ice cream] ∈ [91.0, Inf]
 26 x[hamburger] + 10 x[chicken] + 32 x[hot dog] + 19 x[fries] + 10 x[macaroni] + 12 x[pizza] + 12 x[salad] + 2.5 x[milk] + 10 x[ice cream] ∈ [0.0, 65.0]
 730 x[hamburger] + 1190 x[chicken] + 1800 x[hot dog] + 270 x[fries] + 930 x[macaroni] + 820 x[pizza] + 1230 x[salad] + 125 x[milk] + 180 x[ice cream] ∈ [0.0, 1779.0]
 x[hamburger] ≥ 0.0
 x[chicken] ≥ 0.0
 x[hot dog] ≥ 0.0
 x[fries] ≥ 0.0
 x[macaroni] ≥ 0.0
 x[pizza] ≥ 0.0
 x[salad] ≥ 0.0
 x[milk] ≥ 0.0
 x[ice cream] ≥ 0.0

Solution

optimize!(model)


solution_summary(model)
* Solver : GLPK

* Status
  Termination status : OPTIMAL
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Message from the solver:
  "Solution is optimal"

* Candidate solution
  Objective value      : 11.82886111111111
  Objective bound      : -Inf
  Dual objective value : 11.82886111111111

* Work counters
  Solve time (sec)   : 0.00004

Success! We found an optimal solution. Let's see what the optimal solution is:

for food in foods.name
    println(food, " = ", value(x[food]))
end
hamburger = 0.6045138888888888
chicken = 0.0
hot dog = 0.0
fries = 0.0
macaroni = 0.0
pizza = 0.0
salad = 0.0
milk = 6.9701388888888935
ice cream = 2.591319444444441

That's a lot of milk and ice cream! And sadly, we only get 0.6 of a hamburger.

Problem modification

JuMP makes it easy to take an existing model and modify it by adding extra constraints. Let's see what happens if we add a constraint that we can buy at most 6 units of milk or ice cream combined.

@constraint(model, x["milk"] + x["ice cream"] <= 6)

\[ x_{milk} + x_{ice cream} \leq 6.0 \]

optimize!(model)


solution_summary(model)
* Solver : GLPK

* Status
  Termination status : INFEASIBLE
  Primal status      : NO_SOLUTION
  Dual status        : INFEASIBILITY_CERTIFICATE
  Message from the solver:
  "No feasible primal-dual solution exists."

* Candidate solution
  Objective value      : 11.828861111111111
  Objective bound      : -Inf
  Dual objective value : 3.5614583333333307

* Work counters
  Solve time (sec)   : 0.00004

Uh oh! The exists no feasible solution to our problem. Looks like we're stuck eating ice cream for dinner.


Tip

This tutorial was generated using Literate.jl. View the source .jl file on GitHub.