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

1.主体思想

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

2.数独约束

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

• 9层矩阵的同一位置只可以出现一个1
  $$\sum\nolimits_{k=1}^{9}x(i,j,k)=1$$
• 同理，任何一层的任何一行只能有一个1
  $$\sum\nolimits_{i=1}^{9}x(i,j,k)=1$$
• 同理，任何一层的任何一列只能有一个1
  $$\sum\nolimits_{j=1}^{9}x(i,j,k)=1$$
• 对于$3\times3$的小单位有类似的约束
  $$\sum\nolimits_{i=1}^{3}\sum\nolimits_{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对应位置的系数，公式如下：
  $$Aeq*x=Beq\\ 其中 \min\ f^Tx\ {\rm subject\ to} \left\{\begin{matrix} \\x{\rm (intcon)\ are\ integers} \\A\cdot x\leq b \\Aeq\cdot x=beq \\lb \leq x \leq ub \end{matrix}\right.$$

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.结果图示