# The diet problem

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

## Required packages

This tutorial requires the following packages:

using JuMP
import DataFrames
import HiGHS

## 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 macro-nutrient profile $a_{m,f}$ for each macro-nutrient $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×6 DataFrame
Rownamecostcaloriesproteinfatsodium
AnyAnyAnyAnyAnyAny
1hamburger2.494102426730
2chicken2.8942032101190
3hot dog1.556020321800
4fries1.89380419270
5macaroni2.093201210930
6pizza1.993201512820
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 "macro-nutrient" 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×3 DataFrame
Rownameminmax
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 HiGHS as our optimizer:

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

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

Our objective is to minimize 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)),
);

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[milk] ≥ 0.0
x[ice cream] ≥ 0.0

## Solution

Let's optimize and take a look at the solution:

optimize!(model)
solution_summary(model)
* Solver : HiGHS

* Status
Result count       : 1
Termination status : OPTIMAL
Message from the solver:
"kHighsModelStatusOptimal"

* Candidate solution (result #1)
Primal status      : FEASIBLE_POINT
Dual status        : FEASIBLE_POINT
Objective value    : 1.18289e+01
Objective bound    : 0.00000e+00
Relative gap       : Inf
Dual objective value : 1.18289e+01

* Work counters
Solve time (sec)   : 3.02315e-04
Simplex iterations : 6
Barrier iterations : 0
Node count         : -1


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.6045138888888871
chicken = 0.0
hot dog = 0.0
fries = 0.0
macaroni = 0.0
pizza = 0.0
milk = 6.9701388888888935
ice cream = 2.5913194444444447

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

We can also use the function Containers.rowtable to easily convert the result into a DataFrame:

table = Containers.rowtable(value, x; header = [:food, :quantity])
solution = DataFrames.DataFrame(table)
9×2 DataFrame
Rowfoodquantity
StringFloat64
1hamburger0.604514
2chicken0.0
3hot dog0.0
4fries0.0
5macaroni0.0
6pizza0.0
8milk6.97014
9ice cream2.59132

This makes it easy to perform analyses our solution:

filter!(row -> row.quantity > 0.0, solution)
3×2 DataFrame
Rowfoodquantity
StringFloat64
1hamburger0.604514
2milk6.97014
3ice cream2.59132

## 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)
optimize!(model)
solution_summary(model)
* Solver : HiGHS

* Status
Result count       : 1
Termination status : INFEASIBLE
Message from the solver:
"kHighsModelStatusInfeasible"

* Candidate solution (result #1)
Primal status      : NO_SOLUTION
Dual status        : INFEASIBILITY_CERTIFICATE
Objective value    : 1.18289e+01
Objective bound    : 0.00000e+00
Relative gap       : Inf
Dual objective value : 3.56146e+00

* Work counters
Solve time (sec)   : 5.18560e-04
Simplex iterations : 0
Barrier iterations : 0
Node count         : -1


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

Tip