资源描述
整数计划分支定界法MATLAB程序
1.这种方法绝对能都解出答案, 而且答案正确
function [x,val]=fzdj(n,f,a,b,aeq,beq,lb,ub)
x=zeros(n,1);
x1=zeros(n,1);
m1=2;
m2=1;
[x1,val1]=linprog(f,a,b,aeq,beq,lb,ub);
if (x1==0)
x=x1;
val=val1;
elseif (round(x1)==x1)
x=x1;
val=val1;
else
e1={0,a,b,aeq,beq,lb,ub,x1,val1};
e(1,1)={e1};
zl=0;
zu=-val1;
while (zu~=zl)
for c=1:1:m2
if (m1~=2)
if (cell2mat(e{m1-1,c}(1))==1)
e1={1,[],[],[],[],[],[],[],0};
e(m1,c*2-1)={e1};
e(m1,c*2)={e1};
continue;
end;
end;
x1=cell2mat(e{m1-1,c}(8));
x2=zeros(n,1);
s=0;
s1=1;
s2=1;
lb1=cell2mat(e{m1-1,c}(6));
ub1=cell2mat(e{m1-1,c}(7));
lb2=cell2mat(e{m1-1,c}(6));
ub2=cell2mat(e{m1-1,c}(7));
for d=1:1:n
if (abs((round(x1(d))-x1(d)))>0.0001)&(s==0)
s=1;
lb1(d)=fix(x1(d))+1;
if (a*lb1<=b)
s1=0;
end;
ub2(d)=fix(x1(d));
if (a*lb2<=b)
s2=0;
end;
end;
end;
e1={s1,a,b,aeq,beq,lb1,ub1,[],0};
e2={s2,a,b,aeq,beq,lb2,ub2,[],0};
e(m1,c*2-1)={e1};
e(m1,c*2)={e2};
end;
m1=m1+1;
m2=m2*2;
for c=1:1:m2
if (cell2mat(e{m1-1,c}(1))==0)
[x1,val1]=linprog(f,cell2mat(e{m1-1,c}( 2)),cell2mat(e{m1-1,c}(3)),cell2mat(e{m1-1,c}(4)),cell2mat(e{m1-1,c}(5)),cell2mat(e{m1-1,c}(6)),cell2mat(e{m1-1,c}(7)));
e1={cell2mat(e{m1-1,c}(1)),cell2mat(e{m1-1,c}(2)),cell2mat(e{m1-1,c}(3)),cell2mat(e{m1-1,c}(4)),cell2mat(e{m1-1,c}(5)),cell2mat(e{m1-1,c}(6)),cell2mat(e{m1-1,c}(7)),x1,val1};
e(m1-1,c)={e1};
end;
z=val1;
if ((-z)<(-zl))
e1={1,[],[],[],[],[],[],[],0};
e(m1-1,c)={e1};
elseif (abs(round(x1)-x1)<=0.0001)
zl=z;
end;
end;
for c=1:1:m2
if (cell2mat(e{m1-1,c}(1))==0)
zu=cell2mat(e{m1-1,c}(9));
end;
end;
for c=1:1:m2
if (-cell2mat(e{m1-1,c}(9))>(-zu))
zu=cell2mat(e{m1-1,c}(9));
end;
end;
end;
for c=1:1:m2
if (cell2mat(e{m1-1,c}(1))==0)&(cell2mat(e{m1-1,c}(9))==zu)
x=cell2mat(e{m1-1,c}(8));
end;
end;
val=zu;
end;
2.这种方法是书本上程序, 不过不能解题, 期望高手能将它改善
function [x,y]=IntLp(f,G,h,Geq,heq,lb,ub,x,id,options)
%整数线性计划分支定界法, 可求解全整数线性或混合整数线性计划
%y=min f'*x subject to:G*x<=h Geq*x=heq x为全整数
%数或混合整数列向量
%使用方法
% [x,y]=IntLp(f,G,h)
% [x,y]=IntLp(f,G,h,Geq,heq)
% [x,y]=IntLp(f,G,h,Geq,heq,lb,ub)
% [x,y]=IntLp(f,G,h,Geq,heq,lb,ub,x)
% [x,y]=IntLp(f,G,h,Geq,heq,lb,ub,x,id)
% [x,y]=IntLp(f,G,h,Geq,heq,lb,ub,x,id,options)
%参数说明
% x:最优解列向量; y:目标函数最小值; f:目标函数系数列向量
% G:约束不等式条件系数矩阵; h:约束不等式条件右端列向量
% Gep: 约束等式条件系数矩阵; heq:约束等式条件右端列向量
% lb:解下界列向量(Default:-inf)
% ub:解得上界列向量(Default:inf)
% x:迭代初值列向量;
% id:整数变量指标列向量, 1—整数, 0—实数(default:1)
% options设置请参见 optimset或linprog
%例 min z=x1+4x2
%s.t. 2x1+x2<=8
% x1+2x2>=6
% x1,x2>=0且为整数
%先将 x1+x2>=6 化为 -x1-2x2<=-6
%[x,y]=IntLp([1;4],[2 1;-1 -2],[8;-6],[],[],[0;0])
global upper opt c x0 A b Aeq beq ID options;
if nargin<10,options=optimset({});options.Display='off';
options.LargeScale='off';end
if nargin<9,id=ones(size(f));end
if nargin<8,x=[];end
if nargin<7|isempty(ub),ub=inf*ones(size(f));end
if nargin<6|isempty(lb),lb=zeros(size(f));end
if nargin<5,heq=[];end
if nargin<4,Geq=[];end
upper=inf;c=f;x0=x;A=G;b=h;Aeq=Geq;beq=heq;ID=id;
ftemp=IntLp(lb(:),ub(:));
x=opt;y=upper;
%以下为子函数
function ftemp=IntLp(vlb,vub)
global upper opt c x0 A b Aeq beq ID options;
[x,ftemp,how]=linprog(c,A,b,Aeq,beq,vlb,vub,x0,options);
if how<=0
return;
end;
if ftemp-upper>0.00005%in order to avoid error
return;
end;
if max(abs(x.*ID-round(x.*ID)))<0.00005
if upper-ftemp>0.00005%in order to avoid error
opt=x';upper=ftemp;
return;
else
opt=[opt;x'];
return;
end;
end;
notintx=find(abs(x-round(x))>0.00005);
%in order to avoid error
intx=fix(x);tempvlb=vlb;tempvub=vub;
if vub(notintx(1,1),1)>=intx(notintx(1,1),1)+1
tempvlb(notintx(1,1),1)=intx(notintx(1,1),1)+1;
ftemp=IntLp(tempvlb,vub);
end;
if vlb(notintx(1,1),1)<=intx(notintx(1,1),1)
tempvub(notintx(1,1),1)=intx(notintx(1,1),1);
ftemp=IntLp(vlb,tempvub);
end;
3.这种方法是我网上搜, 也不能解题, 期望高手改善
function [y,fval]=BranchBound(c,A,b,Aeq,beq)
NL=length(c);
UB=inf;
LB=-inf;
FN=[0];
AA(1)={A};
BB(1)={b};
k=0;
flag=0;
while flag==0;
[x,fval,exitFlag]=linprog(c,A,b,Aeq,beq);
if (exitFlag == -2) | (fval >= UB) FN(1)=[];
if isempty(FN)==1
flag=1;
else
k=FN(1);
A=AA{k};
b=BB{k};
end
else
for i=1:NL
if abs(x(i)-round(x(i)))>1e-7
kk=FN(end);
FN=[FN,kk+1,kk+2];
temp_A=zeros(1,NL);
temp_A(i)=1;
temp_A1=[A;temp_A];
AA(kk+1)={temp_A1};
b1=[b;fix(x(i))];
BB(kk+1)={b1};
temp_A2=[A;-temp_A];
AA(kk+2)={temp_A2};
b2=[b;-(fix(x(i))+1)];
BB(kk+2)={b2};
FN(1)=[];
k=FN(1);
A=AA{k};
b=BB{k};
break;
end
end
if (i==NL) & (abs(x(i)-round(x(i)))<=1e-7)
UB=fval;
y=x;
FN(1)=[];
if isempty(FN)==1
flag=1;
else
k=FN(1);
A=AA{k};
b=BB{k};
end
end
end
end
y=round(y);
fval=c*y;
以上程序都是我网上搜到, 还有自己从书本上自己敲, 第一个程序能成功解出答案, 后两个不能, 可能是程序有问题, 可能是我不会调用, 期望高手指点, 帮助更多好友们能够应用这三个程序! ! !
展开阅读全文