# slama.dev

## Linear Programming in Python

YouTube , released 4. 7. 2023, updated 13. 9. 2023, Česká verze článku

### 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 $3$ tons of potato seeds and $4$ tons of carrot seeds. To grow the crops efficiently, you also have $5$ tons of fertilizer, which has to be used when planting in a $1:1$ ratio (i.e. $1$ kilogram of potatoes or carrots requires $1$ kilogram of fertilizer). The profit is $1.2\/\mathrm{kg}$ for potato seeds and $1.7\/\mathrm{kg}$ 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 $n$ items, each with weight $w_i$ and price $p_i$, our task is to maximize the price of the items we take into our backpack without exceeding its carry weight $M$.

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 $k$ such that a graph $G$ is vertex $k$-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 $G$, 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",
[,
[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 $n$ items with weights $w_1, \ldots, w_n$ and an arbitrary number of bins with maximum carry weight $C$, 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 $n$ items with weights $w_1, \ldots, w_n$, 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 $G$, 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 $G$, 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: