Main Content

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

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

您可能看到过数独谜题。谜题是用 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.

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) 中有且仅有一个非零元素。换句话说,对于每个 ij,都满足

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

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

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

二维网格中的每一列 jj 也具有同样的特性:对于每个 jj 和 k,都满足

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

3×3 粗线网格具有类似的约束。对于 1i31j3 的网格元素,对于每层 1k9 都满足:

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

要表示全部九个粗线网格,只需将每个 ij 索引加上 3 或 6:

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

将提示表示为约束

每个初始值(提示)都可以表示为一个约束。假设在 1m9 时,(i,j) 格点处的提示为 m。则 x(i,j,m)=1。约束 k=19x(i,j,k)=1 能确保在 km 时,所有其他 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.6.0: Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
221 rows, 246 cols, 990 nonzeros
45 rows, 30 cols, 196 nonzeros
27 rows, 16 cols, 114 nonzeros
0 rows, 0 cols, 0 nonzeros
Presolve: Optimal

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.00 (total)
                    0.00 (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.
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

相关主题