當(dāng)然整個建模到出結(jié)果花了我8個小時時間.......

看到GM的這期視頻里的玉米拼圖，我發(fā)現(xiàn)這個結(jié)構(gòu)還是比較簡單的，應(yīng)該可以用比較簡單的數(shù)學(xué)模型建模出來然后求解，正好回顧一下去年學(xué)的數(shù)值最優(yōu)化課程。

當(dāng)然我手頭并沒有這個拼圖，所以一共要分為根據(jù)GM的視頻推測出玉米拼圖的具體結(jié)構(gòu)，按照結(jié)構(gòu)來進(jìn)行數(shù)學(xué)建模，求解模型并轉(zhuǎn)換回到拼圖結(jié)果三個部分。

這里略去玉米拼圖的介紹（包括轉(zhuǎn)換為二維拼圖的部分），還請先看完GM的這期視頻后再來看后面的部分。

拼圖結(jié)構(gòu)

視頻并沒有直接呈現(xiàn)出每塊拼圖的結(jié)構(gòu)，但是還是能夠通過GM拼圖的過程（有剪輯），以及其給出的二維的結(jié)果（有一個小錯誤，實(shí)際是完全鋪滿的），來推斷出來每個拼圖塊的具體結(jié)構(gòu)。

這里就不講具體的推測過程了，整體不是很難但是也需要一些邏輯推導(dǎo)，這也是一個puzzle呢

值得注意的是，這里水平的一格對應(yīng)了實(shí)際拼圖中的兩顆玉米粒，這是由于實(shí)際拼圖都是兩個兩個玉米粒在一起的。

建模部分

這里每個玉米片的狀態(tài)是有限的，并且不算特別多，所以會比較好建模。

這里需要讓玉米片平鋪整個平面，自然需要計(jì)算某一個位置的玉米粒的數(shù)目（這里以拼圖中的兩顆玉米粒為一個）。由上，對于給定的一塊玉米片（i），其處于某一個狀態(tài)（j，例如在第1層，旋轉(zhuǎn)了2個單位，正向朝向）時，在位置（α，β）處的玉米粒數(shù)目我們是可以知道的，記為：

將不同狀態(tài)的結(jié)果寫到一起，應(yīng)該滿足這樣的關(guān)系：

而考慮所有的拼圖貢獻(xiàn)，最終處于位置（α，β）處的玉米粒數(shù)目應(yīng)該將其加起來：

這樣全部鋪滿則需要對于所有的（α，β），上結(jié)果都為 1，也就是：

上述狀態(tài)的選擇可以使用變量來實(shí)現(xiàn)，這里使用這樣的方法：

把某個玉米片（i）是否處于一個狀態(tài)（j）記為 Pij，Pij=1 表示其處于這個狀態(tài)，而 Pij=0 表示不處于這個狀態(tài)。這樣上述狀態(tài)選擇就可以使用乘積求和的方式來表示：

其實(shí)就是將其建模成0-1整型規(guī)劃問題，是數(shù)學(xué)建模中非常常見的模型

則上述拼圖問題就是尋找合適的一組 {Pij}，使得其滿足：

當(dāng)然這里 Pij?還需滿足一些約束條件，一個是每個玉米片只能選擇一種狀態(tài)（這是所有0-1整型規(guī)劃問題都要求的），也就是對于某個固定的?i，對于所有的 j，Pij?有且只有一個為1（其余都為0），可以寫成等價的公式形式（其和為1）：

還有一點(diǎn)是對于這個拼圖，每一層都只能有一個玉米片，并且層數(shù)剛好等于玉米片數(shù)（都為14），也就是每層都有且只有一個玉米片，即在同一層的?Pij?中也有且只有一個為1：

其中 l 表示層數(shù)，k 為與層數(shù)無關(guān)的角標(biāo)，這里意思就是對相同層數(shù)內(nèi)的 Pij?的一個求和。

這樣整個模型就構(gòu)造完成了，剩下的是所有狀態(tài)下?Y 的計(jì)算了。

在開始之前先確定所有的狀態(tài)數(shù)，一般來說，對于一個玉米片，其可以處于上下共14層，旋轉(zhuǎn)一共7個單位（以兩顆玉米為一個單位），以及正反2種，共196種狀態(tài)。

再考慮到一共有14個玉米片，所以變量?Pij?的數(shù)目就一共有2744個。這個數(shù)目當(dāng)然不可能窮舉求解，不過現(xiàn)代的整型線性規(guī)劃求解器（這是一個整型線性規(guī)劃問題）求解這種規(guī)模的問題還是很輕松的。

Y 的計(jì)算

這里將 Y 寫成矩陣的形式，這樣矩陣的行和列就能表示位置（α，β），例如上圖中的玉米片7，（在第1層，旋轉(zhuǎn)0單位，正向，其中認(rèn)為數(shù)字編號在第一列是沒有旋轉(zhuǎn)）寫成矩陣形式就有：

其中是倒過來的原因是矩陣行數(shù)是從上往下的，而拼圖中是從下往上的，為了保證數(shù)值是一致的所以展現(xiàn)出來是倒過來的。

整體的角標(biāo) j 也拆成了三個（s, r, l），其中 s?表示朝向?yàn)檎?）或反（2），r 表示旋轉(zhuǎn)的次數(shù)（r-1），l 表示所在層數(shù)（l）。

使用矩陣后，拼圖反向的操作可以通過將矩陣旋轉(zhuǎn)180度（rot90(Y,2)）實(shí)現(xiàn)（此時層數(shù)有一定問題，需要上下平移一定行數(shù)）；拼圖的旋轉(zhuǎn)可以通過矩陣列方向的循環(huán)平移（cricshift(A,r-1,2)）來實(shí)現(xiàn)；而層數(shù)直接通過選取矩陣的一部分來實(shí)現(xiàn)。

這樣在某些層數(shù)時有些拼圖的一部分可能會超出矩陣（例如上第7塊拼圖的在超過11層后，如 Y7,(1,1,12)?），這里直接將這一部分丟掉，由于最終要求完全鋪滿，所以這種丟掉了一部分的狀態(tài)一定不滿足要求，所以不會是最終的解，符合要求。

轉(zhuǎn)換為標(biāo)準(zhǔn)的線性規(guī)劃問題

這里需要將上述模型轉(zhuǎn)換為matlab能夠接受的標(biāo)準(zhǔn)形式：

其中變量 x 就是本文的?Pij，只需將?Pij?重新排列（reshape）成一維向量即可，使用第三個等式約束，則需要構(gòu)造出矩陣 Aeq，beq。

這里略去具體的步驟，僅說明一下這里?Pij?重新排列的順序依次是：玉米片編號（i），正反向編號（s），旋轉(zhuǎn)編號（r），層數(shù)編號（l）；Aeq 的每行依次對應(yīng)的約束是：位置（α，β）處的玉米粒數(shù)（順序不影響結(jié)果，程序先排列的行數(shù)，共7*14=98行），每個玉米片的狀態(tài)選擇約束（14行），每層只有一個玉米片約束（14行）；beq 全為 1。

排列先后順序意思是前面的序號到達(dá)最大后，后面的才會加一（前面的歸一），詳細(xì)的可以參考代碼（https://github.com/CHanzyLazer/Yu-Mi）

最終得到?Aeq 的圖片還挺好看的（？）

指示所有的變量都為整數(shù)，并且上下界為 [0,1]，就可以交給matlab來求解了。

結(jié)果處理

最后得到的結(jié)果可以想象其基本上都是0，然后有幾個1，所以還需要用程序統(tǒng)計(jì)轉(zhuǎn)換為比較好看懂的結(jié)果。

這個步驟就比較簡單了，首先先把輸出重新排列成以 (i, s, r, l) 為下標(biāo)的四維數(shù)組，直接多重循環(huán)遍歷，找到為1的下標(biāo)?(i', s', r', l') ，則第?i' 個拼圖的狀態(tài)就是 (s=s',?r=r', l=l')，最終結(jié)果寫成表格的形式（p_out）：

按照這個編號也可以得到鋪滿的結(jié)果，宣告建模的正確：

和視頻中的結(jié)果對比就會發(fā)現(xiàn)其實(shí)是整體都倒過來了。

總結(jié)

這個寫成代碼后整個運(yùn)行下來確實(shí)在3秒鐘之內(nèi)就能算出來：

不過整個構(gòu)思加上編程，除去從視頻分析玉米片的結(jié)構(gòu)的時間，也不下8個小時了，而且還是我運(yùn)氣比較好，一次性就做對了情況。

當(dāng)然好處是這種玉米粒類型的拼圖都能這樣求解出來了，設(shè)計(jì)師再換種切法那我也能很快算出來。

當(dāng)然最主要的是通過這個建模的過程我鍛煉了我的數(shù)學(xué)建模能力，復(fù)習(xí)了數(shù)值最優(yōu)化里的知識點(diǎn)，不過比較遺憾的是整型規(guī)劃具體操作比較復(fù)雜，這里就直接用了matlab現(xiàn)有的求解器了。

代碼我上傳到github上了（https://github.com/CHanzyLazer/Yu-Mi），感興趣的可以下載自己運(yùn)行看看（代碼寫和github頁面都弄的挺隨意的）。