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® 求解方法。

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

这种方法特别简单，因为您无需给出求解算法。只需将数独的规则以及每条提示表示为对解的约束，然后 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，都满足

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

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

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

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

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

3×3 粗线网格具有类似的约束。对于 1≤i≤3 且 1≤j≤3 的网格元素，对于每层 1≤k≤9 都满足：

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

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

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

将提示表示为约束

每个初始值（提示）都可以表示为一个约束。假设在 1≤m≤9 时，(i,j) 格点处的提示为 m。则 x(i,j,m)=1。约束 ∑k=19x(i,j,k)=1 能确保在 k≠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
LP:                Optimal objective value is 29565.000000.                                         


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance, options.IntegerTolerance = 1e-05.

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

相关主题