整數規划:
clc,clear; c = [-40;-90]; A = [9 7;7 20]; b = [56;70]; lb = zeros(2,1); [x,fval]= intlinprog(c,1:2,A,b,[],[],lb); fval = -fval x
分支定界法或者割平面法求解純或者混合整數線性規划問題;
輸出:
當條件A,B之間不是且關系而是或的時候:
固定成本問題(最優化函數中含有與xi無關的常量,相當於固定成本,優化函數可以寫成總固定成本加上總可變成本之和):
0-1整數規划問題(過濾隱枚舉法,分枝隱枚舉法)
指派問題(0-1規划特殊情形:匈牙利法)
蒙特卡洛法(求解各種類型規划)
下面主要介紹蒙特卡洛法(隨機取樣法):
例題:
如果用顯枚舉法試探,需要計算1010個點,計算量巨大。但是用蒙特卡洛去計算106個點便可以找到滿意解。
前提:整數規划的最優點不是孤立的奇點;
而采集106個點后,我們有很大把握最優值點在106個點之中;
function [f,g] = mengte(x); f = x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)-8*x(1)-2*x(2)-3*x(3)-... x(4)-2*x(5); g = [sum(x)-400 x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800 2*x(1)+x(2)+6*x(3)-200 x(3)+x(4)+5*x(5)-200];
rand('state',sum(clock)); p0 = 0; tic for i = 1:10^6 x = 99*rand(5,1); x1 = floor(x); x2 = ceil(x); [f,g] = mengte(x1); if sum(g<=0)==4 if p0<=f x0 = x1;p0=f; end end [f,g] = mengte(x2); if sum(g<=0)==4 if p0 <= f x0 = x2; p0 = f; end end end x0,p0 toc
輸出:
蒙特卡洛法得到的解為最優解的近似解,10^6個數據已經用了將近7s的時間,所以如果增加十倍,可能得70s時間才能得到結果。
蒙特卡洛法對計算器的計算能力要求很高。
lingo軟件可以求得精確的全局最優解:
程序如下:
model: sets: !初始化變量集合; row/1..4/:b; !b為長度為4的行向量; col/1..5/:c1,c2,x; !c1,c2,x 為長度為5的列向量; link(row,col):a; !a為行列分別為4,5的系數矩陣; endsets data: !實例化變量; c1 = 1,1,3,4,2 c2 = -8,-2,-3,-1,-2; a = 1 1 1 1 1 1 2 2 1 6 2 1 6 0 0 0 0 1 1 5; b = 400,800,200,200; !滿足a*x <= b; enddata max = @sum(col:c1*x^2+c2*x); !優化函數; @for(row(i):@sum(col(j):a(i,j)*x(j))<b(i)); @for(col:@gin(x)); !產生隨機的x; @for(col:@bnd(0,x,99)); !x定義邊界;
整數規划問題的求解可以使用Lingo等專用軟件,對於一般的整數規划問題,無法直接利用matlab的函數;
必須利用Matlab編程實現分支定界解法和割平面解法。
對於指派問題等0-1整數規划問題,可以直接利用Matlab的函數intlinprog求解;
c=[3 8 2 10 3;8 7 2 9 7;6 4 2 7 5 8 4 2 3 5;9 10 6 9 10]; c=c(:); a=zeros(10,25); for i=1:5 a(i,(i-1)*5+(1:5))=1; a(5+i,i:5:25)=1; end b=ones(10,1); lb=zeros(25,1); ub=ones(25,1); [x,y]=intlinprog(c,(1:25),[],[],a,b,lb,ub); x=reshape(x,[5,5]),y
所以 x15 = x23 = x32 = x44 = x51 = 1, 最優值為21
目標函數中含有分段線性函數時候如何求解問題呢?
可以看出最優化函數中含有分段線性函數c(x), 然后寫出約束條件如下;
甲、乙兩種汽油含有原油A的最低比例分別為0.5和0.6;
現在的問題關鍵在於如何處理分段線性函數c(x)。
可以有三種解決方法:
解法一:
非線性規划模型用Lingo軟件更合適;
程序如下:
model: sets: var1/1..4/:y; !這里y(1)=x11,y(2)=x21,y(3)=x12,y(4)=x22; var2/1..3/:x,c; endsets max=4.8*(y(1)+y(2))+5.6*(y(3)+y(4))-@sum(var2:c*x); y(1)+y(3)<@sum(var2:x)+500; y(2)+y(4)<1000; 0.5*(y(1)-y(2))>0; 0.4*y(3)-0.6*y(4)>0; (x(1)-500)*x(2)=0; (x(2)-500)*x(3)=0; @for(var2:@bnd(0,x,500)); data: c=10 8 6; enddata end
!可以用菜單命令“LINGO|Options”在“Global Solver”選項卡上啟動全局優化選 型,並運行上述程序求得全局最有解:購買 1000 噸原油 A ,與庫存的 500 噸原油 A 和 1000 噸原油 B 一起,共生產 2500 噸汽油乙,利潤為 5000(千元)。;
解法二:
通過z1, z2, z3可以把x1, x2, x3限制在0到500之間,且滿足x3 <= 500z3 <= x2 <= 500z2 <= x1 <= 500z1;
這種解法可以用matlab求解;
程序如下:
c = [4.8 5.6 4.8 5.6 0 -10 -8 -6 0 0 0]'; A = [1 1 0 0 -1 zeros(1,6); 0 0 1 1 zeros(1,7); 0 0 0 0 1 zeros(1,6); -1 0 1 zeros(1,8); 0 -2 0 3 zeros(1,7); zeros(1,5) -1 0 0 0 500 0; zeros(1,5) 1 0 0 -500 0 0; zeros(1,6) -1 0 0 0 500; zeros(1,6) 1 0 0 -600 0; zeros(1,7) 1 0 0 -500]; b = [500 1000 1500 0 0 0 0 0 0 0]'; aeq = [zeros(1,4) -1 1 1 1 zeros(1,3)]; beq = [0]; %sum(x(6:8)) = x(5); lb = zeros(11,1); ub = [inf*ones(1,8) 1 1 1]'; [x, y] = intlinprog(-c,(1:11),A,b,aeq,beq,lb,ub);
format short g x = x(1:5),y
輸出為:
與解法一的結果一致。
解法三:將該分段函數直接畫出來,可見:
其中50改成500。不等式(22)可以保證選擇的那一段zi處的w滿足要求。
所以我們在x1 x2 x3 x4 x 上加了w1 w2 w3 w4 z1 z2 z3 其中z1 z2 z3為整數0-1,w1 w2 w3 w4為0~1之間的小數,與解法二相同,
也可以用matlab解出;但是這種方法多加了七個量,相對復雜,但是有一定的普遍性,所有含線性分段函數的整數線性規划問題均
可以用這種方法。