利用Matlab二次规划算法玩通数独游戏

程序员文章站 2022-03-04 20:02:40

...

1.主体思想

将 $9 \times 9$ 的数独矩阵转化为 $9 \times 9 \times 9$ 的三维矩阵，三维矩阵中每一个元素的值为0或1。如原数独矩阵的某一位置数字为7，则三维矩阵第7层的对应位置为1，其它层该位置皆为0.
利用二次规划求解该问题时，需要设定约束和目标。此处只需求得一个可行解，因此无需目标方程，但为了满足整数规划求解需求，此处使用一个非定常目标方程。

2.数独约束

假定 $x$ 代表一个 $9 \times 9 \times 9$ 的矩阵。其有以下约束。

9层矩阵的同一位置只可以出现一个1
$\sum_{k = 1}^{9} x (i, j, k) = 1$
同理，任何一层的任何一行只能有一个1
$\sum_{i = 1}^{9} x (i, j, k) = 1$
同理，任何一层的任何一列只能有一个1
$\sum_{j = 1}^{9} x (i, j, k) = 1$
对于 $3 \times 3$ 的小单位有类似的约束
$\sum_{i = 1}^{3} \sum_{j = 1}^{3} x (i + U, j + V, k) = 1 其中， U, V \subset {0, 3, 6}$

3.求解数独

检查格式，将输入统一为（行，列，线索）模式

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

检查线索矩阵中的数据是否在1~9之间，且为整数

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)])
    error('Entries must be integers from 1 to 9')
end

数独矩阵进行约束
数独矩阵在运算中化为长度为 $(9^{3}, 1)$ 的数组，对于每一条约束，只要将约束矩阵相应位置置为1，即**约束矩阵Aeq对应位置的系数，公式如下：
$A e q * x = B e q 其中 min f^{T} x s u b j e c t t o {\begin{matrix} x (i n t c o n) a r e i n t e g e r s \\ A \cdot x \leq b \\ A e q \cdot x = b e q \\ l b \leq x \leq u b \end{matrix}$

N = 9^3; % x中独立变量的个数
M = 4*9^2; % 约束的个数
Aeq = zeros(M,N); % 等式约束矩阵左侧
beq = ones(M,1); % 等式约束矩阵右侧
f = (1:N)'; % 一个目标方程矩阵，无实际意义
lb = zeros(9,9,9); % 元素下限
ub = lb+1; % 元素上限

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

将有线索的位置下限设为1

for i = 1:size(B,1)
    lb(B(i,1),B(i,2),B(i,3)) = 1;
end

求解整数二次规划

intcon = 1:N;
[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub);

整理结果

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); % 合并第三维
else
    S = [];
end

4.绘制图形

绘图准备

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)

填写数字

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

5.结果图示

利用Matlab二次规划算法玩通数独游戏