Solve Sudokus using MILP
In this vignettes we will solve Sudoku puzzles using MILP. Sudoku in its most popular form is a constraint satisfaction problem and by setting the objective function to \(0\) you transform the optimization problem into a pure constraint satistication problem. In this document we will consider Sudokus in a 9x9 grid with 3x3 sub-matrices.
Of course you can formulate an objective function as well that directs the solver towards solutions maximizing a certain linear function.
The model
The idea is to introduce a binary variable \(x\) with three indexes \(i, j, k\) that is \(1\) if and only if the number \(k\) is in cell \(i, j\).
library(ompr)
library(dplyr)
n <- 9
model <- MIPModel() %>%
# The number k stored in position i,j
add_variable(x[i, j, k], i = 1:n, j = 1:n, k = 1:9, type = "binary") %>%
# no objective
set_objective(0) %>%
# only one number can be assigned per cell
add_constraint(sum_over(x[i, j, k], k = 1:9) == 1, i = 1:n, j = 1:n) %>%
# each number is exactly once in a row
add_constraint(sum_over(x[i, j, k], j = 1:n) == 1, i = 1:n, k = 1:9) %>%
# each number is exactly once in a column
add_constraint(sum_over(x[i, j, k], i = 1:n) == 1, j = 1:n, k = 1:9) %>%
# each 3x3 square must have all numbers
add_constraint(sum_over(x[i, j, k], i = 1:3 + sx, j = 1:3 + sy) == 1,
sx = seq(0, n - 3, 3), sy = seq(0, n - 3, 3), k = 1:9)
model## Mixed integer linear optimization problem
## Variables:
## Continuous: 0
## Integer: 0
## Binary: 729
## Model sense: maximize
## Constraints: 324Solve the model
We will use glpk to solve the above model. Note that we
haven’t fixed any numbers to specific values. That means that the solver
will find a valid sudoku without any prior hints.
library(ompr.roi)
library(ROI.plugin.glpk)
result <- solve_model(model, with_ROI(solver = "glpk", verbose = TRUE))## <SOLVER MSG> ----
## GLPK Simplex Optimizer 5.0
## 324 rows, 729 columns, 2916 non-zeros
## 0: obj = -0.000000000e+00 inf = 3.240e+02 (324)
## 319: obj = -0.000000000e+00 inf = 2.141e-14 (0) 1
## OPTIMAL LP SOLUTION FOUND
## GLPK Integer Optimizer 5.0
## 324 rows, 729 columns, 2916 non-zeros
## 729 integer variables, all of which are binary
## Integer optimization begins...
## Long-step dual simplex will be used
## + 319: mip = not found yet <= +inf (1; 0)
## + 1102: >>>>> 0.000000000e+00 <= 0.000000000e+00 0.0% (46; 0)
## + 1102: mip = 0.000000000e+00 <= tree is empty 0.0% (0; 91)
## INTEGER OPTIMAL SOLUTION FOUND
## <!SOLVER MSG> ----
# the following dplyr statement plots a 9x9 matrix
result %>%
get_solution(x[i,j,k]) %>%
filter(value > 0) %>%
select(i, j, k) %>%
tidyr::spread(j, k) %>%
select(-i)## 1 2 3 4 5 6 7 8 9
## 1 6 8 1 5 9 7 4 3 2
## 2 4 3 2 8 1 6 7 5 9
## 3 5 7 9 3 4 2 8 1 6
## 4 7 6 4 2 5 1 3 9 8
## 5 1 2 8 9 7 3 6 4 5
## 6 3 9 5 4 6 8 1 2 7
## 7 2 1 6 7 3 5 9 8 4
## 8 8 4 7 1 2 9 5 6 3
## 9 9 5 3 6 8 4 2 7 1If you want to solve a concrete sudoku you can fix certain cells to specific values. For example here we solve a sudoku that has the sequence from 1 to 9 in the first 3x3 matrix fixed.
model_fixed <- model %>%
add_constraint(x[1, 1, 1] == 1) %>%
add_constraint(x[1, 2, 2] == 1) %>%
add_constraint(x[1, 3, 3] == 1) %>%
add_constraint(x[2, 1, 4] == 1) %>%
add_constraint(x[2, 2, 5] == 1) %>%
add_constraint(x[2, 3, 6] == 1) %>%
add_constraint(x[3, 1, 7] == 1) %>%
add_constraint(x[3, 2, 8] == 1) %>%
add_constraint(x[3, 3, 9] == 1)
result <- solve_model(model_fixed, with_ROI(solver = "glpk", verbose = TRUE))## <SOLVER MSG> ----
## GLPK Simplex Optimizer 5.0
## 333 rows, 729 columns, 2925 non-zeros
## 0: obj = -0.000000000e+00 inf = 3.330e+02 (333)
## 360: obj = -0.000000000e+00 inf = 2.197e-13 (0) 1
## OPTIMAL LP SOLUTION FOUND
## GLPK Integer Optimizer 5.0
## 333 rows, 729 columns, 2925 non-zeros
## 729 integer variables, all of which are binary
## Integer optimization begins...
## Long-step dual simplex will be used
## + 360: mip = not found yet <= +inf (1; 0)
## + 861: >>>>> 0.000000000e+00 <= 0.000000000e+00 0.0% (30; 0)
## + 861: mip = 0.000000000e+00 <= tree is empty 0.0% (0; 59)
## INTEGER OPTIMAL SOLUTION FOUND
## <!SOLVER MSG> ----
result %>%
get_solution(x[i,j,k]) %>%
filter(value > 0) %>%
select(i, j, k) %>%
tidyr::spread(j, k) %>%
select(-i) ## 1 2 3 4 5 6 7 8 9
## 1 1 2 3 4 8 9 5 7 6
## 2 4 5 6 2 1 7 9 3 8
## 3 7 8 9 5 6 3 1 2 4
## 4 3 1 4 6 7 5 8 9 2
## 5 5 6 8 3 9 2 4 1 7
## 6 9 7 2 1 4 8 6 5 3
## 7 6 3 7 9 5 4 2 8 1
## 8 8 4 5 7 2 1 3 6 9
## 9 2 9 1 8 3 6 7 4 5Feedback
Do you have any questions, ideas, comments? Or did you find a mistake? Let’s discuss on Github.