最初的难题

这是一个数据矩阵Bof clues. The first row,B(1,2,2), means row 1, column 2 has a clue 2. The second row,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)％有关此程序的列表，请参见此示例的结尾。

This puzzle, and an alternative MATLAB® solution technique, was featured inCleve's Cornerin 2009.

有许多方法可以手动解决Sudoku难题以及许多程序化方法。此示例显示了使用二进制整数编程的直接方法。

这种方法特别简单，因为您不提供解决方案算法。只需表达Sudoku的规则，将线索表示为解决方案的约束，然后MATLAB产生解决方案。

Binary Integer Programming Approach

The key idea is to transform a puzzle from a square 9-by-9 grid to a cubic 9-by-9-by-9 array of binary values (0 or 1). Think of the cubic array as being 9 square grids stacked on top of each other, where each layer corresponds to an integer from 1 through 9. The top grid, a square layer of the array, has a 1 wherever the solution or clue has a 1. The second layer has a 1 wherever the solution or clue has a 2. The ninth layer has a 1 wherever the solution or clue has a 9.

This formulation is precisely suited for binary integer programming.

这里不需要目标函数，也可能是一个恒定的项。问题实际上只是为了找到可行的解决方案，这意味着满足所有约束的解决方案。但是，要在整数编程求解器的内部中打破领带，从而提高了解决方案速度，请使用非稳定目标函数。

Express the Rules for Sudoku as Constraints

Suppose a solution $$x$$is represented in a 9-by-9-by-9 binary array. What properties does $$x$$有？首先，2-D网格（I，J）中的每个正方形都有一个值，因此3-D数组条目中恰好有一个非零元素 $$x(i,j,1),...,x(i,j,9)$$. In other words, for every $$i$$和 $$j$$,

Similarly, in each row $$i$$在2-D网格中，每个数字从1到9中都有一个值。换句话说，每个数字 $$i$$和 $$k$$,

And each column $$j$$in the 2-D grid has the same property: for each $$j$$和 $$k$$,

主要的3乘3个网格具有相似的约束。对于网格元素 $$1\le i\le 3$$和 $$1\le j\le 3$$, and for each $$1\le k\le 9$$,

要表示所有九个主要网格，只需向每个添加3或6个 $$i$$和 $$j$$index:

$$\underset{i=1}{\overset{3}{\sum }}\underset{j=1}{\overset{3}{\sum }}x(i+U,j+V,k)=1,$$where $$U,Vϵ\{0,3,6\}.$$

明确线索

每个初始值（线索）可以表示为约束。假设那个 $$(i,j)$$线索是 $$m$$对于一些 $$1\le m\le 9$$. Then $$x(i,j,m)=1$$. The constraint $$\underset{k=1}{\overset{9}{\sum }}x(i,j,k)=1$$确保所有其他 $$x(i,j,k)=0$$for $$k\ne m$$.

Sudoku in Optimization Problem Form

Create an optimization variablexthat is binary and of size 9-by-9-by-9.

x = optimvar('X'，9,9,9，'类型','整数','LowerBound'，0，'UpperBound'，1）;

Create an optimization problem with a rather arbitrary objective function. The objective function can help the solver by destroying the inherent symmetry of the problem.

sudpuzzle = optimproblem; mul = ones(1,1,9); mul = cumsum(mul,3); sudpuzzle.Objective = sum(sum(sum(x,1),2).*mul);

表示约束的总和x在每个坐标方向上是一个。

sudpuzzle.constraints.consx = sum（x，1）== 1;sudpuzzle.constraints.consy = sum（x，2）== 1;sudpuzzle.constraints.consz = sum（x，3）== 1;

创建一个约束，即主要网格的总和也达到一个。

majorg = optimconstr(3,3,9);foru = 1:3forv = 1：3 arr = x（3*（u-1）+1：3*（u-1）+3,3*（v-1）+1：3*（v-1）+3，：：）；majorg（u，v，:) = sum（sum（arr，1），2）==一个（1,1,9）;endendsudpuzzle.constraints.majorg = majorg;

Include the initial clues by setting lower bounds of 1 at the clue entries. This setting fixes the value of the corresponding entry to be 1, and so sets the solution at each clued value to be the clue entry.

foru = 1：size（b，1）x.LowerBound（b（u，1），b（u，2），b（u，3））= 1;end

Solve the Sudoku puzzle.

sudsoln =solve(sudpuzzle)

使用IntlinProg解决问题。LP：最佳目标值为405.000000。找到最佳解决方案。IntlinProg停止在根节点上，因为目标值在最佳值的间隙公差之内。INTCON变量是公差，options.integertolerance = 1E-05（默认值）的整数。

sudsoln =带有字段的结构：x: [9x9x9 double]

围绕解决方案确保所有条目都是整数，并显示解决方案。

sudsoln.x = round(sudsoln.x); y = ones(size(sudsoln.x));fork = 2：9 y（：，：，：，k）= k;% multiplier for each depth kendS = sudsoln.x.*y;% multiply each entry by its depths = sum（s，3）;% S is 9-by-9 and holds the solved puzzle绘制

您可以轻松检查解决方案是否正确。

功能绘制Sudoku难题

类型drawSudoku

function drawSudoku(B) % Function for drawing the Sudoku board % Copyright 2014 The MathWorks, Inc. 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) % Fill in the clues % % The rows of B are of the form (i,j,k) where i is the row counting from % the top, j is the column, and k is the clue. To place the entries in the % boxes, j is the horizontal distance, 10-i is the vertical distance, and % we subtract 0.5 to center the clue in the box. % % If B is a 9-by-9 matrix, convert it to 3 columns first 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 hold off end