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

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

您可能看到过数独谜题。谜题是用 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® 求解方法。

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

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

二元整数规划方法

关键思路是将谜题从 9×9 正方形网格转换为一个由二元值（0 或 1）组成的 9×9×9 立方数组。将这个立方数组想象成由 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$ 。

编写数独的规则

尽管数独规则可以很方便地用 9×9×9 解数组 x 表示，但线性约束是用向量解矩阵 x(:) 给出的。因此，当您编写数独程序时，必须使用从 9×9×9 初始数组派生的约束矩阵。

以下是一种建立数独规则的方法，也包括作为约束的提示。您的软件随附 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.7.1: Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+00, 7e+02]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 1e+00]
Presolving model
223 rows, 248 cols, 996 nonzeros  0s
44 rows, 30 cols, 185 nonzeros  0s
38 rows, 24 cols, 186 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
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)

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

另请参阅

主题

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