The Stigler Diet Linear Programming Problem in Python

The Stigler Diet problem is similar to the diet problem we solved previously. The difference is that we now have much more data to work with.

The objective is to determine how many dollar to be spend of each foodstuff at a minimum white satisfying the minimum nutrition requirement. The table below provides the list of foodstuff(commodity) along with their cost per unit and the nutrient content.

For example, Macaroni cost 14.1 cents and contains 418g of Protein, 0.7g of Calcium etc.

We are also provided with the minimum daily requirements for the 8 nutrients. This is as given in the table below:

Since the nutrients have all been normalized by the price, the goal is now to minimize the sum of the foods, therefore,

• the variables will be the food items (Commodity column of the table)
• cost function will be the sum of all the foods, each food with coefficient of 1

Note: the code for the above two tables can be found in here. You can also watch the video and see the link in the video description.

Create the variables and define the cost function

The code below creates the variables and stores them in a list called food. Then we create the objective function using the variables and assign the coefficient of each variable the value 1

food = [[]] * len(data) # list to hold the variables

objective = solver.Objective() # define the objective
for in range(0 to len(data)):
    food[i] = solver.NumVar(0, solver.infinity(), data[i][0])
    objective.SetCoefficient(food[i], 1)
objective.SetMinimization()

# food[i] is the normalized amount to be spent on foodstuff i

Define the Constraints

The constraints require that the total amount of nutrients provided by all the foods be at least the minimum daily nutrient requirement for each nutrient. The model the constraints, we consider the following:

the amount of nutrient i provided by food j per dollar is data[j][i+3] (note that it’s j+3 because the nutrient data begins in the 4th column of the list)

since the amount of to be spent on food j is food[j], therefore, the amount of nutrient provided by food j is data[j][i+3] * food[j]

the total nutrient provided by all the foods is then given by the sum of all the nutrients

we can then write the constraint as

The following code creates the constraints:

# Setup the constraints

constraints = [0] * len(nutrients) 
for i in range(0, len(nutrients)):
    constraints[i] = solver.Constraint(nutrients[i][1], solver.infinity())
    for j in range(0, len(data)):
        constraints[i].SetCoefficient(food[j], data[j][i+3])
        
# the outer loop i, moves across the columns, creating the constraints
# the inner loop j, moves down the rows, assigning one coefficient to each constraint

The following code calls the solver and displays the values of the solution:

# Call the solver and print the results
status = solver.Solve()

if status == solver.OPTIMAL:
    price = 0
    num_nutrients = len(data[i]) - 3
    nutrients = [0] * (len(data[i]) - 3)
    
    # loop through the rows of data
    for i in range(0, len(data)):
        price += food[i].solution_value()
        
        # loop through the nutrients
        for nutrient in range(0, num_nutrients):
            nutrients[nutrient] += data[i][nutrient+3] * food[i].solution_value()
        
        # print the cost of food if the value is not 0
        if food[i].solution_value() > 0:
            print('%s = %f' %(data[i][0], food[i].solution_value()))
            
    print('Optimal yearly price: $%.2f' %(365 * price))
else:
    if status == solver.FEASIBLE:
        print('A potentially suboptimal solution was found.')
    else:
        print('The solver could not find a solution to the problem')