QUBO Models入門(mén)資料推薦以及編程求解
Quadratic unconstrained binary optimization，QUBO中文名是二次無(wú)約束二元優(yōu)化，它是在二次規(guī)劃(QP, Quadratic Programming)的基礎(chǔ)上添加了兩個(gè)限制條件：（1）只有目標(biāo)函數(shù)，沒(méi)有約束條件，例如等式約束、不等式約束等；（2）決策變量的取值只能是0和1。
下面給出它的標(biāo)準(zhǔn)形式：
圖片來(lái)源：參考文獻(xiàn)1，見(jiàn)文章最后
從上圖可知：
（1）x是決策變量，只有兩種取值情況，通常用0和1表示。Q是一個(gè)常數(shù)組成的方陣！
（2）QUBO有兩種形式，取決于矩陣Q是一個(gè)對(duì)稱(chēng)的矩陣還是一個(gè)上三角矩陣。兩種形式是等價(jià)的，對(duì)稱(chēng)的形式看起來(lái)更加直觀。
下面來(lái)看一個(gè)實(shí)際的例子：
大家可以思考：
（1）有幾個(gè)決策變量？
（2）矩陣Q的大小是多少？
（3）xi的平方和xi有什么關(guān)系？
（4）你能試著寫(xiě)出Q矩陣嗎？
下面給出答案：
如果大家學(xué)過(guò)線(xiàn)性代數(shù)的話(huà)，應(yīng)該會(huì)想到二次型變換：
再次強(qiáng)調(diào)：除了決策變量是01約束之外，不能加入其他的約束條件，所有的約束條件都要想辦法轉(zhuǎn)換到Q矩陣中。
此外，許多其他看起來(lái)與QUBO問(wèn)題無(wú)關(guān)的問(wèn)題可以重新表述為QUBO模型。下面看具體的例子。
Number Partitioning Problem 問(wèn)題，也稱(chēng)為數(shù)字劃分問(wèn)題
給定一個(gè)集合S = {3,1,1,2,2,1}
如何將其元素劃分給兩個(gè)集合，使得兩個(gè)集合中的元素的和相等。
例如集合1：S(1) = {1,1,1,2} 、集合2：S(2) = {2,3}
此時(shí)每個(gè)集合中的元素相加都為5，這就是問(wèn)題的一個(gè)解。
然而，有時(shí)候無(wú)論如何劃分兩個(gè)集合的和也不會(huì)相等，
例如：S = {8,1,1,2,2,1}，此時(shí)我們需要找到兩個(gè)集合中的元素的和盡可能接近的一組解。
下面我們來(lái)看建模過(guò)程：
假設(shè)S中有m個(gè)元素，分別是s1,s2,...,sm，如果第j個(gè)元素sj屬于第一個(gè)集合，那么sj就等于1（否則sj屬于第二個(gè)集合，此時(shí)sj等于0）。
那么，第一個(gè)集合中的元素和就是sum1，第二個(gè)集合中的元素和就可以用S中所有的元素和減去sum1,得到的結(jié)果是sum2。
兩個(gè)集合的元素和的差異是：
c是S中所有元素的和，因此是一個(gè)常數(shù)。
我們的目標(biāo)函數(shù)是使得diff盡可能接近于0，因此有兩種處理方法：
（1）min abs(diff)? 求diff的絕對(duì)值，絕對(duì)值越小越好。
（2）min (diff)^2 ??求diff的平方，平方越小越好。
那么，哪個(gè)目標(biāo)函數(shù)和QUBO的關(guān)系更接近呢？
很明顯，第二種方法和QUBO的關(guān)系更大！
注意：這里要將平方項(xiàng)展開(kāi)（具體的展開(kāi)形式論文中省略了，大家有時(shí)間可以自己推導(dǎo)看看，應(yīng)該有點(diǎn)復(fù)雜），然后得到Q矩陣，Q矩陣是一個(gè)m階的對(duì)稱(chēng)矩陣！里面所有元素都是已知的。
下面給了一個(gè)數(shù)值計(jì)算的例子：
對(duì)于上面這個(gè)優(yōu)化問(wèn)題，有很多種求解方法：
（1）使用MATLAB求解
由于MATLAB沒(méi)有求解QUBO的內(nèi)置函數(shù)，因此需要借助第三方求解器。第三方求解器有很多，這里推薦一個(gè)免費(fèi)的工具箱：OPTI Toolbox，這個(gè)工具箱中集成了很多第三方的求解器。官網(wǎng)下載地址：
https://www.controlengineering.co.nz/Wikis/OPTI/pmwiki.php/DL/DownloadOPTI
OPTI Toolbox集成了一系列優(yōu)化求解器，安裝好后在Matlab命令窗口輸入optiSolver，即可看到自帶的求解器。
然而，對(duì)于QUBO優(yōu)化問(wèn)題，工具箱中自帶的求解器很弱，容易求出局部最優(yōu)解，因此我們我們需要補(bǔ)充安裝一個(gè)更強(qiáng)大的SCIP求解器。
學(xué)術(shù)版本的SCIP可以免費(fèi)使用，只需要提交一個(gè)郵箱即可下載scip.mexw64，下載好后需要把文件放到OPTI的Solvers文件夾。
接下來(lái)以我們這個(gè)問(wèn)題為例，試試工具箱的計(jì)算。
OPTI官網(wǎng)給出了支持的標(biāo)準(zhǔn)的MIQP（混合整數(shù)二次規(guī)劃問(wèn)題）形式：
https://www.controlengineering.co.nz/Wikis/OPTI/pmwiki.php/Probs/MIQP
我們要求解的QUBO可以看成MIQP問(wèn)題的特例！
第一步：得到Q矩陣：
s=[25,7,13,31,42,17,21,10]; ?% 8個(gè)元素
c = sum(s);
m = numel(s); ?
Q = zeros(m,m);
for ii =1:m
? ?for jj = ii:m
? ? ? ?if ii == jj
? ? ? ? ? ?Q(ii,jj) = s(ii)*(s(ii)-c);
? ? ? ?else
? ? ? ? ? ? Q(ii,jj) = s(ii)*s(jj);
? ? ? ? ? ? ?Q(jj,ii) = s(ii)*s(jj);
? ? ? ?end
? ?end
end
Q
第二步：調(diào)用求解器
H = Q;
f = zeros(m,1);
xtype = char(join( repmat("B",1,m),""));
% xtype = 'BBBBBBBB' ? ?8個(gè)B表示有8個(gè)二值約束
Opt = opti('qp',H,f,'xtype',xtype)
% Solve the MIQP problem
[x,fval,exitflag,info] = solve(Opt)
這個(gè)警告的意思是，目標(biāo)函數(shù)可能是非凸的，也就是說(shuō)，它可能有多個(gè)局部最優(yōu)解，而不是一個(gè)全局最優(yōu)解。這樣的話(huà)，你得到的解可能不是全局最優(yōu)解，而是一個(gè)局部最優(yōu)解。
從倒數(shù)第二段可以看出，使用的是SCIP工具箱，算法是空間分枝定界法。
最優(yōu)解對(duì)應(yīng)的目標(biāo)函數(shù)是-3444.5（實(shí)際上要乘以2等于-6889，工具箱的標(biāo)準(zhǔn)型上將目標(biāo)函數(shù)乘以了1/2）。
兩個(gè)集合包含的元素：
這個(gè)答案和參考資料里面給的答案不同：
但是最優(yōu)結(jié)果是一樣的，這說(shuō)明我們這個(gè)問(wèn)題有多個(gè)解。
（2）使用Python求解
這里推薦一個(gè)很強(qiáng)大的第三方庫(kù)：PyQUBO。
官網(wǎng)：https://pyqubo.readthedocs.io/en/latest/getting_started.html
它可以將組合優(yōu)化問(wèn)題映射為QUBO形式，剛好找到了一篇文章介紹這個(gè)庫(kù)的使用：
https://cloud.tencent.com/developer/article/2111550
如果你將QUBO看成MIQP問(wèn)題的特例，那么可以求解的第三方包就更多了，例如CVXPY、pyMIQP等。
（3）使用lingo求解
LINGO是Linear Interactive and General Optimizer的縮寫(xiě)，即“交互式的線(xiàn)性和通用優(yōu)化求解器”，由美國(guó)LINDO系統(tǒng)公司（Lindo System Inc.）推出的，可以用于求解非線(xiàn)性規(guī)劃，也可以用于一些線(xiàn)性和非線(xiàn)性方程組的求解等，功能十分強(qiáng)大，是求解優(yōu)化模型的最佳選擇。其特色在于內(nèi)置建模語(yǔ)言，提供了許多內(nèi)部函數(shù)，可以允許決策變量是整數(shù)（即整數(shù)規(guī)劃，包括 0-1 整數(shù)規(guī)劃），方便靈活，而且執(zhí)行速度非?？?。
關(guān)于lingo求解二次整數(shù)規(guī)劃的問(wèn)題，大家可以自己查閱相關(guān)資料。
有約束條件時(shí)如何轉(zhuǎn)換？
注意，上面我們的例子中除了要求x是01變量外，是沒(méi)有其他約束的，然而在實(shí)際例子中，往往有其他的約束限制，我們需要對(duì)這些約束進(jìn)行轉(zhuǎn)換，轉(zhuǎn)換的思想很簡(jiǎn)單，假設(shè)我們的目的是計(jì)算目標(biāo)函數(shù)的最小值，那么我們可以在目標(biāo)函數(shù)上增加一個(gè)非負(fù)的懲罰函數(shù)（設(shè)計(jì)為關(guān)于決策變量的二次函數(shù)形式），該函數(shù)滿(mǎn)足下面兩點(diǎn)要求：
（1）如果當(dāng)前解滿(mǎn)足約束，那么懲罰函數(shù)的值等于0；
（2）如果當(dāng)前解不滿(mǎn)足約束，那么懲罰函數(shù)是一個(gè)很大的正數(shù)。
我們來(lái)看參考文獻(xiàn)的介紹：
P是一個(gè)很大的正數(shù)，如果違反約束的話(huà)，目標(biāo)函數(shù)會(huì)受到一個(gè)很大的懲罰，以至于我們最終得到的解會(huì)傾向于符合約束。
我們還是以之前的數(shù)字劃分問(wèn)題為例：
假如我們要求第一個(gè)元素25和第四個(gè)元素31要分布在不同的集合中，即要求：x1+x4=1
那么我們可以根據(jù)上表在原目標(biāo)函數(shù)的基礎(chǔ)上要添加懲罰：
不妨假設(shè)P=10000，原目標(biāo)函數(shù)中的Q需要進(jìn)行更改：
那么Q中第一行第一列的元素以及第四行第四列的元素都需要減去10000，第一行第四列的元素以及第四行第一列的元素則需要增加10000。
P*1是一個(gè)常數(shù)可以暫時(shí)忽略，不影響最優(yōu)解。
P = 10000;
Q(1,1) = Q(1,1) - P;
Q(4,4) = Q(4,4) - P;
Q(1,4) = Q(1,4) + P;
Q(4,1) = Q(4,1) + P;
H = Q;
f = zeros(m,1);
xtype = char(join( repmat("B",1,m),""));
Opt = opti('qp',H,f,'xtype',xtype);
[x,fval,exitflag,info] = solve(Opt)
fval * 2 ?+ P最優(yōu)的方案發(fā)生了變化，但是最優(yōu)的目標(biāo)函數(shù)沒(méi)有變化：
大家可以思考下：P值的選擇有什么要求嗎？
（1）P值過(guò)小可能有什么問(wèn)題？
（2）P值過(guò)大可能有什么問(wèn)題？
P值的選取技巧：
英文不好的同學(xué)可以看下面的機(jī)器翻譯：
正如我們所指出的，許多問(wèn)題轉(zhuǎn)換為QUBO形式需要引入一個(gè)懲罰參數(shù)P，而且P必須給出一個(gè)具體的數(shù)值。懲罰參數(shù)P值并不是唯一的，可以選擇許多不同的P值。對(duì)于一個(gè)特定的問(wèn)題，一個(gè)可行的P值通常是根據(jù)領(lǐng)域知識(shí)和需要完成的目標(biāo)來(lái)設(shè)定的。通常，我們對(duì)所有的約束使用相同的懲罰，但是如果有充分的理由區(qū)別對(duì)待不同的約束，那么對(duì)不同的約束使用不同的懲罰也沒(méi)有問(wèn)題。如果一個(gè)約束必須絕對(duì)滿(mǎn)足，即“硬”約束，那么P必須足夠大，以防止違反。然而，有些約束是“軟”的，意味著滿(mǎn)足它們是可取的，但是可以容忍輕微的違反。對(duì)于這樣的情況，一個(gè)更適度的懲罰值就足夠了。?
一個(gè)太大的懲罰值會(huì)阻礙解決過(guò)程，因?yàn)閼土P項(xiàng)壓倒了原始目標(biāo)函數(shù)的信息，使得難以區(qū)分一個(gè)解和另一個(gè)解的質(zhì)量。另一方面，一個(gè)太小的懲罰值會(huì)危及尋找可行解的過(guò)程。一般來(lái)說(shuō)，有一個(gè)相當(dāng)大的“Goldilocks region”，包含了可以很好工作的懲罰值。對(duì)模型進(jìn)行一些初步思考可以得到原始目標(biāo)函數(shù)值的一個(gè)大致估計(jì)。將P取為這個(gè)估計(jì)的一定百分比（75%到150%）通常是一個(gè)好的起點(diǎn)。最后，總是可以檢查生成的解是否可行，從而根據(jù)需要改變懲罰值并進(jìn)行進(jìn)一步的解決過(guò)程，以找到一個(gè)可接受的解。
簡(jiǎn)單點(diǎn)說(shuō)就是提前估計(jì)一個(gè)最優(yōu)的目標(biāo)函數(shù)，然后將P值設(shè)置在這個(gè)值的周?chē)?，如果最后求出?lái)的解違反了你的約束條件，那么你可以嘗試增加P的大小來(lái)提高違反約束的懲罰程度。
思考題：將下面的約束問(wèn)題轉(zhuǎn)換為QUBO問(wèn)題，式中xi均是變量。
有兩個(gè)約束，均需要轉(zhuǎn)換為懲罰函數(shù)：
轉(zhuǎn)換后：
存在線(xiàn)性約束怎么轉(zhuǎn)換？
一樣的加入懲罰項(xiàng)：
例如：
下面給出具體的求解過(guò)程：
如果有線(xiàn)性不等式約束應(yīng)該怎么處理？
這里面涉及到運(yùn)籌學(xué)里面的概念：松弛變量。
它的作用是把不等式約束轉(zhuǎn)換成等式約束，當(dāng)約束條件為“≤”（或者“≥”）類(lèi)型時(shí)的不等式約束時(shí)，可在不等式左邊加上（或者減去）一個(gè)非負(fù)的新變量s，即可將這個(gè)不等式約束轉(zhuǎn)換成等式約束。這個(gè)新增的非負(fù)變量稱(chēng)為松弛變量（或剩余變量），也可統(tǒng)稱(chēng)為松弛變量。在目標(biāo)函數(shù)中一般認(rèn)為新增的松弛變量的系數(shù)為0。
下面我們來(lái)看參考資料中給的例子：
顯然，第一個(gè)約束和第三個(gè)約束都是不等式約束，我們需要進(jìn)行轉(zhuǎn)換：
對(duì)于第一個(gè)約束條件，可以在不等式左邊加上一個(gè)非負(fù)的新變量s1，即可化為等式：2x1+2x2+4x3+3x4+2x5+s1=7，其中s1≥0。
對(duì)于第三個(gè)約束條件，可以在不等式左邊減去一個(gè)非負(fù)的新變量s3，即可化為等式：3x1+3x2+2x3+4x4+4x5-s3=5，其中s3≥0。
由于我們要求解QUBO問(wèn)題，因此這里的s1和s3不能直接寫(xiě)成這種形式，需要將其進(jìn)行轉(zhuǎn)換：轉(zhuǎn)換成由01變量構(gòu)成。
為此我們需要先估計(jì)出s1和s2的上界：
為啥s1的上界是3呢？我們可以將2x1+2x2+4x3+3x4+2x5+s1=7減去第二個(gè)等式約束x1+2x2+2x3+x4+2x5=4，得到：x1+2x3+2x4+s1 = 3，即s1 =3-(x1+2x3+2x4)
又因?yàn)閤1+2x3+2x4理論上的最小值是0，因此可以估計(jì)出s1的上界是3。
為啥s3的上界是6呢？我們將3x1+3x2+2x3+4x4+4x5-s3=5減去兩倍的第二個(gè)等式約束，得到：x1-x2-2x3+2x4-s3 =?-3，那么s3?=x1-x2-2x3+2x4+3，當(dāng)x1和x4為1，x2和x3為0時(shí)，可以估計(jì)出s3的上界是6。
（注意：這里估計(jì)的上界可以更大，比如s3的上界我們還可以這樣估計(jì)：
由3x1+3x2+2x3+4x4+4x5-s3=5得s3 =?3x1+3x2+2x3+4x4+4x5-5，當(dāng)xi全取1時(shí)s3為11，那么為什么我們不將s3的上界估計(jì)成11呢？這是因?yàn)閟3的上界我們給的越大，后續(xù)求解時(shí)你的解空間中也可能遇到更多的情況，會(huì)降低搜索的效率，因此這里松弛變量的上界我們盡可能去找一個(gè)小的滿(mǎn)足原來(lái)約束的，這里需要很多的技巧?。?/p>
接下來(lái)對(duì)每個(gè)松弛變量進(jìn)行了變形，讓他們變成01變量的和：
注意：為啥s1和s3要這樣拆分呢？能不能直接讓s1取為3x6呢？注意，如果你直接讓s1=3x6，那么就限制了s1只能為0或者3，但實(shí)際上s1是可能取為1和2的，因此讓s1=3x6不可??；那能否讓s1=x6+x7+x8呢？答案是肯定的，但這樣多引入了三個(gè)變量，而直接讓s1=x6+2x7既可以保證s1取到0 1 2 3，又可以只引入兩個(gè)變量，可以節(jié)省計(jì)算的開(kāi)銷(xiāo)。
類(lèi)似的，s3 = x8+2x9+4x10也能保證s3取到0 1 2 3 4 5 6，又是引入新的變量最少的一種方案。
后續(xù)的計(jì)算步驟就很簡(jiǎn)單了，引入懲罰項(xiàng)轉(zhuǎn)換為標(biāo)準(zhǔn)的QUBO形式：
二次分配（指派）問(wèn)題
很明顯這里面存在n的平方個(gè)01決策變量，因此當(dāng)n很大時(shí)問(wèn)題的求解規(guī)模很大。下面給出n=3的例子：
這里將雙下標(biāo)轉(zhuǎn)換為單下標(biāo)的形式，然后進(jìn)行建模：
添加懲罰項(xiàng)：
以上介紹來(lái)自參考論文：
GLOVER, FRED, KOCHENBERGER, GARY, DU, YU. Quantum Bridge Analytics I: a tutorial on formulating and using QUBO models[J]. 4OR: Quarterly Journal of the Belgian, French and Italian Operations Research Societies,2019,17(4):335-371. DOI:10.1007/s10288-019-00424-y.
上面這篇文章的內(nèi)容對(duì)于初學(xué)者而言是一個(gè)很好的資料，下面這篇博客是一個(gè)學(xué)者整理的QUBO在經(jīng)典優(yōu)化問(wèn)題上的應(yīng)用，有100多個(gè)經(jīng)典優(yōu)化問(wèn)題，包括最大割、圖著色、旅行商問(wèn)題、背包問(wèn)題等：
https://blog.xa0.de/post/List-of-QUBO-formulations/
里面比較有名的問(wèn)題有：
最大割問(wèn)題（Max-Cut Problem）
圖著色問(wèn)題（Graph Coloring Problem）
旅行商問(wèn)題（Travelling Salesman Problem）
背包問(wèn)題（Knapsack Problem）
作業(yè)調(diào)度問(wèn)題（Job Shop Scheduling Problem）
車(chē)間調(diào)度問(wèn)題（Flow Shop Scheduling Problem）
二次指派問(wèn)題（Quadratic Assignment Problem）
集合覆蓋問(wèn)題（Set Covering Problem）
博客后面給出了這些問(wèn)題的參考文獻(xiàn)：
