通过整数规划求解数独谜题：基于问题

此示例说明如何使用二元整数规划来求解数独谜题。有关基于求解器的方法，请参阅通过整数规划求解数独谜题：基于求解器。

您可能看到过数独谜题。谜题是用 1 到 9 的整数填充 9×9 网格，要求这九个数字中的每个数字在每一行、每一列和每个九宫格（3×3 粗线正方形）中都只出现一次，不能重复。网格已根据提示进行了部分填充，您的任务是填充网格的其余部分。

初始谜题

这是一个带有提示的数据矩阵 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.

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

2009 年发布的 Cleve's Corner 文章中介绍了此谜题以及一种备选的 MATLAB® 求解方法。

求解数独谜题有许多手动方法，也有许多编程方法。此示例说明一种使用二元整数规划的简单方法。

这种方法特别简单，因为您无需给出求解算法。只需将数独的规则以及每条提示表示为对解的约束，然后 MATLAB 就会生成解。

二元整数规划方法

关键思路是将谜题从 9×9 正方形网格转换为一个由二元值（0 或 1）组成的 9×9×9 立方数组。将这个立方数组想象成由 9 个正方形网络上下堆叠在一起，其中每一层对应于 1 到 9 中的一个整数。顶部网格是数组的第一个正方形层，只要解或提示包含整数 1，则该层对应位置就为 1。只要解或提示包含整数 2，第二层对应位置就为 1。只要解或提示包含整数 9，第九层对应位置就为 1。

这种表示非常适合二元整数规划。

此处不需要目标函数，目标函数也可能是常数项 0。问题实际上只是求可行解，也就是满足所有约束的解。但是，为在整数规划求解器内部打破平衡，同时保证求解速度，这里使用一个非常数目标函数。

将数独规则表示为约束

假设解 $x$ 用一个 9×9×9 的二元数组表示。那 $x$ 应具备哪些特性呢？首先，由于二维网格 (i,j) 的每个方格中都有且仅有一个值，因此三维数组项 $x (i, j, 1), . . ., x (i, j, 9)$ 中有且仅有一个非零元素。换句话说，对于每个 $i$ 和 $j$ ，都满足

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

类似地，在二维网格的每一行 $i$ 中，1 到 9 中的每个数字都只能出现一次。换句话说，对于每个 $i$ 和 $k$ ，都满足

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

二维网格中的每一列 $j$ j 也具有同样的特性：对于每个 $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 .$

要表示全部九个粗线网格，只需将每个 $i$ 和 $j$ 索引加上 3 或 6：

$\sum_{i = 1}^{3} \sum_{j = 1}^{3} x (i + U, j + V, k) = 1,$ 其中 $U, V ϵ {0, 3, 6} .$

将提示表示为约束

每个初始值（提示）都可以表示为一个约束。假设在 $1 \leq m \leq 9$ 时， $(i, j)$ 格点处的提示为 $m$ 。则 $x (i, j, m) = 1$ 。约束 $\sum_{k = 1}^{9} x (i, j, k) = 1$ 能确保在 $k \neq m$ 时，所有其他 $x (i, j, k) = 0$ 。

以优化问题的形式求解数独

创建一个优化变量 x，它是一个大小为 9×9×9 的二元变量。

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.7.1: Copyright (c) 2024 HiGHS under MIT licence terms
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
42 rows, 28 cols, 192 nonzeros  0s
22 rows, 18 cols, 86 nonzeros  0s
Objective function is integral with scale 1

Solving MIP model with:
   22 rows
   18 cols (18 binary, 0 integer, 0 implied int., 0 continuous)
   86 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   372             inf                  inf        0      0      0         0     0.0s
         0       0         0   0.00%   405             inf                  inf        0      0      2         7     0.0s

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

您可以轻松检查解是否正确。

用于绘制数独谜题的函数

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