{"id":862,"date":"2014-03-14T11:53:48","date_gmt":"2014-03-14T16:53:48","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=862"},"modified":"2014-03-14T11:53:48","modified_gmt":"2014-03-14T16:53:48","slug":"integer-programming-and-hyper-sudoku","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2014\/03\/14\/integer-programming-and-hyper-sudoku\/","title":{"rendered":"Integer Programming and Hyper Sudoku"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>I'd like to introduce this week's guest blogger <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/answers\/contributors\/1033975-alan-weiss\">Alan Weiss<\/a>. Alan writes documentation for mathematical toolboxes here at MathWorks.<\/p><p>Hi, folks. While I have not written a blog post for Loren before, if you use Optimization Toolbox&#8482; or Global Optimization Toolbox then you might have read my work.<\/p><p>I am excited to describe how to use a new solver. Beginning with Release 2014a, <a href=\"https:\/\/www.mathworks.com\/help\/optim\/index.html\">Optimization Toolbox<\/a> has mixed-integer linear programming. The <a href=\"https:\/\/www.mathworks.com\/help\/optim\/ug\/intlinprog.html\"><tt>intlinprog<\/tt><\/a> solver attempts to solve problems with linear objective functions, linear constraints, and (this is the new part) the constraint that some variables must be integer-valued.<\/p><p>The addition of integer constraints opens a surprisingly wide field of problems that are now solvable with Optimization Toolbox. There are examples in the documentation of a complicated factory production problem and solution, a travelling salesman problem, and <a href=\"https:\/\/www.mathworks.com\/help\/optim\/examples\/solve-sudoku-puzzles-via-integer-programming.html\">a Sudoku puzzle solver<\/a>.<\/p><p>This article starts with an explanation of how the Sudoku puzzle solver works. Then it extends the method to solving Hyper Sudoku puzzles. You might not be familiar with Hyper Sudoku; I will describe it presently.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#907e2487-d9d2-4fc5-92b0-bb892fed8d99\">Sudoku<\/a><\/li><li><a href=\"#5a14d408-f200-45fb-bc9f-35ab802c2333\">How to Model Sudoku As a Binary Integer Program<\/a><\/li><li><a href=\"#e4e9e932-9e8b-47b3-b1d4-2b53ab3064a2\">Put the Rules in <tt>intlinprog<\/tt> Form<\/a><\/li><li><a href=\"#a29a2632-4274-4445-94a6-68fd8bb1e889\">Represent the Sudoku Clues<\/a><\/li><li><a href=\"#dc594d32-b936-4841-a88b-732c437c260f\">Solve the Sudoku Puzzle<\/a><\/li><li><a href=\"#9d3e2632-7c2d-4644-80c0-c5764ea190d5\">Hyper Sudoku<\/a><\/li><li><a href=\"#2e1fc9b1-d5bd-4ca7-81d3-5f310d8eff7b\">Include and Solve a Hyper Sudoku Puzzle<\/a><\/li><li><a href=\"#bc9b6377-0594-4237-9d78-6b7f80e0a173\">Further experiments<\/a><\/li><\/ul><\/div><h4>Sudoku<a name=\"907e2487-d9d2-4fc5-92b0-bb892fed8d99\"><\/a><\/h4><p>Let's look at the Sudoku puzzle solver. As you probably know, a Sudoku puzzle is to fill a 9-by-9 grid with integers from 1 through 9 so that each integer appears only once in each row, column, and major 3-by-3 square, as in the following diagram, which includes puzzle clues.<\/p><p>The clues are encoded as follows: each row of a clue matrix has the form [i,j,k]. This means row i, column j, has clue k. If the first row is [1,2,2], then the clue at row 1, column 2 is equal to 2. Both the <tt>sudokuEngine<\/tt> and <tt>drawSudoku<\/tt> functions accept clues in this form, or the more conventional 9-by-9 matrix form.<\/p><pre class=\"codeinput\">B = [1,2,2;\r\n    1,5,3;\r\n    1,8,4;\r\n    2,1,6;\r\n    2,9,3;\r\n    3,3,4;\r\n    3,7,5;\r\n    4,4,8;\r\n    4,6,6;\r\n    5,1,8;\r\n    5,5,1;\r\n    5,9,6;\r\n    6,4,7;\r\n    6,6,5;\r\n    7,3,7;\r\n    7,7,6;\r\n    8,1,4;\r\n    8,9,8;\r\n    9,2,3;\r\n    9,5,4;\r\n    9,8,2];\r\n\r\ndrawSudoku(B)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2014\/hypersudoku_AlanW_01.png\" alt=\"\"> <p><a href=\"https:\/\/www.mathworks.com\/help\/optim\/examples\/solve-sudoku-puzzles-via-integer-programming.html\">The documentation<\/a> shows how to solve this Sudoku puzzle using the <tt>sudokuEngine<\/tt> function:<\/p><pre class=\"codeinput\">drawSudoku(sudokuEngine(B));\r\n<\/pre><pre class=\"codeoutput\">LP:                Optimal objective value is 29565.000000.                                         \r\n\r\nCut Generation:    Applied 1 strong CG cut,                                                         \r\n                   and 2 zero-half cuts.                                                            \r\n                   Lower bound is 29565.000000.                                                     \r\n                   Relative gap is 0.00%.                                                          \r\n\r\n\r\nOptimal solution found.\r\n\r\nIntlinprog stopped at the root node because the objective value is within a gap\r\ntolerance of the optimal value; options.TolGapAbs = 0 (the default value). The\r\nintcon variables are integer within tolerance, options.TolInteger = 1e-05 (the\r\ndefault value).\r\n\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2014\/hypersudoku_AlanW_02.png\" alt=\"\"> <p>How does <tt>sudokuEngine<\/tt> solve Sudoku puzzles? It simply encodes the rules and clues in a binary integer programming problem. The <tt>intlinprog<\/tt> function does the work. Isn't that great? You don't have to figure out a solution algorithm. Compare this with several Sudoku solvers in the File Exchange and <a href=\"https:\/\/www.mathworks.com\/company\/newsletters\/articles\/solving-sudoku-with-matlab.html\">Cleve's article about how to program a Sudoku solver<\/a>. Those all use ingenuity to develop a Sudoku-solving algorithm. But we don't have to figure out an algorithm. We just have to figure out how to represent the rules of Sudoku as an integer linear programming problem. <tt>sudokuEngine<\/tt> uses a binary integer programming formulation.<\/p><h4>How to Model Sudoku As a Binary Integer Program<a name=\"5a14d408-f200-45fb-bc9f-35ab802c2333\"><\/a><\/h4><p>The key trick in turning the rules of Sudoku into a binary integer program is to turn the 9-by-9 grid of clues and answers into a 9-by-9-by-9 cube of binary numbers. Think of the cube as a stack of 9-by-9 grids, one at level 1, one at level 2, up to a grid at level 9. A clue of, say, 5 in the original (3,4) grid becomes a 1 at level 5 in the cube at grid coordinate (3,4). In other words, the (i,j) grid can be viewed as a stack of 9 levels, and a value k in the (i,j) grid becomes a 1 at stack level k.<\/p><pre class=\"codeinput\">RGB = imread(<span class=\"string\">'sudokuGrid3D.jpg'<\/span>);\r\nimage(RGB);\r\naxis <span class=\"string\">image<\/span>;\r\naxis <span class=\"string\">off<\/span>;\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2014\/hypersudoku_AlanW_03.png\" alt=\"\"> <p>Of course, for this to make sense, the solution needs to have exactly one 1 in the stack at each grid coordinate. So if <tt>x(i,j,k)<\/tt> represents the value of the <tt>(i,j)<\/tt> grid at level <tt>k<\/tt>, then we must have<\/p><p>$$\\sum_{k=1}^9 x(i,j,k) = 1$$<\/p><p>at a solution.<\/p><p>What are the other rules of Sudoku? Each row has exactly one number of each value from 1 to 9. In equations, for each $j$ and $k$, we have<\/p><p>$$\\sum_{i=1}^9 x(i,j,k) = 1$$<\/p><p>Each column has exactly one number of each value, too. So for each $i$ and $k$,<\/p><p>$$\\sum_{j=1}^9 x(i,j,k) = 1$$<\/p><p>Each major 3-by-3 grid has exactly one number of each type (1 through 9). For the grid elements $1\\le i\\le 3$ and $1\\le j\\le 3$, and for each $1\\le k\\le 9$,<\/p><p>$$\\sum_{i=1}^3\\sum_{j=1}^3 x(i,j,k) = 1.$$<\/p><p>To represent all nine major grids, just add 3 or 6 to each $i$ and $j$ index:<\/p><p>$$\\sum_{i=1}^3\\sum_{j=1}^3 x(i+U,j+V,k) = 1,$$<\/p><p>where $U,V~\\epsilon~\\{0,3,6\\}.$<\/p><h4>Put the Rules in <tt>intlinprog<\/tt> Form<a name=\"e4e9e932-9e8b-47b3-b1d4-2b53ab3064a2\"><\/a><\/h4><p>The <tt>intlinprog<\/tt> function has a similar syntax to most Optimization Toolbox functions. The constraints on the solution <tt>x<\/tt> are bounds and linear equalities. All the entries in <tt>x<\/tt> are binary numbers, so all the lower bounds are <tt>0<\/tt> and all the upper bounds are <tt>1<\/tt>.<\/p><pre class=\"codeinput\">lb = zeros(9,9,9);\r\nub = 1 + lb;\r\n<\/pre><p>The other constraints are all linear equalities. Optimization Toolbox solvers use the syntax<\/p><p>$$Aeq\\cdot x = beq$$<\/p><p>to represent linear equalities, where <tt>x<\/tt> is regarded as a column vector.<\/p><p><tt>x<\/tt> is naturally a 9-by-9-by-9 matrix. However, you can transform it to a column vector using the colon operator <tt>xcolumn = x(:)<\/tt>. You can recover the 9-by-9-by-9 matrix representation using the <tt>reshape<\/tt> function: <tt>x = reshape(xcolumn,9,9,9)<\/tt>.<\/p><p>Let's see how to express the first equation,<\/p><p>$$\\sum_{k=1}^9 x(i,j,k) = 1,$$<\/p><p>as a matrix <tt>Aeq * x = beq<\/tt>. Each row of the <tt>Aeq<\/tt> matrix is matrix multiplied by the <tt>x<\/tt> matrix. This is the same as saying that entry <tt>i<\/tt> of <tt>beq<\/tt> is equal to the dot product of row <tt>i<\/tt> of <tt>Aeq<\/tt> with the vector <tt>x<\/tt>. So to determine the entries of row <tt>i<\/tt> of <tt>Aeq<\/tt>, make a 9-by-9-by-9 array of the multipliers of <tt>x<\/tt> in the equation, then turn the array into a row vector. For <tt>i<\/tt> and <tt>j<\/tt> ranging from <tt>1<\/tt> to <tt>9<\/tt>, we have<\/p><pre>temp = lb; % A zero array, 9-by-9-by-9\r\ntemp(i,j,:) = 1; % Set the coefficients of x to 1\r\n% Then take temp(:)' as a row of Aeq<\/pre><pre class=\"codeinput\">N = 9^3; <span class=\"comment\">% number of independent variables in x, a 9-by-9-by-9 array<\/span>\r\nM = 4*9^2; <span class=\"comment\">% number of constraints, see the construction of Aeq<\/span>\r\nAeq = zeros(M,N); <span class=\"comment\">% allocate equality constraint matrix Aeq*x = beq<\/span>\r\nbeq = ones(M,1); <span class=\"comment\">% allocate constant vector of ones, beq<\/span>\r\n\r\ncounter = 1;\r\n<span class=\"keyword\">for<\/span> i = 1:9 <span class=\"comment\">% one in each depth<\/span>\r\n    <span class=\"keyword\">for<\/span> j = 1:9\r\n        temp = lb; <span class=\"comment\">% Clear temp<\/span>\r\n        temp(i,j,1:end) = 1; <span class=\"comment\">% Set each depth entry to 1<\/span>\r\n        Aeq(counter,:) = temp(:)'; <span class=\"comment\">% Change the temp array to a row, put it in Aeq<\/span>\r\n        counter = counter + 1;\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>Now for the equation representing each row,<\/p><p>$$\\sum_{i=1}^9 x(i,j,k) = 1$$<\/p><p>do a similar mapping from array to row of <tt>Aeq<\/tt><\/p><pre class=\"codeinput\"><span class=\"keyword\">for<\/span> j = 1:9 <span class=\"comment\">% one in each row<\/span>\r\n    <span class=\"keyword\">for<\/span> k = 1:9\r\n        temp = lb; <span class=\"comment\">% clear temp<\/span>\r\n        temp(1:end,j,k) = 1; <span class=\"comment\">% one row in Aeq*x = beq<\/span>\r\n        Aeq(counter,:) = temp(:)'; <span class=\"comment\">% put temp in a row of Aeq<\/span>\r\n        counter = counter + 1;\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>Similarly, put in the column constraints.<\/p><pre class=\"codeinput\"><span class=\"keyword\">for<\/span> i = 1:9 <span class=\"comment\">% one in each column<\/span>\r\n    <span class=\"keyword\">for<\/span> k = 1:9\r\n        temp = lb;\r\n        temp(i,1:end,k) = 1;\r\n        Aeq(counter,:) = temp(:)';\r\n        counter = counter + 1;\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>Finally, set the coefficients for the major squares.<\/p><pre class=\"codeinput\"><span class=\"keyword\">for<\/span> U = 0:3:6 <span class=\"comment\">% one in each square<\/span>\r\n    <span class=\"keyword\">for<\/span> V = 0:3:6\r\n        <span class=\"keyword\">for<\/span> k = 1:9\r\n            temp = lb;\r\n            temp(U+(1:3),V+(1:3),k) = 1;\r\n            Aeq(counter,:) = temp(:)';\r\n            counter = counter + 1;\r\n        <span class=\"keyword\">end<\/span>\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4>Represent the Sudoku Clues<a name=\"a29a2632-4274-4445-94a6-68fd8bb1e889\"><\/a><\/h4><p>Include the initial clues in the <tt>lb<\/tt> array by setting corresponding entries to 1. This forces the solution to have <tt>x(i,j,k) = 1<\/tt>.<\/p><pre class=\"codeinput\"><span class=\"keyword\">for<\/span> i = 1:size(B,1)\r\n    lb(B(i,1),B(i,2),B(i,3)) = 1;\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4>Solve the Sudoku Puzzle<a name=\"dc594d32-b936-4841-a88b-732c437c260f\"><\/a><\/h4><p>The <tt>Aeq<\/tt> matrix, <tt>beq<\/tt> vector, and bounds <tt>lb<\/tt> and <tt>ub<\/tt> represent the rules of Sudoku and the clues. The objective function vector <tt>f<\/tt> is not important in this case; we just want a solution to the puzzle. We have to tell the solver that all variables are integers.<\/p><pre class=\"codeinput\">opts = optimoptions(<span class=\"string\">'intlinprog'<\/span>,<span class=\"string\">'Display'<\/span>,<span class=\"string\">'off'<\/span>); <span class=\"comment\">% no exit messages<\/span>\r\nintcon = 1:N; <span class=\"comment\">% all variables are integers<\/span>\r\nf = []; <span class=\"comment\">% no objective function, we just want a feasible solution<\/span>\r\n<\/pre><p>Now call <tt>intlinprog<\/tt> to solve the puzzle.<\/p><pre class=\"codeinput\">[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub,opts);\r\n<\/pre><p>To convert the solution to usable form, simply add up the numbers at each $(i,j)$ entry, multiplied by the depth at which the numbers appear:<\/p><pre class=\"codeinput\"><span class=\"keyword\">if<\/span> eflag &gt; 0 <span class=\"comment\">% good solution<\/span>\r\n    x = reshape(x,9,9,9); <span class=\"comment\">% change back to a 9-by-9-by-9 array<\/span>\r\n    x = round(x); <span class=\"comment\">% clean up non-integer solutions<\/span>\r\n    y = ones(size(x));\r\n    <span class=\"keyword\">for<\/span> k = 2:9\r\n        y(:,:,k) = k; <span class=\"comment\">% multiplier for each depth k<\/span>\r\n    <span class=\"keyword\">end<\/span>\r\n\r\n    S = x.*y; <span class=\"comment\">% multiply each entry by its depth<\/span>\r\n    S = sum(S,3); <span class=\"comment\">% S is 9-by-9 and holds the solved puzzle<\/span>\r\n<span class=\"keyword\">else<\/span>\r\n    S = [];\r\n<span class=\"keyword\">end<\/span>\r\n\r\ndrawSudoku(S)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2014\/hypersudoku_AlanW_04.png\" alt=\"\"> <h4>Hyper Sudoku<a name=\"9d3e2632-7c2d-4644-80c0-c5764ea190d5\"><\/a><\/h4><p>Hyper Sudoku differs from Sudoku by having more constraints. There are four 3-by-3 squares in addition to the major 3-by-3 squares that also require exactly one entry of each numeral from 1 through 9. The following image shows <a href=\"http:\/\/en.wikipedia.org\/wiki\/Sudoku#Hypersudoku\">a Hyper Sudoku puzzle taken from Wikipedia<\/a><\/p><pre class=\"codeinput\">RGB = imread(<span class=\"string\">'hyperBoard.png'<\/span>); <span class=\"comment\">% Colored squares are the new constraint regions<\/span>\r\nimage(RGB);\r\naxis <span class=\"string\">image<\/span>;\r\naxis <span class=\"string\">off<\/span>;\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2014\/hypersudoku_AlanW_05.png\" alt=\"\"> <p>It is very easy to write the rules for Hyper Sudoku based on our work to date. There are four new constraint regions, each requiring 9 equations. So allocate a new matrix <tt>Anew<\/tt>, and fill it in.<\/p><pre class=\"codeinput\">Anew = zeros(4*9,N); <span class=\"comment\">% allocate 4*9 new rows<\/span>\r\nbeq = [beq;ones(4*9,1)]; <span class=\"comment\">% extended beq matrix<\/span>\r\n<\/pre><p>The constraints involve the 3-by-3 matrices with row and column indices 2 through 4, and also indices 6 through 8.<\/p><pre class=\"codeinput\">counter = 1; <span class=\"comment\">% reset counter<\/span>\r\nlb = zeros(9,9,9); <span class=\"comment\">% reset lb to all zeros<\/span>\r\n<span class=\"keyword\">for<\/span> U = [2,6] <span class=\"comment\">% start at 2 and 6<\/span>\r\n    <span class=\"keyword\">for<\/span> V = [2,6] <span class=\"comment\">% for both rows and columns<\/span>\r\n        <span class=\"keyword\">for<\/span> k = 1:9\r\n            temp = lb;\r\n            temp(U+(0:2),V+(0:2),k) = 1;\r\n            Anew(counter,:) = temp(:)';\r\n            counter = counter + 1;\r\n        <span class=\"keyword\">end<\/span>\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n\r\nAeq = [Aeq;Anew]; <span class=\"comment\">% Append Anew to Aeq<\/span>\r\n<\/pre><h4>Include and Solve a Hyper Sudoku Puzzle<a name=\"2e1fc9b1-d5bd-4ca7-81d3-5f310d8eff7b\"><\/a><\/h4><p>The Hyper Sudoku puzzle pictured above can be represented by the following matrix.<\/p><pre class=\"codeinput\">B2 = [1     8     1\r\n     2     3     2\r\n     2     8     3\r\n     2     9     4\r\n     3     5     5\r\n     3     6     1\r\n     4     6     6\r\n     4     7     5\r\n     5     2     7\r\n     5     4     3\r\n     5     8     8\r\n     6     3     3\r\n     7     5     8\r\n     8     1     5\r\n     8     2     8\r\n     8     7     9\r\n     9     1     6\r\n     9     2     9];\r\n<\/pre><p>Include the clues in the lower bound.<\/p><pre class=\"codeinput\"> <span class=\"keyword\">for<\/span> i = 1:size(B2,1)\r\n    lb(B2(i,1),B2(i,2),B2(i,3)) = 1;\r\n <span class=\"keyword\">end<\/span>\r\n<\/pre><p>Solve the puzzle.<\/p><pre class=\"codeinput\">[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub,opts);\r\n\r\n<span class=\"keyword\">if<\/span> eflag &gt; 0 <span class=\"comment\">% good solution<\/span>\r\n    x = reshape(x,9,9,9); <span class=\"comment\">% change back to a 9-by-9-by-9 array<\/span>\r\n    x = round(x); <span class=\"comment\">% clean up non-integer solutions<\/span>\r\n    y = ones(size(x));\r\n    <span class=\"keyword\">for<\/span> k = 2:9\r\n        y(:,:,k) = k; <span class=\"comment\">% multiplier for each depth k<\/span>\r\n    <span class=\"keyword\">end<\/span>\r\n\r\n    S = x.*y; <span class=\"comment\">% multiply each entry by its depth<\/span>\r\n    S = sum(S,3); <span class=\"comment\">% S is 9-by-9 and holds the solved puzzle<\/span>\r\n<span class=\"keyword\">else<\/span>\r\n    S = [];\r\n<span class=\"keyword\">end<\/span>\r\n\r\n<span class=\"comment\">% drawSudoku(S) % this would draw the board. But for a Hyper Sudoku figure:<\/span>\r\nRGB = imread(<span class=\"string\">'hyperBoardSolved.png'<\/span>);\r\nimage(RGB);\r\naxis <span class=\"string\">image<\/span>;\r\naxis <span class=\"string\">off<\/span>;\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2014\/hypersudoku_AlanW_06.png\" alt=\"\"> <h4>Further experiments<a name=\"bc9b6377-0594-4237-9d78-6b7f80e0a173\"><\/a><\/h4><p>The <tt>drawSudoku<\/tt> and <tt>sudokuEngine<\/tt> utilities ship with R2014a, so you can solve your own puzzles easily. However, perhaps not as easily as you might like. I also submitted an <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/45791-visual-sudoku\">entry to the File Exchange called <tt>visualSudoku<\/tt><\/a> that is a visual editor for Sudoku puzzles. <tt>visualSudoku<\/tt> has no real code, it relies on <tt>sudokuEngine<\/tt> to solve puzzles, but it makes it easy to solve, save, and import Sudoku puzzles.<\/p><p>I also uploaded a <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/45792-visual-hyper-sudoku\">visual Hyper Sudoku solver<\/a>, and an engine for solving Hyper Sudoku puzzles.<\/p><p>The constraint matrices for Sudoku and Hyper Sudoku are quite sparse. I did not use that sparsity in the formulation. For examples using sparse constraint matrices, see the <a href=\"https:\/\/www.mathworks.com\/help\/optim\/examples\/factory-warehouse-sales-allocation-model.html\">factory example<\/a> or the <a href=\"https:\/\/www.mathworks.com\/help\/optim\/examples\/travelling-salesman-problem.html\">travelling salesman problem<\/a> in the documentation.<\/p><p>Did you enjoy seeing how to solve Sudoku or Hyper Sudoku puzzles with integer programming? Do you have any other types of puzzles you solve this way? Tell us about it <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=862#respond\">here<\/a>.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_d624770e53c64c3ba88478fd5772d649() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='d624770e53c64c3ba88478fd5772d649 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' d624770e53c64c3ba88478fd5772d649';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2014 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_d624770e53c64c3ba88478fd5772d649()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2014a<br><\/p><p class=\"footer\"><br>\r\n      Published with MATLAB&reg; R2014a<br><\/p><\/div><!--\r\nd624770e53c64c3ba88478fd5772d649 ##### SOURCE BEGIN #####\r\n%% Integer Programming and Hyper Sudoku\r\n% I'd like to introduce this week's guest blogger\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/answers\/contributors\/1033975-alan-weiss Alan Weiss>.\r\n% Alan writes documentation for mathematical toolboxes here at MathWorks. \r\n%\r\n% Hi, folks. While I have not written a blog post for Loren before, if you\r\n% use Optimization Toolbox(TM) or Global Optimization Toolbox then you\r\n% might have read my work.\r\n%\r\n% I am excited to describe how to use a new solver. Beginning with Release\r\n% 2014a, <https:\/\/www.mathworks.com\/help\/optim\/index.html Optimization\r\n% Toolbox> has mixed-integer linear programming. The\r\n% <https:\/\/www.mathworks.com\/help\/optim\/ug\/intlinprog.html |intlinprog|>\r\n% solver attempts to solve problems with linear objective functions, linear\r\n% constraints, and (this is the new part) the constraint that some\r\n% variables must be integer-valued.\r\n%\r\n% The addition of integer constraints opens a surprisingly wide field of\r\n% problems that are now solvable with Optimization Toolbox. There are\r\n% examples in the documentation of a complicated factory production problem\r\n% and solution, a travelling salesman problem, and\r\n% <https:\/\/www.mathworks.com\/help\/optim\/examples\/solve-sudoku-puzzles-via-integer-programming.html\r\n% a Sudoku puzzle solver>.\r\n%\r\n% This article starts with an explanation of how the Sudoku puzzle solver\r\n% works. Then it extends the method to solving Hyper Sudoku puzzles. You\r\n% might not be familiar with Hyper Sudoku; I will describe it presently.\r\n\r\n%% Sudoku\r\n% Let's look at the Sudoku puzzle solver. As you probably know, a Sudoku\r\n% puzzle is to fill a 9-by-9 grid with integers from 1 through 9 so that\r\n% each integer appears only once in each row, column, and major 3-by-3\r\n% square, as in the following diagram, which includes puzzle clues.\r\n%\r\n% The clues are encoded as follows: each row of a clue matrix has the form\r\n% [i,j,k]. This means row i, column j, has clue k. If the first row is\r\n% [1,2,2], then the clue at row 1, column 2 is equal to 2. Both the\r\n% |sudokuEngine| and |drawSudoku| functions accept clues in this form, or\r\n% the more conventional 9-by-9 matrix form.\r\n\r\nB = [1,2,2;\r\n    1,5,3;\r\n    1,8,4;\r\n    2,1,6;\r\n    2,9,3;\r\n    3,3,4;\r\n    3,7,5;\r\n    4,4,8;\r\n    4,6,6;\r\n    5,1,8;\r\n    5,5,1;\r\n    5,9,6;\r\n    6,4,7;\r\n    6,6,5;\r\n    7,3,7;\r\n    7,7,6;\r\n    8,1,4;\r\n    8,9,8;\r\n    9,2,3;\r\n    9,5,4;\r\n    9,8,2];\r\n\r\ndrawSudoku(B)\r\n\r\n%%\r\n% <https:\/\/www.mathworks.com\/help\/optim\/examples\/solve-sudoku-puzzles-via-integer-programming.html\r\n% The documentation> shows how to solve this Sudoku puzzle using the\r\n% |sudokuEngine| function:\r\n\r\ndrawSudoku(sudokuEngine(B));\r\n\r\n%%\r\n% How does |sudokuEngine| solve Sudoku puzzles? It simply encodes the rules\r\n% and clues in a binary integer programming problem. The |intlinprog|\r\n% function does the work. Isn't that great? You don't have to figure out a\r\n% solution algorithm. Compare this with several Sudoku solvers in the File\r\n% Exchange and\r\n% <https:\/\/www.mathworks.com\/company\/newsletters\/articles\/solving-sudoku-with-matlab.html\r\n% Cleve's article about how to program a Sudoku solver>. Those all use\r\n% ingenuity to develop a Sudoku-solving algorithm. But we don't have to\r\n% figure out an algorithm. We just have to figure out how to represent the\r\n% rules of Sudoku as an integer linear programming problem. |sudokuEngine|\r\n% uses a binary integer programming formulation.\r\n\r\n%% How to Model Sudoku As a Binary Integer Program\r\n%\r\n% The key trick in turning the rules of Sudoku into a binary integer\r\n% program is to turn the 9-by-9 grid of clues and answers into a\r\n% 9-by-9-by-9 cube of binary numbers. Think of the cube as a stack of\r\n% 9-by-9 grids, one at level 1, one at level 2, up to a grid at level 9. A\r\n% clue of, say, 5 in the original (3,4) grid becomes a 1 at level 5 in the\r\n% cube at grid coordinate (3,4). In other words, the (i,j) grid can be\r\n% viewed as a stack of 9 levels, and a value k in the (i,j) grid becomes a\r\n% 1 at stack level k.\r\nRGB = imread('sudokuGrid3D.jpg');\r\nimage(RGB);\r\naxis image;\r\naxis off;\r\n%%\r\n% Of course, for this to make sense, the solution needs to have exactly one\r\n% 1 in the stack at each grid coordinate. So if |x(i,j,k)| represents the\r\n% value of the |(i,j)| grid at level |k|, then we must have\r\n%\r\n% $$\\sum_{k=1}^9 x(i,j,k) = 1$$\r\n%\r\n% at a solution.\r\n%\r\n% What are the other rules of Sudoku? Each row has exactly one number of\r\n% each value from 1 to 9. In equations, for each $j$ and $k$, we have\r\n%\r\n% $$\\sum_{i=1}^9 x(i,j,k) = 1$$\r\n%\r\n% Each column has exactly one number of each value, too. So for each $i$\r\n% and $k$,\r\n%\r\n% $$\\sum_{j=1}^9 x(i,j,k) = 1$$\r\n%\r\n% Each major 3-by-3 grid has exactly one number of each type (1 through 9).\r\n% For the grid elements\r\n% $1\\le i\\le 3$ and $1\\le j\\le 3$, and for each $1\\le k\\le 9$,\r\n%\r\n% $$\\sum_{i=1}^3\\sum_{j=1}^3 x(i,j,k) = 1.$$\r\n%\r\n% To represent all nine major grids, just add 3 or 6 to each $i$ and $j$\r\n% index:\r\n%\r\n% $$\\sum_{i=1}^3\\sum_{j=1}^3 x(i+U,j+V,k) = 1,$$\r\n%\r\n% where $U,V~\\epsilon~\\{0,3,6\\}.$\r\n%% Put the Rules in |intlinprog| Form\r\n% The |intlinprog| function has a similar syntax to most Optimization\r\n% Toolbox functions. The constraints on the solution |x| are bounds and\r\n% linear equalities. All the entries in |x| are binary numbers, so all the\r\n% lower bounds are |0| and all the upper bounds are |1|.\r\nlb = zeros(9,9,9);\r\nub = 1 + lb;\r\n%%\r\n% The other constraints are all linear equalities. Optimization Toolbox\r\n% solvers use the syntax\r\n%\r\n% $$Aeq\\cdot x = beq$$\r\n%\r\n% to represent linear equalities, where |x| is regarded as a column vector.\r\n% \r\n% |x| is naturally a 9-by-9-by-9 matrix. However, you can transform it to a\r\n% column vector using the colon operator |xcolumn = x(:)|. You can recover\r\n% the 9-by-9-by-9 matrix representation using the |reshape| function: |x =\r\n% reshape(xcolumn,9,9,9)|.\r\n%\r\n% Let's see how to express the first equation,\r\n%\r\n% $$\\sum_{k=1}^9 x(i,j,k) = 1,$$\r\n%\r\n% as a matrix |Aeq * x = beq|. Each row of the |Aeq| matrix is matrix\r\n% multiplied by the |x| matrix. This is the same as saying that entry |i|\r\n% of |beq| is equal to the dot product of row |i| of |Aeq| with the vector\r\n% |x|. So to determine the entries of row |i| of |Aeq|, make a 9-by-9-by-9\r\n% array of the multipliers of |x| in the equation, then turn the array into\r\n% a row vector. For |i| and |j| ranging from |1| to |9|, we have\r\n%\r\n%  temp = lb; % A zero array, 9-by-9-by-9\r\n%  temp(i,j,:) = 1; % Set the coefficients of x to 1\r\n%  % Then take temp(:)' as a row of Aeq\r\nN = 9^3; % number of independent variables in x, a 9-by-9-by-9 array\r\nM = 4*9^2; % number of constraints, see the construction of Aeq\r\nAeq = zeros(M,N); % allocate equality constraint matrix Aeq*x = beq\r\nbeq = ones(M,1); % allocate constant vector of ones, beq\r\n\r\ncounter = 1;\r\nfor i = 1:9 % one in each depth\r\n    for j = 1:9\r\n        temp = lb; % Clear temp\r\n        temp(i,j,1:end) = 1; % Set each depth entry to 1 \r\n        Aeq(counter,:) = temp(:)'; % Change the temp array to a row, put it in Aeq\r\n        counter = counter + 1;\r\n    end\r\nend\r\n\r\n%%\r\n% Now for the equation representing each row,\r\n%\r\n% $$\\sum_{i=1}^9 x(i,j,k) = 1$$\r\n%\r\n% do a similar mapping from array to row of |Aeq|\r\nfor j = 1:9 % one in each row\r\n    for k = 1:9\r\n        temp = lb; % clear temp\r\n        temp(1:end,j,k) = 1; % one row in Aeq*x = beq\r\n        Aeq(counter,:) = temp(:)'; % put temp in a row of Aeq\r\n        counter = counter + 1;\r\n    end\r\nend\r\n%%\r\n% Similarly, put in the column constraints.\r\nfor i = 1:9 % one in each column\r\n    for k = 1:9\r\n        temp = lb;\r\n        temp(i,1:end,k) = 1;\r\n        Aeq(counter,:) = temp(:)';\r\n        counter = counter + 1;\r\n    end\r\nend\r\n%%\r\n% Finally, set the coefficients for the major squares.\r\nfor U = 0:3:6 % one in each square\r\n    for V = 0:3:6\r\n        for k = 1:9\r\n            temp = lb;\r\n            temp(U+(1:3),V+(1:3),k) = 1;\r\n            Aeq(counter,:) = temp(:)';\r\n            counter = counter + 1;\r\n        end\r\n    end\r\nend\r\n%% Represent the Sudoku Clues\r\n% Include the initial clues in the |lb| array by setting corresponding\r\n% entries to 1. This forces the solution to have |x(i,j,k) = 1|.\r\n\r\nfor i = 1:size(B,1)\r\n    lb(B(i,1),B(i,2),B(i,3)) = 1;\r\nend\r\n%% Solve the Sudoku Puzzle\r\n% The |Aeq| matrix, |beq| vector, and bounds |lb| and |ub| represent the\r\n% rules of Sudoku and the clues. The objective function vector |f| is not\r\n% important in this case; we just want a solution to the puzzle. We have to\r\n% tell the solver that all variables are integers.\r\nopts = optimoptions('intlinprog','Display','off'); % no exit messages\r\nintcon = 1:N; % all variables are integers\r\nf = []; % no objective function, we just want a feasible solution\r\n%%\r\n% Now call |intlinprog| to solve the puzzle.\r\n[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub,opts);\r\n%%\r\n% To convert the solution to usable form, simply add up the numbers at each\r\n% $(i,j)$ entry, multiplied by the depth at which the numbers appear:\r\n\r\nif eflag > 0 % good solution\r\n    x = reshape(x,9,9,9); % change back to a 9-by-9-by-9 array\r\n    x = round(x); % clean up non-integer solutions\r\n    y = ones(size(x));\r\n    for k = 2:9\r\n        y(:,:,k) = k; % multiplier for each depth k\r\n    end\r\n\r\n    S = x.*y; % multiply each entry by its depth\r\n    S = sum(S,3); % S is 9-by-9 and holds the solved puzzle\r\nelse\r\n    S = [];\r\nend\r\n\r\ndrawSudoku(S)\r\n\r\n%% Hyper Sudoku\r\n% Hyper Sudoku differs from Sudoku by having more constraints. There are\r\n% four 3-by-3 squares in addition to the major 3-by-3 squares that also\r\n% require exactly one entry of each numeral from 1 through 9. The following\r\n% image shows <http:\/\/en.wikipedia.org\/wiki\/Sudoku#Hypersudoku a Hyper\r\n% Sudoku puzzle taken from Wikipedia>\r\n\r\nRGB = imread('hyperBoard.png'); % Colored squares are the new constraint regions\r\nimage(RGB);\r\naxis image;\r\naxis off;\r\n\r\n%%\r\n% It is very easy to write the rules for Hyper Sudoku based on our work to\r\n% date. There are four new constraint regions, each requiring 9 equations.\r\n% So allocate a new matrix |Anew|, and fill it in.\r\nAnew = zeros(4*9,N); % allocate 4*9 new rows\r\nbeq = [beq;ones(4*9,1)]; % extended beq matrix\r\n%%\r\n% The constraints involve the 3-by-3 matrices with row and column indices 2\r\n% through 4, and also indices 6 through 8.\r\ncounter = 1; % reset counter\r\nlb = zeros(9,9,9); % reset lb to all zeros\r\nfor U = [2,6] % start at 2 and 6\r\n    for V = [2,6] % for both rows and columns\r\n        for k = 1:9\r\n            temp = lb;\r\n            temp(U+(0:2),V+(0:2),k) = 1;\r\n            Anew(counter,:) = temp(:)';\r\n            counter = counter + 1;\r\n        end\r\n    end\r\nend\r\n\r\nAeq = [Aeq;Anew]; % Append Anew to Aeq\r\n%% Include and Solve a Hyper Sudoku Puzzle\r\n% The Hyper Sudoku puzzle pictured above can be represented by the\r\n% following matrix.\r\nB2 = [1     8     1\r\n     2     3     2\r\n     2     8     3\r\n     2     9     4\r\n     3     5     5\r\n     3     6     1\r\n     4     6     6\r\n     4     7     5\r\n     5     2     7\r\n     5     4     3\r\n     5     8     8\r\n     6     3     3\r\n     7     5     8\r\n     8     1     5\r\n     8     2     8\r\n     8     7     9\r\n     9     1     6\r\n     9     2     9];\r\n %%\r\n % Include the clues in the lower bound.\r\n for i = 1:size(B2,1)\r\n    lb(B2(i,1),B2(i,2),B2(i,3)) = 1;\r\n end\r\n%%\r\n% Solve the puzzle.\r\n[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub,opts);\r\n\r\nif eflag > 0 % good solution\r\n    x = reshape(x,9,9,9); % change back to a 9-by-9-by-9 array\r\n    x = round(x); % clean up non-integer solutions\r\n    y = ones(size(x));\r\n    for k = 2:9\r\n        y(:,:,k) = k; % multiplier for each depth k\r\n    end\r\n\r\n    S = x.*y; % multiply each entry by its depth\r\n    S = sum(S,3); % S is 9-by-9 and holds the solved puzzle\r\nelse\r\n    S = [];\r\nend\r\n\r\n% drawSudoku(S) % this would draw the board. But for a Hyper Sudoku figure:\r\nRGB = imread('hyperBoardSolved.png');\r\nimage(RGB);\r\naxis image;\r\naxis off;\r\n\r\n%% Further experiments\r\n% The |drawSudoku| and |sudokuEngine| utilities ship with R2014a, so you\r\n% can solve your own puzzles easily. However, perhaps not as easily as you\r\n% might like. I also submitted an\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/45791-visual-sudoku entry to the\r\n% File Exchange called |visualSudoku|> that is a visual editor for Sudoku\r\n% puzzles. |visualSudoku| has no real code, it relies on |sudokuEngine| to\r\n% solve puzzles, but it makes it easy to solve, save, and import Sudoku\r\n% puzzles.\r\n%\r\n% I also uploaded a\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/45792-visual-hyper-sudoku  visual Hyper\r\n% Sudoku solver>, and an engine for solving Hyper Sudoku puzzles.\r\n%\r\n% The constraint matrices for Sudoku and Hyper Sudoku are quite sparse. I\r\n% did not use that sparsity in the formulation. For examples using sparse\r\n% constraint matrices, see the\r\n% <https:\/\/www.mathworks.com\/help\/optim\/examples\/factory-warehouse-sales-allocation-model.html\r\n% factory example> or the\r\n% <https:\/\/www.mathworks.com\/help\/optim\/examples\/travelling-salesman-problem.html\r\n% travelling salesman problem> in the documentation.\r\n%\r\n% Did you enjoy seeing how to solve Sudoku or Hyper Sudoku puzzles with\r\n% integer programming? Do you have any other types of puzzles you solve\r\n% this way? Tell us about it <https:\/\/blogs.mathworks.com\/loren\/?p=862#respond here>.\r\n\r\n\r\n##### SOURCE END ##### d624770e53c64c3ba88478fd5772d649\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2014\/hypersudoku_AlanW_06.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>I'd like to introduce this week's guest blogger <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/answers\/contributors\/1033975-alan-weiss\">Alan Weiss<\/a>. Alan writes documentation for mathematical toolboxes here at MathWorks.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2014\/03\/14\/integer-programming-and-hyper-sudoku\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[33,6,29],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/862"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=862"}],"version-history":[{"count":6,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/862\/revisions"}],"predecessor-version":[{"id":869,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/862\/revisions\/869"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=862"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=862"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=862"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}