Solving Sudoku puzzles
This example demonstrates how to find a solution to a Sudoku puzzle using the methods of mixed Integer Linear Programming (MILP) and Constraint Programming (CP).
We will solve this problem using the [problem-oriented optimization] workflow (https://engee.com/community/catalogs/projects/algoritm-resheniia-zadach-optimizatsii ).
In this example, we will use the functions of the [JuMP.jl] library(https://engee.com/helpcenter/stable/ru/julia/JuMP/index.html ) to formulate an optimization problem using the linear optimization library HiGHS.jl, the random number generation library Random.jl and the [Plots.jl] library(https://engee.com/helpcenter/stable/ru/julia/juliaplots/index.html ) to visualize the results.
Installing Libraries
If the latest version of the package is not installed in your environment JuMP, uncomment and run the cell below:
Pkg.add(["JuMP", "HiGHS"])
Task description
Sudoku is a logical combinatorial number placement puzzle. The goal is to fill a 9x9 grid with digits so that each column, each row, and each of the nine 3x3 subgrid contains all the digits from 1 to 9.
Rules:
- Each line must contain numbers from 1 to 9 without repetition.
- Each column must contain numbers from 1 to 9 without repetition.
- Each 3x3 subnet must contain numbers from 1 to 9 without repetition.
In order for the puzzle to be solved, the initial hints are provided.:
.png)
An example of a solved Sudoku:
.png)
Sudoku puzzles range from simple to extremely difficult, with difficulty often determined by the number of initial clues provided and the complexity of the strategies required to solve them. A well-designed Sudoku puzzle should have only one solution.
Connecting libraries
Connect the library JuMP:
using JuMP;
Connect the solver library HiGHS:
using HiGHS;
Connect the charting library Plots:
using Plots;
Connect the random number generation library Random:
using Random;
Creating a Sudoku visualization function
Create a function plot_sudoku(), visualizing 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 visualize Sudoku in this example.
Sudoku Puzzle Generation
In this section, we will write a program that will generate Sudoku puzzles to solve.
Create a function generate_full_grid() 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!() filling a 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(), checking whether the placement of a given number is acceptable (num) in a specific 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 feature is crucial for following all three rules in the process of completing the Sudoku grid.
Create a function remove_numbers!() 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 difficulty levels, 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 the removal of numbers from the grid (in the function remove_numbers!) did not lead to multiple solutions.
Create a function generate_sudoku_puzzle(), which generates a Sudoku with a set number of hints:
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 single—solution Sudoku puzzle is 17.
The program we wrote can stably generate Sudoku with a single solution, where the minimum number of hints is somewhere between 20 and 30, which is more than sufficient for this example.
Generate a Sudoku with 25 hints:
sudoku_puzzle = generate_sudoku_puzzle(25);
Visualize Sudoku with hints:
plot_sudoku(sudoku_puzzle)
Formulation of the mixed integer linear programming problem
Let's 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 optimization task using the function Model() and specify the name of the solver in parentheses.:
sudoku_milp = Model(HiGHS.Optimizer)
Create a variable x, which is a three - dimensional array , where i, j and k are in the range from 1 to 9. The argument Bin indicates that the variable x it is binary.
@variable(sudoku_milp, x[i = 1:9, j = 1:9, k = 1:9], Bin);
Array x consists of:
i(Row index): Represents the row number in the Sudoku grid in the range from 1 to 9j(Column index): Represents the column number in the Sudoku grid, also in the range from 1 to 9k(digit index): Represents a potential digit (from 1 to 9) that can be placed in a cell (i,j)
So the binary representation is x means that:
-
— the number 8 is placed in the cell located in row 3, column 3
-
= 0 — the digit 8 is not placed in this cell.
Using cycles for to iterate through the lines i and columns j create the following condition for the optimization problem:
This condition ensures that each of the cells in the Sudoku grid contains exactly one number from 1 to 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 cycles for to iterate through the lines i and columns j, create the following conditions for the optimization task:
This constraint ensures that each digit k appears exactly once in each line 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 upstream cell adds constraints to the model 162:
- 81 row limits (9 rows × 9 digits)
- 81 column limits (9 columns × 9 digits)
With these two conditions, we have fulfilled the observance of the first two rules of Sudoku.:
- Each line must contain numbers from 1 to 9 without repetition.
- Each column must contain numbers from 1 to 9 without repetition.
Using cycles for create the following condition for the optimization problem:
The outermost loop for i at 1:3:7, it iterates through rows 1, 4, and 7, which are the starting rows of each 3x3 subnet.
The average cycle for j at 1:3:7, it iterates through columns 1, 4, and 7, which are the starting columns of each 3x3 subnet.
The innermost loop for k at 1:9, he goes through all possible digits (from 1 to 9).
Variables r and c allows you to define constraints for each 3x3 subnet in Sudoku.
r- represents row indexes in a 3x3 subnet and takes values fromibeforei+2c- represents column indexes in a 3x3 subnet and takes values fromjbeforej+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 upstream cell adds 81 constraints to the model (9 digits for each of the 9 subnets), ensuring that each 3x3 subnet contains all the digits from 1 to 9 without repetition. This is how we fulfill the third rule of Sudoku puzzles.
Using cycles for to iterate through the cells, set the values to 1 (true) in a binary variable x in the cells where the initial hints are (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 fulfillment of all three rules for solving puzzles and added hints to the variable x. Now you can solve the optimization problem.:
optimize!(sudoku_milp)
Display the status of the solver result:
game_status = termination_status(sudoku_milp)
println("Статус: ", game_status)
if game_status == MOI.OPTIMAL || game_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, a position in the grid was assigned to a variable x_val.
Create a new 9x9 matrix filled with zeros, which will be used to convert the results from binary format to the standard Sudoku grid format. Argument Int indicates that the values of the matrix are integers.
sol_milp = zeros(Int, 9, 9);
Using cycles for to iterate through all the cells and possible numbers in the Sudoku grid, transform the values of the three-dimensional array x_val into 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
Now the variable sol_milp It contains the actual solution to the puzzle in the format of a two-dimensional array of integers.
Visualize the results:
plot_sudoku(sol_milp)
We have demonstrated the solution of the problem using mixed integer linear programming. The solver finds a possible number configuration that conforms to the Sudoku solving rules formulated as row, column, and subnet constraints.
Formulation of the programming problem in constraints
You can also simulate this task using constraint programming. Let's apply the condition AllDifferent()("all are different"), which states that no two elements of a vector can take the same value.
Create an optimization task using the function Model() and specify the name of the solver in parentheses.:
sudoku_cp = Model(HiGHS.Optimizer)
Create a variable sudoku_cp, containing a two-dimensional array . Each variable ranges from 1 to 9. The argument is Int specifies that the values of the variables are integers.
@variable(sudoku_cp , 1 <= x[1:9, 1:9] <= 9, Int);
Specify the condition applicable to each row i sudoku grids. MOI.AllDifferent(9) ensures that all the elements in the row are i they are different.
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));
Specify the condition applicable to each column j sudoku grids. MOI.AllDifferent(9) ensures that all the items in the column j they are different.
This means that each 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 applicable to each 3x3 subnet in Sudoku. MOI.AllDifferent(9) ensures that all the elements in each of the subnets are different. This means that each number from 1 to 9 appears exactly once in each 3x3 subnet.
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 [i, j] transfer the original hints from sudoku_puzzle to a 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 conditions for the fulfillment of all three rules for solving puzzles and added hints to the variable x. Now you can solve the optimization problem.:
optimize!(sudoku_cp)
Display the status of the solver result:
game_status = termination_status(sudoku_cp)
println("Статус: ", game_status)
if game_status == MOI.OPTIMAL || game_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 values received are of the type Float64. You can use the function round along with the argument Int to convert values to integers.
sol_cp = round.(Int, value.(x));
Visualize the resulting Sudoku solution:
plot_sudoku(sol_cp)
You can compare the solutions obtained by using mixed integer linear programming and by programming in constraints.:
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 modeling and problem solving. As part of the solution to the problem demonstrated in this example, we can identify some 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 in 1:9) == 1)
A larger model with more variables, but with simpler constraints
- It is generally more resilient to small model changes and scales better with increasing task size
Programming in constraints:
- Uses integer variables
@variable(sudoku_cp, 1 <= x[1:9, 1:9] <= 9, Int)
- Uses global constraints such as
AllDifferent, to ensure compliance with the rules of Sudoku
@constraint(sudoku_cp, [i= 1:9], x[i, :] in MOI.AllDifferent(9))
- A smaller model with fewer variables, but with more complex constraints
- May be more sensitive to the size and structure of the task.
In practice, the choice between mixed integer linear programming and constraint programming often depends on the specific characteristics of the task, and for some tasks one method is better suited than the other. In some cases, hybrid approaches combining both methods may be effective.