割平面法Gommoy

割平面法Gommoy算法步骤:①求解原整数规划对应的线性规划minf(x)=cx,.为整数xixbAtsx.0.,设最优解为x*.②如果最优解的分量均为整数,则x*为原整数规划的最优解:否则任选一个x*中不是整数的分量,设其对应的基变量为xp,,定义包含这个基变量的切割约束方程jcomjijbxrpx,其中xj为非基变量.③令][bbb],[rcomcomcomijijijrr,其中[]为高斯函数符号,表示不大于某数的最大整数.将切割约束方程变换为jjijjjijpxrxrxcomcomb][b][,由于10,1r0comijb,所以有jijcomxrb1,因为自变量为整数,则jijcomxrb也为整数,所以进一步有jijcomxrb=0.④将切割方程加入约束方程中,用对偶单纯算法求解线性规划00bbAx.,)(mincomxxrtscxxfjjij,转.算法MATLAB实现代码:function[intx,intf]=Gomory(A,c,b,base)%约束矩阵:A;%目标函数系数向量:c%约束右端向量:b%初始基向量base%目标函数取最小化时的自变量值:x%目标函数的最小值:minfsz=size(A);nVia=sz(2);n=sz(1);xx=1:nVia;iflength(base)~=ndisp('基变量的个数要与约束矩阵的行数相等!');mx=NaN;mf=NaN;return;endM=0;sigma=-[transpose(c)zeros(1,(nVia-length(c)))];xb=b;while1[maxs,ind]=max(sigma);ifmaxs=0vr=find(c~=0,1,'last');forl=1:vrele=find(base==l,1);if(isempty(ele))mx(l)=0;elsemx(l)=xb(ele);endendifmax(abs(round(mx)-mx))1.0e-7intx=mx;intf=mx*c;return;elsesz=size(A);sr=sz(1);sc=sz(2);[max_x,index_x]=max(ads(round(mx)-mx));[isB,num]=find(index_x==base);fi=xb(num)-floor(xb(num));fori=1:(index_x-1)Atmp(1,i)=A(num,i)-floor(A(num,i));endfori=(infex_x+1):scAtmp(1,i)=A(num,i)-floor(A(num,i));end%构建对单纯形法的初始表格Atmp(1,index_x)=0;A=[Azeros(sr,1);-Atmp(1,:)1];xb=[xb;-fi];base=[basesc+1];sigma=[sigma0];%对偶单纯形法迭代过程while1if(xb)=0ifmax(abs(round(xb)-xb))1.0e-7%用对偶单纯形法求得了整数解vr=find(c~=0,1,'last');forl=1:vrele=find(base==l,1);if(isempty(ele))mx_1(l)=0;elsemx_1(l)=xb(ele);endendintx=mx_1;intf=mx_1*c;return;elsesz=size(A);sr=sz(1);sc=sz(2);[max_x,index_x]=max(abs(round(mx_1)-mx_1));[isB,num]=find(index_x==base);fi=xb(num)-flooor(xb(num));fori=1:(index_x-1)Atmp(1,i)=a(num,i)-floor(A(num,i));endfori=(index_x+1):scAtmp(1,i)=a(num,i)-floor(a(num,i));end%下一次对偶单纯形法迭代的初始表格Atmp(1,index_x)=0;A=[Azeros(sr,1);-Atmp(1,:)1];xb=[xb;-fi];base=[basesc+1];sigma=[sigma0];continue;end%对偶单纯形法的换基变量过程elseminb_1=inf;chagB_1=inf;sA=size(A);[br,idb]=min(xb);forj=1:sA(2)ifA(idb,j)0bm=sigma(j)/A(idb,j);ifbmminb_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);fori=1:sA(1)ifi~=idbxb(i)=xb(i)-A(i,chagB_1)*xb(idb);A(i,:)=A(i,:)-A(i,chagB_1)*A(idb,:);endendbase=chagB_1;endendendelseminb=inf;chagB=inf;forj=1:nifA(j,ind)0bz=xb(j)/A(j,ind);ifbzminbminb=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);fori=1:nifi~=chagBxb(i)=xb(i)-A(i,ind)*xb(chagB);A(i,:)=A(i,:)-A(i,ind)*A(chagB,:);endendbase(chagB)=ind;endM=M+1;if(M==1000000)disp('找不到最优解!');mx=NaN;minf=NaN;return;endend算法举例为整数且2121212121,,0,4222.,)(infmxxxxxxxxtsxxx解:首先引入人工变量x3,x4,将约束条件化为等式形式为整数且21432142132121,,0,,,422x2.,)(infmxxxxxxxxxxxtsxxx在MATLAB命令行输入下列命令:A=[-1210;2101];C=[1;-1];B=[2;4];[intx,intf]=Gomory(A,c,b,[34])结果为:intx=01Intf=-1