Optimization of Daily Diet Using Pyomo

Optimization of Daily Diet Using Pyomo#

Prepared by: Landi Wang (lwang25@nd.edu, 2024) and Yi Liu (yliu66@nd.edu, 2024) based on prior work by Prof. Alexander Dowling

Learning Objective#

Review syntax for Pyomo package
Demonstrate the power of Pyomo in solving and optimizing complex dynamic process systems
Learn how to formulate linear optimization problem
Gain more knowledge on how changes in one attribute of the system change the optimization results

Primary Source#

This notebook was adapted from an earlier version Continuous Optimization: Linear Programming

Install Packages#

# This code cell installs packages on Colab

import sys
if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper
    helper.easy_install()
else:
    sys.path.insert(0, '../')
    import helper
helper.set_plotting_style()

--2024-10-07 22:02:11--  https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 6493 (6.3K) [text/plain]
Saving to: ‘helper.py’

helper.py             0%[                    ]       0  --.-KB/s               
helper.py           100%[===================>]   6.34K  --.-KB/s    in 0s      

2024-10-07 22:02:11 (59.7 MB/s) - ‘helper.py’ saved [6493/6493]

Installing idaes via pip...
idaes was successfully installed
idaes, version 2.6.0

Running idaes get-extensions to install Ipopt, k_aug, and more...
Checking solver versions:
Ipopt 3.13.2 (x86_64-pc-linux-gnu), ASL(20190605)

[K_AUG] 0.1.0, Part of the IDAES PSE framework
Please visit https://idaes.org/ (x86_64-pc-linux-gnu), ASL(20190605)

Couenne 0.5.8 -- an Open-Source solver for Mixed Integer Nonlinear Optimization
Mailing list: couenne@list.coin-or.org
Instructions: http://www.coin-or.org/Couenne
couenne (x86_64-pc-linux-gnu), ASL(20190605)

Bonmin 1.8.8 using Cbc 2.10.10 and Ipopt 3.13.2
bonmin (x86_64-pc-linux-gnu), ASL(20190605)

Ipopt 3.13.3 (x86_64-pc-linux-gnu), ASL(20190605)

1.0 dot_1' (x86_64-pc-linux-gnu), ASL(20190605)

import pandas as pd
import pyomo.environ as pyo
import numpy as np
import matplotlib.pyplot as plt

You want to save money eating while remaining healthy. A healthy diet requires at least 70 grams of protein, 250 grams of carbohydrates, 70 grams of fats, 35 grams of fiber per day. Now we want to add one more constraint of food volume that one is able to intake, as well as more food options to more accurately simulate real life conditions. We has set the daily intake volume to be around 90 ounces, assuming 30 ounces per meal and 3 meals per day. The following table describes food attributes retrived directly from the source.

Model Without Water#

Model Description and Formulation#

We formulated a model without water involved (as it is in the dataframe). There are four types of food: broccoli, bread, steak, and Popeyes (four piece combo). We chose these items because they are the most common, meanwhile each has its unique attributes: broccoli is only high in fiber, bread is only high in carbs, steak is high in protein and fat, and Popeyes is high in protein, fat, and decent in carbs. Among these options, Popeyes is the most well-rounded food item, providing most nutrients in abundance. Before we solve the model, we believe that the best combination would be Popeyes with some bread as well as broccoli that is enough to satisfy the fiber needs, whereas steak would be ignored because of Popeyes’ nutritional value (though sounds odd).

model = pyo.ConcreteModel()

# Define the food items
food_items = ['broccoli', 'bread', 'steak', 'Popeyes']

# Define the nutrition and price data
nutrition_data = {
    'Protein': [2.6, 3.5, 26, 76],
    'Carbohydrates': [6, 19.5, 0, 36],
    'Fat': [0.3, 0.84, 7.6, 64],
    'Fiber': [2, 1.16, 0, 0],
    'Volume': [2.6, 1.41, 3, 25],
    'Price': [2, 1.5, 6, 25]
}

# Define the bounds for each constraint
bounds = {
    'Protein': 70,
    'Carbohydrates': 250,
    'Fat': 70,
    'Fiber': 35,
    'Volume': 90
}

# Define the decision variables for the amount of each food item
model.food_quantities = pyo.Var(food_items, domain=pyo.NonNegativeReals)

# Define the objective function to minimize the total price
model.obj = pyo.Objective(expr=sum(model.food_quantities[i] * nutrition_data['Price'][j]
                                   for j, i in enumerate(food_items)), sense=pyo.minimize)

# Add constraints for total intake of protein, carbohydrates, fat, fiber, and volume
model.protein_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Protein'][j]
                                                   for j, i in enumerate(food_items)) >= bounds['Protein'])
model.carbohydrates_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Carbohydrates'][j]
                                                         for j, i in enumerate(food_items)) >= bounds['Carbohydrates'])
model.fat_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fat'][j]
                                               for j, i in enumerate(food_items)) >= bounds['Fat'])
model.fiber_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fiber'][j]
                                                 for j, i in enumerate(food_items)) >= bounds['Fiber'])
model.volume_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Volume'][j]
                                                  for j, i in enumerate(food_items)) <= bounds['Volume'])

# Solve the problem using IPOPT solver
solver = pyo.SolverFactory('ipopt')
results = solver.solve(model, tee=True)

# Print the results
print("Status:", results.solver.termination_condition)
print("Optimal solution:")
for i in food_items:
    print(f"{i}: {pyo.value(model.food_quantities[i]):.2f} units")

# Print the total price
total_price = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Price'][j]
                  for j, i in enumerate(food_items))
print(f"Total price: ${total_price:.2f}")

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:       17
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        4
                     variables with only lower bounds:        4
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        5
        inequality constraints with only lower bounds:        4
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        1

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  3.4499966e-01 2.49e+02 1.03e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  5.1239594e-01 2.49e+02 1.95e+00  -1.0 7.00e+01    -  1.67e-03 4.77e-03h  1
   2  1.0829719e+00 2.45e+02 4.47e+00  -1.0 3.81e+01    -  7.04e-03 1.54e-02h  1
   3  1.6283397e+00 2.42e+02 2.80e+00  -1.0 6.28e+01    -  2.19e-02 9.87e-03h  1
   4  6.0978594e+01 0.00e+00 1.07e+01  -1.0 6.66e+01    -  3.02e-02 1.00e+00h  1
   5  6.1024713e+01 0.00e+00 8.47e-02  -1.0 1.33e-01    -  9.94e-01 1.00e+00f  1
   6  6.0925244e+01 0.00e+00 2.00e-07  -1.7 8.10e-01    -  1.00e+00 1.00e+00f  1
   7  6.0848635e+01 0.00e+00 3.80e-04  -3.8 3.21e+00    -  9.99e-01 9.82e-01f  1
   8  6.0847087e+01 0.00e+00 1.84e-11  -5.7 4.87e-03    -  1.00e+00 1.00e+00f  1
   9  6.0847080e+01 0.00e+00 2.93e-14  -8.6 3.61e-04    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 9

                                   (scaled)                 (unscaled)
Objective...............:   6.0847079778165650e+01    6.0847079778165650e+01
Dual infeasibility......:   2.9317359951248336e-14    2.9317359951248336e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5118925268508237e-09    2.5118925268508237e-09
Overall NLP error.......:   2.5118925268508237e-09    2.5118925268508237e-09


Number of objective function evaluations             = 10
Number of objective gradient evaluations             = 10
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 10
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 10
Number of Lagrangian Hessian evaluations             = 9
Total CPU secs in IPOPT (w/o function evaluations)   =      0.002
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
Status: optimal
Optimal solution:
broccoli: 13.47 units
bread: 6.94 units
steak: 0.00 units
Popeyes: 0.94 units
Total price: $60.85

From the model output, we can see that the results are pretty similar to what we imagined at first: 0.94 units of Popeyes, 6.94 units of bread, 13.47 units of broccoli, and no steak. These units translated into normal weights dimensions are: almost one serving of the Popeyes four-piece chicken combo, 2 pounds of broccoli, and 0.6 pound of bread. Steak is not selected at all primarily because it is not a economically efficient option for its nutritional values, as the two of its competitive edges (protein and fat) are both covered by Popeyes.

Degree of Freedom Analysis#

The degree of freedom of this model is determined to be 0: there are 5 variables, one for each type of food, and five constraints, one for each type of nutrient.

Model With Water#

Model Description and Formulation#

Now we want to formulate a model with water involved. Water has no nutritional value but volume. We also added a constraint specifying the amount of water to be consumed per day: 3 bottles of 10 oz size, in order to imitate real-life situations more accurately.

model = pyo.ConcreteModel()

# Update food_items to include 'water'
food_items = ['broccoli', 'bread', 'steak', 'Popeyes', 'water']

# Update the nutrition and price data to include 'water'
nutrition_data = {
    'Protein': [2.6, 3.5, 26, 76, 0],
    'Carbohydrates': [6, 19.5, 0, 36, 0],
    'Fat': [0.3, 0.84, 7.6, 64, 0],
    'Fiber': [2, 1.16, 0, 0, 0],
    'Volume': [2.6, 1.41, 3, 25, 10],
    'Price': [2, 1.5, 6, 25, 2.5]
}

# Define the bounds for each constraint
bounds = {
    'Protein': 70,
    'Carbohydrates': 250,
    'Fat': 70,
    'Fiber': 35,
    'Volume': 90
}

# Define the decision variables for the amount of each food item
model.food_quantities = pyo.Var(food_items, domain=pyo.NonNegativeReals)

# Define the objective function to minimize the total price
model.obj = pyo.Objective(expr=sum(model.food_quantities[i] * nutrition_data['Price'][j]
                                   for j, i in enumerate(food_items)), sense=pyo.minimize)

# Add constraints for total intake of protein, carbohydrates, fat, fiber, and volume
model.protein_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Protein'][j]
                                                   for j, i in enumerate(food_items)) >= bounds['Protein'])
model.carbohydrates_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Carbohydrates'][j]
                                                         for j, i in enumerate(food_items)) >= bounds['Carbohydrates'])
model.fat_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fat'][j]
                                               for j, i in enumerate(food_items)) >= bounds['Fat'])
model.fiber_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fiber'][j]
                                                 for j, i in enumerate(food_items)) >= bounds['Fiber'])
model.volume_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Volume'][j]
                                                  for j, i in enumerate(food_items)) <= bounds['Volume'])
model.popeyes_constraint = pyo.Constraint(expr=model.food_quantities['water'] >= 3)

# Solve the problem using IPOPT solver
solver = pyo.SolverFactory('ipopt')
results = solver.solve(model, tee=True)

# Print the results
print("Status:", results.solver.termination_condition)
print("Optimal solution:")
for i in food_items:
    print(f"{i}: {pyo.value(model.food_quantities[i]):.2f} units")

# Print the total price
total_price = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Price'][j]
                  for j, i in enumerate(food_items))
print(f"Total price: ${total_price:.2f}")

def check_active_constraint(lhs_value, bound, operator):
    if operator == ">=":
        return abs(lhs_value - bound) < 1e-5
    elif operator == "<=":
        return abs(lhs_value - bound) < 1e-5
    return False

# Calculate the total intake for each nutrient
protein_intake = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Protein'][j] for j, i in enumerate(food_items))
carbohydrates_intake = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Carbohydrates'][j] for j, i in enumerate(food_items))
fat_intake = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Fat'][j] for j, i in enumerate(food_items))
fiber_intake = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Fiber'][j] for j, i in enumerate(food_items))
volume_intake = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Volume'][j] for j, i in enumerate(food_items))

# Check the active constraints
print("\nActive Constraints:")
if check_active_constraint(protein_intake, bounds['Protein'], '>='):
    print("Protein constraint is active.")
else:
    print("Protein constraint is not active.")

if check_active_constraint(carbohydrates_intake, bounds['Carbohydrates'], '>='):
    print("Carbohydrates constraint is active.")
else:
    print("Carbohydrates constraint is not active.")

if check_active_constraint(fat_intake, bounds['Fat'], '>='):
    print("Fat constraint is active.")
else:
    print("Fat constraint is not active.")

if check_active_constraint(fiber_intake, bounds['Fiber'], '>='):
    print("Fiber constraint is active.")
else:
    print("Fiber constraint is not active.")

if check_active_constraint(volume_intake, bounds['Volume'] ,'<='):
    print("Volume constraint is active.")
else:
    print("Volume constraint is not active.")

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:       19
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        5
                     variables with only lower bounds:        5
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        6
        inequality constraints with only lower bounds:        5
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        1

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  3.6999963e-01 2.49e+02 1.04e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  5.4691760e-01 2.49e+02 1.79e+00  -1.0 7.00e+01    -  1.67e-03 4.77e-03h  1
   2  1.0535935e+00 2.45e+02 3.16e+00  -1.0 5.27e+01    -  7.06e-03 1.24e-02h  1
   3  1.7059280e+00 2.43e+02 2.42e+00  -1.0 9.11e+01    -  1.18e-02 1.06e-02h  1
   4  8.0202224e+00 2.19e+02 1.96e+01  -1.0 9.52e+01    -  2.45e-02 9.62e-02h  1
   5  6.2016951e+01 2.16e+01 7.91e+00  -1.0 8.71e+01    -  1.99e-01 8.93e-01h  1
   6  6.2348101e+01 2.02e+01 5.18e+01  -1.0 6.10e+00    -  8.15e-01 5.93e-02h  1
   7  6.3479329e+01 5.20e+00 1.01e+01  -1.0 2.52e+02    -  8.47e-03 1.57e-01h  1
   8  7.0202576e+01 5.12e-02 3.94e+01  -1.0 3.24e+02    -  3.90e-01 9.83e-01h  1
   9  7.0270853e+01 9.87e-03 8.17e+00  -1.0 7.88e+00    -  1.00e+00 4.01e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  7.0271205e+01 9.84e-03 3.71e+01  -1.0 2.75e+00    -  1.00e+00 4.63e-01h  1
  11  7.0290969e+01 0.00e+00 4.24e+01  -1.0 2.07e+00    -  1.00e+00 6.30e-01h  1
  12  7.0283057e+01 0.00e+00 1.00e-06  -1.0 2.55e-01    -  1.00e+00 1.00e+00h  1
  13  7.0224950e+01 0.00e+00 7.39e+00  -2.5 3.26e-01    -  9.73e-01 1.00e+00f  1
  14  7.0218543e+01 0.00e+00 7.14e+00  -2.5 3.60e-01    -  1.00e+00 4.00e-01f  2
  15  7.0207487e+01 0.00e+00 2.83e-08  -2.5 6.23e-02    -  1.00e+00 1.00e+00f  1
  16  7.0202817e+01 0.00e+00 2.13e+00  -3.8 8.98e-02    -  7.48e-01 1.00e+00f  1
  17  7.0199636e+01 0.00e+00 1.50e-09  -3.8 2.82e-01    -  1.00e+00 1.00e+00f  1
  18  7.0199069e+01 0.00e+00 1.11e-02  -5.7 3.21e-02    -  9.54e-01 1.00e+00f  1
  19  7.0199053e+01 0.00e+00 1.84e-11  -5.7 2.02e-03    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  7.0199043e+01 0.00e+00 2.68e-14  -8.6 6.82e-04    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 20

                                   (scaled)                 (unscaled)
Objective...............:   7.0199043466935052e+01    7.0199043466935052e+01
Dual infeasibility......:   2.6815287305497230e-14    2.6815287305497230e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   3.6353123001924004e-09    3.6353123001924004e-09
Overall NLP error.......:   3.6353123001924004e-09    3.6353123001924004e-09


Number of objective function evaluations             = 23
Number of objective gradient evaluations             = 21
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 23
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 21
Number of Lagrangian Hessian evaluations             = 20
Total CPU secs in IPOPT (w/o function evaluations)   =      0.007
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
Status: optimal
Optimal solution:
broccoli: 0.02 units
bread: 30.13 units
steak: 0.00 units
Popeyes: 0.70 units
water: 3.00 units
Total price: $70.20

Active Constraints:
Protein constraint is not active.
Carbohydrates constraint is not active.
Fat constraint is active.
Fiber constraint is active.
Volume constraint is active.

In this model with an extra item, water, and an extra constraint on it that limits the amount of water consumed per day to be at least 3 units (30 oz in total). Now the solution turned into almost no broccoli, 30.13 units of bread, 0.7 unit of Popeyes, and 3 bottles of water. Translating into common units: 3 pieces of Popeyes chicken, 2.655 pounds of bread, and 3 bottles of water. This seems very unrealistic and unhealthy but at the same time satisfies daily needs, if we only look at these included variables. In this model, we also found that two of the five constraints are not active (protein and carbs), so we made the following 3-D visualization on the feasible area of this model:

Degree of Freedom Analysis#

The degree of freedom of the model is 0 as well: 6 variables, one for each type of food, and 6 constraints, one for each type of nutrient and an extra one on the amount of water consumed.

Steak Nutritional Value Modifications#

Protein Modifications#

In the following section, we want to investigate under what conditions we would choose steak as a part of our diet. We first wanted to loop through values from 26 to 61 for grams of protein in a piece of 3oz steak.

# Initialize the range of protein values for steak
protein_values = np.linspace(26, 61, 50)
selected_steak_quantities = []

# Define a function to solve the model with a given protein value for steak
def solve_model(protein_value):
    model = pyo.ConcreteModel()
    food_items = ['broccoli', 'bread', 'steak', 'Popeyes', 'water']
    nutrition_data = {
        'Protein': [2.6, 3.5, protein_value, 76, 0],  # Update protein value for steak
        'Carbohydrates': [6, 19.5, 0, 36, 0],
        'Fat': [0.3, 0.84, 7.6, 64, 0],
        'Fiber': [2, 1.16, 0, 0, 0],
        'Volume': [2.6, 1.41, 3, 25, 10],
        'Price': [2, 1.5, 6, 25, 2.5]
    }

    bounds = {
        'Protein': 70,
        'Carbohydrates': 250,
        'Fat': 70,
        'Fiber': 35,
        'Volume': 90
    }

    model.food_quantities = pyo.Var(food_items, domain=pyo.NonNegativeReals)

    model.obj = pyo.Objective(expr=sum(model.food_quantities[i] * nutrition_data['Price'][j]
                                       for j, i in enumerate(food_items)), sense=pyo.minimize)

    model.protein_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Protein'][j]
                                                       for j, i in enumerate(food_items)) >= bounds['Protein'])
    model.carbohydrates_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Carbohydrates'][j]
                                                             for j, i in enumerate(food_items)) >= bounds['Carbohydrates'])
    model.fat_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fat'][j]
                                                   for j, i in enumerate(food_items)) >= bounds['Fat'])
    model.fiber_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fiber'][j]
                                                     for j, i in enumerate(food_items)) >= bounds['Fiber'])
    model.volume_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Volume'][j]
                                                      for j, i in enumerate(food_items)) <= bounds['Volume'])
    model.water_constraint = pyo.Constraint(expr=model.food_quantities['water'] == 3)

    solver = pyo.SolverFactory('ipopt')
    results = solver.solve(model)

    # Return the selected quantity of steak
    return pyo.value(model.food_quantities['steak'])

# Loop over the range of protein values and solve the model for each
for protein_value in protein_values:
    selected_steak_quantities.append(solve_model(protein_value))

# Plot the results
plt.plot(protein_values, selected_steak_quantities, marker='o')
plt.title('Effect of Protein Value of Steak on Amount Selected')
plt.xlabel('Protein Value of Steak')
plt.ylabel('Amount of Steak Selected (units)')
plt.grid(True)
plt.show()

../../_images/f9d82e6d30b714deef2f74a67c4f93e17d0aef25fe4054f1a94e74b00499b73b.png

In this graph, we can see that even if the steak contains 61 grams of protein per 3oz, it would still not be selected as a part of the diet. Thus, protein might not be a determining factor. Next, we tried to change the amount of fat contained in a piece of 3oz steak, from 5 grams to 20 grams per 3oz.

Fat Modifications#

In this section, we looped through fat values for a 3oz steak from 5 to 20, and figure out how does the amount of steak chosen changes over the loop.

# Initialize the range of fat values for steak
fat_values = np.linspace(5, 20, 50)
selected_steak_quantities = []

# Define a function to solve the model with a given fat value for steak
def solve_model(fat_value):
    model = pyo.ConcreteModel()
    food_items = ['broccoli', 'bread', 'steak', 'Popeyes', 'water']

    nutrition_data = {
        'Protein': [2.6, 3.5, 26, 76, 0],
        'Carbohydrates': [6, 19.5, 0, 36, 0],
        'Fat': [0.3, 0.84, fat_value, 64, 0],  # Update fat value for steak
        'Fiber': [2, 1.16, 0, 0, 0],
        'Volume': [2.6, 1.41, 3, 25, 10],
        'Price': [2, 1.5, 6, 25, 2.5]
    }

    bounds = {
        'Protein': 70,
        'Carbohydrates': 250,
        'Fat': 70,
        'Fiber': 35,
        'Volume': 90
    }

    model.food_quantities = pyo.Var(food_items, domain=pyo.NonNegativeReals)

    model.obj = pyo.Objective(expr=sum(model.food_quantities[i] * nutrition_data['Price'][j]
                                       for j, i in enumerate(food_items)), sense=pyo.minimize)

    model.protein_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Protein'][j]
                                                       for j, i in enumerate(food_items)) >= bounds['Protein'])
    model.carbohydrates_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Carbohydrates'][j]
                                                             for j, i in enumerate(food_items)) >= bounds['Carbohydrates'])
    model.fat_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fat'][j]
                                                   for j, i in enumerate(food_items)) >= bounds['Fat'])
    model.fiber_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fiber'][j]
                                                     for j, i in enumerate(food_items)) >= bounds['Fiber'])
    model.volume_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Volume'][j]
                                                      for j, i in enumerate(food_items)) <= bounds['Volume'])
    model.water_constraint = pyo.Constraint(expr=model.food_quantities['water'] == 3)

    solver = pyo.SolverFactory('ipopt')
    results = solver.solve(model)

    # Return the selected quantity of stea
    return pyo.value(model.food_quantities['steak'])

# Loop over the range of fat values and solve the model for each
for fat_value in fat_values:
    selected_steak_quantities.append(solve_model(fat_value))

# Plot the results
plt.plot(fat_values, selected_steak_quantities, marker='o')
plt.title('Effect of Fat Value of Steak on Amount Selected')
plt.xlabel('Fat Value of Steak')
plt.ylabel('Amount of Steak Selected (units)')
plt.grid(True)
plt.show()

../../_images/5ab11aef64f271e31bf32862a775d7dd9aaf3e716981ed869615013010d21ca4.png

Now, in the above graph, we can see that the amount of steak selected dramatically increased to 3 when the fat value reaches around 13, and then started to decrease as the steak contains more fat but still not enough for the model to exclude Popeyes. As the fat value increases more to around 15.5, the amount of steak selected increased again and reached around 3.7. We speculate that this is because that the steak’s relative nutritional value has exceeded that of Popeyes.

Turning Point 1: Fat per 3oz steak = 14#

model = pyo.ConcreteModel()
food_items = ['broccoli', 'bread', 'steak', 'Popeyes', 'water']

nutrition_data = {
    'Protein': [2.6, 3.5, 26, 76, 0],
    'Carbohydrates': [6, 19.5, 0, 36, 0],
    'Fat': [0.3, 0.84, 14, 64, 0],
    'Fiber': [2, 1.16, 0, 0, 0],
    'Volume': [2.6, 1.41, 3, 25, 10],
    'Price': [2, 1.5, 6, 25, 2.5]
}

bounds = {
    'Protein': 70,
    'Carbohydrates': 250,
    'Fat': 70,
    'Fiber': 35,
    'Volume': 90
}

model.food_quantities = pyo.Var(food_items, domain=pyo.NonNegativeReals)

model.obj = pyo.Objective(expr=sum(model.food_quantities[i] * nutrition_data['Price'][j]
                                   for j, i in enumerate(food_items)), sense=pyo.minimize)

model.protein_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Protein'][j]
                                                   for j, i in enumerate(food_items)) >= bounds['Protein'])
model.carbohydrates_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Carbohydrates'][j]
                                                         for j, i in enumerate(food_items)) >= bounds['Carbohydrates'])
model.fat_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fat'][j]
                                               for j, i in enumerate(food_items)) >= bounds['Fat'])
model.fiber_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fiber'][j]
                                                 for j, i in enumerate(food_items)) >= bounds['Fiber'])
model.volume_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Volume'][j]
                                                  for j, i in enumerate(food_items)) <= bounds['Volume'])
model.popeyes_constraint = pyo.Constraint(expr=model.food_quantities['water'] >= 3)

solver = pyo.SolverFactory('ipopt')
results = solver.solve(model, tee=True)

# Print the results
print("Status:", results.solver.termination_condition)
print("Optimal solution:")
for i in food_items:
    print(f"{i}: {pyo.value(model.food_quantities[i]):.2f} units")

# Print the total price
total_price = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Price'][j]
                  for j, i in enumerate(food_items))
print(f"Total price: ${total_price:.2f}")

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:       19
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        5
                     variables with only lower bounds:        5
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        6
        inequality constraints with only lower bounds:        5
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        1

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  3.6999963e-01 2.49e+02 1.04e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  5.1425900e-01 2.49e+02 1.06e+00  -1.0 6.70e+01    -  1.76e-03 3.66e-03h  1
   2  1.0303972e+00 2.45e+02 1.91e+00  -1.0 5.05e+01    -  7.23e-03 1.29e-02h  1
   3  2.0012814e+00 2.42e+02 3.41e+00  -1.0 9.15e+01    -  1.12e-02 1.57e-02h  1
   4  9.0768387e+00 2.16e+02 1.25e+01  -1.0 9.51e+01    -  2.84e-02 1.08e-01h  1
   5  6.2047976e+01 2.14e+01 7.64e+00  -1.0 8.55e+01    -  2.21e-01 8.92e-01h  1
   6  6.2385147e+01 2.00e+01 8.78e+00  -1.0 6.35e+00    -  2.52e-01 5.98e-02h  1
   7  6.8333553e+01 2.68e+00 1.08e+01  -1.0 3.62e+01    -  7.33e-03 7.64e-01h  1
   8  7.0254966e+01 0.00e+00 2.22e-01  -1.0 9.42e+00    -  9.90e-01 1.00e+00h  1
   9  7.0219835e+01 0.00e+00 2.00e-07  -1.7 3.58e+00    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  7.0181446e+01 0.00e+00 3.99e-02  -3.8 1.32e+01    -  9.70e-01 4.75e-01f  1
  11  7.0139909e+01 0.00e+00 1.50e-09  -3.8 1.03e+00    -  1.00e+00 1.00e+00f  1
  12  7.0139128e+01 0.00e+00 1.84e-11  -5.7 1.02e+00    -  1.00e+00 1.00e+00f  1
  13  7.0139120e+01 0.00e+00 2.57e-14  -8.6 6.71e-03    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 13

                                   (scaled)                 (unscaled)
Objective...............:   7.0139119670067416e+01    7.0139119670067416e+01
Dual infeasibility......:   2.5719026771974613e-14    2.5719026771974613e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.6389668491914908e-09    2.6389668491914908e-09
Overall NLP error.......:   2.6389668491914908e-09    2.6389668491914908e-09


Number of objective function evaluations             = 14
Number of objective gradient evaluations             = 14
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 14
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 14
Number of Lagrangian Hessian evaluations             = 13
Total CPU secs in IPOPT (w/o function evaluations)   =      0.005
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
Status: optimal
Optimal solution:
broccoli: 12.56 units
bread: 8.52 units
steak: 3.14 units
Popeyes: 0.24 units
water: 3.00 units
Total price: $70.14

As seen in the above model, the fat value per 3oz of steak is now 14 grams. In the model output, Popeyes is only selected 0.24 units, which is significantly different from the baseline model output. The amount of broccoli and bread also changed significantly, leading to a more balanced diet.

Turning Point 2: Fat per 3oz steak = 15.5#

model = pyo.ConcreteModel()
food_items = ['broccoli', 'bread', 'steak', 'Popeyes', 'water']

nutrition_data = {
    'Protein': [2.6, 3.5, 26, 76, 0],
    'Carbohydrates': [6, 19.5, 0, 36, 0],
    'Fat': [0.3, 0.84, 15.5, 64, 0],
    'Fiber': [2, 1.16, 0, 0, 0],
    'Volume': [2.6, 1.41, 3, 25, 10],
    'Price': [2, 1.5, 6, 25, 2.5]
}

bounds = {
    'Protein': 70,
    'Carbohydrates': 250,
    'Fat': 70,
    'Fiber': 35,
    'Volume': 90
}

model.food_quantities = pyo.Var(food_items, domain=pyo.NonNegativeReals)

model.obj = pyo.Objective(expr=sum(model.food_quantities[i] * nutrition_data['Price'][j]
                                   for j, i in enumerate(food_items)), sense=pyo.minimize)

model.protein_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Protein'][j]
                                                   for j, i in enumerate(food_items)) >= bounds['Protein'])
model.carbohydrates_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Carbohydrates'][j]
                                                         for j, i in enumerate(food_items)) >= bounds['Carbohydrates'])
model.fat_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fat'][j]
                                               for j, i in enumerate(food_items)) >= bounds['Fat'])
model.fiber_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fiber'][j]
                                                 for j, i in enumerate(food_items)) >= bounds['Fiber'])
model.volume_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Volume'][j]
                                                  for j, i in enumerate(food_items)) <= bounds['Volume'])
model.popeyes_constraint = pyo.Constraint(expr=model.food_quantities['water'] >= 3)

solver = pyo.SolverFactory('ipopt')
results = solver.solve(model, tee=True)

# Print the results
print("Status:", results.solver.termination_condition)
print("Optimal solution:")
for i in food_items:
    print(f"{i}: {pyo.value(model.food_quantities[i]):.2f} units")

# Print the total price
total_price = sum(pyo.value(model.food_quantities[i]) * nutrition_data['Price'][j]
                  for j, i in enumerate(food_items))
print(f"Total price: ${total_price:.2f}")

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:       19
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        5
                     variables with only lower bounds:        5
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        6
        inequality constraints with only lower bounds:        5
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        1

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  3.6999963e-01 2.49e+02 1.12e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  5.1089735e-01 2.49e+02 1.14e+00  -1.0 6.66e+01    -  1.79e-03 3.51e-03h  1
   2  1.0291764e+00 2.45e+02 1.58e+00  -1.0 5.03e+01    -  7.24e-03 1.30e-02h  1
   3  2.2285306e+00 2.41e+02 3.82e+00  -1.0 9.17e+01    -  1.11e-02 1.94e-02h  1
   4  1.0056887e+01 2.12e+02 1.10e+01  -1.0 9.50e+01    -  3.10e-02 1.19e-01h  1
   5  6.2071409e+01 2.13e+01 7.36e+00  -1.0 8.41e+01    -  2.46e-01 8.91e-01h  1
   6  6.2417799e+01 1.98e+01 2.85e+00  -1.0 6.52e+00    -  1.59e-01 6.11e-02h  1
   7  6.8410733e+01 0.00e+00 8.31e+00  -1.0 2.73e+01    -  2.06e-02 1.00e+00h  1
   8  6.8450525e+01 0.00e+00 1.02e-01  -1.0 3.69e-01    -  9.94e-01 1.00e+00f  1
   9  6.8338609e+01 0.00e+00 3.79e-02  -2.5 1.19e+00    -  8.91e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  6.8323287e+01 0.00e+00 2.11e-02  -2.5 1.57e+01    -  6.66e-01 5.68e-01f  1
  11  6.8309284e+01 0.00e+00 1.50e-09  -3.8 8.71e-01    -  1.00e+00 1.00e+00f  1
  12  6.8308476e+01 0.00e+00 1.84e-11  -5.7 6.77e-02    -  1.00e+00 1.00e+00f  1
  13  6.8308465e+01 0.00e+00 2.55e-14  -8.6 2.97e-03    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 13

                                   (scaled)                 (unscaled)
Objective...............:   6.8308464554317595e+01    6.8308464554317595e+01
Dual infeasibility......:   2.5491926439145658e-14    2.5491926439145658e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.8716956147740828e-09    2.8716956147740828e-09
Overall NLP error.......:   2.8716956147740828e-09    2.8716956147740828e-09


Number of objective function evaluations             = 14
Number of objective gradient evaluations             = 14
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 14
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 14
Number of Lagrangian Hessian evaluations             = 13
Total CPU secs in IPOPT (w/o function evaluations)   =      0.006
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
Status: optimal
Optimal solution:
broccoli: 12.25 units
bread: 9.05 units
steak: 3.79 units
Popeyes: 0.00 units
water: 3.00 units
Total price: $68.31

Now as the fat value of a 3oz steak increase to 15.5, the amount chosen has increased to 3.79 units. Popeyes is completely excluded from the final result as well, having 0 units selected.

Changes in Amount of Other Types of Food#

The following graph shows how each of the items selected change as we increase steak’s fat nutritional value:

# Initialize the range of fat values for steak
fat_values = np.linspace(5, 20, 50)

food_items = ['broccoli', 'bread', 'steak', 'Popeyes', 'water']
quantities = {item: [] for item in food_items}  # To store quantities of all food items

# Define a function to solve the model with a given fat value for steak
def solve_model(fat_value):
    model = pyo.ConcreteModel()
    nutrition_data = {
        'Protein': [2.6, 3.5, 26, 76, 0],
        'Carbohydrates': [6, 19.5, 0, 36, 0],
        'Fat': [0.3, 0.84, fat_value, 64, 0],  # Update fat value for steak
        'Fiber': [2, 1.16, 0, 0, 0],
        'Volume': [2.6, 1.41, 3, 25, 10],
        'Price': [2, 1.5, 6, 25, 2.5]
    }

    bounds = {
        'Protein': 70,
        'Carbohydrates': 250,
        'Fat': 70,
        'Fiber': 35,
        'Volume': 90
    }

    model.food_quantities = pyo.Var(food_items, domain=pyo.NonNegativeReals)

    model.obj = pyo.Objective(expr=sum(model.food_quantities[i] * nutrition_data['Price'][j]
                                       for j, i in enumerate(food_items)), sense=pyo.minimize)

    model.protein_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Protein'][j]
                                                       for j, i in enumerate(food_items)) >= bounds['Protein'])
    model.carbohydrates_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Carbohydrates'][j]
                                                             for j, i in enumerate(food_items)) >= bounds['Carbohydrates'])
    model.fat_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fat'][j]
                                                   for j, i in enumerate(food_items)) >= bounds['Fat'])
    model.fiber_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Fiber'][j]
                                                     for j, i in enumerate(food_items)) >= bounds['Fiber'])
    model.volume_constraint = pyo.Constraint(expr=sum(model.food_quantities[i] * nutrition_data['Volume'][j]
                                                      for j, i in enumerate(food_items)) <= bounds['Volume'])
    model.water_constraint = pyo.Constraint(expr=model.food_quantities['water'] == 3)

    solver = pyo.SolverFactory('ipopt')
    results = solver.solve(model)

    # Return the selected quantity of each food item
    return {item: pyo.value(model.food_quantities[item]) for item in food_items}

# Loop over the range of fat values and solve the model for each
for fat_value in fat_values:
    food_quantities = solve_model(fat_value)

    # Store the quantity for each food item at the current fat value
    for item in food_items:
        quantities[item].append(food_quantities[item])

# Plot the results
plt.figure(figsize=(10, 6))
for item in food_items:
    plt.plot(fat_values, quantities[item], label=item)

plt.title('Change in Quantities of Food Items as Fat Value of Steak Changes')
plt.xlabel('Fat Value of Steak')
plt.ylabel('Quantity Selected (units)')
plt.legend()
plt.grid(True)
plt.show()

../../_images/109dc8d5141d6e0e256f5fab0ac72f8a86453de5d0cd1ee352c467e801a967c8.png

In this graph, we can observe the changes in the unit quantity selected of each food type while looping through the fat value of steak from 5 to 20 with 50 even intermediate values. The first turning point occurs when steak’s fat value is increased to around 14, where the quantity of bread selected sharply decreases from 30 to around 8, while the quantity selected for broccoli sharply increases to 12.5, for steak from 0 to 3. Then the composition stays stable, while the quantity selected of steak slowly decreases as the fat value per unit increases. The second turning point occurs when the steak’s fat value is increased to around 15.5, where smaller changes compared the first turning point occurs. Now steak is selected around 4 units, and Popeyes is not selected at all. The quantity of steak selected continues to slowly decrease after the second turning point as the fat value continues to increase. The quantity of water always stay unchanged at 3, primarily due to the set up of the constraint limiting the quantity of water consumed to be 3.

References#

Pyomo-optimization modeling in python (Vol. 67, p. 277). Berlin: Springer.
Continuous Optimization: Linear Programming
https://plk-use1-prod.sites.rbictg.com/nutrition/PLK_Nutrition.pdf
https://www.verywellfit.com/nutrition-4157081

Optimization of Daily Diet Using Pyomo

Contents

Optimization of Daily Diet Using Pyomo#

Learning Objective#

Primary Source#

Install Packages#

Model Without Water#

Model Description and Formulation#

Degree of Freedom Analysis#

Model With Water#

Model Description and Formulation#

Degree of Freedom Analysis#

Steak Nutritional Value Modifications#

Protein Modifications#

Fat Modifications#

Turning Point 1: Fat per 3oz steak = 14#

Turning Point 2: Fat per 3oz steak = 15.5#

Changes in Amount of Other Types of Food#

References#