我正在用 Gekko 编写一个数独解算器（主要是为了好玩，我知道有更好的工具可以实现这一点。） 我想表达每行、每列和正方形中的变量的约束最不同于 彼此。这是代码：
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)
目前，解算器尊重我已填充的预填充方块，但用 填充其他所有方块
1.0m.Equation(board[row,j] != board[row,i])!=gekko1011这是一个用 PuLP 编写的 MILP 版本，并使用 MILP CBC 求解器进行求解（归功于 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)
使用 CBC 求解器，求解时间非常快。
+ ----------- + ----------- + ----------- +
| 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 |
+ ----------- + ----------- + ----------- +
使用两种建模语言的相似语法，将这个等效问题转换为 gekko 很容易：
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)
但是，APOPT 求解器是非线性混合整数规划 (MINLP) 求解器，比线性混合整数问题的专用 MILP 求解器慢得多。有时使用这些求解器选项会有所帮助：
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 1000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1']
但是，在大多数情况下，最好使用专用的 MILP 求解器。