利用Matlab优化工具箱解数独问题

前一阵使用Matlab的优化工具箱，发现可以求解数独问题，有意思！实际上，Matlab优化工具箱有两个：Optimization Toolbox和Global Optimization Toolbox，前者求解常见的优化问题，后者求解全局最优问题，如常见的方法包括遗传算法、粒子群算法、模拟退火算法等。

    而数独问题则是借助Optimization Toolbox工具箱的 整数线性规划函数intlinprog完成求解。Matlab帮助中有详细介绍，搜索Solve Sudoku Puzzles Via Integer Programming可以看到，这里对帮助做下翻译。

     数独问题就是将9x9的网格填充上1~9的整数值，同时保证每个整数在每一行、每一列和每个3x3方格只出现一次。网格中分布有一些线索，你的任务就是填充剩余的网格数字，如下图。对于此问题，可以使用二元整数优化求解。

    Matlab提供了两个文件SudokuExample.m和sudokuEngine.m，前者用于发布讲解文件，后者为真正求解数独的函数调用，这里将两个文件中的英文注释翻译为中文，如下

    SudokuExample.m文件：

%% 利用整数规划求解数独问题

% 这个例子展示了使用二元整数规划来求解数独问题。

% 数独问题就是将9x9的网格填充上1~9的整数值，

% 同时保证每个整数在每一行、每一列和每个3x3方格只出现一次。

% 网格中分布有一些线索，你的任务就是填充剩余的网格数字。

% 

% 数独有很多手工和自动解法，如

% <http://www.mathworks.com/company/newsletters/articles/solving-sudoku-with-matlab.html Cleve's Corner(2009)>。

% 这个例子展示了一种基于二元整数规划的直接解法。

% 这种解法特别简单，因为不需要复杂的求解流程，而只需要表达清数独的规则，

% 并把线索表示为约束，然后调用函数intlinprog完成求解。

 

%% 初始化网格

clc

clear

close all

 

% 这是一个包含有线索的数据矩阵B，矩阵第一行的(1,2,2)表示在第一行第二列有线索2；

% 第二行的(1,5,3)表示在第一行第五列有线索3。

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)

 

%% 二元整数规划解法

% 求解的关键是将数独问题从9x9方格转化为九层9x9二元矩阵，即9个矩阵的元素值只能为0或1，

% 而求解结果中为1的地方第一层相应位置设为1，为2的地方第二层相应位置为1，依次类推，为9的地方第九层相应位置为1。

% 如此，问题可以通过二元整数规划求解。

%

% 这里不需要目标函数，也可以设为0，因为问题实际上就是找到一个满足所有约束的可行解。

% 然而，为了满足整数规划内部求解器的要求，同时保证求解速度，这里使用一个非常数目标函数。

 

%% 将数独的规则表示为约束

% 假设求解结果 $x$为一个9x9x9的二元数组，那么 $x$应该具有哪些性质呢？

% 首先，数独格中每个格点位置 $(i,j)$都要有一个值，所以在三维9层数组的9个元素 $x(i,j,1),...,x(i,j,9)$中必有一个非零。

% 换句话说，对每个 $i$和 $j$，有

% 

% $\sum_{k=1}^9 x(i,j,k) = 1.$

%

% 类似地，数独格的每一行，只能有一个1~9间的数字，即

% 

% $\sum_{j=1}^9 x(i,j,k) = 1.$

%

% 数独格的每一列有相同的性质，即

% 

% $\sum_{i=1}^9 x(i,j,k) = 1.$

%

% 对3x3的大格也有相似的约束，即对元素 $1\le i\le 3$和 $1\le j\le 3$，对各层 $k$都有

% 

% $\sum_{i=1}^3\sum_{j=1}^3 x(i,j,k) = 1.$

%

% 为了表达所有9个大格的这种属性，可以如下表示

% 

% $\sum_{i=1}^3\sum_{j=1}^3 x(i+U,j+V,k) = 1$

% 

% 其中

% 

% $U,V~\epsilon~\{0,3,6\}.$

 

%% 将线索表达为约束

% 每个初始值，即线索，也可以表达为约束。假设 $(i,j)$位置的线索为 $m$，其中 $1\le m\le 9$，

% 则有 $x(i,j,m) = 1$。而约束 $\sum_{k=1}^9 x(i,j,k) = 1$能保证对其他 $k\neq m$， $x(i,j,k) = 0$。

 

%% 数独规则的数学表达

% 尽管数独的规则可以方便地由9x9x9的x表达，但是线性约束通常需要表达为x(:)的向量形式。

% 因此，当你表达一个数独问题时，需要由9x9x9的初始数组来得到约束矩阵。

% 函数sudokuEngine中提供了一种配置数独规则和线索的的方式。

 

%% 数独求解

tic

S = sudokuEngine(B);

toc

drawSudoku(S)

    sudokuEngine.m文件：

function [S,eflag] = sudokuEngine(B)

% 求解数独问题。函数读取矩阵B中的初始数据，然后调用intlinprog求解，最后返回结果到S中。

% 矩阵B为3列、至少17行的矩阵（因为数独问题至少要有17个初始元素才能保证求解的唯一性）。

% 其中，前两列是线索的i、j坐标，第三列为线索值。

% 如果B是一个9x9矩阵，函数会首先将其转化为3列矩阵的形式。

 

% 判断线索的存储格式

if isequal(size(B),[9,9]) % 如果是9x9格式

    % 转化为81x3

    [SM,SN] = meshgrid(1:9);

    B = [SN(:),SM(:),B(:)];

    

    % 删除线索为0的行

    [rrem,~] = find(B(:,3) == 0);

    B(rrem,:) = [];

end

 

% 判断B是否满足所需格式要求

if size(B,2) ~= 3 || length(size(B)) > 2

    error('输入矩阵必须为Nx3或9x9')

end

 

% 判断元素是否为1~9的整数

if sum([any(B ~= round(B)),any(B < 1),any(B > 9)])

    error('元素必须为1~9的整数')

end

 

%% 数独的规则

N = 9^3; % 求解结果x中独立变量的个数，x为一个9^3长度的列向量

M = 4*9^2; % 约束的个数

Aeq = zeros(M,N); % 等式约束矩阵Aeq*x = beq

beq = ones(M,1); % 常向量beq

f = (1:N)'; % 目标函数可以是任意值，但是使用非常数可以加快求解速度

lb = zeros(9,9,9); % 各元素下限

ub = lb+1; % 各元素上限，上下限如此配置可以保证结果为二元规划

 

counter = 1; % 指向约束次序的指针

% 约束：第k层第j列的和为1

for j = 1:9

    for k = 1:9

        Astuff = lb; % Astuff清零

        Astuff(1:end,j,k) = 1;

        Aeq(counter,:) = Astuff(:)'; % 将Astuff存为Aeq的一行

        counter = counter + 1;

    end

end

% 约束：第k层第i行的和为1

for i = 1:9

    for k = 1:9

        Astuff = lb;

        Astuff(i,1:end,k) = 1;

        Aeq(counter,:) = Astuff(:)';

        counter = counter + 1;

    end

end

% 约束：第k层每个3x3格的和为1

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

% 约束：第i行第j列的各层和为1

for i = 1:9

    for j = 1:9

        Astuff = lb;

        Astuff(i,j,1:end) = 1;

        Aeq(counter,:) = Astuff(:)';

        counter = counter + 1;

    end

end

 

%% 将数独线索表达为约束

% 将初始线索植入lb数组，即有线索的地方下限设为1，这就保证了结果满足x(i,j,k) = 1。

for i = 1:size(B,1)

    lb(B(i,1),B(i,2),B(i,3)) = 1;

end

 

%% 数独求解

% 数独问题的数学表达：规则由Aeq和beq矩阵表达，线索由lb中的1表达。然后调用intlinprog函数求解。

% 为了保证问题的解为二元解，设置intcon参数为1:N（即所有变量为整数），同时他们的上下限为0和1。

intcon = 1:N;

[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub); % 虽然lb，ub是三层数组的形式，但intlinprog第一步会将它们转化为lb(:)和ub(:)的向量格式

 

%% 将结果转化为可用格式

% 从x得到数独网格只需叠加$(i,j)$位置元素乘其所在的层：

if eflag > 0 % 线性规划顺利完成

    x = reshape(x,9,9,9); % 将x转化为9x9x9的数组

    x = round(x); % 四舍五入非整数解

    y = ones(size(x));

    for k = 2:9

        y(:,:,k) = k; % 保存每层层数

    end

    

    S = x.*y; % 各元素乘以各层层数

    S = sum(S,3); % 累加得到9x9的结果S

else

    S = [];

end

    对于SudokuExample.m，运行Publish，可以得到如下html文件结果：

from: http://chunqiu.blog.ustc.edu.cn/?p=820