Solving the Stigler Diet problem with Gurobi, CPlex, and Glop

First some of the motivation behind this this small project:

• A while ago I did some optimization work and this project is to ‘awaken some knowledge’ so I can eventually work on some more complicated optimization projects
• When I did do some optimization work, I had a lot of difficulties with different solvers on the different machines I had access to. Some machines were on an academic domain, others weren’t, and others were behind weird firewalls. Hence, licensing was difficult. At the time, I used a mix of some horrible code that made me want to cry. I always wanted to make a wrapper class for the different solvers and thought that now is a great time to do so
• When I did my optimization work I wanted to try out Glop. So this project allowed me to play around a bit with it.
• Also, I like the diet problem. I plan to do some experiments on myself with it.

Stigler Diet Problem…

The Stigler Diet problem is perhaps one of the more famous Linear programming problems out there. The problem was made up in the early 1940s by George Stigler where the goal is find the minimum costing set of foods that meet the dietary requirements. Wikipedia describes it as:

For a moderately active man weighing 154 pounds, how much of each of 77 foods should be eaten on a daily basis so that the man’s intake of nine nutrients will be at least equal to the recommended dietary allowances (RDAs) suggested by the National Research Council in 1943, with the cost of the diet being minimal?

Mathematically the problem can be described as:

Solving the Diet Problem

With a larger set of foods (that is decision variable), the diet problem can take ages and ages to solve by hand on paper. Either one needs to write a program to solve the problem or use a pre-existing tool. Since the techniques to solve optimization problems are fairly well known (for example the simplex algorithm, branch and bound, and the cutting plane method), there are a fair bit of open source and commercial tools available including: Gurobi, CPlex, Glop, SciPy.Optimize, and PuLP. With this project Gurobi, CPlex and Glop were used.

The main pattern when using any of these solvers is:

• Build the Model: create decision variables, specify constraints, and the objective function
• Solve the Model: let the solver do what it was designed for
• Report Results: output details if an optimal solution was found including decision variable and objective function values

Gurobi Solution

Gurobi is a commercial product originally created by Zonghao Gu, Edward Rothberg, and Robert Bixby and is currently sold by Gurobi Optimization LLC. They claim that they are the fastest solver out there.

Note that Gurobi needs a licence to run. It can either be purchased or used for free in an academic environment. Also, a temporary 2 month licence can be obtained for free which is nice.

Building the Model

The following is required to build the model, that is create the model, decision variables, constraints, and objective function. Note that all code snippets are in Python. 

model = gp.Model ('DietProblem')

# Create the decision variables... one by one
amountsToBuy = {}
for food in foodNames:
    amountsToBuy [food] = model.addVar (name=food)

# make and assign the objective function
expression = gp.LinExpr ()
for food in self._foodNames:
    expression.addTerms (foodCosts [food]['unitCost'], amountsToBuy [food])
model.setObjective (expression, GRB.MINIMIZE)

# Add constraints...
for nutrient in nutrientNames:
    expression = gp.LinExpr ()
    for food in foodNames:
        expression.addTerms (nutrientData[food][nutrient], amountsToBuy[food])
        model.addRange (expression, \
                        dietaryReqs[nutrient]['min'], \
                        dietaryReqs[nutrient]['max'], \
                        nutrient) 

Solving the Model

If the model was successfully built, then a call to optimize will attempt to solve it, that is:

model.optimize ()

Reporting the Results

If the model was successfully solved, then the values of the decision variables can be outputted and used for various calculations:

modelVars = model.getVars ()

# Display the amounts to purchase of each food.
for var in modelVars:
    # note both decision vars and constraints are variables
    if var.varName in foodNames and var.x > 0.0001:
        print (f'{var.varName:25} --> {var.x:8.3f} units  $ {(var.x * foodCosts [var.varName]["unitCost"]):.3f}')

print ()
print (f'Objective function results: {model.objVal}')
print ()

# Nutrient amounts... 
for nutrient in self._nutrientNames:
    total = 0
    for var in modelVars:
        if var.varName in foodNames and var.x > 0.0001:
            total += var.x * nutrientData [var.varName][nutrient]

    print (f'   {nutrient:15}  {total:7.2f}')

CPlex Solution

CPlex is currently an IBM product and it was originally made by some of the same people who created Gurobi. The free/demo version is limited to 1000 variables/constraints which makes any complicated example impossible to run (such as the online CPlex examples supplied by IBM). However, this limitation is fine for the diet problem. 

Note that there is a bit of a history with CPlex. Once the original project reached a certain level of maturity, it was purchased by ILOG (which is a french company, and the name is a combination of “Intelligence” and “Logiciel”). Eventually IBM purchased a portion of ILOG and CPlex became IBM ILOG CPlex. Then in 2020, IBM finalized the purchase of ILOG and now in the code CPlex is referred to as DOCPlex where the DO stands for Decision Optimization.

Building the Model

The code is similar to the Gurobi code however there are some differences. The two main ones are: 1) IBM’s documentation isn’t as nice a Gurobi’s and one needs to spend a fair bit of time with google; and 2) it seems to make life a bit easier for the users, CPlex has a lot of extra helper functions (like with computing KPIs).

Model = Model (name='DietProblem')

# Decision variables, limited to a range 
amountToBuy = {}
for food in foodNames:
    amountToBuy [food] = Model.continuous_var (lb = self._foodCosts[food]['minQuantity'], \
                                               ub = self._foodCosts[food]['maxQuantity'], \
                                               name = (food))

# adding constraints 
for nutrient in nutrientNames:
    categoryAmount = Model.sum (amountToBuy [food] * nutrientData [food][nutrient] for food in foodNames )
    Model.add_range (dietaryReqs[nutrient]['min'], \ 
                     categoryAmount, \ 
                     dietaryReqs[nutrient]['max'] )

# adding objective function 
objectiveFunction = Model.sum (amountToBuy [food] * foodCosts [food]['unitCost'] for food in foodNames )
Model.minimize (objectiveFunction)

Solving the Model

One funny thing is that the model.solve function does not return a value; rather state variables need to be inspected afterwards.

Model.solve ()
solveStatus = Model.solve_details.status_code

Reporting the Results

Again, CPlex has some extra functions which can make it easier in reporting KPIs. Also, it has a few functions which output the state of the model including:

Model.print_information ()
Model.print_solution ()

The reporting can be performed in a manner similar to Gurobi…

if solveStatus == True:

    totalCost = 0
    print (f"{'Food':22}   Amount  Cost")
    for food in foodNames:
        decVar = Model.get_var_by_name (food)
        cost = foodCosts [food]['unitCost'] * decVar.solution_value
        totalCost += cost
        if decVar.solution_value > 0.001:
            print (f'{food:22}  {decVar.solution_value:6.3f}   {cost:.2f}')

    print ()
    print (f'Total Cost: {totalCost:.2f}')
    print ()

    print ('Daily Requirement Categories')
    for nutrient in nutrientNames:
        categoryAmount = 0
        for food in foodNames:
            decVar = Model.get_var_by_name (food)
            categoryAmount += ( decVar.solution_value * nutrientData [food][nutrient] )

        print (f'{nutrient:15}  {categoryAmount:7.2f}')

Glop Solution

Glop is a linear optimization/programming solver made by Google. It is not as polished and mature as either Gurobi or CPlex, however it is still good enough to solve a most of optimization problems. The code is fairly easy to use and much more familiar that both CPlex or Gurobi.

Building the Model

The model building code is mostly similar to Gurobi and CPlex, however there are a few differences.

model = pywraplp.Solver ('DietProblem', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING) 

        # make an empty objective function to be filled later
        objective = self._model.Objective()

        # make decision variables and complete the objective function
        self._amountToBuy = {}
        for food in self._foodNames:
            self._amountToBuy [food] = self._model.NumVar (self._foodCosts[food]['minQuantity'], \
                                                           self._foodCosts[food]['maxQuantity'], \
                                                           food)

            # sets the coefficient for the objective function
            objective.SetCoefficient (self._amountToBuy[food], self._foodCosts[food]['unitCost'])

        objective.SetMinimization()

        # make the constraints
        constraints = {}

        for nutrient in self._nutrientNames:
            # Make a constraint for each nutrient... min is defined, max is infinity
            constraints [nutrient] = self._model.Constraint (self._dietaryReqs[nutrient]['min'], \
                                                             self._dietaryReqs[nutrient]['max'], \
                                                             nutrient)

            # for each food item
            for food in self._foodNames:
                constraints [nutrient].SetCoefficient (self._amountToBuy[food], self._nutrientData[food][nutrient])

Solving the Model

This seems to be a mix between Gurobi and CPlex…

solveStatus = model.Solve ()

Reporting the Results

And reporting the results is slightly more work than the previous solutions but still not too bad:

if solveStatus == model.OPTIMAL:

    cost = 0
    for food in foodNames:
        cost += amountToBuy[food].solution_value() * foodCosts[food]['unitCost'];

        if amountToBuy[food].solution_value() > 0.001:
            print (f'{food:25} --> {amountToBuy[food].solution_value():7.3f} units')
            print (f'          --> {amountToBuy[food].solution_value()*foodCosts[food]["unitCost"]:.2f} dollars')

    totalNutrients = {}
    for nutrient in nutrientNames:
        totalNutrients [nutrient] = 0
        for food in foodNames:
            totalNutrients [nutrient] += amountToBuy[food].solution_value() * nutrientData[food][nutrient]

    print ()
    print (f'Optimal daily cost:  {cost:7.2f}')
    print (f'Optimal annual cost: {(365 * cost):7.2f}')
    print ()

    for nutrient in totalNutrients:
        print (f'{nutrient:15} --> {totalNutrients[nutrient]:7.2f}')

Solution Wrapper

Back when I did some pretty basic optimization research, I really wanted to make some class wrappers to make interacting with the different solvers a bit easier. That way I could go between different machines where different solvers were installed (due to licensing issues because of the different environments the machines belonged to). One of the higher ups was really against this idea since ‘it was a huge waste of time’ and ‘won’t produce any papers’. Since I am doing this project for myself, I am going to do it the way I want to so I’m adding some class wrappers.

The main solver class is as follows. Note that parameters and data are placeholders for the actual lists and dictionary values.

class DietProblemSolver ():

    def __init__ (self, parameters, data):

        self._parameters  = parameters
        self._data        = data
        self._model       = None
        self._solveStatus = -1

    def GetSolverName (self): 
        pass

    def BuildModel (self):
        pass

    def SolveModel (self):
        pass

    def PrintResults (self):
        pass

Next each different solver is derived from DietProblemSolver. For example, the CPlex version is as follows:

from DietProblemSolver import DietProblemSolver

from docplex.mp.model import Model
from docplex.util.environment import get_environment
import cplex


class DietProblemSolver_CPlex (DietProblemSolver):

  def GetSolverName (self): 
    return 'CPlex ' + cplex.Cplex().get_version()

  def BuildModel (self):

    # Make the model
    self._Model = Model (name='DietProblem')

    # Decision variables, limited to a range 
    self._amountToBuy = {}
    for food in self._foodNames:
      self._amountToBuy [food] = self._Model.continuous_var (lb = self._foodCosts[food]['minQuantity'], \
                                                             ub = self._foodCosts[food]['maxQuantity'], \                                                                                              
                                                             name = (food))

    # adding constraints 
    for nutrient in self._nutrientNames:
      categoryAmount = self._Model.sum ( self._amountToBuy [food] * self._nutrientData [food][nutrient] for food in self._foodNames )
      self._Model.add_range ( self._dietaryReqs[nutrient]['min'], \
                              categoryAmount, \ 
                              self._dietaryReqs[nutrient]['max'] )

    # adding objective function 
    objectiveFunction = self._Model.sum ( self._amountToBuy [food] * self._foodCosts [food]['unitCost'] for food in self._foodNames )
    self._Model.minimize (objectiveFunction)

  def SolveModel (self):

    if self._Model != None:
      self._Model.solve ()
      self._solveStatus = self._Model.solve_details.status_code

  def PrintResults (self):

    if self._solveStatus == True:

      totalCost = 0
      print (f"{'Food':22}   Amount  Cost")
      for food in self._foodNames:
        decVar = self._Model.get_var_by_name (food)
        cost = self._foodCosts [food]['unitCost'] * decVar.solution_value
        totalCost += cost
        if decVar.solution_value > 0.001:
          print (f'{food:22}  {decVar.solution_value:6.3f}   {cost:.2f}')

        print (f'Total Cost: {totalCost:.2f}')

        print ('Daily Requirement Categories')
        for nutrient in self._nutrientNames:
          categoryAmount = 0
          for food in self._foodNames:
            decVar = self._Model.get_var_by_name (food)
            categoryAmount += ( decVar.solution_value * self._nutrientData [food][nutrient] )

            print (f'{nutrient:15}  {categoryAmount:7.2f}')

    else:  
      print ('no solution')

Finally, the main code would build and use the required solver in the following manner…

DietProblemSolver = None
if solverName == 'Glop':
    DietProblemSolver = DietProblemSolver_Glop (params...)
elif solverName == 'Gurobi':
    DietProblemSolver = DietProblemSolver_Gurobi (params...) 
elif solverName == 'Cplex':
    DietProblemSolver = DietProblemSolver_CPlex (params...)

if DietProblemSolver:
    DietProblemSolver.BuildModel ()
      .
      .
      .

Code

All the code is available here. Again it is all written in python.

Final Note

In the code some basic datasets were used (and are provided in the above link). All of the datasets are fairly simple and not too realistic; either the cases are too simple or out of date (Stigler’s original 1939 dataset). This project will be revisited in the future with a much larger dataset and perhaps even the ability to obtain values from a database.