Matlab學(xué)習(xí)系列26.整數(shù)規(guī)劃

26.整數(shù)規(guī)劃全部變量限制為整數(shù)的規(guī)劃問題，稱為純整數(shù)規(guī)劃；部分變量限制為整數(shù)的規(guī)劃問題，稱為混合整數(shù)規(guī)劃；變量只取0或1的規(guī)劃問題，稱為0-1整數(shù)規(guī)劃。整數(shù)規(guī)劃問題，建議使用Lingo軟件求解。常用的整數(shù)規(guī)劃問題解法有：（1）分枝定界法：可求純或混合整數(shù)線性規(guī)劃；（2）割平面法：可求純或混合整數(shù)線性規(guī)劃；（3）隱枚舉法：用于求解0-1整數(shù)規(guī)劃，有過濾法和分枝法；（4）匈牙利法：解決指派問題（0-1規(guī)劃特殊情形）；（5）蒙特卡羅法：求解各種類型規(guī)劃。一、分枝定界法分支定界法的基本思想是：設(shè)有最大化的整數(shù)規(guī)劃問題A，先解與之相應(yīng)的線性規(guī)劃問題B，若B的最優(yōu)解不符合A的整數(shù)條件，那么B的最優(yōu)目標函數(shù)必是A的最優(yōu)目標函數(shù)z*的上界，記作z2,而A的任意可行解的目標函數(shù)值將是z*的一個下界z1,分支定界法就是將B的可行域分成子區(qū)域（稱為分支）的方法，逐步減小z2和增大z1,最終求到z*。例1分枝定界法原理示例：用Lingo軟件求解：代碼：max5x1+8x2stx1+x2<=65x1+9x2<=45endgin2運行結(jié)果：Globaloptimalsolutionfound.Objectivevalue:40.00000Objectivebound:40.00000Infeasibilities:0.000000Extendedsolversteps:0Totalsolveriterations:0VariableValueReducedCostX10.000000-5.000000X25.000000-8.000000RowSlackorSurplusDualPrice140.000001.00000021.0000000.00000030.0000000.000000二、0-1整數(shù)規(guī)劃變量xi只能取值0,1，該約束條件可表示為：0≤xi≤1,xi∈N或xi(1-xi)=01.隱枚舉法例2求解下列0-1規(guī)劃問題：求解思路：（1）先試探性地求一個可行解，易看出(x1,x2,x3)=(1,0,0)滿足約束條件，故是一個可行解，相應(yīng)的目標函數(shù)值為z=3.（2）由于是求極大值，故目標值z<3的解，不必檢驗是否滿足約束條件即可刪除，于是可增加一個約束條件（稱為過濾條件）：(e)（3）用全部枚舉法，3個變量共23=8種可能的組合，用過濾條件（并計算目標函數(shù)值，不斷改進過濾條件）篩選每個可能的組合，最終得到問題的最優(yōu)解。(x1,x2,x3)目標值約束條件過濾條件(a)(b)(c)(d)(e)(0,0,0)0×(1,0,0)3√√√√√3x1-2x2+x3≥3(0,1,0)-2×(0,0,1)5√√√√√3x1-2x2+x3≥5(1,1,0)1×(1,0,1)8√√√√√3x1-2x2+x3≥8(1,1,1)6×(0,1,1)3×從而得到最優(yōu)解為(1,0,1),最優(yōu)值為z*=8.用Lingo軟件求解：代碼：max3x1-2x2+5x3stx1+2x2-x3<=2x1+4x2+x3<=4x1+x2<=34x2+x3<=6endint3運行結(jié)果：Globaloptimalsolutionfound.Objectivevalue:8.000000Objectivebound:8.000000Infeasibilities:0.000000Extendedsolversteps:0Totalsolveriterations:0VariableValueReducedCostX11.000000-3.000000X20.0000002.000000X31.000000-5.000000RowSlackorSurplusDualPrice18.0000001.00000022.0000000.00000032.0000000.00000042.0000000.00000055.0000000.0000002.指派問題某公司準備分派n個工人X1,…,Xn做n件工作Y1,…,Yn,但由于任務(wù)性質(zhì)和各人專長不同，因此各人完成每件工作的效率也就不同，于是產(chǎn)生一個問題：應(yīng)指派哪個人去完成哪件工作，使得工作效率最高？用效率矩陣表示：其中，元素cij≥0（i,j=1,…,n）表示指派第i人去完成第j項任務(wù)時的效率（或時間、成本等）。令則指派問題的一般數(shù)學(xué)模型為：對于上述問題，其實通俗來說，就是n×n矩陣中，選取n個元素，每行每列各1個元素，使得和最小。根據(jù)指派問題的特點，可以使用匈牙利算法來進行求解。匈牙利算法原理（匈牙利數(shù)學(xué)家狄·康尼格證明的兩個定理）：定理1如果從指派問題效率矩陣[cij]的每一行元素中分別減去（或加上）一個常數(shù)ui（稱為該行的位勢），從每一列分別減去（或加上）一個常數(shù)vj（稱為該列的位勢），得到一個新的效率矩陣[bij],若其中bij=cij-ui-vj,則[bij]的最優(yōu)解的結(jié)果等價于[cij]的最優(yōu)解的結(jié)構(gòu)。因此，指派問題的最優(yōu)解具有性質(zhì)：若從矩陣的一行（列）各元素中分別減去該行（列）的最小元素，得到規(guī)約矩陣，其最優(yōu)解和原矩陣的最優(yōu)解相同。經(jīng)過減去行列最小元素的轉(zhuǎn)化之后的效率矩陣[bij]中，把位于不同行不同列的0元素，簡稱為獨立0元素。目的便是找出位于不同行不同列的n個獨立的0元素，從而也就得到了原問題的最優(yōu)解。定理2若效率矩陣C的元素看可分成“0”和“非0”兩部分，則覆蓋所有“0”元素的最少直線數(shù)=獨立的“0”元素的最多個數(shù)。例3.設(shè)指派問題的效率矩陣為Lingo代碼：model:sets:var/1..5/;link(var,var):c,x;endsetsdata:c=3821038729764275842359106910;enddatamin=@sum(link:c*x);@for(var(i):@sum(var(j):x(i,j))=1);@for(var(j):@sum(var(i):x(i,j))=1);@for(link:@bin(x));end運行結(jié)果（部分）：Globaloptimalsolutionfound.Objectivevalue:21.00000Objectivebound:21.00000Infeasibilities:0.000000Extendedsolversteps:0Totalsolveriterations:0VariableValueReducedCostX(1,1)0.0000003.000000X(1,2)0.0000008.000000X(1,3)0.0000002.000000X(1,4)0.00000010.00000X(1,5)1.0000003.000000X(2,1)0.0000008.000000X(2,2)0.0000007.000000X(2,3)1.0000002.000000X(2,4)0.0000009.000000X(2,5)0.0000007.000000X(3,1)0.0000006.000000X(3,2)1.0000004.000000X(3,3)0.0000002.000000X(3,4)0.0000007.000000X(3,5)0.0000005.000000X(4,1)0.0000008.000000X(4,2)0.0000004.000000X(4,3)0.0000002.000000X(4,4)1.0000003.000000X(4,5)0.0000005.000000X(5,1)1.0000009.000000X(5,2)0.00000010.00000X(5,3)0.0000006.000000X(5,4)0.0000009.000000X(5,5)0.00000010.00000RowSlackorSurplusDualPrice121.00000-1.00000020.0000000.000000三、蒙特卡羅法（隨機取樣法）前面的方法，主要是針對線性整數(shù)規(guī)劃而言，對于非線性整數(shù)規(guī)劃沒有通用的有效解法。整數(shù)規(guī)劃由于限制變量是整數(shù)，增加了求解難度，但整數(shù)解是有限個，所以可以采用枚舉法。當枚舉個數(shù)很多時，顯性枚舉是不現(xiàn)實的，但利用蒙特卡羅隨機取樣法，在一定的計算量下是可以得到滿意解的。例4.求解如下非線性整數(shù)規(guī)劃問題：若用顯枚舉法，共需1005=1010個點，計算量非常大。但用蒙特卡羅法隨機計算106個點，便可找到滿意解。代碼：目標函數(shù)：function[f,g]=fun1(x)f=x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)^2-8*x(1)-2*x(2)-3*x(3)-x(4)-2*x(5);g(1)=sum(x)-400;g(2)=x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800;g(3)=2*x(1)+x(2)+6*x(3)-200;g(4)=x(3)+x(4)+5*x(5)-200;主程序：rand('state',sum(clock));p0=0;ticfori=1:10^6x=99*rand(5,1);x1=floor(x);x2=ceil(x);[f,g]=fun1(x1);ifsum(g<=0)==4ifp0<=fx0=x1;p0=f;endend[f,g]=fun1(x2);ifsum(g<=0)==4ifp0<=fx0=x2;p0=f;endendendx0,p0toc運行結(jié)果：x0=43941995p0=49298Elapsedtimeis45.494952seconds.注意：由于是隨機取樣，故每次的運行結(jié)果可能不一樣。利用lingo軟件可以求得全局最優(yōu)解：代碼：model:sets:row/1..4/:b;col/1..5/:c1,c2,x;link(row,col):a;endsetsdata:c1=1,1,3,4,2;c2=-8,-2,-3,-1,-2;a=11111122162160000115;b=400,800,200,200;enddatamax=@sum(col:c1*x^2+c2*x);@for(row(i):@sum(col(j):a(i,j)*x(j))<b(i));@for(col:@gin(x));@for(col:@bnd(0,x,99));End運行結(jié)果（部分）：Localoptimalsolutionfound.Objectivevalue:51568.00Objectivebound:51568.00Infeasibilities:0.000000Extendedsolversteps:0Totalsolveriterations:51VariableValueReducedCostX(1)50.00000-92.00000X(2)99.00000-196.0000X(3)0.0000003.000000X(4)99.00000-791.0000X(5)20.00000-78.00000RowSlackorSurplusDualPrice151568.001.0000002132.00000.0000003333.00000.00000041.0000000.00000051.0000000.000000四、混合整數(shù)規(guī)劃（配送選址問題）圖1三級配送的網(wǎng)絡(luò)結(jié)構(gòu)物流配送中心選址問題，是在給定某一地區(qū)所有備選點的地址集合中選出一定數(shù)目的地址建立配送中心，從而建立一系列的配送區(qū)域，以實現(xiàn)選出點建立的配送中心與各需求點和工廠(供貨點)形成的配送系統(tǒng)總物流費用最小?？傎M用=運輸費用+配送費用+倉儲費用+固定成本設(shè)有m個供應(yīng)工廠，n個客戶，l個備選配送中心；pk表示第k個工廠的產(chǎn)量(k=1,…,m);dj表示第j個客戶的需求量(j=1,…,n);gi表示第i個配送中心的單位管理成本(i=1,…,l);fi表示建立第i個配送中心的固定成本；zi指示是否選中第i個配送中心，1表示選中，0表示未選中；zmax表示配送中心的最大限制數(shù)目；cki表示第k個供應(yīng)工廠到第i個配送中心的運輸單價；hij表示第i個配送中心到第j個客戶的配送單價；wki表示第k個供應(yīng)工廠到第i個配送中心的運量；xij表示第i個配送中心到第j個客戶的運量。則配送選址問題數(shù)學(xué)模型的一般形式如下：例5設(shè)有4個備選物流配送中心地址，6個工廠為其供貨，6個客戶需要產(chǎn)品，最多設(shè)置3個物流配送中心。工廠到物流配送中心的運輸價格見表1，物流配送中心到客戶的運輸價格見表2，工廠的總生產(chǎn)能力見表3，物流配送中心的固定成本、單位管理成本，及容量見表4，客戶的需求量見表5：表1工廠到配送中心的運輸價格w1w2w3w4p16542p22349p36875p47423p54251p63417表2配送中心到客戶的運輸價格c1c2c3c4c5c6w1327475w2614253w3245368w4563746表3工廠的總生產(chǎn)能力工廠p1p2p3p4p5p6總生產(chǎn)能力(p)40,00050,00060,00070,00060,00040,000表4備選物流配送中心的固定成本，單位管理成本，容量物流配送中心w1w2w3w4固定成本(f)500,000300,000400,000400,000單位管理成本(g)3254倉庫容量(a)10,00060,00070,00050,000表5客戶的需求量顧客c1c2c3c4c5c6需求(d)10,00020,00010,00020,00030,00010,000利用Lingo軟件求解以上混合整數(shù)規(guī)劃。代碼：model:sets:factory/p1..p6/:p;warhouse/w1..w4/:a,f,g;customer/c1..c6/:d;tr/tr1..tr4/:z;link1(factory,warhouse):c,w;link2(warhouse,customer):h,x;endsetsdata:p=40000,50000,60000,70000,60000,40000;a=70000,60000,70000,50000;f=500000,300000,400000,400000;g=3,2,5,4;d=10000,20000,10000,20000,30000,10000;c=654223496875742342513417;H=327475614253245368563746;enddatamin=@sum(link1(k,i):c(k,i)*w(k,i))+@sum(link2(i,j):h(i,j)*x(i,j))+@sum(link1(k,i):g(i)*w(k,i))+@sum(warhouse(i):f(i)*z(i));@for(factory(k):@sum(link1(k,i):w(k,i))<=p(k));@for(warhouse(i):@sum(link2(i,j):x(i,j))=@sum(link1(k,i):w(k,i)));@for(customer(j):@sum(link2(i,j):x(i,j))>=d(j));@for(warhouse(i):@sum(link1(k,i):w(k,i))<=(a(i)*z(i)));@sum(tr(i):z(i))<=3;@for(tr(i):@bin(z));End運行結(jié)果（部分）：Globaloptimalsolutionfound.Objectivevalue:1480000.Objectivebound:1480000.Infeasibilities:0.000000Extendedsolversteps:4Totalsolveriterations:50VariableValueReducedCostZ(TR1)0.000000290000.0Z(TR2)1.000000300000.0Z(TR3)0.000000190000.0Z(TR4)1.000000400000.0C(P1,W1)6.0000000.000000C(P1,W2)5.0000000.000000C(P1,W3)4.0000000.000000C(P1,W4)2.0000000.000000C(P2,W1)2.0000000.000000C(P2,W2)3.0000000.000000C(P2,W3)4.0000000.000000C(P2,W4)9.0000000.000000C(P3,W1)6.0000000.000000C(P3,W2)8.0000000.000000C(P3,W3)7.0000000.000000C(P3,W4)5.0000000.000000C(P4,W1)7.0000000.000000C(P4,W2)4.0000000.000000C(P4,W3)2.0000000.000000C(P4,W4)3.0000000.000000C(P5,W1)4.0000000.000000C(P5,W2)2.0000000.000000C(P5,W3)5.0000000.000000C(P5,W4)1.0000000.000000C(P6,W1)3.0000000.000000C(P6,W2)4.0000000.000000C(P6,W3)1.0000000.000000C(P6,W4)7.0000000.000000W(P1,W1)0.0000004.000000W(P1,W2)0.0000002.000000W(P1,W3)0.0000003.000000W(P1,W4)40000.000.000000W(P2,W1)0.0000000.000000W(P2,W2)0.0000000.000000W(P2,W3)0.0000003.000000W(P2,W4)0.0000007.000000W(P3,W1)0.0000004.000000W(P3,W2)0.0000005.000000W(P3,W3)0.0000006.000000W(P3,W4)0.0000003.000000W(P4,W1)0.0000005.000000W(P4,W2)0.0000001.000000W(P4,W3)0.0000001.000000W(P4,W4)0.0000001.000000W(P5,W1)0.0000003.000000W(P5,W2)60000.000.000000W(P5,W3)0.0000005.000000W(P5,W4)0.0000000.000000W(P6,W1)0.0000001.000000W(P6,W2)0.0000001.000000W(P6,W3)0.0000000.000000W(P6,W4)0.0000005.000000H(W1,C1)3.0000000.000000H(W1,C2)2.0000000.000000H(W1,C3)7.0000000.000000H(W1,C4)4.0000000.000000H(W1,C5)7.0000000.000000H(W1,C6)5.0000000.000000H(W2,C1)6.0000000.000000H(W2,C2)1.0000000.000000H(W2,C3)4.0000000.000000H(W2,C4)2.0000000.000000H(W2,C5)5.0000000.000000H(W2,C6)3.0000000.000000H(W3,C1)2.0000000.000000H(W3,C2)4.0000000.000000H(W3,C3)5.0000000.000000H(W3,C4)3.0000000.000000H(W3,C5)6.0000000.000000H(W3,C6)8.0000000.000000H(W4,C1)5.0000000.000000H(W4,C2)6.0000000.000000H(W4,C3)3.0000000.000000H(W4,C4)7.0000000.000000H(W4,C5)4.0000000.000000H(W4,C6)6.0000000.000000X(W1,C1)