跳转至

数独

解决方法

传统的解决方法是暴力搜索，即以 BFS 和 DFS 两种典型算法为主的解决方法。这在网上有很多参考方案。 用 MATLAB 还可以使用内置的 intlinprog 函数利用整数线性规划的思路求解问题，现依次介绍如下。

线性规划

即将数独问题转化为线性规划问题。

Note

我们将数独的9宫格用一3维张量\(N(x,y,z)\)表示 例如: \(N(1,2,3)\) 表示原9宫格的第1行第2列填有数字3

首先翻译约束条件:

原始条件 翻译条件 数学表示
每列1-9有 且仅有1个 每层矩阵的 每列只有一个1 \(\(\sum_{j=1}^{9}N(i,j,k)=1 \space i,k=1,...,9\)\)
每行1-9有 且仅有1个 每层矩阵的 每行只有一个1 \(\(\sum_{i=1}^{9}N(i,j,k)=1\\ j,k=1,...,9\)\)
小9宫格1-9有 且仅有1个 每层矩阵的每一个 9宫格中只出现一个1 $$\sum_{i=1}^{3} \sum_{j=1}^{3}N(i+M,j+S,k)=1 \ k=1,...,9,M,S\in{3,6,9} $$
9层矩阵中同 一个位置仅有一个1 \(\(\sum_{k=1}^{9}N(i,j,k)=1\\ i,j=1,...,9\)\)

而 MATLAB 内置的函数 intlinprog 求解的问题模型如下所示:

\[ \begin{gathered}minf^Tx\\\text{subject to}\left\{\begin{array}{l}A\cdot x<=b\\A_{eq}\cdot x=b_{eq}\\lb<=x<=ub\\x(\text{intcon) are integers}\end{array}\right.\end{gathered} \]

重点求解条件矩阵$ A_{eq} $,

对于条件1,2,3,4分别写出即可:

%% 约束方程
numOfConstraint = 9 * 9 + 3 * 9 * 9;
% 每格只能取一个值 + 每行/列/块每个数字只能出现一次
Aeq = zeros(numOfConstraint, N);
Constraint = 1;

% 每格只能取一个值与每行/列每个数字只能出现一次
for m = 1:9
    for n = 1:9
        % 每格
        layer = 0:81:8 * 81;
        location = m + (n - 1) * 9;
        Aeq(Constraint, location + layer) = 1; 
        Constraint = Constraint + 1;
        % 每行
        y = 1:9;      
        location = (n - 1) * 81 + m;
        Aeq(Constraint, location + (y - 1) * 9) = 1; 
        Constraint = Constraint + 1;
        % 每列
        x = 1:9;
        location = (n - 1) * 81 + (m - 1) * 9;
        Aeq(Constraint, location + x) = 1; 
        Constraint = Constraint + 1;
    end
end
% 每格每个数字只能出现一次 
for x = 1:3
    for y = 1:3
        for z = 1:9
            location = (z - 1) * 81 + (x - 1) * 3 + (y - 1) * 3 * 9;
            shift = [1 2 3 10 11 12 19 20 21];
            Aeq(Constraint, location + shift) = 1;
            Constraint = Constraint + 1;
        end
    end
end
beq = ones(numOfConstraint, 1);

利用 gurobi 求解器求解

如果利用前人写好的求解器求解，则代码会大大简化。 本次介绍 MATLAB+Gurobi+Yalmip

前置

安装 Yalmip 和 Gurobi YALMIP是MATLAB下的一个工具箱(app), 安装之后即可通过它使用Gurobi之类的求解器(solver),求解许多问题。尽管MATLAB内置的intlinprog等函数已经非常好用，但是对于大型模型的求解仍然需要使用Gurobi这样的商业求解器。而YALMIP则使代码更贴近于数学语言，即不必写出条件矩阵，而是直接输入数学语言即可求解，大大方便了模型建立。

详细的安装教程在此:MATLAB+Gurobi+Yalmip安装和使用教程

接下来看一些运用YALMIP求解优化问题的例子

用YALMIP解问题的格式如下

  • 创建决策变量 = 目标函数 z
  • 约束条件设置 C
  • 参数设置
  ops = sdpsetting('solver','Cplex','verbose',0);% verbose:显示冗余度 0为只显示结果
  • 求解result = solvesdp(C,z,ops)

Note

有 \(n\) 项任务，由 \(n\) 个人来完成，每个人只能做一件，第 \(i\) 个人完成第 \(j\) 项任务要 \(c_{ij}\) 小时 ，如何合理安排时间才能使总用时最小？

用 \(x_{ij}\) 表示第 \(i\) 个人完成第 \(j\) 项工作所需要的资源数，称之为价值系数。因此指派问题的数学模型是:

\[ \min z=\sum_{i=1}^n \sum_{j=1}^n c_{ij}x_{ij} \]
\[ s.t= \begin{cases} \sum_{i=1}^n x_{ij}=1,\quad i=1,2,\cdots,n\\ \sum_{j=1}^n x_{ij}=1,\quad j=1,2,\cdots,n\\ x_{ij}=0\text{或}1,\quad i,j=1,2,\cdots,n \end{cases} \]

例如

人任务耗时 A B C D E
甲 12 7 9 7 9
乙 8 9 6 6 6
丙 7 17 12 14 9
丁 15 14 6 6 10
戊 4 10 7 10 9 
clear;
clc;
close all;
c=load('data');
%决策变量
x = binvar(5,5,'full');

%目标函数
z = sum(sum(c.*x));

%添加约束
C = [];

C = [C;sum(x,1)==1];   %  1 横向相加
C = [C;sum(x,2)==1];   %  2 纵向相加 

ops=sdpsettings('verbose',0);
%求解
result = optimize(C,z,ops);
if result.problem == 0    %求解成功
    x_star = double(x)
    z_star = double(z)
    else
    disp('求解过程中出错');
end

回到正题 Gurobi官方也提供了数独的求解过程:Sudoku 总的来说，利用优化问题的思路求解要比原有的暴力搜索来得简单。