How to optimize a linear problem?

PuLP is a simple Linear Problem Solver. It solves systems with many linear equations, e.g. Sudokus, scheduling trains and many more. For complex LES problems you will need commercial software (like Gurobi), but PuLP helps you get the taste of it.

Install it with:

pip install pulp

The example below solves a magic square (insert the numbers 1..9 in to a 3x3 square so that all rows and columns sum up to 15):

from pulp import LpProblem, LpMinimize, LpVariable, LpInteger, lpSum, LpStatus, value

prob = LpProblem("PULPTEST", LpMinimize)

# model variables
XCOORD = [0, 1, 2]
YCOORD = [0, 1, 2]
NUMBERS = [1, 2, 3, 4, 5, 6, 7, 8, 9]

# variable is a 3 x 3 x 9 matrix of binary values
allocation = LpVariable.dicts("square", (XCOORD, YCOORD, NUMBERS), 0, 1, LpInteger)

# target function
prob += 0, "Arbitrary Objective Function"

# constraint: sum over rows
for x in XCOORD:
    prob += lpSum([n * allocation[x][y][n] for y in YCOORD for n in NUMBERS]) == 15

# constraint: sum over columns
for y in YCOORD:
    prob += lpSum([n * allocation[x][y][n] for x in XCOORD for n in NUMBERS]) == 15

# constraint: each number only once
for n in NUMBERS:
    prob += lpSum([allocation[x][y][n] for x in XCOORD for y in YCOORD]) == 1

# constraint: three numbers per column
for x in XCOORD:
    prob += lpSum([allocation[x][y][n] for y in YCOORD for n in NUMBERS]) == 3

# constraint: three numbers per row
for y in YCOORD:
    prob += lpSum([allocation[x][y][n] for x in XCOORD for n in NUMBERS]) == 3

# constraint: 9 numbers set
prob += lpSum([allocation[x][y][n] for x in XCOORD for y in YCOORD for n in NUMBERS]) == 9

# run the solver
prob.solve()

print("Status:", LpStatus[prob.status])

# print the numbers that have been found
for y in YCOORD:
    for x in XCOORD:
        for n in NUMBERS:
            if value(allocation[x][y][n]) == 1:
                print(n, end=' ')
    print()