前一陣使用Matlab的優(yōu)化工具箱，發(fā)現(xiàn)可以求解數(shù)獨(dú)問題，有意思！實(shí)際上，Matlab優(yōu)化工具箱有兩個：Optimization Toolbox和Global Optimization Toolbox，前者求解常見的優(yōu)化問題，后者求解全局最優(yōu)問題，如常見的方法包括遺傳算法、粒子群算法、模擬退火算法等。

而數(shù)獨(dú)問題則是借助Optimization Toolbox工具箱的整數(shù)線性規(guī)劃函數(shù)intlinprog完成求解。Matlab幫助中有詳細(xì)介紹，搜索Solve Sudoku Puzzles Via Integer Programming可以看到，這里對幫助做下翻譯。 數(shù)獨(dú)問題就是將9x9的網(wǎng)格填充上1~9的整數(shù)值，同時保證每個整數(shù)在每一行、每一列和每個3x3方格只出現(xiàn)一次。網(wǎng)格中分布有一些線索，你的任務(wù)就是填充剩余的網(wǎng)格數(shù)字，如下圖。對于此問題，可以使用二元整數(shù)優(yōu)化求解。
? ? Matlab提供了兩個文件SudokuExample.m和sudokuEngine.m，前者用于發(fā)布講解文件，后者為真正求解數(shù)獨(dú)的函數(shù)調(diào)用，這里將兩個文件中的英文注釋翻譯為中文，如下 ? ??SudokuExample.m文件： %%?利用整數(shù)規(guī)劃求解數(shù)獨(dú)問題 %?這個例子展示了使用二元整數(shù)規(guī)劃來求解數(shù)獨(dú)問題。 %?數(shù)獨(dú)問題就是將9x9的網(wǎng)格填充上1~9的整數(shù)值， %?同時保證每個整數(shù)在每一行、每一列和每個3x3方格只出現(xiàn)一次。 %?網(wǎng)格中分布有一些線索，你的任務(wù)就是填充剩余的網(wǎng)格數(shù)字。 %? %?數(shù)獨(dú)有很多手工和自動解法，如 %?<http://www.mathworks.com/company/newsletters/articles/solving-sudoku-with-matlab.html?Cleve's?Corner(2009)>。 %?這個例子展示了一種基于二元整數(shù)規(guī)劃的直接解法。 %?這種解法特別簡單，因?yàn)椴恍枰獜?fù)雜的求解流程，而只需要表達(dá)清數(shù)獨(dú)的規(guī)則， %?并把線索表示為約束，然后調(diào)用函數(shù)intlinprog完成求解。 ? %%?初始化網(wǎng)格 clc clear close?all ? %?這是一個包含有線索的數(shù)據(jù)矩陣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]; ? %?顯示數(shù)獨(dú)圖 drawSudoku(B) ? %%?二元整數(shù)規(guī)劃解法 %?求解的關(guān)鍵是將數(shù)獨(dú)問題從9x9方格轉(zhuǎn)化為九層9x9二元矩陣，即9個矩陣的元素值只能為0或1， %?而求解結(jié)果中為1的地方第一層相應(yīng)位置設(shè)為1，為2的地方第二層相應(yīng)位置為1，依次類推，為9的地方第九層相應(yīng)位置為1。 %?如此，問題可以通過二元整數(shù)規(guī)劃求解。 % %?這里不需要目標(biāo)函數(shù)，也可以設(shè)為0，因?yàn)閱栴}實(shí)際上就是找到一個滿足所有約束的可行解。 %?然而，為了滿足整數(shù)規(guī)劃內(nèi)部求解器的要求，同時保證求解速度，這里使用一個非常數(shù)目標(biāo)函數(shù)。 ? %%?將數(shù)獨(dú)的規(guī)則表示為約束 %?假設(shè)求解結(jié)果?$x$為一個9x9x9的二元數(shù)組，那么?$x$應(yīng)該具有哪些性質(zhì)呢？ %?首先，數(shù)獨(dú)格中每個格點(diǎn)位置?$(i,j)$都要有一個值，所以在三維9層數(shù)組的9個元素?$x(i,j,1),...,x(i,j,9)$中必有一個非零。 %?換句話說，對每個?$i$和?$j$，有 %? %?$\sum_{k=1}^9?x(i,j,k)?=?1.$ % %?類似地，數(shù)獨(dú)格的每一行，只能有一個1~9間的數(shù)字，即 %? %?$\sum_{j=1}^9?x(i,j,k)?=?1.$ % %?數(shù)獨(dú)格的每一列有相同的性質(zhì)，即 %? %?$\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.$ % %?為了表達(dá)所有9個大格的這種屬性，可以如下表示 %? %?$\sum_{i=1}^3\sum_{j=1}^3?x(i+U,j+V,k)?=?1$ %? %?其中 %? %?$U,V~\epsilon~\{0,3,6\}.$ ? %%?將線索表達(dá)為約束 %?每個初始值，即線索，也可以表達(dá)為約束。假設(shè)?$(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$。 ? %%?數(shù)獨(dú)規(guī)則的數(shù)學(xué)表達(dá) %?盡管數(shù)獨(dú)的規(guī)則可以方便地由9x9x9的x表達(dá)，但是線性約束通常需要表達(dá)為x(:)的向量形式。 %?因此，當(dāng)你表達(dá)一個數(shù)獨(dú)問題時，需要由9x9x9的初始數(shù)組來得到約束矩陣。 %?函數(shù)sudokuEngine中提供了一種配置數(shù)獨(dú)規(guī)則和線索的的方式。 ? %%?數(shù)獨(dú)求解 tic S?=?sudokuEngine(B); toc drawSudoku(S) ? ??sudokuEngine.m文件： function?[S,eflag]?=?sudokuEngine(B) %?求解數(shù)獨(dú)問題。函數(shù)讀取矩陣B中的初始數(shù)據(jù)，然后調(diào)用intlinprog求解，最后返回結(jié)果到S中。 %?矩陣B為3列、至少17行的矩陣（因?yàn)閿?shù)獨(dú)問題至少要有17個初始元素才能保證求解的唯一性）。 %?其中，前兩列是線索的i、j坐標(biāo)，第三列為線索值。 %?如果B是一個9x9矩陣，函數(shù)會首先將其轉(zhuǎn)化為3列矩陣的形式。 ? %?判斷線索的存儲格式 if?isequal(size(B),[9,9])?%?如果是9x9格式 ????%?轉(zhuǎn)化為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的整數(shù) if?sum([any(B?~=?round(B)),any(B?<?1),any(B?>?9)]) ????error('元素必須為1~9的整數(shù)') end ? %%?數(shù)獨(dú)的規(guī)則 N?=?9^3;?%?求解結(jié)果x中獨(dú)立變量的個數(shù)，x為一個9^3長度的列向量 M?=?4*9^2;?%?約束的個數(shù) Aeq?=?zeros(M,N);?%?等式約束矩陣Aeq*x?=?beq beq?=?ones(M,1);?%?常向量beq f?=?(1:N)';?%?目標(biāo)函數(shù)可以是任意值，但是使用非常數(shù)可以加快求解速度 lb?=?zeros(9,9,9);?%?各元素下限 ub?=?lb+1;?%?各元素上限，上下限如此配置可以保證結(jié)果為二元規(guī)劃 ? 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 ? %%?將數(shù)獨(dú)線索表達(dá)為約束 %?將初始線索植入lb數(shù)組，即有線索的地方下限設(shè)為1，這就保證了結(jié)果滿足x(i,j,k)?=?1。 for?i?=?1:size(B,1) ????lb(B(i,1),B(i,2),B(i,3))?=?1; end ? %%?數(shù)獨(dú)求解 %?數(shù)獨(dú)問題的數(shù)學(xué)表達(dá)：規(guī)則由Aeq和beq矩陣表達(dá)，線索由lb中的1表達(dá)。然后調(diào)用intlinprog函數(shù)求解。 %?為了保證問題的解為二元解，設(shè)置intcon參數(shù)為1:N（即所有變量為整數(shù)），同時他們的上下限為0和1。 intcon?=?1:N; [x,~,eflag]?=?intlinprog(f,intcon,[],[],Aeq,beq,lb,ub);?%?雖然lb，ub是三層數(shù)組的形式，但intlinprog第一步會將它們轉(zhuǎn)化為lb(:)和ub(:)的向量格式 ? %%?將結(jié)果轉(zhuǎn)化為可用格式 %?從x得到數(shù)獨(dú)網(wǎng)格只需疊加$(i,j)$位置元素乘其所在的層： if?eflag?>?0?%?線性規(guī)劃順利完成 ????x?=?reshape(x,9,9,9);?%?將x轉(zhuǎn)化為9x9x9的數(shù)組 ????x?=?round(x);?%?四舍五入非整數(shù)解 ????y?=?ones(size(x)); ????for?k?=?2:9 ????????y(:,:,k)?=?k;?%?保存每層層數(shù) ????end ???? ????S?=?x.*y;?%?各元素乘以各層層數(shù) ????S?=?sum(S,3);?%?累加得到9x9的結(jié)果S else ????S?=?[]; end ? ? 對于SudokuExample.m，運(yùn)行Publish，可以得到如下html文件結(jié)果：

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

總結(jié)

以上是生活随笔為你收集整理的利用Matlab优化工具箱解数独问题的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。