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\).

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 0 
##   Integer: 0 
##   Binary: 729 
## Model sense: maximize 
## Constraints: 324

Solve 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.

## <SOLVER MSG>  ----
## GLPK Simplex Optimizer, v4.63
## 324 rows, 729 columns, 2916 non-zeros
##       0: obj =  -0.000000000e+00 inf =   3.240e+02 (324)
##     340: obj =  -0.000000000e+00 inf =   4.087e-14 (0) 1
## OPTIMAL LP SOLUTION FOUND
## GLPK Integer Optimizer, v4.63
## 324 rows, 729 columns, 2916 non-zeros
## 729 integer variables, all of which are binary
## Integer optimization begins...
## +   340: mip =     not found yet <=              +inf        (1; 0)
## +  1075: >>>>>   0.000000000e+00 <=   0.000000000e+00   0.0% (43; 0)
## +  1075: mip =   0.000000000e+00 <=     tree is empty   0.0% (0; 85)
## INTEGER OPTIMAL SOLUTION FOUND
## <!SOLVER MSG> ----
##   1 2 3 4 5 6 7 8 9
## 1 4 2 3 5 7 8 1 9 6
## 2 1 7 5 9 3 6 4 8 2
## 3 9 8 6 2 1 4 7 3 5
## 4 8 6 4 1 5 9 3 2 7
## 5 3 5 9 6 2 7 8 1 4
## 6 2 1 7 8 4 3 5 6 9
## 7 7 9 1 3 6 5 2 4 8
## 8 6 4 2 7 8 1 9 5 3
## 9 5 3 8 4 9 2 6 7 1

If 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, v4.63
## 333 rows, 729 columns, 2925 non-zeros
##       0: obj =  -0.000000000e+00 inf =   3.330e+02 (333)
##     354: obj =  -0.000000000e+00 inf =   1.250e-13 (0) 1
## OPTIMAL LP SOLUTION FOUND
## GLPK Integer Optimizer, v4.63
## 333 rows, 729 columns, 2925 non-zeros
## 729 integer variables, all of which are binary
## Integer optimization begins...
## +   354: mip =     not found yet <=              +inf        (1; 0)
## +   992: >>>>>   0.000000000e+00 <=   0.000000000e+00   0.0% (44; 0)
## +   992: mip =   0.000000000e+00 <=     tree is empty   0.0% (0; 87)
## 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 5 8 4 9 7 6
## 2 4 5 6 9 7 3 8 2 1
## 3 7 8 9 2 1 6 3 5 4
## 4 3 1 5 6 4 7 2 9 8
## 5 9 7 4 1 2 8 5 6 3
## 6 8 6 2 3 9 5 4 1 7
## 7 2 4 8 7 5 1 6 3 9
## 8 5 3 1 8 6 9 7 4 2
## 9 6 9 7 4 3 2 1 8 5

Feedback

Do you have any questions, ideas, comments? Or did you find a mistake? Let’s discuss on Github.