基本思想
松弛問題:線性規劃
割掉一塊全部都是小數的區域(這一部分取不到整數)
案例
1)橫坐標x1,縱坐標x2
2)藍色小三角形的區域:x2:(1,7/4) x1:(0,3/4)
這塊區域,x1與x2完全取不到整數,所以直接切去
所以,此時取值范圍變化了:
x2<=1把此約束條件帶入,得到x1=1,x2=1,z=2
3)能夠取到整數的區域就不能切掉
引入松弛變量:(解出x1=1,x2=1,z=2的過程)
1)松弛變量:引入之后的效果與原先是一致的
如:-x1+x2<=1 引入x3>=0之后 得到 -x1+x2+x3=1 則此時-x1+x2仍然<=1,所以不影響結果
2)把式子4與5的系數與單個數字拆分為(整數+小數,小數>=0)
即:x1=(1+0)x1 -1/4x3=(-1+3/4)x3 1/4x4=(0+1/4)x4 3/4=(0+3/4)
然后再把整數部分(系數為整數與單個整數)放在左邊,小數部分放在右邊(系數為小數+單個小數)
所以現在變為了
3/4-正數=一個整數
而且0=<3/4<=1 x3,x4>=0
所以 ,3/4-正數<=0
即 3x3+x4>=3
4x2+3x3+x4>=7
基本步驟
引入松弛變量,變不等式為等式
aikxk 松弛變量
aik=[aik]+fik 松弛變量的系數化為正數部分和小數部分
[aik] xk正數部分匯合
fik xk小數部分匯合
[aik]xk -[bi]整數部分放在左側
[bi]+fi 小數部分放在右側
切割平面法流程
案例
解答:
引入松弛變量:
matlab中只有min,所以求最大值要加上負號
matlab代碼
DividePlane.m
function [intx,intf] = DividePlane(A,c,b,baseVector)
%功能：用割平面法求解整數規劃
%調用格式：[intx,intf]=DividePlane(A,c,b,baseVector)
%其中，A：約束矩陣；
% c：目標函數系數向量；
% b：約束右端向量；
% baseVector：初始基向量；
% intx：目標函數取最小值時的自變量值；
% intf：目標函數的最小值；
sz = size(A);
nVia = sz(2);%獲取有多少決策變量
n = sz(1);%獲取有多少約束條件
xx = 1:nVia;if length(baseVector) ~= ndisp('基變量的個數要與約束矩陣的行數相等！');mx = NaN;mf = NaN;return;
endM = 0;
sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
xb = b;%首先用單純形法求出最優解
while 1 [maxs,ind] = max(sigma);
%--------------------用單純形法求最優解--------------------------------------if maxs <= 0 %當檢驗數均小于0時，求得最優解。 vr = find(c~=0 ,1,'last');for l=1:vrele = find(baseVector == l,1);if(isempty(ele))mx(l) = 0;elsemx(l)=xb(ele);endendif max(abs(round(mx) - mx))<1.0e-7 %判斷最優解是否為整數解，如果是整數解。intx = mx;intf = mx*c;return;else %如果最優解不是整數解時，構建切割方程sz = size(A);sr = sz(1);sc = sz(2);[max_x, index_x] = max(abs(round(mx) - mx));[isB, num] = find(index_x == baseVector);fi = xb(num) - floor(xb(num));for i=1:(index_x-1)Atmp(1,i) = A(num,i) - floor(A(num,i));endfor i=(index_x+1):scAtmp(1,i) = A(num,i) - floor(A(num,i));endAtmp(1,index_x) = 0; %構建對偶單純形法的初始表格A = [A zeros(sr,1);-Atmp(1,:) 1];xb = [xb;-fi];baseVector = [baseVector sc+1];sigma = [sigma 0];%-------------------對偶單純形法的迭代過程----------------------while 1%----------------------------------------------------------if xb >= 0 %判斷如果右端向量均大于0，求得最優解if max(abs(round(xb) - xb))<1.0e-7 %如果用對偶單純形法求得了整數解，則返回最優整數解vr = find(c~=0 ,1,'last');for l=1:vrele = find(baseVector == l,1);if(isempty(ele))mx_1(l) = 0;elsemx_1(l)=xb(ele);endendintx = mx_1;intf = mx_1*c;return;else %如果對偶單純形法求得的最優解不是整數解，繼續添加切割方程sz = size(A);sr = sz(1);sc = sz(2);[max_x, index_x] = max(abs(round(mx_1) - mx_1));[isB, num] = find(index_x == baseVector);fi = xb(num) - floor(xb(num));for i=1:(index_x-1)Atmp(1,i) = A(num,i) - floor(A(num,i));endfor i=(index_x+1):scAtmp(1,i) = A(num,i) - floor(A(num,i));endAtmp(1,index_x) = 0; %下一次對偶單純形迭代的初始表格A = [A zeros(sr,1);-Atmp(1,:) 1];xb = [xb;-fi];baseVector = [baseVector sc+1];sigma = [sigma 0];continue;endelse %如果右端向量不全大于0，則進行對偶單純形法的換基變量過程minb_1 = inf;chagB_1 = inf;sA = size(A);[br,idb] = min(xb);for j=1:sA(2)if A(idb,j)<0bm = sigma(j)/A(idb,j);if bm<minb_1minb_1 = bm;chagB_1 = j;endendendsigma = sigma -A(idb,:)*minb_1;xb(idb) = xb(idb)/A(idb,chagB_1);A(idb,:) = A(idb,:)/A(idb,chagB_1);for i =1:sA(1)if i ~= idbxb(i) = xb(i)-A(i,chagB_1)*xb(idb);A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);endendbaseVector(idb) = chagB_1;end%------------------------------------------------------------end %--------------------對偶單純形法的迭代過程--------------------- end else %如果檢驗數有不小于0的，則進行單純形算法的迭代過程minb = inf;chagB = inf;for j=1:nif A(j,ind)>0bz = xb(j)/A(j,ind);if bz<minbminb = bz;chagB = j;endendendsigma = sigma -A(chagB,:)*maxs/A(chagB,ind);xb(chagB) = xb(chagB)/A(chagB,ind);A(chagB,:) = A(chagB,:)/A(chagB,ind);for i =1:nif i ~= chagBxb(i) = xb(i)-A(i,ind)*xb(chagB);A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);endendbaseVector(chagB) = ind;endM = M + 1;if (M == 1000000)disp('找不到最優解！');mx = NaN; minf = NaN;return;end
end
