stigler_diet2023.py

gistfile1.txt

"""The Stigler diet problem.

A description of the problem can be found here:
https://en.wikipedia.org/wiki/Stigler_diet.
Forked from ortools tutorial: https://developers.google.com/optimization/lp/stigler_diet
"""
from ortools.linear_solver import pywraplp
import pandas as pd

def main():
    """Entry point of the program."""
    # Instantiate the data problem.
    # Nutrient minimums.
    # https://www.nutrition.org.uk/media/1z2ekndj/nutrition-requirements-update.pdf
    # http://www.mydailyintake.net/daily-intake-levels/
    nutrients = [
        ["Calories (cal)", 3000],
        ["Protein (g)", 70],
        # ["Calcium (g)", 0.8],
        # ["Iron (mg)", 12],
        ["Vitamin A (IU)", 3000],
        # ["Vitamin B1 (mg)", 1.8],
        # ["Vitamin B2 (mg)", 2.7],
        # ["Niacin (mg)", 18],
        ["Vitamin C (mg)", 90],
        ["Fat (g)", 78],
        ["Carbohydrate (g)", 275],
        ["Fibre (g)", 28],
        ["Salt (g)", 2.3],       
    ]

    sheet_url = "https://docs.google.com/spreadsheets/d/1xUptnXwkS5fUioOz0YTzHApADOWPRyzF6LtmnLgt7ag/export?gid=0&format=csv"
    data = pd.read_csv(sheet_url).fillna(0.)
    # Instantiate a Glop solver and naming it.
    solver = pywraplp.Solver.CreateSolver("GLOP")
    if not solver:
        return

    # Declare an array to hold our variables.
    foods = [solver.NumVar(0.0, solver.infinity(), item) for item in data.Product]

    print("Number of variables =", solver.NumVariables())

    # Create the constraints, one per nutrient.
    constraints = []
    for i, nutrient in enumerate(nutrients):
        constraints.append(solver.Constraint(nutrient[1], solver.infinity()))
        for j, (_, item) in enumerate(data.iterrows()):
            constraints[i].SetCoefficient(foods[j], item[nutrient[0]])

    print("Number of constraints =", solver.NumConstraints())

    # Objective function: Minimize the sum of (price-normalized) foods.
    objective = solver.Objective()
    for food in foods:
        objective.SetCoefficient(food, 1)
    objective.SetMinimization()

    status = solver.Solve()

    # Check that the problem has an optimal solution.
    if status != solver.OPTIMAL:
        print("The problem does not have an optimal solution!")
        if status == solver.FEASIBLE:
            print("A potentially suboptimal solution was found.")
        else:
            print("The solver could not solve the problem.")
            exit(1)

    # Display the amounts (in dollars) to purchase of each food.
    nutrients_result = [0] * len(nutrients)
    print("\nAnnual Foods:")
    for i, food in enumerate(foods):
        if food.solution_value() > 0.0:
            print("{}: €{}".format(data.iloc[i][0], 365.0 * food.solution_value()))
            for j, (name, _) in enumerate(nutrients):
                nutrients_result[j] += data.iloc[i][name] * food.solution_value()
    print("\nOptimal annual price: €{:.4f}".format(365.0 * objective.Value()))

    print("\nNutrients per day:")
    for i, nutrient in enumerate(nutrients):
        print(
            "{}: {:.2f} (min {})".format(nutrient[0], nutrients_result[i], nutrient[1])
        )

if __name__ == "__main__":
    main()