最初的难题

这是一个数据矩阵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)有关此程序的列表，请参见本例的末尾。

这个难题，和替代MATLAB®解决技术，在克里夫的角落在2009年。

有许多手动解决数独谜题的方法，以及许多程序化方法。这个例子展示了一种使用二进制整数编程的简单方法。

这种方法特别简单，因为不需要给出求解算法。只要把数独的规则表达出来，把线索表达成解的约束，然后MATLAB就能生成解。

二进制整数规划方法

关键思想是将一个谜题从一个9 × 9的正方形网格转换为一个9 × 9 × 9的二进制值(0或1)的立方数组。将立方数组想象成9个相互堆叠的正方形网格，其中每一层对应从1到9的整数。最上面的网格，数组的一个方形层，在解决方案或线索有1的地方都有1。第二层是1，只要答案或线索是2。第九层是1，只要答案或线索是9。

这个公式精确地适用于二进制整数规划。

这里不需要目标函数，它也可以是常数项0。问题其实就是找到一个可行的解，也就是满足所有约束条件的解。然而，对于整数规划求解器内部的打结问题，使用非常数目标函数来提高求解速度。

将数独游戏规则表达为约束条件

假设一个解决方案 $$x$$用9 × 9 × 9二进制数组表示。什么属性 $$x$$有什么?首先，二维网格(i,j)中的每个正方形只有一个值，因此在三维数组条目中只有一个非零元素 $$x（我，j，1），．．．，x（我，j，9）$$．换句话说，对于每一个 $$我$$而且 $$j$$，

类似地，在每一行中 $$我$$在二维网格中，从1到9的每一个数字中只有一个值。换句话说，对于每一个 $$我$$而且 $$k$$，

每一列 $$j$$在二维网格中具有相同的属性:对于每一个 $$j$$而且 $$k$$，

主要的3 × 3网格也有类似的限制。对于网格元素 $$1\le 我\le 3.$$而且 $$1\le j\le 3.$$，对于每个 $$1\le k\le 9$$，

要表示所有9个主要网格，只需在每个网格上添加3或6 $$我$$而且 $$j$$指数:

$$\underset{我＝1}{\overset{3.}{\sum }}\underset{j＝1}{\overset{3.}{\sum }}x（我+U，j+V，k）＝1，$$在哪里 $$U，Vϵ｛0，3.，6｝．$$

表达的线索

每个初始值(线索)都可以表示为一个约束。假设 $$（我，j）$$线索是 $$米$$对于一些 $$1\le 米\le 9$$．然后 $$x（我，j，米）＝1$$．约束 $$\underset{k＝1}{\overset{9}{\sum }}x（我，j，k）＝1$$确保所有其他 $$x（我，j，k）＝0$$为 $$k\ne 米$$．

优化问题形式的数独

创建优化变量x它是二进制的，大小为9 × 9 × 9。

X = optimvar(“x”9日,9日,9日“类型”，“整数”，下界的0,“UpperBound”1);

创建一个具有任意目标函数的优化问题。目标函数可以通过破坏问题的固有对称性来帮助求解者。

Sudpuzzle =优化问题;Mul = ones(1,1,9);Mul = cumsum(Mul,3);sudpuzzle。目标= sum(sum(sum(x,1)，2).*mul);

表示约束条件x每个坐标方向为1。

sudpuzzle. constraints . conx = sum(x,1) == 1;sudpuzzle.Constraints.consy = sum(x,2) == 1;sudpuzzle. constraints . conz = sum(x,3) == 1;

创建约束条件，使主要网格的和也等于1。

Majorg = optimconstr(3,3,9);为U = 1:3为v = 1:3 arr = x (3 * (u-1) + 1:3 * (u-1) + 3, 3 *(它们)+ 1:3 *(它们)+ 3:);maorg (u,v，:) = sum(sum(arr,1)，2) == ones(1,1,9);结束结束sudpuzzle.Constraints.majorg = majorg;

通过在线索条目处设置下界1来包含初始线索。该设置将对应条目的值固定为1，因此将每个线索值处的解设置为线索条目。

为u = 1:尺寸(B, 1) x.LowerBound (B (u, 1), B (u, 2), B (u, 3)) = 1;结束

解决数独难题。

Sudsoln = solve(sudpuzzle)

使用intlinprog解决问题。LP:最优目标值为405.000000。找到最优解。Intlinprog停止在根节点，因为目标值在最佳值的间隙公差范围内，选项。AbsoluteGapTolerance = 0(默认值)。intcon变量是公差范围内的整数，选项。IntegerTolerance = 1e-05(默认值)。

sudsoln =带字段的结构:X: [9x9x9 double]

四舍五入解决方案，以确保所有条目都是整数，并显示解决方案。

sudsoln。x＝round(sudsoln.x); y = ones(size(sudsoln.x));为K = 2:9 y(:，:， K) = K;每个深度k的%乘数结束S = sudsoln.x.*y;%将每个条目乘以它的深度S = sum(S,3);% S是9乘9的，保存着已解决的谜题drawSudoku (S)

您可以很容易地检查解决方案是否正确。

绘制数独谜题的函数

类型drawSudoku

function drawSudoku(B) %用于绘制数独板的函数%版权所有公司图;坚持;轴,轴平等%准备画矩形(“位置”,[0 0 9 9],“线宽”,3,“剪裁”,“关闭”)%外边界矩形(“位置”,[3 0 3 9],“线宽”,2)%沉重的竖线矩形(“位置”,[3]0、3、9日,“线宽”,2)%沉重的水平线矩形(“位置”,[1]0 1 9日,“线宽”,1)%小横线矩形(“位置”,[1]0、4、9日,“线宽”,1)矩形(“位置”,[1]0、7、9日,“线宽”,1)矩形(“位置”,[1,0,1,9],“线宽”,1)%小垂直直线矩形('Position'，[4,0,1,9]，'LineWidth'，1)矩形('Position'，[7,0,1,9]，'LineWidth'，1) %填充线索% % B的行为(i,j,k)的形式，其中i是从%开始计数的行，j是列，k是线索。为了将条目放在%框中，j是水平距离，10-i是垂直距离，%我们减去0.5以使线索在框中居中。如果B是一个9乘9的矩阵，如果size(B,2) == 9 % 9列[SM,SN] = meshgrid(1:9)，则先将其转换为3列;% make i,j entries B = [SN(:)，SM(:)，B(:)];% i,j,k行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