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

この例では、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.

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

このパズルと代替の 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$ について次のようになります。

$\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$ になります。

最適化問題形式の数独

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.
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, 9e+00]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 1e+00]
Presolving model
221 rows, 246 cols, 990 nonzeros  0s
50 rows, 34 cols, 182 nonzeros  0s
29 rows, 14 cols, 108 nonzeros  0s
22 rows, 11 cols, 91 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%   405             405                0.00%        0      0      0         0     0.0s

Solving report
  Status            Optimal
  Primal bound      405
  Dual bound        405
  Gap               0% (tolerance: 0.01%)
  P-D integral      0
  Solution status   feasible
                    405 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.00 (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.

sudsoln = struct with fields:
    x: [9×9×9 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)

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 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