Linear Programming in Python

YouTube · 4. 7. 2023 · Česká verze článku · [edit]

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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

1
2
3
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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="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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
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("set size =", int(model.objective.value()))

    for i in range(n):
        print(int(variables[i].value()), end=" ")
    print()
    print()

Output

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
"""
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("set size =", int(model.objective.value()))

    for i in range(n):
        print(int(variables[i].value()), end=" ")
    print()
    print()

Output

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
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: