I am writing a Sudoku solver in Gekko (mostly for fun, I'm aware there are better tools for this.) I would like to express the constraints the variables in each row, column, and square most be different from each other. Here is the code:
import itertools
from gekko import GEKKO
SQUARES = [
( (0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2) ),
( (0, 3), (0, 4), (0, 5), (1, 3), (1, 4), (1, 5), (2, 3), (2, 4), (2, 5) ),
( (0, 6), (0, 7), (0, 8), (1, 6), (1, 7), (1, 8), (2, 6), (2, 7), (2, 8) ),
( (3, 0), (3, 1), (3, 2), (4, 0), (4, 1), (4, 2), (5, 0), (5, 1), (5, 2) ),
( (3, 3), (3, 4), (3, 5), (4, 3), (4, 4), (4, 5), (5, 3), (5, 4), (5, 5) ),
( (3, 6), (3, 7), (3, 8), (4, 6), (4, 7), (4, 8), (5, 6), (5, 7), (5, 8) ),
( (6, 0), (6, 1), (6, 2), (7, 0), (7, 1), (7, 2), (8, 0), (8, 1), (8, 2) ),
( (6, 3), (6, 4), (6, 5), (7, 3), (7, 4), (7, 5), (8, 3), (8, 4), (8, 5) ),
( (6, 6), (6, 7), (6, 8), (7, 6), (7, 7), (7, 8), (8, 6), (8, 7), (8, 8) )
]
BOARD_DIMS = (9, 9)
PRE_FILLED_SQUARES = [
(0, 1, 2), (0, 3, 5), (0, 5, 1),
(0, 7, 9), (1, 0, 8), (1, 3, 2),
(1, 3, 2), (1, 5, 3), (1, 8, 6),
(2, 1, 3), (2, 4, 6), (2, 7, 7),
(3, 2, 1), (3, 6, 6), (4, 0, 5),
(4, 1, 4), (4, 7, 1), (4, 8, 9),
(5, 2, 2), (5, 6, 7), (6, 1, 9),
(6, 4, 3), (6, 7, 8), (7, 0, 2),
(7, 3, 8), (7, 5, 4), (7, 8, 7),
(8, 1, 1), (8, 3, 9), (8, 5, 7),
(8, 7, 6)
]
m = GEKKO()
m.options.SOLVER = 1
board = m.Array(m.Var, BOARD_DIMS, lb=1, ub=9, integer=True, value=1)
for i, j, val in PRE_FILLED_SQUARES:
board[i,j].value = val
m.Equation(board[i,j]==val)
for row in range(BOARD_DIMS[0]):
for i, j in itertools.combinations(range(BOARD_DIMS[1]), r=2):
diff = m.if3(board[row, i] - board[row, j], -1, 1)
diff_prime = m.if3(board[row, j] - board[row, i], -1, 1)
m.Equation(diff + diff_prime == 0)
for col in range(BOARD_DIMS[1]):
for i, j in itertools.combinations(range(BOARD_DIMS[0]), r=2):
diff = m.if3(board[i, col] - board[j, col], -1, 1)
diff_prime = m.if3(board[j, col] - board[i, col], -1, 1)
m.Equation(diff + diff_prime == 0)
for square in SQUARES:
for (y, x), (y_prime, x_prime) in itertools.combinations(square, r=2):
diff = m.if3(board[y, x] - board[y_prime, x_prime], -1, 1)
diff_prime = m.if3(board[y_prime, x_prime] - board[y, x], -1, 1)
m.Equation(diff + diff_prime == 0)
m.solve(disp=True)
print(board)
Currently, the solvers respects the pre-filled squares I've filled in but fills every other square with
1.0. The current expression of inequality, with the double subtraction, is because trying something along the lines of m.Equation(board[row,j] != board[row,i]) gave an error to the effect that objects of type "int" don't have a length.
There is no != (not equal) constraint in gekko. You can set the sum of each row to be 1 where each value can be 0 or 1. If the value is equal to 1 then that corresponding value is filled in that cell.
Here is an MILP version that is written in PuLP and solved with the MILP CBC solver (credit to TowardsDataScience)
import pulp as plp
def add_default_sudoku_constraints(prob, grid_vars, rows, cols, grids, values):
# Constraint to ensure only one value is filled for a cell
for row in rows:
for col in cols:
prob.addConstraint(plp.LpConstraint(e=plp.lpSum([grid_vars[row][col][value] for value in values]),
sense=plp.LpConstraintEQ, rhs=1, name=f"constraint_sum_{row}_{col}"))
# Constraint to ensure that values from 1 to 9 is filled only once in a row
for row in rows:
for value in values:
prob.addConstraint(plp.LpConstraint(e=plp.lpSum([grid_vars[row][col][value]*value for col in cols]),
sense=plp.LpConstraintEQ, rhs=value, name=f"constraint_uniq_row_{row}_{value}"))
# Constraint to ensure that values from 1 to 9 is filled only once in a column
for col in cols:
for value in values:
prob.addConstraint(plp.LpConstraint(e=plp.lpSum([grid_vars[row][col][value]*value for row in rows]),
sense=plp.LpConstraintEQ, rhs=value, name=f"constraint_uniq_col_{col}_{value}"))
# Constraint to ensure that values from 1 to 9 is filled only once in the 3x3 grid
for grid in grids:
grid_row = int(grid/3)
grid_col = int(grid%3)
for value in values:
prob.addConstraint(plp.LpConstraint(e=plp.lpSum([grid_vars[grid_row*3+row][grid_col*3+col][value]*value for col in range(0,3) for row in range(0,3)]),
sense=plp.LpConstraintEQ, rhs=value, name=f"constraint_uniq_grid_{grid}_{value}"))
def add_diagonal_sudoku_constraints(prob, grid_vars, rows, cols, values):
# Constraint from top-left to bottom-right - numbers 1 - 9 should not repeat
for value in values:
prob.addConstraint(plp.LpConstraint(e=plp.lpSum([grid_vars[row][row][value]*value for row in rows]),
sense=plp.LpConstraintEQ, rhs=value, name=f"constraint_uniq_diag1_{value}"))
# Constraint from top-right to bottom-left - numbers 1 - 9 should not repeat
for value in values:
prob.addConstraint(plp.LpConstraint(e=plp.lpSum([grid_vars[row][len(rows)-row-1][value]*value for row in rows]),
sense=plp.LpConstraintEQ, rhs=value, name=f"constraint_uniq_diag2_{value}"))
def add_prefilled_constraints(prob, input_sudoku, grid_vars, rows, cols, values):
for row in rows:
for col in cols:
if(input_sudoku[row][col] != 0):
prob.addConstraint(plp.LpConstraint(e=plp.lpSum([grid_vars[row][col][value]*value for value in values]),
sense=plp.LpConstraintEQ,
rhs=input_sudoku[row][col],
name=f"constraint_prefilled_{row}_{col}"))
def extract_solution(grid_vars, rows, cols, values):
solution = [[0 for col in cols] for row in rows]
grid_list = []
for row in rows:
for col in cols:
for value in values:
if plp.value(grid_vars[row][col][value]):
solution[row][col] = value
return solution
def print_solution(solution, rows,cols):
# Print the final result
print(f"\nFinal result:")
print("\n\n+ ----------- + ----------- + ----------- +",end="")
for row in rows:
print("\n",end="\n| ")
for col in cols:
num_end = " | " if ((col+1)%3 == 0) else " "
print(solution[row][col],end=num_end)
if ((row+1)%3 == 0):
print("\n\n+ ----------- + ----------- + ----------- +",end="")
def solve_sudoku(input_sudoku, diagonal = False ):
# Create the linear programming problem
prob = plp.LpProblem("Sudoku_Solver")
rows = range(0,9)
cols = range(0,9)
grids = range(0,9)
values = range(1,10)
# Decision Variable/Target variable
grid_vars = plp.LpVariable.dicts("grid_value", (rows,cols,values), cat='Binary')
# Set the objective function
# Sudoku works only on the constraints - feasibility problem
# There is no objective function that we are trying maximize or minimize.
# Set a dummy objective
objective = plp.lpSum(0)
prob.setObjective(objective)
# Create the default constraints to solve sudoku
add_default_sudoku_constraints(prob, grid_vars, rows, cols, grids, values)
# Add the diagonal constraints if flag is set
if (diagonal):
add_diagonal_sudoku_constraints(prob, grid_vars, rows, cols, values)
# Fill the prefilled values from input sudoku as constraints
add_prefilled_constraints(prob, input_sudoku, grid_vars, rows, cols, values)
# Solve the problem
prob.solve()
# Print the status of the solution
solution_status = plp.LpStatus[prob.status]
print(f'Solution Status = {plp.LpStatus[prob.status]}')
# Extract the solution if an optimal solution has been identified
if solution_status == 'Optimal':
solution = extract_solution(grid_vars, rows, cols, values)
print_solution(solution, rows,cols)
normal_sudoku = [
[3,0,0,8,0,0,0,0,1],
[0,0,0,0,0,2,0,0,0],
[0,4,1,5,0,0,8,3,0],
[0,2,0,0,0,1,0,0,0],
[8,5,0,4,0,3,0,1,7],
[0,0,0,7,0,0,0,2,0],
[0,8,5,0,0,9,7,4,0],
[0,0,0,1,0,0,0,0,0],
[9,0,0,0,0,7,0,0,6]
]
solve_sudoku(input_sudoku=normal_sudoku, diagonal=False)
diagonal_sudoku = [
[0,3,0,2,7,0,0,0,0],
[0,0,0,0,0,0,0,0,0],
[8,0,0,0,0,0,0,0,0],
[5,1,0,0,0,0,0,8,4],
[4,0,0,5,9,0,0,7,0],
[2,9,0,0,0,0,0,1,0],
[0,0,0,0,0,0,1,0,5],
[0,0,6,3,0,8,0,0,7],
[0,0,0,0,0,0,3,0,0]
]
solve_sudoku(input_sudoku=diagonal_sudoku, diagonal=True)
The solution time is very fast with the CBC solver.
+ ----------- + ----------- + ----------- +
| 6 3 9 | 2 7 4 | 8 5 1 |
| 1 4 2 | 6 8 5 | 7 3 9 |
| 8 7 5 | 1 3 9 | 6 4 2 |
+ ----------- + ----------- + ----------- +
| 5 1 3 | 7 6 2 | 9 8 4 |
| 4 6 8 | 5 9 1 | 2 7 3 |
| 2 9 7 | 8 4 3 | 5 1 6 |
+ ----------- + ----------- + ----------- +
| 3 8 4 | 9 2 7 | 1 6 5 |
| 9 5 6 | 3 1 8 | 4 2 7 |
| 7 2 1 | 4 5 6 | 3 9 8 |
+ ----------- + ----------- + ----------- +
Translating this equivalent problem to gekko is easy with the similar syntax of the two modeling languages:
from gekko import GEKKO
def add_default_sudoku_constraints(m, grid_vars, rows, cols, grids, values):
for row in rows:
for col in cols:
m.Equation(sum([grid_vars[row][col][value-1] \
for value in values]) == 1)
for row in rows:
for value in values:
m.Equation(sum([grid_vars[row][col][value-1]*value \
for col in cols]) == value)
for col in cols:
for value in values:
m.Equation(sum([grid_vars[row][col][value-1]*value \
for row in rows]) == value)
for grid in grids:
grid_row = int(grid/3)
grid_col = int(grid%3)
for value in values:
m.Equation(sum([grid_vars[grid_row*3+row][grid_col*3+col][value-1]*value \
for col in range(0,3) for row in range(0,3)]) == value)
def add_diagonal_sudoku_constraints(m, grid_vars, rows, cols, values):
for value in values:
m.Equation(sum([grid_vars[row][row][value-1]*value for row in rows]) == value)
m.Equation(sum([grid_vars[row][len(rows)-row-1][value-1]*value \
for row in rows]) == value)
def add_prefilled_constraints(m, input_sudoku, grid_vars, rows, cols, values):
for row in rows:
for col in cols:
if(input_sudoku[row][col] != 0):
m.Equation(sum([grid_vars[row][col][value-1]*value \
for value in values]) == input_sudoku[row][col])
def extract_solution(grid_vars, rows, cols, values):
solution = [[0 for col in cols] for row in rows]
for row in rows:
for col in cols:
for v, value in enumerate(values):
if int(grid_vars[row][col][v].value[0]) == 1:
solution[row][col] = value
return solution
def print_solution(solution, rows, cols):
print(f"\nFinal result:")
print("\n\n+ ----------- + ----------- + ----------- +",end="")
for row in rows:
print("\n",end="\n| ")
for col in cols:
num_end = " | " if ((col+1)%3 == 0) else " "
print(solution[row][col],end=num_end)
if ((row+1)%3 == 0):
print("\n\n+ ----------- + ----------- + ----------- +",end="")
def solve_sudoku(input_sudoku, diagonal=False):
m = GEKKO(remote=True)
rows = range(0,9)
cols = range(0,9)
grids = range(0,9)
values = range(1,10)
grid_vars = [[[m.Var(lb=0, ub=1, integer=True) \
for _ in values] for _ in cols] for _ in rows]
add_default_sudoku_constraints(m, grid_vars, rows, cols, grids, values)
if diagonal:
add_diagonal_sudoku_constraints(m, grid_vars, rows, cols, values)
add_prefilled_constraints(m, input_sudoku, grid_vars, rows, cols, values)
m.options.SOLVER = 3
m.solve(disp=True)
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 1000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1']
m.options.SOLVER = 1
m.solve(disp=True)
solution = extract_solution(grid_vars, rows, cols, values)
print_solution(solution, rows, cols)
normal_sudoku = [
[3,0,0,8,0,0,0,0,1],
[0,0,0,0,0,2,0,0,0],
[0,4,1,5,0,0,8,3,0],
[0,2,0,0,0,1,0,0,0],
[8,5,0,4,0,3,0,1,7],
[0,0,0,7,0,0,0,2,0],
[0,8,5,0,0,9,7,4,0],
[0,0,0,1,0,0,0,0,0],
[9,0,0,0,0,7,0,0,6]
]
solve_sudoku(input_sudoku=normal_sudoku, diagonal=False)
diagonal_sudoku = [
[0,3,0,2,7,0,0,0,0],
[0,0,0,0,0,0,0,0,0],
[8,0,0,0,0,0,0,0,0],
[5,1,0,0,0,0,0,8,4],
[4,0,0,5,9,0,0,7,0],
[2,9,0,0,0,0,0,1,0],
[0,0,0,0,0,0,1,0,5],
[0,0,6,3,0,8,0,0,7],
[0,0,0,0,0,0,3,0,0]
]
solve_sudoku(input_sudoku=diagonal_sudoku, diagonal=True)
However, the APOPT solver is a Nonlinear Mixed Integer Programming (MINLP) solver and is much slower than dedicated MILP solvers for linear mixed integer problems. Sometimes using these solver options can help:
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 1000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1']
However, in most cases it is better to use a dedicated MILP solver.