割平面法Gommoy分析和总结.docxVIP

割平面法Gommoy

算法步骤:

?A ?b

① 求解原整数规划对应的线性规划minf(x)=cx,.s.t?

?

x

x?0

.xi为整数,设最优解为

x*.

② 如果最优解的分量均为整数,则x*为原整数规划的最优解:否则任选一个x*中不是

整数的分量,设其对应的基变量为x ,定义包含这个基变量的切割约束方程

x ?? rx

p ij j

j

?b

com

,其中x

p,

j为非基变量.

j

③ 令r

ij

?r?[r

ij ij

],b

com

?b

com

?[b

com

],其中[]为高斯函数符号,表示不大于某数的最

大整数. 将切割约束方程变换为x

p

?

??[r]x

ij j

j

?[b

com

]?b

com

?rx

ij j

?j

,由于

0?r

ij

?1,0?b

com

?1,所以有b

?

com

rx

ij j

1,因为自变量为整数,则b

?

com

rx

ij j

也为整

数,所以进一步有b

com

?rx

ij j

=0.

④ 将切割方程加入约束方程中,用对偶单纯算法求解线性规划

? Ax?b

minf(x)?

?

cx,s.t?b

?

??

com

? rx

jij j

j

x?0

?0,转 .

算法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)~=n

disp(基变量的个数要与约束矩阵的行数相等!);mx=NaN;

mf=NaN;return;

end

M=0;

sigma=-[transpose(c)zeros(1,(nVia-length(c)))];xb=b;

while1

[maxs,ind]=max(sigma);ifmaxs=0

vr=find(c~=0,1,last);forl=1:vr

ele=find(base==l,1);if(isempty(ele))

mx(l)=0;

end

elseend

mx(l)=xb(ele);

ifmax(abs(round(mx)-mx))1.0e-7intx=mx;

intf=mx*c;return;

else

sz=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));

end

fori=(infex_x+1):sc

Atmp(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];

%对偶单纯形法迭代过程while1

if(xb)=0

ifmax(abs(round(xb)-xb))1.0e-7

%用对偶单纯形法求得了整数解vr=find(c~=0,l1a,st);

forl=1:vr

ele=find(base==l,1);if(isempty(ele))

mx_1(l)=0;

end

elseend

mx_1(l)=xb(ele);

intx=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==base);

fi=xb(num)-flooor(xb(num));fori=1:(index_x-1)

Atmp(1,i)=a(num,i)-floor(A(num,i));

end

fori=(index_x+1):sc

Atmp(1,i)=a(num,i)-floor(a(num