Sudoku puzzle solving
This example demonstrates how to find a solution to a Sudoku puzzle using Mixed Integer Linear Programming (MILP) and Constraint Programming (CP) methods.
We will solve this problem using the [problem-oriented optimisation] workflow(https://engee.com/community/catalogs/projects/algoritm-resheniia-zadach-optimizatsii).
Installing the libraries
If your environment does not have the latest version of the JuMP package installed , uncomment and run the box below:
Pkg.add(["JuMP", "HiGHS"])
#Pkg.add("JuMP");
To launch a new version of the library after the installation is complete, click on the "My Account" button:
Then click on the "Stop" button:
Restart the session by pressing the "Start Engee" button:
Task Description
Sudoku is a logic combinatorial number placement puzzle game. The goal is to fill a 9x9 grid with numbers so that each column, each row and each of the nine 3x3 sub-grids contains all numbers from 1 to 9.
Rules:
- Each line must contain digits from 1 to 9 without repetition.
- Each column must contain digits from 1 to 9 without repetition.
- Each 3x3 subnet must contain digits from 1 to 9 without repetition.
Initial hints are provided to make solving the puzzle possible:
.png)
An example of a solved sudoku:
.png)
Sudoku puzzles range from simple to extremely complex, with the difficulty often determined by the number of initial clues provided and the complexity of the strategies required to solve them. A properly constructed sudoku puzzle should have only one solution.
Connecting libraries
Connect the library JuMP:
using JuMP;
Connect the solver library HiGHS:
using HiGHS;
Connect the graphing library Plots:
using Plots;
Connect the random number generation library Random:
using Random;
Creating a Sudoku visualisation function
Create a function plot_sudoku(), which visualises sudoku:
function plot_sudoku(grid)
p = plot(aspect_ratio=:equal, legend=false, axis=false, ticks=false)
for i in 0:9
plot!(p, [i, i], [0, 9], color=:black, linewidth=(i % 3 == 0 ? 2 : 1))
plot!(p, [0, 9], [i, i], color=:black, linewidth=(i % 3 == 0 ? 2 : 1))
end
for i in 1:9, j in 1:9
if grid[i,j] != 0
annotate!(p, j-0.5, 9.5-i, text(string(grid[i,j]), 14))
end
end
return p
end
We will use this function to visualise sudoku in this example.
Generating Sudoku puzzles
In this section, we will write a programme that will generate sudoku puzzles to solve.
Create a function generate_full_grid(), which is responsible for building the sudoku grid:
function generate_full_grid()
grid = zeros(Int, 9, 9)
fill_grid!(grid, 1)
return grid
end
Create a recursive function fill_grid!(), which fills the sudoku grid with numbers:
function fill_grid!(grid, pos)
if pos > 81
return true
end
row, col = divrem(pos - 1, 9) .+ 1
nums = shuffle(1:9)
for num in nums
if is_valid(grid, row, col, num)
grid[row, col] = num
if fill_grid!(grid, pos + 1)
return true
end
grid[row, col] = 0
end
end
return false
end
This function implements a backtracking algorithm to fill the sudoku grid with valid numbers.
Create a function is_valid(), which checks if a given number (num) is allowed to be placed in a certain cell (row, column) of the sudoku grid:
function is_valid(grid, row, col, num)
for i in 1:9
if grid[row, i] == num || grid[i, col] == num
return false
end
end
box_row, box_col = 3 * ((row - 1) ÷ 3) + 1, 3 * ((col - 1) ÷ 3) + 1
for i in 0:2, j in 0:2
if grid[box_row + i, box_col + j] == num
return false
end
end
return true
end
This function is crucial for following all three rules in the sudoku grid filling process.
Create a function remove_numbers!(), which is responsible for removing numbers from a filled sudoku grid, to create a puzzle:
function remove_numbers!(grid, num_to_remove)
positions = [(i, j) for i in 1:9 for j in 1:9]
shuffle!(positions)
for _ in 1:num_to_remove
for (row, col) in positions
if grid[row, col] != 0
temp = grid[row, col]
grid[row, col] = 0
solution_count = count_solutions(grid)
if solution_count != 1
grid[row, col] = temp
else
break
end
end
end
end
end
This is a key part of creating sudoku puzzles of varying levels of difficulty, as the number of cells removed can affect the difficulty of the puzzle.
Create a recursive function count_solutions(), which counts the number of valid solutions for a given sudoku grid:
function count_solutions(grid, limit=2)
empty_cell = findfirst(iszero, grid)
if empty_cell === nothing
return 1
end
row, col = Tuple(empty_cell)
count = 0
for num in 1:9
if is_valid(grid, row, col, num)
grid[row, col] = num
count += count_solutions(grid, limit - count)
grid[row, col] = 0
if count >= limit
return count
end
end
end
return count
end
This function is crucial to ensure that removing numbers from the grid (in the function remove_numbers!) does not result in multiple solutions.
Create a function generate_sudoku_puzzle() that generates Sudoku with a given number of clues:
function generate_sudoku_puzzle(num_clues)
full_grid = generate_full_grid()
num_to_remove = 81 - num_clues
remove_numbers!(full_grid, num_to_remove)
return full_grid
end
The theoretical minimum number of hints for a sudoku puzzle with a single solution is 17.
The programme we have written can consistently generate Sudoku with a single solution, where the minimum number of hints is somewhere between 20 and 30, which is more than enough for this example.
Generate Sudoku with 25 clues:
sudoku_puzzle = generate_sudoku_puzzle(25);
Visualise Sudoku with clues:
plot_sudoku(sudoku_puzzle)
Formulation of a mixed integer linear programming problem
We will solve the problem using mixed integer linear programming. The goal is to model the puzzle as a feasibility problem where each cell must satisfy certain constraints.
Create an optimisation problem using the function Model() and give the name of the solver in brackets:
sudoku_milp = Model(HiGHS.Optimizer)
Create a variable x, representing a three-dimensional array , where i, j, and k are in the range 1 to 9. The argument Bin indicates that the variable x is binary.
@variable(sudoku_milp, x[i = 1:9, j = 1:9, k = 1:9], Bin);
The array x consists of:
i(row index): represents the row number in the sudoku grid, in the range 1 to 9j(column index): represents the number of a column in the sudoku grid, also in the range 1 to 9k(digit index): represents a potential digit (from 1 to 9) that can be placed in a cell (i,j)
Thus, the binary representation of x means that:
-
- the digit 8 is placed in the cell located in row 3, column 3
-
= 0 - the digit 8 is not placed in this cell
Using the loop for to loop through the rows i and columns j, create the following condition for the optimisation problem:
This condition ensures that each of the cells in the sudoku grid contains exactly one number between 1 and 9.
for i in 1:9
for j in 1:9
@constraint(sudoku_milp, sum(x[i, j, k] for k in 1:9) == 1)
end
end
Using loops for to loop through rows i and columns j, create the following conditions for the optimisation problem:
This constraint ensures that each digit of k appears exactly once in each row of ind.
This constraint ensures that each digit k appears exactly once in each column ind.
for ind in 1:9
for k in 1:9
@constraint(sudoku_milp, sum(x[ind, j, k] for j in 1:9) == 1)
@constraint(sudoku_milp, sum(x[i, ind, k] for i in 1:9) == 1)
end
end
The code in the above cell adds 162 constraints to the model:
- 81 constraints for rows (9 rows × 9 digits)
- 81 column constraints (9 columns × 9 digits)
With these two conditions, we have fulfilled the fulfilment of the first two sudoku rules:
- Each row must contain digits from 1 to 9 without repetition
- Each column must contain digits from 1 to 9 without repetition
Using loops for, create the following condition for the optimisation problem:
The outermost loop for i at 1:3:7 searches rows 1, 4, and 7, which are the initial rows of each 3x3 subnet.
The middle loop for j at 1:3:7 loops through columns 1, 4, and 7, which are the initial columns of each 3x3 subnet.
The innermost loop for k at 1:9 searches all possible digits (from 1 to 9).
The variables r and c allow you to define constraints for each 3x3 subnet in sudoku.
r- represents the row indices in the 3x3 subgrid and takes values fromitoi+2c- represents column indices in a 3x3 subnet and takes values fromjtoj+2
for i in 1:3:7
for j in 1:3:7
for k in 1:9
@constraint(sudoku_milp, sum(x[r, c, k] for r in i:(i+2), c in j:(j+2)) == 1)
end
end
end
The code in the above cell adds 81 constraints (9 digits for each of the 9 subgrids) to the model, ensuring that each 3x3 subgrid contains all digits from 1 to 9 without repetition. This is how we fulfil the third rule of Sudoku puzzles.
Using the loops for to loop through the cells, put the values 1 (true) in the binary variable x in the cells that have the initial clues (sudoku_puzzle) using the function fix():
for i in 1:9
for j in 1:9
if sudoku_puzzle[i, j] != 0
fix(x[i, j, sudoku_puzzle[i, j]], 1; force = true)
end
end
end
You have created conditions for the fulfilment of all three puzzle solving rules and entered the clues into the variable x. Now you can solve the optimisation problem:
optimize!(sudoku_milp)
Output the status of the solver's result:
status = termination_status(sudoku_milp)
println("Статус: ", status)
if status == MOI.OPTIMAL || status == MOI.FEASIBLE_POINT
println("Оптимальное решение найдено")
else
println("Оптимальное решение не найдено")
end
The OPTIMAL status means that the solver has found the global optimal solution to the problem.
Save the values of the binary variable x in the variable x_val:
x_val = value.(x)
You can make sure that each value k from the solution has been assigned a grid position in the variable x_val.
Create a new 9x9 matrix filled with zeros that will be used to convert the results from binary format to the standard sudoku grid format. The argument Int indicates that the matrix values are integers.
sol_milp = zeros(Int, 9, 9);
Using for loops to loop through all the cells and possible digits in the sudoku grid, convert the values of the three-dimensional array x_val to a two-dimensional array of integers sol_milp.
sol_milp = zeros(Int, 9, 9)
for i in 1:9
for j in 1:9
for k in 1:9
if round(Int, x_val[i, j, k]) == 1
sol_milp[i, j] = k
end
end
end
end
The variable sol_milp now contains the actual solution to the puzzle in the format of a two-dimensional array of integers.
Visualise the results:
plot_sudoku(sol_milp)
We have demonstrated solving the problem using mixed integer linear programming. The solver finds a possible configuration of numbers that satisfies the rules for solving sudoku, formulated as row, column, and subgrid constraints.
Formulation of the programming problem in constraints
You can also model this problem using constraint programming. Let's apply the condition AllDifferent()("all different"), which states that no two elements of the vector can take the same value.
Create an optimisation problem using the function Model() and specify the name of the solver in brackets:
sudoku_cp = Model(HiGHS.Optimizer)
Create a variable sudoku_cp, containing a two-dimensional array . Each variable ranges from 1 to 9. The argument Int indicates that the values of the variables are integers.
@variable(sudoku_cp , 1 <= x[1:9, 1:9] <= 9, Int);
Set the condition that applies to each row i of the sudoku grid. MOI.AllDifferent(9) ensures that all elements in row i are distinct.
This means that each number from 1 to 9 appears exactly once in each row.
@constraint(sudoku_cp , [i = 1:9], x[i, :] in MOI.AllDifferent(9));
Set the condition that applies to each column j of the sudoku grid. MOI.AllDifferent(9) guarantees that all elements in column j are distinct.
This means that every number from 1 to 9 appears exactly once in each column.
@constraint(sudoku_cp, [j = 1:9], x[:, j] in MOI.AllDifferent(9));
Set the condition that applies to each 3x3 subgrid in sudoku. MOI.AllDifferent(9) ensures that all elements in each subgrid are distinct. This means that every number from 1 to 9 appears exactly once in each 3x3 subgrid.
for i in 1:3:7, j in 1:3:7
@constraint(sudoku_cp, vec(x[i:i+2, j:j+2]) in MOI.AllDifferent(9))
end
For each of the cells in [i, j], transfer the original clues from sudoku_puzzle to the variable x using the function fix:
for i in 1:9, j in 1:9
if sudoku_puzzle[i, j] != 0
fix(x[i, j], sudoku_puzzle[i, j]; force = true)
end
end
You have created the fulfilment conditions for all three puzzle solving rules and put the clues into the variable x. Now you can solve the optimisation problem:
optimize!(sudoku_cp)
Output the status of the solver's result:
status = termination_status(sudoku_cp)
println("Статус: ", status)
if status == MOI.OPTIMAL || status == MOI.FEASIBLE_POINT
println("Оптимальное решение найдено")
else
println("Оптимальное решение не найдено")
end
The OPTIMAL status means that the solver has found the global optimal solution to the problem.
Save the results in a variable. The obtained values are of type Float64. You can use the function round together with the argument Int, to convert the values to integers.
sol_cp = round.(Int, value.(x));
Visualise the resulting sudoku solution:
plot_sudoku(sol_cp)
You can compare the solutions obtained by using mixed integer linear programming and constraint programming:
sol_milp == sol_cp
The results of both methods of solving the puzzle are the same.
Conclusion
In this example, we solved the Sudoku puzzle using two different approaches - Mixed Integer Linear Programming (MILP) and Constraint Programming (CP). Both approaches can solve the puzzle, but they formulate and solve the problem in fundamentally different ways, demonstrating different characteristics in modelling and problem solving. In the context of solving the problem demonstrated in this example, we can highlight some of the differences between these approaches.
Mixed integer linear programming:
** Uses binary variables
@variable(sudoku_milp, x[i = 1:9, j = 1:9, k = 1:9], Bin)
- Uses linear constraints to enforce sudoku rules
@constraint(sudoku_milp, sum(x[i, j, k] for k at 1:9) == 1)
- Larger model with more variables, but with simpler constraints
` * Generally more robust to small model changes and scales better with increasing problem size
** Constraint programming:**
- Uses integer variables
@variable(sudoku_cp, 1 <= x[1:9, 1:9] <= 9, Int)
- Uses global constraints such as
AllDifferent, to enforce sudoku rules
@constraint(sudoku_cp, [i = 1:9], x[i, :] in MOI.AllDifferent(9))
- Smaller model with fewer variables but more complex constraints
` * May be more sensitive to problem size and structure
In practice, the choice between mixed integer linear programming and constraint programming often depends on the particular characteristics of the problem, and for some problems one method is better suited than the other. In some cases, hybrid approaches combining both methods may be effective.