Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

整数計画法により数独パズルを解く: 問題ベース

この例では、0-1 整数計画法を使用して数独パズルを解く方法を説明します。ソルバーベースのアプローチについては、整数計画法により数独パズルを解く: ソルバーベースを参照してください。

おそらく数独パズルを見たことがあるでしょう。各整数が各行、列および主な 3 行 3 列の正方形に 1 回だけ表示されるように 1 ~ 9 の整数を使って 9 行 9 列のグリッドを埋めるパズルです。グリッドには部分的にヒントが表示されており、タスクは残りのグリッドを埋めることです。

初期パズル

ヒントのデータ行列 B を以下に示します。最初の行 B(1,2,2) は 1 行 2 列目にヒント 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 で特集されていました。

数独パズルを手作業で解く多数の方法があるのと同じように、プログラムによって解く多数の方法があります。この例では、0-1 整数計画法を使用した簡単な方法を説明します。

この方法は解のアルゴリズムを指定しないので特に単純です。単に数独のルールを表現し、ヒントを解の制約として記述して、MATLAB が解を求めます。

0-1 整数計画アプローチ

基本となる考え方は、正方形の 9 行 9 列のグリッドを 2 進数値 (0 または 1) の 3 次元 9 x 9 x 9 の配列に変換することです。3 次元配列を相互に積み重ねられた 9 個の正方形グリッドとして考えます。各レイヤーが整数の 1 ~ 9 に対応します。配列の正方形レイヤーである最上位グリッドは、解またはヒントとして 1 が表示されている位置に 1 が表示されています。2 番目のレイヤーは、解またはヒントとして 2 が表示されている位置に 1 が表示されています。9 番目のレイヤーは、解またはヒントとして 9 が表示されている位置に 1 が表示されています。

この定式は 0-1 整数計画法に適しています。

ここでは目的関数は不要であり、定数項 0 である可能性もあります。問題は実行可能解、つまりすべての制約を満たす解を実際に求めることです。ただし、整数計画ソルバーの内部をタイ ブレーキングし、解を求める速度を向上させるために非定数目的関数を使用します。

数独のルールを制約として記述する

x が 9 x 9 x 9 のバイナリ配列で表されると仮定します。x にはどのような特徴があるでしょうか。まず、2 次元グリッド (i,j) 内の各四角形には厳密に 1 つの値が存在するので、3 次元配列のエントリ x(i,j,1),...,x(i,j,9) には厳密に 1 つの非ゼロ要素が存在します。つまり、すべての i および j について次のようになります。

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

同様に、2 次元グリッドの各行 i には、1 ~ 9 の各数値から厳密に 1 つの値が含まれます。つまり、ik のそれぞれについて次のようになります。

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

さらに 2 次元グリッドの各列 j は同じ特徴になり、jk のそれぞれについて次のようになります。

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

主な 3 行 3 列のグリッドは同様の制約をもちます。グリッド要素 1i31j3 および 1k9 のそれぞれについて次のようになります。

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) は、1m9 に対して m であると仮定します。すると、x(i,j,m)=1 になります。制約 k=19x(i,j,k)=1 によって、km に対して他はすべて x(i,j,k)=0 になります。

最適化問題形式の数独

2 値でサイズが 9 x 9 x 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 intcon variables are integer within tolerance, options.IntegerTolerance = 1e-05.
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

関連するトピック