이 페이지의 최신 내용은 아직 번역되지 않았습니다. 최신 내용은 영문으로 볼 수 있습니다.

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

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

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

초기 퍼즐

다음은 단서로 구성된 데이터 행렬 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에서 해를 구해줍니다.

이진 정수 계획법

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

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

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

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

x가 9x9x9 이진 배열로 표현되어 있다고 가정하겠습니다. 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.

3x3 주 그리드는 유사한 제약 조건을 가집니다. 그리드 요소 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이 되도록 합니다.

최적화 문제 형식의 스도쿠

이진 변수이며 크기가 9x9x9인 최적화 변수 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)
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