# The diet problem

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

## Required 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 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
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 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)
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[milk]
x[ice cream]

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

## Solution

optimize!(model)

solution_summary(model)
* Solver : HiGHS

* Status
Termination status : OPTIMAL
Primal status      : FEASIBLE_POINT
Dual status        : FEASIBLE_POINT
Message from the solver:
"kHighsModelStatusOptimal"

* Candidate solution
Objective value      : 1.18289e+01
Objective bound      : 0.00000e+00
Relative gap         : Inf
Dual objective value : 1.18289e+01

* Work counters
Solve time (sec)   : 4.21286e-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.6045138888888888
chicken = 0.0
hot dog = 0.0
fries = 0.0
macaroni = 0.0
pizza = 0.0
milk = 6.970138888888889
ice cream = 2.5913194444444447

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 : HiGHS

* Status
Termination status : INFEASIBLE
Primal status      : NO_SOLUTION
Dual status        : INFEASIBILITY_CERTIFICATE
Message from the solver:
"kHighsModelStatusInfeasible"

* Candidate solution
Objective value      : 1.18289e+01
Objective bound      : 0.00000e+00
Relative gap         : Inf
Dual objective value : 3.56146e+00

* Work counters
Solve time (sec)   : 6.37054e-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