Main Content

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

이 예제에서는 이진 정수 계획법을 사용하여 스도쿠 퍼즐을 푸는 방법을 보여줍니다. 솔버 기반 접근법은 Solve Sudoku Puzzles Via Integer Programming: Solver-Based 항목을 참조하십시오.

스도쿠 퍼즐을 본 적이 있을 것입니다. 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.

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

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

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

이진 정수 계획법

핵심 발상은 정사각 9×9 그리드에서 이진 값(0 또는 1)으로 구성된 3차원 9×9×9 배열로 퍼즐을 변환하는 것입니다. 이 3차원 배열을 서로 쌓여있는 9개의 정사각 그리드로 생각하십시오. 여기서 각 층은 1부터 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이 아닌 요소가 정확히 하나 있습니다. 즉, 모든 ij에 대해 다음이 충족됩니다.

k=19x(i,j,k)=1.

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

j=19x(i,j,k)=1.

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

i=19x(i,j,k)=1.

3×3 주 그리드는 유사한 제약 조건을 가집니다. 그리드 요소 1i31j31k9 각각에 대해 다음이 충족됩니다.

i=13j=13x(i,j,k)=1.

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

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

단서 표현하기

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

최적화 문제 형식의 스도쿠

이진 변수이며 크기가 9×9×9인 최적화 변수 x를 생성합니다.

x = optimvar('x',9,9,9,'Type','integer','LowerBound',0,'UpperBound',1);

임의의 목적 함수를 사용하여 최적화 문제를 생성합니다. 목적 함수는 문제의 내재된 대칭성을 제거하여 솔버에 도움을 줄 수 있습니다.

sudpuzzle = optimproblem;
mul = ones(1,1,9);
mul = cumsum(mul,3);
sudpuzzle.Objective = sum(sum(sum(x,1),2).*mul);

각 좌표 방향에서 x의 합이 1이 되도록 제약 조건을 나타냅니다.

sudpuzzle.Constraints.consx = sum(x,1) == 1;
sudpuzzle.Constraints.consy = sum(x,2) == 1;
sudpuzzle.Constraints.consz = sum(x,3) == 1;

주 그리드의 합도 1이 되도록 제약 조건을 생성합니다.

majorg = optimconstr(3,3,9);

for u = 1:3
    for v = 1:3
        arr = x(3*(u-1)+1:3*(u-1)+3,3*(v-1)+1:3*(v-1)+3,:);
        majorg(u,v,:) = sum(sum(arr,1),2) == ones(1,1,9);
    end
end
sudpuzzle.Constraints.majorg = majorg;

단서 요소의 하한을 1로 설정하여 초기 단서를 포함시킵니다. 이렇게 설정하면 대응되는 요소의 값이 1로 고정되며 단서로 지정된 각 값의 해가 단서 항목으로 설정됩니다.

for u = 1:size(B,1)
    x.LowerBound(B(u,1),B(u,2),B(u,3)) = 1;
end

스도쿠 퍼즐을 풉니다.

sudsoln = solve(sudpuzzle)
Solving problem using intlinprog.
LP:                Optimal objective value is 405.000000.                                           


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap
tolerance of the optimal value, options.AbsoluteGapTolerance = 0 (the default
value). The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05 (the default value).
sudsoln = struct with fields:
    x: [9x9x9 double]

모든 요소가 정수가 되도록 해를 반올림한 후 해를 표시합니다.

sudsoln.x = round(sudsoln.x);

y = ones(size(sudsoln.x));
for k = 2:9
    y(:,:,k) = k; % multiplier for each depth k
end
S = sudsoln.x.*y; % multiply each entry by its depth
S = sum(S,3); % S is 9-by-9 and holds the solved puzzle
drawSudoku(S)

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

스도쿠 퍼즐을 그리는 함수

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

%   Copyright 2014 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

관련 항목