정수 계획법을 통해 스도쿠 퍼즐 풀기: 솔버 기반

이 예제에서는 이진 정수 계획법을 사용하여 스도쿠 퍼즐을 푸는 방법을 보여줍니다. 문제 기반 접근법에 대해서는 정수 계획법을 통해 스도쿠 퍼즐 풀기: 문제 기반 항목을 참조하십시오.

스도쿠 퍼즐을 본 적이 있을 것입니다. 1부터 9까지의 정수 중 각 정수가 각 행과 열, 그리고 3×3 주 정사각형에서 한 번만 나오도록 정수를 9×9 그리드에 채우는 퍼즐입니다. 그리드에는 단서가 부분적으로 채워져 있으며 그리드의 나머지 부분을 채우는 것이 여러분의 몫입니다.

초기 퍼즐

다음은 단서로 구성된 데이터 행렬 B입니다. 첫 번째 행 B(1,2,2)는 1행, 2열에 단서 2가 있음을 의미합니다. 두 번째 행 B(1,5,3)은 1행, 5열에 단서 3이 있음을 의미합니다. 전체 행렬 B는 다음과 같습니다.

B = [1,2,2;
    1,5,3;
    1,8,4;
    2,1,6;
    2,9,3;
    3,3,4;
    3,7,5;
    4,4,8;
    4,6,6;
    5,1,8;
    5,5,1;
    5,9,6;
    6,4,7;
    6,6,5;
    7,3,7;
    7,7,6;
    8,1,4;
    8,9,8;
    9,2,3;
    9,5,4;
    9,8,2];

drawSudoku(B) % For the listing of this program, see the end of this example.

Figure contains an axes object. The hidden axes object contains 30 objects of type rectangle, text.

이 퍼즐과 대체 MATLAB® 풀이 기법은 2009년도 게시된 Cleve's Corner에서 다루었습니다.

스도쿠 퍼즐을 수동으로 푸는 접근법으로는 여러 가지가 있으며 계획법을 사용한 접근법도 많이 있습니다. 이 예제에서는 이진 정수 계획법을 사용하는 간단한 접근법을 보여줍니다.

해 알고리즘을 제공하지 않으므로 이 접근법은 특히 간단합니다. 스도쿠 규칙을 표현하고 단서를 해에 대한 제약 조건으로 표현하기만 하면 intlinprog가 해를 구해줍니다.

이진 정수 계획법

핵심 발상은 정사각 9×9 그리드에서 이진 값(0 또는 1)으로 구성된 3차원 9×9×9 배열로 퍼즐을 변환하는 것입니다. 이 3차원 배열을 서로 쌓여있는 9개의 정사각 그리드로 생각하십시오. 배열의 정사각 층인 상단 그리드는 해나 단서가 1을 가지는 경우 항상 1을 가집니다. 두 번째 층은 해나 단서가 2를 가지는 경우 항상 1을 가집니다. 9번째 층은 해나 단서가 9를 가지는 경우 항상 1을 가집니다.

이 문제 정식화는 이진 정수 계획법에 적합합니다.

목적 함수가 굳이 필요치 않으며 원하면 0으로 놓아도 됩니다. 이 문제의 목적은 실현 가능한 해, 즉 모든 제약 조건을 충족하는 하나의 해를 구하는 것뿐입니다. 그러나, 정수 계획법 솔버 내부에서 동점에 대한 우선 순위를 결정하고 풀이 속도를 높이려면 상수가 아닌 목적 함수를 사용하십시오.

스도쿠 규칙을 제약 조건으로 표현하기

해 $x$ 가 9×9×9 이진 배열로 표현되어 있다고 가정하겠습니다. $x$ 가 어떤 속성을 가질까요? 먼저, 2차원 그리드 (i,j)의 각 정사각형은 정확히 하나의 값을 가지므로 3차원 배열 요소 $x (i, j, 1), . . ., x (i, j, 9)$ 중에는 0이 아닌 요소가 정확히 하나 있습니다. 즉, 모든 $i$ 와 $j$ 에 대해 다음이 충족됩니다.

$\sum_{k = 1}^{9} x (i, j, k) = 1 .$

마찬가지로, 2차원 그리드의 각 행 $i$ 에는 1부터 9까지의 숫자 중에서 정확히 하나의 값이 있습니다. 다시 말해서 $i$ 와 $k$ 각각에 대해 다음이 충족됩니다.

$\sum_{j = 1}^{9} x (i, j, k) = 1 .$

2차원 그리드의 각 열 $j$ 는 동일한 속성을 가집니다. 즉, $j$ 와 $k$ 각각에 대해 다음이 충족됩니다.

$\sum_{i = 1}^{9} x (i, j, k) = 1 .$

3×3 주 그리드는 유사한 제약 조건을 가집니다. 그리드 요소 $1 \leq i \leq 3$ 및 $1 \leq j \leq 3$ 과 $1 \leq k \leq 9$ 각각에 대해 다음이 충족됩니다.

$\sum_{i = 1}^{3} \sum_{j = 1}^{3} x (i, j, k) = 1 .$

9개의 주 그리드를 모두 표현하려면 $i$ 인덱스와 $j$ 인덱스 각각에 3 또는 6을 추가하기만 하면 됩니다.

$\sum_{i = 1}^{3} \sum_{j = 1}^{3} x (i + U, j + V, k) = 1,$ 여기서 $U, V ϵ {0, 3, 6} .$

단서 표현하기

각각의 초기값(단서)은 제약 조건으로 나타낼 수 있습니다. $(i, j)$ 단서가 $m$ ( $1 \leq m \leq 9$ )이라고 가정하겠습니다. 그러면 $x (i, j, m) = 1$ 입니다. 제약 조건 $\sum_{k = 1}^{9} x (i, j, k) = 1$ 은 $k \neq m$ 라면, 모든 나머지 경우에 대해 $x (i, j, k) = 0$ 이 되도록 합니다.

스도쿠 규칙 만들기

스도쿠 규칙은 9×9×9 해 배열 x로 편리하게 표현되지만, 선형 제약 조건이 벡터 해 행렬 x(:)로 주어집니다. 따라서 스도쿠 프로그램을 작성하는 경우 9×9×9 초기 배열에서 파생된 제약 조건 행렬을 사용해야 합니다.

여기에 스도쿠 규칙을 설정하는 한 가지 접근법이 있으며, 단서가 제약 조건으로 포함됩니다.

function [S,eflag] = sudokuEngine(B)
% This function sets up the rules for Sudoku. It reads in the puzzle
% expressed in matrix B, calls intlinprog to solve the puzzle, and returns
% the solution in matrix S.
%
% The matrix B should have 3 columns and at least 17 rows (because a Sudoku
% puzzle needs at least 17 entries to be uniquely solvable). The first two
% elements in each row are the i,j coordinates of a clue, and the third
% element is the value of the clue, an integer from 1 to 9. If B is a
% 9-by-9 matrix, the function first converts it to 3-column form.

%   Copyright 2014 The MathWorks, Inc. 

if isequal(size(B),[9,9]) % 9-by-9 clues
    % Convert to 81-by-3
    [SM,SN] = meshgrid(1:9); % make i,j entries
    B = [SN(:),SM(:),B(:)]; % i,j,k rows
    % Now delete zero rows
    [rrem,~] = find(B(:,3) == 0);
    B(rrem,:) = [];
end

if size(B,2) ~= 3 || length(size(B)) > 2
    error("The input matrix must be N-by-3 or 9-by-9")
end

if sum([any(B ~= round(B)),any(B < 1),any(B > 9)]) % enforces entries 1-9
    error("Entries must be integers from 1 to 9")
end

%% The rules of Sudoku:
N = 9^3; % number of independent variables in x, a 9-by-9-by-9 array
M = 4*9^2; % number of constraints, see the construction of Aeq
Aeq = zeros(M,N); % allocate equality constraint matrix Aeq*x = beq
beq = ones(M,1); % allocate constant vector beq
f = (1:N)'; % the objective can be anything, but having nonconstant f can speed the solver
lb = zeros(9,9,9); % an initial zero array
ub = lb+1; % upper bound array to give binary variables

counter = 1;
for j = 1:9 % one in each row
    for k = 1:9
        Astuff = lb; % clear Astuff
        Astuff(1:end,j,k) = 1; % one row in Aeq*x = beq
        Aeq(counter,:) = Astuff(:)'; % put Astuff in a row of Aeq
        counter = counter + 1;
    end
end

for i = 1:9 % one in each column
    for k = 1:9
        Astuff = lb;
        Astuff(i,1:end,k) = 1;
        Aeq(counter,:) = Astuff(:)';
        counter = counter + 1;
    end
end

for U = 0:3:6 % one in each square
    for V = 0:3:6
        for k = 1:9
            Astuff = lb;
            Astuff(U+(1:3),V+(1:3),k) = 1;
            Aeq(counter,:) = Astuff(:)';
            counter = counter + 1;
        end
    end
end

for i = 1:9 % one in each depth
    for j = 1:9
        Astuff = lb;
        Astuff(i,j,1:end) = 1;
        Aeq(counter,:) = Astuff(:)';
        counter = counter + 1;
    end
end

%% Put the particular puzzle in the constraints
% Include the initial clues in the |lb| array by setting corresponding
% entries to 1. This forces the solution to have |x(i,j,k) = 1|.

for i = 1:size(B,1)
    lb(B(i,1),B(i,2),B(i,3)) = 1;
end

%% Solve the Puzzle
% The Sudoku problem is complete: the rules are represented in the |Aeq|
% and |beq| matrices, and the clues are ones in the |lb| array. Solve the
% problem by calling |intlinprog|. Ensure that the integer program has all
% binary variables by setting the intcon argument to |1:N|, with lower and
% upper bounds of 0 and 1.

intcon = 1:N;

[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub);

%% Convert the Solution to a Usable Form
% To go from the solution x to a Sudoku grid, simply add up the numbers at
% each $(i,j)$ entry, multiplied by the depth at which the numbers appear:

if eflag > 0 % good solution
    x = reshape(x,9,9,9); % change back to a 9-by-9-by-9 array
    x = round(x); % clean up non-integer solutions
    y = ones(size(x));
    for k = 2:9
        y(:,:,k) = k; % multiplier for each depth k
    end

    S = x.*y; % multiply each entry by its depth
    S = sum(S,3); % S is 9-by-9 and holds the solved puzzle
else
    S = [];
end
end

스도쿠 솔버 호출하기

S = sudokuEngine(B); % Solves the puzzle pictured at the start

Running HiGHS 1.11.0: Copyright (c) 2025 HiGHS under MIT licence terms
MIP  has 324 rows; 729 cols; 2916 nonzeros; 729 integer variables (708 binary)
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+00, 7e+02]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 1e+00]
Presolving model
223 rows, 248 cols, 996 nonzeros  0s
46 rows, 31 cols, 192 nonzeros  0s
38 rows, 25 cols, 163 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve: Optimal

Src: B => Branching; C => Central rounding; F => Feasibility pump; J => Feasibility jump;
     H => Heuristic; L => Sub-MIP; P => Empty MIP; R => Randomized rounding; Z => ZI Round;
     I => Shifting; S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point

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

         0       0         0   0.00%   29565           29565              0.00%        0      0      0         0     0.0s

Solving report
  Status            Optimal
  Primal bound      29565
  Dual bound        29565
  Gap               0% (tolerance: 0.01%)
  P-D integral      0
  Solution status   feasible
                    29565 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.01 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 0
  Nodes             0
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     0 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)

Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 1e-06. The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.

drawSudoku(S)

Figure contains an axes object. The hidden axes object contains 90 objects of type rectangle, text.

해가 올바르다는 것을 손쉽게 확인할 수 있습니다.

스도쿠 퍼즐을 그리는 함수

다음 코드는 drawSudoku 헬퍼 함수를 생성합니다.

function drawSudoku(B)
% Function for drawing the Sudoku board

%   Copyright 2014-2026 The MathWorks, Inc. 


figure;hold on;axis off;axis equal % prepare to draw
rectangle(Position=[0 0 9 9],LineWidth=3,Clipping="off") % outside border
rectangle(Position=[3,0,3,9],LineWidth=2) % heavy vertical lines
rectangle(Position=[0,3,9,3],LineWidth=2) % heavy horizontal lines
rectangle(Position=[0,1,9,1],LineWidth=1) % minor horizontal lines
rectangle(Position=[0,4,9,1],LineWidth=1)
rectangle(Position=[0,7,9,1],LineWidth=1)
rectangle(Position=[1,0,1,9],LineWidth=1) % minor vertical lines
rectangle(Position=[4,0,1,9],LineWidth=1)
rectangle(Position=[7,0,1,9],LineWidth=1)

% Fill in the clues
%
% The rows of B are of the form (i,j,k) where i is the row counting from
% the top, j is the column, and k is the clue. To place the entries in the
% boxes, j is the horizontal distance, 10-i is the vertical distance, and
% we subtract 0.5 to center the clue in the box.
%
% If B is a 9-by-9 matrix, convert it to 3 columns first

if size(B,2) == 9 % 9 columns
    [SM,SN] = meshgrid(1:9); % make i,j entries
    B = [SN(:),SM(:),B(:)]; % i,j,k rows
end

for ii = 1:size(B,1)
    text(B(ii,2)-0.5,9.5-B(ii,1),num2str(B(ii,3)))
end

hold off

end

참고 항목

도움말 항목

정수 계획법을 통해 스도쿠 퍼즐 풀기: 문제 기반