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

In this example, we will use the functions of the JuMP.jl library to formulate the optimisation problem, the linear optimisation library HiGHS.jl, the random number generation library Random.jl and the Plots.jl library to visualise the results.

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:

1. Each line must contain digits from 1 to 9 without repetition.
2. Each column must contain digits from 1 to 9 without repetition.
3. Each 3x3 subnet must contain digits from 1 to 9 without repetition.

Initial hints are provided to make solving the puzzle possible:

An example of a solved sudoku:

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

plot_sudoku (generic function with 1 method)

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

generate_full_grid (generic function with 1 method)

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

fill_grid! (generic function with 1 method)

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

is_valid (generic function with 1 method)

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

remove_numbers! (generic function with 1 method)

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

count_solutions (generic function with 2 methods)

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

generate_sudoku_puzzle (generic function with 1 method)

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)

A JuMP Model
├ solver: HiGHS
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

Create a variable x, representing a three-dimensional array $x(i, j, k)$, 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 9
• j (column index): represents the number of a column in the sudoku grid, also in the range 1 to 9
• k (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:

• $x(3, 3, 8) = 1$ - the digit 8 is placed in the cell located in row 3, column 3

• $x(3, 3, 8)$ = 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:

$$\sum_{k=1}^{9} x_{i, j, k} = 1$$

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: $$\sum_{j=1}^{9} x_{ind, j, k} = 1$$

This constraint ensures that each digit of k appears exactly once in each row of ind.

$$\sum_{i=1}^{9} x_{i, ind, k} = 1$$

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:

$$\sum_{r=i}^{i+2} \sum_{c=j}^{j+2} x_{r, c, k} = 1$$

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 from i to i+2
• c - represents column indices in a 3x3 subnet and takes values from j to j+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)

Running HiGHS 1.7.2 (git hash: 5ce7a2753): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [0e+00, 0e+00]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 1e+00]
Presolving model
200 rows, 184 cols, 736 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve: Optimal

Solving report
  Status            Optimal
  Primal bound      0
  Dual bound        0
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.00 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             0
  LP iterations     0 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)

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

Статус: OPTIMAL
Оптимальное решение найдено

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)

9×9×9 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0

[:, :, 2] =
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0

[:, :, 3] =
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0

[:, :, 4] =
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0

[:, :, 5] =
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0

[:, :, 6] =
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0

[:, :, 7] =
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

[:, :, 8] =
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0

[:, :, 9] =
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

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)

A JuMP Model
├ solver: HiGHS
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

Create a variable sudoku_cp, containing a two-dimensional array $x(i, j)$. 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)

Running HiGHS 1.7.2 (git hash: 5ce7a2753): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 9e+00]
  Cost   [0e+00, 0e+00]
  Bound  [1e+00, 9e+00]
  RHS    [1e+00, 1e+00]
Presolving model
747 rows, 1568 cols, 5788 nonzeros  0s
504 rows, 1568 cols, 4108 nonzeros  0s
Objective function is integral with scale 1

Solving MIP model with:
   504 rows
   1568 cols (1512 binary, 56 integer, 0 implied int., 0 continuous)
   4108 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0               inf                  inf        0      0      0         0     0.0s
         0       0         0   0.00%   0               inf                  inf        0      0      3       831     0.0s

10.4% inactive integer columns, restarting
Model after restart has 302 rows, 570 cols (514 bin., 56 int., 0 impl., 0 cont.), and 2281 nonzeros

         0       0         0   0.00%   0               inf                  inf       13      0      0      8462     1.7s

Solving report
  Status            Optimal
  Primal bound      0
  Dual bound        0
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0 (objective)
                    0 (bound viol.)
                    1.34559030585e-13 (int. viol.)
                    0 (row viol.)
  Timing            1.77 (total)
                    0.08 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     8783 (total)
                    0 (strong br.)
                    1996 (separation)
                    5635 (heuristics)

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

Статус: OPTIMAL
Оптимальное решение найдено

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

true

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 $x(i, j, k)$

@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 $x(i, j)$

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