Solving MILP Problems in Python - Part 2¶
The Diet Example¶
So far we have written models with embedded data. However, when building an optimization model, it is typical to separate the optimization model itself from the data used to create an instance of the model. These two model ingredients are often stored in completely different files.
There are alternate approaches to providing data to the optimization model: they can be embedded in the source file, read from an SQL database (using the Python sqlite3 package), or read them from an Excel spreadsheet (using the Python xlrd package) and more.
Diet Problem¶
Bob wants to plan a nutritious diet, but he is on a limited budget, so he wants to spend as little money as possible. His nutritional requirements are as follows:
----------------
2000 Kcal
55 g protein
800 mg calcium
----------------Bob is considering the following foods with corresponding nutritional values
Serving Size Price per serving Energy (Kcal) Protein (g) Calcium (mg)
----------- -------------- ------------------- --------------- ------------- --------------
Oatmeal 28 g 0.3 110 4 2
Chicken 100 g 2.4 205 32 12
Eggs 2 large 1.3 160 13 54
Milk 237 cc 0.9 160 8 285
Apple Pie 170 g 2 420 4 22
Pork 260 g 1.9 260 14 80With the help of Gurobi Python, find the amount of servings of each type of food in the diet.
We specify the model independently from the data. We could put the model in a file, eg, dietmodel.py and the data in another file, eg, diet1.py. Make
sure you understand the model and read about the Gurobi global function
quicksum.
%pip install -i https://pypi.gurobi.com gurobipy
import gurobipy as gp
def solve(categories, minNutrition, maxNutrition, foods, cost, nutritionValues):
# Model
m = gp.Model("diet")
# Create decision variables for the nutrition information,
# which we limit via bounds
nutrition = {}
for c in categories:
nutrition[c] = m.addVar(lb=minNutrition[c], ub=maxNutrition[c], name=c)
# Create decision variables for the foods to buy
buy = {}
for f in foods:
buy[f] = m.addVar(obj=cost[f], name=f)
# The objective is to minimize the costs
m.modelSense=gp.GRB.MINIMIZE
# Nutrition constraints
for c in categories:
m.addConstr(gp.quicksum(nutritionValues[f, c] * buy[f] for f in foods) == nutrition[c],c)
def printSolution():
if m.status == gp.GRB.status.OPTIMAL:
print('\nCost: %g' % m.ObjVal)
print('\nBuy:')
for f in foods:
if buy[f].X > 0.0001:
print('%s %g' % (f, buy[f].X))
print('\nNutrition:')
for c in categories:
print('%s %g' % (c, nutrition[c].X))
else:
print('No solution')
# Solve
m.update()
m.display()
m.write("diet.lp")
m.optimize()
printSolution()
To arrange the data in Python data structures we use another global function from Gurobi:
multidict. Here is an example of what it
does:
keys, dict1, dict2 = gp.multidict( {
'key1': [1, 2],
'key2': [1, 3],
'key3': [1, 4] } )
print(keys, dict1, dict2)
categories, minNutrition, maxNutrition = gp.multidict({
'Calories': [1800, 2200],
'Protein': [91, gp.GRB.INFINITY],
'Calcium': [0, 1779] })
foods, cost = gp.multidict({
'Oatmeal': 0.30,
'Chicken': 2.40,
'Eggs': 1.30,
'Milk': 0.90,
'Apple Pie': 2.00,
'Pork': 1.90});
# Nutrition values for the foods
nutritionValues = {
('Oatmeal', 'Calories' ): 110,
('Oatmeal', 'Protein' ): 4,
('Oatmeal', 'Calcium' ): 2,
('Chicken', 'Calories' ): 205,
('Chicken', 'Protein' ): 32,
('Chicken', 'Calcium' ): 12,
('Eggs', 'Calories' ): 160,
('Eggs', 'Protein' ): 13,
('Eggs', 'Calcium' ): 54,
('Milk', 'Calories' ): 160,
('Milk', 'Protein' ): 8,
('Milk', 'Calcium' ): 285,
('Apple Pie', 'Calories' ): 420,
('Apple Pie', 'Protein' ): 4,
('Apple Pie', 'Calcium' ): 22,
('Pork', 'Calories' ): 260,
('Pork', 'Protein' ): 14,
('Pork', 'Calcium' ): 80 };
If we put the code in the two files, then we could solve from diet1.py as follows:
import dietmodel
dietmodel.solve(categories, minNutrition, maxNutrition, foods, cost, nutritionValues)
Here it suffices calling:
solve(categories, minNutrition, maxNutrition, foods, cost, nutritionValues)
Your Task¶
A pill salesman offers Bob Calories, Protein, and Calcium pills to fulfill his nutritional needs. He needs to estimate the prices of units of serving, that is, the cost of 1 kcal, the cost of 1 g of protein, the cost of 1 mg of calcium. He wants to make as much money as possible, given Bob's constraints. He knows that Bob wants 2200 kcal, 55 g protein, and 1779 mg calcium. How can we help him in guaranteeing that he does not make a bad deal?
Solution:
The dual seeks to maximize the profit of the salesman. Let $y_i \geq 0$, $i\in N$ be the prices of the pills.
$$\begin{aligned} \text{min}\quad &\sum_{j \in F}c_jx_j\\ \sum_{j\in F} a_{ij}x_j&\geq N_{min,i}, \qquad\forall i \in N\\ %\sum_{j\in F}a_{ij}x_{j}&\leq N_{max,i}, \qquad\forall i \in N\\ x_j&\geq 0, \qquad \forall j\in F\\ %x_j&\leq F_{max,j}, \qquad \forall j \in F\end{aligned}$$$$\begin{aligned} \text{max}\quad &\sum_{i \in F}N_{min,i}y_i\\ \sum_{i\in N} a_{ji}y_i&\leq c_{j}, \qquad\forall j \in F\\ y_i&\geq 0, \qquad \forall i\in N\\\end{aligned}$$However the values of the dual variables can be determined by the last
tableau of the solution to the primal problem by printing the Pi
attribute of the constraints.
Particular Cases¶
The two following LP problems lead to two particular cases when solved by the simplex algorithm. Identify these cases and characterize them, that is, give indication of which conditions generate them in general. Then, implement the models in Gurobi Python and observe the behaviour.
$$\begin{array}{rllllllllll} \mbox{maximize} & 2x_1& +& x_2&\\ \mbox{subject to} &&& x_2& \leq &5&\\ & - x_1& +& x_2 &\leq &1& \\ & x_1, &&x_2&\geq &0& \\ \end{array}$$$$\begin{array}{rlllllllllll} \mbox{maximize}& x_1 &+& x_2&\\ \mbox{subject to} & 5 x_1& +& 10 x_2& \leq& 60&\\ & 4 x_1 &+ &4 x_2 &\leq &40& \\ & x_1, &&x_2&\geq &0& \\ \end{array}$$Pathological Cases¶
This exercise asks you to check the behavior of the solvers on the two pathological cases:
$$\begin{array}{rlllllllllll} \mbox{maximize}&&& 4x_2&\\ \mbox{subject to}&&& 2x_2 &\geq& 0\\ & -3x_1& +&4 x_2& \geq &1\\ & x_1, &&x_2&\geq &0& \\ \end{array}$$$$\begin{array}{rlllllllllll} \mbox{maximize}\ \ \ 10x_1 -57 x_2 -9x_3-24x_4&\\ \mbox{subject to}\ \ -0.5x_1+5.5x_2+2.5x_3-9x_4 &\leq 0\\ -0.5 x_1 + 1.5 x_2 +0.5 x_3 -x_4& \leq 0\\ x_1 &\leq 1&\\ x_1, x_2, x_3, x_4&\geq 0& \\ \end{array}$$What happens with the solver? Can you detect which pathological cases are from the output of the solver? How?
Shortest Path¶
Model the shortest path problem as an LP problem. Write the model in Python using the skeleton below and the data available from http://www.imada.sdu.dk/~marco/DM559/Files/SP/20points.txt In this data the source is node 1 and the target is node 20.
Model the problem in LP and solve it with Gurobi Python. Check the correctness of your solution with the help of the visualization in the template below.
It may be worth looking at the examples netflow and tsp from the
Gurobi documentation to get inspiration for the model. For the
implementation it may be helpful using the Gurobi tuplelist class
tuplelist.
import sys
import math
from gurobipy import *
import matplotlib.pyplot as plt
if len(sys.argv) < 2:
print('Usage: sp.py file')
exit(1)
f = open(sys.argv[1])
lines=f.readlines();
f.close();
N=0
V=set()
points={}
for l in lines:
if l[0]!='#':
elem = l.split("\t")
if len(elem)==1:
N=int(elem[0])
elif len(elem)==3:
V.add(elem[0])
points[elem[0]]=(float(elem[1]),float(elem[2]))
source = '1'
target = str(len(points))
# calculate Euclidean distance and round-towards-zero (truncate)
def distance(points, i, j):
dx = points[i][0] - points[j][0]
dy = points[i][1] - points[j][1]
return math.floor(math.sqrt(dx*dx + dy*dy))
m=Model()
##############################################################################
# Begin: Change here
##############################################################################
# Create variables
##############################################################################
m.update()
##############################################################################
# Post the constraints
##############################################################################
# End: Change here
##############################################################################
# Optimize model
m.write("sp.lp")
m.display()
m.optimize()
solution = m.getAttr('x', vars)
selected = [(i,j) for i,j in arcs if solution[i,j] > 0.5]
print(filter(lambda x: x[1]>0, solution.iteritems()))
print( map(lambda i, j: distance(points, i, j), arcs))
print('')
print('Optimal path: %s' % str(selected))
print('Optimal cost: %g' % m.objVal)
print('')
def plot_path(points, path, style='bo-'):
"Plot lines to connect a series of points."
plt.plot(map(lambda x: x[1][0], points.iteritems()), map(lambda x: x[1][1], points.iteritems()), 'bo')
target = str(len(points))
plt.plot([points['1'][0],points[target][0]],[points['1'][1],points[target][1]], 'rs')
plt.plot(map(lambda x: (points[x[0]][0],points[x[1]][0]), path),
map(lambda x: (points[x[0]][1],points[x[1]][1]), path), 'b-')
plt.axis('scaled'); plt.axis('off')
plt.show()
print(points, selected)
plot_path(points, selected)