目錄

• • • 一、問題描述
    • 二、模型建立
    • 三、求解思路
    • 四、源程序
    • • 1. 主程序 shudu.m
      • 2. 三維矩陣的轉(zhuǎn)置函數(shù) trans.m
      • 3. 顯示矩陣的函數(shù) showShudu.m
    • 五、結(jié)果展示

一、問題描述

??數(shù)獨(dú)是源自18世紀(jì)瑞士的一種數(shù)學(xué)游戲，盤面是個(gè)九宮，每一宮又分為九個(gè)小格。在這八十一格中根據(jù)已知數(shù)字和解題條件，利用邏輯和推理，在其他的空格上填入1-9的數(shù)字，使1-9中的每個(gè)數(shù)字在每一行、每一列和每一宮中都只出現(xiàn)一次。

二、模型建立

??將 $9\times 9$ 的數(shù)獨(dú)矩陣轉(zhuǎn)化為 $9\times 9\times 9$ 的三維矩陣 ，三維矩陣中每一個(gè)元素的值為0或1。例如，若原數(shù)獨(dú)矩陣的第一行第二列的數(shù)字為7，則三維矩陣第7層的對應(yīng)位置為1，其它層該位置皆為0，總共個(gè) $9\times 9\times 9$ 變量，用 $X(i,j,k)$ 表示第 $k$ 層矩陣的第 $i$ 行第 $j$ 列的元素。

??約束條件：

9層矩陣的同一位置只可以出現(xiàn)一個(gè)1：
$$\sum _{k=1}^{9}X(i,j,k)=1,i=1,2,?,9,j=1,2,?,9$$

任何一層的任何一行只能有一個(gè)1：
$$\sum _{j=1}^{9}X(i,j,k)=1,i=1,2,?,9,k=1,2,?,9$$

任何一層的任何一列只能有一個(gè)1：
$$\sum _{i=1}^{9}X(i,j,k)=1,j=1,2,?,9,k=1,2,?,9$$

每層的九宮格中只能有一個(gè)1：
$$\sum _{i=1}^3\sum _{j=1}^{3}X(i+U,j+V,k)=1,k=1,2,?,9,$$$$U,V\in \{0,3,6\}$$

??總共 $4\times 9^2$ 個(gè)約束條件。
??利用整數(shù)規(guī)劃求解，MATLAB整數(shù)規(guī)劃命令：
??[x,fval,exitflag]=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)
??建立起來的數(shù)學(xué)模型：
$$min{f}^{T}x$$$$\text{subject?to}?\begin{array}{l}A?x<=b\\ Aeq?x=beq\\ lb<=x<=ub\\ \text{x(intcon)?are?integers}\end{array}$$

三、求解思路

??聊之前，先說說系數(shù)矩陣 $Aeq$ 的構(gòu)造，從簡單說起：
??如果要用矩陣表示${\sum }_{i=1}^{n}c{}_{i}x{}_i$，利用矩陣的乘法，有兩種方式：
??1.
$$\left(\begin{array}{c}{c}_{{}_1},c{}_2,?,c{}_n\end{array}\right)?\begin{array}{c}x{}_1\\ x{}_2\\ ?\\ x{}_n\end{array}?$$
??2.
$$\left(\begin{array}{c}{x}_{{}_1},x{}_2,?,x{}_n\end{array}\right)?\begin{array}{c}c{}_1\\ c{}_2\\ ?\\ c{}_n\end{array}?$$
??我們采用第一種，利用系數(shù)向量的轉(zhuǎn)置乘以變量向量：
??將系數(shù)放在相應(yīng)變量的位置，然后對整個(gè)矩陣或者是向量進(jìn)行轉(zhuǎn)置，有幾個(gè)約束條件，與之對應(yīng)地，系數(shù)矩陣 $Aeq$ 就有多少行，具體做法見下文代碼部分。

??##求解思路##：

四、源程序

1. 主程序 shudu.m

clear;clc; N=9^3; % 變量個(gè)數(shù) M=4*(9^2); % 約束條件個(gè)數(shù)% 總體思路是利用整數(shù)規(guī)劃進(jìn)行求解 % MATLAB整數(shù)規(guī)劃命令： % [x,fval,exitflag]=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub) % 分別為目標(biāo)函數(shù)、整型約束、不等式約束，等式約束和解的上下界% 下面分為幾步進(jìn)行: % 一、設(shè)置參數(shù); % 二、求解; % 三、結(jié)果顯示;% 一、設(shè)置參數(shù) % 1.目標(biāo)函數(shù)，并無實(shí)際意義 f=zeros(1,N); % 2.intcon:設(shè)置某些變量的取值為整數(shù), % 這里設(shè)置全部變量 intcon=1:N; % 3.設(shè)置A,b A=[]; b=[];% 4.設(shè)置上下限 % (1).上限 ub=ones(N,1); % (2).確定下限 X0=zeros(9,9,9); % 設(shè)置數(shù)獨(dú)矩陣B,9x9矩陣 B=[0,0,0,0,0,0,2,0,00,1,2,0,0,6,0,9,00,3,0,2,5,0,0,0,10,0,8,1,0,0,0,4,00,0,5,0,7,0,3,0,00,4,0,0,0,2,9,0,02,0,0,0,3,4,0,7,00,7,0,9,0,0,6,5,00,0,3,0,0,0,0,0,0]; % 顯示原始數(shù)獨(dú)矩陣 showShudu(B); % 根據(jù)數(shù)獨(dú)矩陣確定某些變量的下限為1 % 這些變量為已知條件，上限也為2，在求解過程中不會變動 for i=1:9for j=1:9if B(i,j)~=0X0(i,j,B(i,j))=1;endend end % trans:自定義函數(shù)，對參數(shù)中的每一個(gè)二維矩陣進(jìn)行轉(zhuǎn)置操作 % 這樣做是因?yàn)閞eshape命令是按照列來提取矩陣中的元素 % 在本問題中，變量是按行排列的 % 所以為了獲取等式約束中的系數(shù)，需要進(jìn)行轉(zhuǎn)置操作 X0=trans(X0); lb=reshape(X0,N,1);% 5.設(shè)置beq beq=ones(M,1); % 6.設(shè)置Aeq Aeq=zeros(M,N); count=1; % 用于遍歷Aeq的列數(shù) % i,j,k分別表示行數(shù)、列數(shù)、層數(shù) % Amid用于向Aeq中輸入變量前的系數(shù) % 為了與變量相對應(yīng)，同樣需要進(jìn)行轉(zhuǎn)置操作% (1).任何一層的任何一行只能有一個(gè)1 for k=1:9for i=1:9Amid=zeros(9,9,9); % 清空AmidAmid(i,:,k)=1; % 某一層的某一行設(shè)為1 Aeq(count,:)=reshape(trans(Amid),1,N);count=count+1;end end% (2).任何一層的任何一列只能有一個(gè)1 for k=1:9for j=1:9Amid=zeros(9,9,9); % 清空AmidAmid(:,j,k)=1; % 某一層的某一列設(shè)為1Aeq(count,:)=reshape(trans(Amid),1,N);count=count+1;end end% (3).9層矩陣的同一位置只可以出現(xiàn)一個(gè)1 for i=1:9for j=1:9Amid=zeros(9,9,9); % 清空AmidAmid(i,j,:)=1; % 所有層的某一位置設(shè)為1Aeq(count,:)=reshape(trans(Amid),1,N);count=count+1;end end% (4).九宮格約束 for k=1:9for u=0:3:6for v=0:3:6Amid=zeros(9,9,9); % 清空AmidAmid(u+1:u+3,v+1:v+3,k)=1;% 小九宮格位置設(shè)為1Aeq(count,:)=reshape(trans(Amid),1,N);count=count+1;endend end% 二、求解 % 整數(shù)規(guī)劃范式： % x=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub) [x,~,exitflag]=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub);% 三、結(jié)果顯示 if exitflag<=0fprintf('解不存在.\n'); else% exitflag>0,說明存在解X1=reshape(x,9,9,9);% 利用intlinprog得到的解不一定是十分準(zhǔn)確的整數(shù)，在此整數(shù)化X1=round(trans(X1)); % 每層的矩陣都乘以對應(yīng)的層數(shù)，從而得到該位置應(yīng)填入的數(shù)字for k=1:9 X1(:,:,k)=k*X1(:,:,k);end% 合并第三維(對應(yīng)位置元素相加)，得到完成后的數(shù)獨(dú)矩陣C=round(sum(X1,3)); % 結(jié)果顯示showShudu(C); end

2. 三維矩陣的轉(zhuǎn)置函數(shù) trans.m

function B=trans(A) % 實(shí)現(xiàn)三維數(shù)組的轉(zhuǎn)置 [~,~,p]=size(A); B=zeros(size(A)); for k=1:pB(:,:,k)=A(:,:,k)'; end

3. 顯示矩陣的函數(shù) showShudu.m

function showShudu(A) % 以九宮格的形式顯示矩陣,A:9x9矩陣figure;hold on; axis off;axis equal; % 畫線 rectangle('Position',[0.2,0,9,9],'LineWidth',1); % outside border rectangle('Position',[0.2,3,9,3],'LineWidth',1); % horizontal lines rectangle('Position',[3.2,0,3,9],'LineWidth',1); % vertical linesfor i=1:9for j=1:9text(j-0.5,9.5-i,num2str(A(i,j)),'FontSize',18);end end

五、結(jié)果展示

?......再來一組

感謝欣賞，下期再見！

總結(jié)

以上是生活随笔為你收集整理的数独游戏——数模处理(MATLAB整数规划)的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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