Nelioton Ruudukko MIP

n = 8
p = MixedIntegerLinearProgram(maximization=True, solver = "GLPK")
w = p.new_variable(binary=True)

for i in range(n):
    for j in range(n):
        for d in range(1, min(n-i, n-j)):
            p.add_constraint(w[i, j] + w[i+d, j] + w[i+d, j+d] + w[i, j+d] <=3)
p.set_objective(sum(w[i,j] for i in range(n) for j in range(n)))

show(p)

solVal = p.solve()
grid = [[int(p.get_values(w[i,j])) for j in range(n)] for i in range(n)]

print (solVal)
show (matrix_plot(grid))