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

この例では、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 整数計画法を使用した簡単な方法を説明します。

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

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

基本となる考え方は、正方形の 9 行 9 列のグリッドを 2 進数値 (0 または 1) の 3 次元 9 x 9 x 9 の配列に変換することです。3 次元配列を相互に積み重ねられた 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$ について次のようになります。

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

同様に、2 次元グリッドの各行 $i$ には、1 ～ 9 の各数値から厳密に 1 つの値が含まれます。つまり、 $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)$ は、 $1 \leq m \leq 9$ に対して $m$ であると仮定します。すると、 $x (i, j, m) = 1$ になります。制約 $\sum_{k = 1}^{9} x (i, j, k) = 1$ によって、 $k \neq m$ に対して他はすべて $x (i, j, k) = 0$ になります。

数独のルールを書き込む

数独ルールは 9 x 9 x 9 解の配列 x を用いて簡単に記述できますが、線形制約はベクトル解行列 x(:) として与えられます。したがって、数独プログラムを書き込む際は、9 x 9 x 9 初期配列から派生した制約を使用する必要があります。

以下に、制約としてヒントも含む数独ルールを設定する 1 つの方法を示します。sudokuEngine ファイルはソフトウェアに付属しています。

type sudokuEngine

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

数独ソルバーの呼び出し

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

Running HiGHS 1.6.0: Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
223 rows, 248 cols, 996 nonzeros
46 rows, 31 cols, 182 nonzeros
25 rows, 14 cols, 125 nonzeros
4 rows, 6 cols, 15 nonzeros
3 rows, 5 cols, 11 nonzeros
0 rows, 0 cols, 0 nonzeros
Presolve: Optimal

Solving report
  Status            Optimal
  Primal bound      29565
  Dual bound        29565
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    29565 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.01 (total)
                    0.01 (presolve)
                    0.00 (postsolve)
  Nodes             0
  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)

解が正しいかどうかを容易に確認できます。

数独パズルを描画する関数

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