# 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

name | cost | calories | protein | fat | sodium | |
---|---|---|---|---|---|---|

Any | Any | Any | Any | Any | Any | |

1 | hamburger | 2.49 | 410 | 24 | 26 | 730 |

2 | chicken | 2.89 | 420 | 32 | 10 | 1190 |

3 | hot dog | 1.5 | 560 | 20 | 32 | 1800 |

4 | fries | 1.89 | 380 | 4 | 19 | 270 |

5 | macaroni | 2.09 | 320 | 12 | 10 | 930 |

6 | pizza | 1.99 | 320 | 15 | 12 | 820 |

7 | salad | 2.49 | 320 | 31 | 12 | 1230 |

8 | milk | 0.89 | 100 | 8 | 2.5 | 125 |

9 | ice cream | 1.59 | 330 | 8 | 10 | 180 |

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

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

name | min | max | |
---|---|---|---|

Any | Any | Any | |

1 | calories | 1800 | 2200 |

2 | protein | 91 | Inf |

3 | fat | 0 | 65 |

4 | sodium | 0 | 1779 |

## 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[salad]
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[salad] ≥ 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
salad = 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.

This tutorial was generated using Literate.jl. View the source `.jl`

file on GitHub.