Linear Programming in Python
Introduction
This post contains additional materials to my newly released video about linear programming, namely a number of practical examples of how it can be used to solve a variety of problems using Python and its pulp package.
Practical examples
Farmer’s problem
With the planting season steadily approaching, your farmer friend presents you with the following problem.
You have tons of potato seeds and tons of carrot seeds. To grow the crops efficiently, you also have tons of fertilizer, which has to be used when planting in a ratio (i.e. kilogram of potatoes or carrots requires kilogram of fertilizer). The profit is for potato seeds and for carrot seeds.
How much potatoes and carrots should you plant to maximize your profit this season?
Source code
from pulp import *
# creating the model
model = LpProblem(sense=LpMaximize)
# variables - how much do we plant?
x_p = LpVariable(name="potatoes", lowBound=0)
x_c = LpVariable(name="carrots", lowBound=0)
# inequalities - don't use more than we have
model += x_p <= 3000 # potatoes
model += x_c <= 4000 # carrots
model += x_p + x_c <= 5000 # fertilizer
# objective function - maximize the profit
model += x_p * 1.2 + x_c * 1.7
# solve (ignoring debug messages)
status = model.solve(PULP_CBC_CMD(msg=False))
print("potatoes:", x_p.value())
print("carrots:", x_c.value())
print("profit:", model.objective.value())
Output
potatoes: 1000.0
carrots: 4000.0
profit: 8000.0
Knapsack problem
Given items, each with weight and price , our task is to maximize the price of the items we take into our backpack without exceeding its carry weight .
Source code
from pulp import *
# three datasets (small/medium/large)
# each contains weights / prices / knapsack carry weight
data = [
[
"Small",
[12, 7, 11, 8, 9],
[24, 13, 23, 15, 16],
26,
],
[
"Medium",
[23, 31, 29, 44, 53, 38, 63, 85, 89, 82],
[92, 57, 49, 68, 60, 43, 67, 84, 87, 72],
165,
],
[
"Large",
[382745, 799601, 909247, 729069, 467902, 44328, 34610, 698150, 823460, 903959, 853665, 551830, 610856, 670702, 488960, 951111, 323046, 446298, 931161, 31385, 496951, 264724, 224916, 169684],
[825594, 1677009, 1676628, 1523970, 943972, 97426, 69666, 1296457, 1679693, 1902996, 1844992, 1049289, 1252836, 1319836, 953277, 2067538, 675367, 853655, 1826027, 65731, 901489, 577243, 466257, 369261],
6404180,
],
]
for name, weights, prices, carry_weight in data:
n = len(weights)
model = LpProblem(sense=LpMaximize)
# variables - binary based on if we take the item or not
variables = [LpVariable(name=f"x_{i}", cat=LpBinary) for i in range(n)]
# inequalities - don't exceed the carry weight
model += lpDot(weights, variables) <= carry_weight
# objective function - price of the items we take
model += lpDot(prices, variables)
status = model.solve(PULP_CBC_CMD(msg=False))
print(name + "\n" + "-" * len(name))
print("price:", model.objective.value())
print("take:", *[variables[i].value() for i in range(n)])
print()
Output
Small
-----
price: 51.0
take: 0.0 1.0 1.0 1.0 0.0
Medium
------
price: 309.0
take: 1.0 1.0 1.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0
Large
-----
price: 13549094.0
take: 1.0 1.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0
Graph coloring
We want to find a minimal such that a graph is vertex -colorable.
Source code
from pulp import *
# three graphs, given as an edge list
# we're assuming that vertices are index 0 to n-1
graphs = [
[
"Complete graph on 5 vertices",
[0, 1, 2, 3, 4],
[(0, 1), (2, 1), (1, 3), (2, 3), (4, 3), (2, 4), (0, 4), (1, 4), (2, 0), (0, 3)],
],
[
"Path of 6 vertices",
[0, 1, 2, 3, 4, 5],
[(0, 1), (2, 1), (2, 3), (3, 4), (4, 5)],
],
[
"Petersen graph",
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
[(0, 1), (0, 4), (0, 5), (1, 2), (1, 6), (2, 3), (2, 7), (3, 4), (3, 8), (4, 9), (5, 7), (5, 8), (6, 8), (6, 9), (7, 9)],
]
]
for name, vertices, edges in graphs:
n = len(vertices)
# note: this is a minimization problem, so we use LpMinimize!
model = LpProblem(sense=LpMinimize)
# n^2 variables, where x[i][k] signifies if the i-th vertex is colored with the k-th color
variables = [
[LpVariable(name=f"x_{i}_{j}", cat=LpBinary) for j in range(n)]
for i in range(n)
]
# one more variable - the chromatic number itself
chromatic_number = LpVariable(name="chromatic number", cat='Integer')
# inequalities
# each vertex has exactly one color
for u in range(n):
model += lpSum(variables[u]) == 1
# adjacent edge colors are different
for u, v in edges:
for k in range(n):
model += variables[u][k] + variables[v][k] <= 1
# we also restrict the chromatic number to be the number of the highest used color
for u in range(n):
for k in range(n):
model += chromatic_number >= (k + 1) * variables[u][k]
# objective function - minimize the chromatic number
model += chromatic_number
status = model.solve(PULP_CBC_CMD(msg=False))
print(name + "\n" + "-" * len(name))
print("k =", int(chromatic_number.value()))
for u in range(n):
for k in range(n):
if variables[u][k].value() != 0:
print(f"v_{u} has color {k + 1}")
print()
Output
Complete graph on 5 vertices
----------------------------
k = 5
v_0 has color 5
v_1 has color 2
v_2 has color 1
v_3 has color 3
v_4 has color 4
Path of 6 vertices
------------------
k = 2
v_0 has color 2
v_1 has color 1
v_2 has color 2
v_3 has color 1
v_4 has color 2
v_5 has color 1
Petersen graph
--------------
k = 3
v_0 has color 1
v_1 has color 2
v_2 has color 3
v_3 has color 2
v_4 has color 3
v_5 has color 3
v_6 has color 3
v_7 has color 2
v_8 has color 1
v_9 has color 1
Traveling salesman problem
Given a weighted oriented graph , we want to find the longest Hamiltonian cycle.
Source code
from pulp import *
from itertools import *
# graphs as an adjacency list - how far is it from vertex i to vertex j
# taken from https://people.sc.fsu.edu/~jburkardt/datasets/tsp/tsp.html
graphs = [
[
"FIVE",
[[0],
[3, 0],
[4, 4, 0],
[2, 6, 5, 0],
[7, 3, 8, 6, 0]],
],
[
"P01",
[[ 0],
[29, 0],
[82, 55, 0],
[46, 46, 68, 0],
[68, 42, 46, 82, 0],
[52, 43, 55, 15, 74, 0],
[72, 43, 23, 72, 23, 61, 0],
[42, 23, 43, 31, 52, 23, 42, 0],
[51, 23, 41, 62, 21, 55, 23, 33, 0],
[55, 31, 29, 42, 46, 31, 31, 15, 29, 0],
[29, 41, 79, 21, 82, 33, 77, 37, 62, 51, 0],
[74, 51, 21, 51, 58, 37, 37, 33, 46, 21, 65, 0],
[23, 11, 64, 51, 46, 51, 51, 33, 29, 41, 42, 61, 0],
[72, 52, 31, 43, 65, 29, 46, 31, 51, 23, 59, 11, 62, 0],
[46, 21, 51, 64, 23, 59, 33, 37, 11, 37, 61, 55, 23, 59, 0]],
],
[
"GR17",
[[ 0],
[633, 0],
[257, 390, 0],
[ 91, 661, 228, 0],
[412, 227, 169, 383, 0],
[150, 488, 112, 120, 267, 0],
[ 80, 572, 196, 77, 351, 63, 0],
[134, 530, 154, 105, 309, 34, 29, 0],
[259, 555, 372, 175, 338, 264, 232, 249, 0],
[505, 289, 262, 476, 196, 360, 444, 402, 495, 0],
[353, 282, 110, 324, 61, 208, 292, 250, 352, 154, 0],
[324, 638, 437, 240, 421, 329, 297, 314, 95, 578, 435, 0],
[ 70, 567, 191, 27, 346, 83, 47, 68, 189, 439, 287, 254, 0],
[211, 466, 74, 182, 243, 105, 150, 108, 326, 336, 184, 391, 145, 0],
[268, 420, 53, 239, 199, 123, 207, 165, 383, 240, 140, 448, 202, 57, 0],
[246, 745, 472, 237, 528, 364, 332, 349, 202, 685, 542, 157, 289, 426, 483, 0],
[121, 518, 142, 84, 297, 35, 29, 36, 236, 390, 238, 301, 55, 96, 153, 336, 0]],
]
]
for name, distances in graphs:
n = len(distances)
model = LpProblem(sense=LpMinimize)
# variables - binary x_{i, j} based on if the edge (i, j) is on the cycle
# they're symmetric so we only care about i > j
# variables x_{i, i} will always be zero, they're there for nicer code
variables = [
[LpVariable(name=f"x_{i}_{j}", cat=LpBinary) for j in range(i + 1)]
for i in range(n)
]
# inequalities
## each edge has one incoming and one outgoing vertex
for i in range(n):
model += lpSum([variables[max(i, j)][min(i, j)] for j in range(n)]) == 2
model += variables[i][i] == 0
## each subset of vertices has an outgoing edge (no smaller loops!)
for size in range(1, n - 1):
for subset in combinations(range(n), r = size):
model += lpSum([variables[max(i, j)][min(i, j)]
for i in subset
for j in range(n)
if j not in subset]) >= 1
# objective function - minimize the length of the cycle
model += lpSum([variables[i][j] * distances[i][j]
for i, j in product(range(n), repeat=2)
if i > j])
status = model.solve(PULP_CBC_CMD(msg=False))
print(name + "\n" + "-" * len(name))
print("minimal tour: ", int(model.objective.value()))
for i in range(n):
for j in range(n):
print(int(variables[max(i, j)][min(i, j)].value()), " ", end="")
print()
print()
Output
FIVE
----
minimal tour: 19
0 0 1 1 0
0 0 1 0 1
1 1 0 0 0
1 0 0 0 1
0 1 0 1 0
P01
---
minimal tour: 291
0 0 0 0 0 0 0 0 0 0 1 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
0 0 0 0 0 0 1 0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 1 0 1 0 0 0 0 0 0
0 0 0 1 0 0 0 1 0 0 0 0 0 0 0
0 0 1 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 1 0 0 0 0 0 1 0
1 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 1 0
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 1 0 0 0
0 1 0 0 0 0 0 0 1 0 0 0 0 0 0
GR17
----
minimal tour: 2085
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0
1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0
0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0
0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0
1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0
Bin packing
Given items with weights and an arbitrary number of bins with maximum carry weight , determine the lowest number of bins that can contain all the items without exceeding their carry weight.
Source code
from pulp import *
from itertools import *
bin_sets = [
("Set 1", 100, [70, 60, 50, 33, 33, 33, 11, 7, 3]),
("Set 2", 100, [99, 94, 79, 64, 50, 46, 43, 37, 32, 19, 18, 7, 6, 3]),
("Set 3", 100, [49, 41, 34, 33, 29, 26, 26, 22, 20, 19]),
("Set 4", 524, [442, 252, 252, 252, 252, 252, 252, 252, 127, 127, 127, 127, 127, 106, 106, 106, 106, 85, 84, 46, 37, 37, 12, 12, 12, 10, 10, 10, 10, 10, 10, 9, 9]),
]
for name, maximum_bin_load, weights in bin_sets:
n = len(weights)
model = LpProblem(sense=LpMinimize)
# variables - binary x[i][j] based on if the i-th item is in the j-th bin
variables = [
[LpVariable(name=f"x_{i}_{j}", cat=LpBinary) for j in range(n)]
for i in range(n)
]
# ... and also the number of bins
bin_count = LpVariable(name="bin count", cat=LpInteger)
# inequalities
## each bin doesn't contain more than its maximum load
for j in range(n):
model += lpSum([variables[i][j] * weights[i] for i in range(n)]) <= maximum_bin_load
## each item must be placed in exactly one bin
for i in range(n):
model += lpSum([variables[i][j] for j in range(n)]) == 1
## we also want the bin count to be the highest bin number
for i in range(n):
for j in range(n):
model += bin_count >= (j + 1) * variables[i][j]
# objective function - number of bins
model += bin_count
status = model.solve(PULP_CBC_CMD(msg=False))
print(name + "\n" + "-" * len(name))
print("bin count:", int(bin_count.value()))
for i in range(n):
for j in range(n):
if variables[i][j].value() != 0:
print(f"item {i} (weight {weights[i]}): bin {j}")
print()
Output
Set 1
-----
bin count: 4
item 0 (weight 70): bin 2
item 1 (weight 60): bin 3
item 2 (weight 50): bin 0
item 3 (weight 33): bin 0
item 4 (weight 33): bin 1
item 5 (weight 33): bin 1
item 6 (weight 11): bin 2
item 7 (weight 7): bin 2
item 8 (weight 3): bin 2
Set 2
-----
bin count: 7
item 0 (weight 99): bin 4
item 1 (weight 94): bin 5
item 2 (weight 79): bin 1
item 3 (weight 64): bin 0
item 4 (weight 50): bin 2
item 5 (weight 46): bin 2
item 6 (weight 43): bin 6
item 7 (weight 37): bin 3
item 8 (weight 32): bin 0
item 9 (weight 19): bin 3
item 10 (weight 18): bin 1
item 11 (weight 7): bin 3
item 12 (weight 6): bin 5
item 13 (weight 3): bin 2
Set 3
-----
bin count: 3
item 0 (weight 49): bin 1
item 1 (weight 41): bin 0
item 2 (weight 34): bin 2
item 3 (weight 33): bin 0
item 4 (weight 29): bin 1
item 5 (weight 26): bin 2
item 6 (weight 26): bin 0
item 7 (weight 22): bin 1
item 8 (weight 20): bin 2
item 9 (weight 19): bin 2
Set 4
-----
bin count: 7
item 0 (weight 442): bin 0
item 1 (weight 252): bin 3
item 2 (weight 252): bin 6
item 3 (weight 252): bin 4
item 4 (weight 252): bin 6
item 5 (weight 252): bin 3
item 6 (weight 252): bin 4
item 7 (weight 252): bin 1
item 8 (weight 127): bin 1
item 9 (weight 127): bin 5
item 10 (weight 127): bin 5
item 11 (weight 127): bin 1
item 12 (weight 127): bin 5
item 13 (weight 106): bin 2
item 14 (weight 106): bin 5
item 15 (weight 106): bin 2
item 16 (weight 106): bin 2
item 17 (weight 85): bin 2
item 18 (weight 84): bin 2
item 19 (weight 46): bin 0
item 20 (weight 37): bin 2
item 21 (weight 37): bin 5
item 22 (weight 12): bin 0
item 23 (weight 12): bin 0
item 24 (weight 12): bin 0
item 25 (weight 10): bin 6
item 26 (weight 10): bin 3
item 27 (weight 10): bin 6
item 28 (weight 10): bin 4
item 29 (weight 10): bin 4
item 30 (weight 10): bin 3
item 31 (weight 9): bin 1
item 32 (weight 9): bin 1
Partition problem
Given items with weights , split them into two parts such that the difference in their weights is minimized.
Source code
from pulp import *
from itertools import *
weight_sets = [
("Set 1", [19, 17, 13, 9, 6]),
("Set 2", [3, 4, 3, 1, 3, 2, 3, 2, 1, 5]),
("Set 3", [2, 10, 3, 8, 5, 7, 9, 5, 3, 2]),
("Set 4", [771, 121, 281, 854, 885, 734, 486, 1003, 83, 62]),
("Set 5", [484, 114, 205, 288, 506, 125, 503, 201, 127, 410, 312, 132, 312, 542]),
]
for name, weights in weight_sets:
n = len(weights)
model = LpProblem(sense=LpMinimize)
# variables - binary x[i] based on if the item belongs to the first partition
variables = [LpVariable(name=f"x_{i}", cat=LpBinary) for i in range(n)]
# ... and also the difference between the left and right side
difference = LpVariable(name=f"difference", lowBound=0, cat=LpInteger)
# inequalities
## the difference of the parts equals their difference (variable)
model += difference == lpSum([variables[i] * weights[i] for i in range(n)]) \
- lpSum([(1 - variables[i]) * weights[i] for i in range(n)])
# objective function - minimize the difference between the sides
model += difference
status = model.solve(PULP_CBC_CMD(msg=False))
print(name + "\n" + "-" * len(name))
print("difference:", difference.value())
for i in range(n):
print(int(variables[i].value()), end=" ")
print()
print()
Output
Set 1
-----
difference: 0.0
1 0 1 0 0
Set 2
-----
difference: 1.0
1 1 1 1 0 1 0 0 1 0
Set 3
-----
difference: 0.0
1 1 1 0 0 1 0 0 1 1
Set 4
-----
difference: 0.0
0 1 0 0 1 0 1 1 1 1
Set 5
-----
difference: 1.0
1 1 1 1 0 0 1 0 1 1 0 0 0 0
Maximum independent set
Given a graph , find the largest set of vertices such that no two share an edge.
Source code
from pulp import *
graphs = [
[
"Complete graph on 5 vertices",
[0, 1, 2, 3, 4],
[(0, 1), (2, 1), (1, 3), (2, 3), (4, 3), (2, 4), (0, 4), (1, 4), (2, 0), (0, 3)],
],
[
"Path of 6 vertices",
[0, 1, 2, 3, 4, 5],
[(0, 1), (2, 1), (2, 3), (3, 4), (4, 5)],
],
[
"Petersen graph",
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
[(0, 1), (0, 4), (0, 5), (1, 2), (1, 6), (2, 3), (2, 7), (3, 4), (3, 8), (4, 9), (5, 7), (5, 8), (6, 8), (6, 9), (7, 9)],
]
]
for name, vertices, edges in graphs:
n = len(vertices)
model = LpProblem(sense=LpMaximize)
# variables - x[i] binary based on if the i-th vertex is in the set
variables = [LpVariable(name=f"x_{i}", cat=LpBinary) for i in range(n)]
# inequalities
## each edge must have at most one vertex from the set
for u, v in edges:
model += variables[u - 1] + variables[v - 1] <= 1
# minimize the number of selected edges
model += lpSum(variables)
status = model.solve(PULP_CBC_CMD(msg=False))
print(name + "\n" + "-" * len(name))
print(f"set size =", int(model.objective.value()))
for i in range(n):
print(int(variables[i].value()), end=" ")
print()
print()
Output
Complete graph on 5 vertices
----------------------------
set size = 1
1 0 0 0 0
Path of 6 vertices
------------------
set size = 3
0 1 0 0 1 1
Petersen graph
--------------
set size = 4
1 0 1 0 1 0 0 0 1 0
Minimum vertex cover
Given a graph , find the smallest set of vertices that cover all edges.
Source code
"""
NOTE: This code is the same as max-independent-set.py (the problems are complimentary),
with two changes:
- line 30: we minimize instead of maximize
- line 39: we want >= instead of <=
"""
from pulp import *
graphs = [
[
"Complete graph on 5 vertices",
[0, 1, 2, 3, 4],
[(0, 1), (2, 1), (1, 3), (2, 3), (4, 3), (2, 4), (0, 4), (1, 4), (2, 0), (0, 3)],
],
[
"Path of 6 vertices",
[0, 1, 2, 3, 4, 5],
[(0, 1), (2, 1), (2, 3), (3, 4), (4, 5)],
],
[
"Petersen graph",
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
[(0, 1), (0, 4), (0, 5), (1, 2), (1, 6), (2, 3), (2, 7), (3, 4), (3, 8), (4, 9), (5, 7), (5, 8), (6, 8), (6, 9), (7, 9)],
]
]
for name, vertices, edges in graphs:
n = len(vertices)
model = LpProblem(sense=LpMinimize)
# variables - x[i] binary based on if the i-th vertex is in the set
variables = [LpVariable(name=f"x_{i}", cat=LpBinary) for i in range(n)]
# inequalities
## each edge must have at most one vertex from the set
for u, v in edges:
model += variables[u - 1] + variables[v - 1] >= 1
# minimize the number of selected edges
model += lpSum(variables)
status = model.solve(PULP_CBC_CMD(msg=False))
print(name + "\n" + "-" * len(name))
print(f"set size =", int(model.objective.value()))
for i in range(n):
print(int(variables[i].value()), end=" ")
print()
print()
Output
Complete graph on 5 vertices
----------------------------
set size = 4
1 0 1 1 1
Path of 6 vertices
------------------
set size = 3
0 1 0 1 0 1
Petersen graph
--------------
set size = 6
0 1 0 1 0 1 1 1 0 1
Resources
Here are all the resources I used for data and documentation: