分支定界法的Matlab实现.doc

一个调用例子：
ifint=[0 1];
f=[10 9];
a=[1 0;0 1;-5 -3];
b=[8 10 -45];
[x,fval,exitflag] = linprogdis(ifint,f,a,b,[],[],[],[],[],[])

function r=checkint(x)
%算法：如果x(i)是整数，则返回r(i)=1；否则返回r(i)=0
function r=ifrowinmat(arow,amat)
%输入参数：
% arow 向量，
% amat 矩阵
%
%设计：如果 arow与矩阵amat中的某一行相等，则返回1，如果找不到现等的一行，则返回0

可以使用ismember(arrow,amat,’rows’)替换ifrowinmat的调用，2005-10-28标注
使用时，将下面的代码存入文件：linprogdis.m
function [x,fval,exitflag,output,lambda]=...
linprogdis(ifint,f,A,b,Aeq,beq,lb,ub,x0,options)
%Title:
% 分支定届法求解混合整数线性规划模型
%
%初步完成：2002年12月
%最新修订: 2004-03-06
%最新注释：2004-11-20

%数据处理
[t1,t2] = size(b);
if t2~=1,
b=b';%将b转置为列向量
end

%调用线性规划求解
[x,fval,exitflag,output,lambda] = linprog(f,A,b,Aeq,beq,lb,ub,x0,options);
if exitflag<=0,%如果线性规划失败，则本求解也失败
return
end

%得到有整数约束的决策变量的序号
v1=find(ifint==1);%整数变量的index

tmp=x(v1);%【整数约束之决策变量】的当前值
if isempty(tmp),
%无整数约束，则是一般的线性规划，直接返回即可
return
end

v2=find(checkint(tmp)==0);%寻找不是整数的index
if isempty(v2),
%如果整数约束决策变量确实均为整数，则调用结束
return
end

%第k个决策变量还不是整数解
%注意先处理第1个不满足整数约束的决策变量
k=v1(v2(1));

%分支1：左分支
tmp1=zeros(1,length(f));%线性约束之系数向量
tmp1(k)=1;
low=floor(x(k));
%thisA 分支后实际调用线性规划的不等式约束的系数矩阵A
%thisb 分支后实际调用线性规划的不等式约束向量b
if ifrowinmat([tmp1,low],[A,b])==1
%如果分支的约束已经存在旧的A，b中，则不改变约束
thisA= A;
thisb= b;
else
thisA=[A;tmp1];
thisb=b;
thisb(end+1)=low;
end
%disp('fenzhi1'),thisA,thisb

%递归调用
[x1,fval1,exitflag1,output1,lambda1] = ...
linprogdis(ifint,f,thisA,thisb,Aeq,beq,lb,ub,x0,options);

%分支2：右分支
tmp2=zeros(1,length(f));%tmp1;
tmp2(k)=-1;
high= - ceil(x(k));
if ifrowinmat([tmp2,high],[A,b])==1
thisA= A;
thisb= b;
else
thisA=[A;tmp2];
thisb=b;
thisb(end+1)=high;
end
%disp('fenzhi2'),thisA,thisb

[x2,fval2,exitflag2,output2,lambda2] = ...
linprogdis(ifint,f,thisA,thisb,Aeq,beq,lb,ub,x0,options);

if isempty(v2) & ((exitflag1>0 & exitflag2<=0 & fval<=fval1 ) ...
| (exitflag2>0 & exitflag1<=0 & fval<=fval2 )...
| (exitflag1>0 & exitflag2>0 & fval<=fval1 & fval<=fval2 )),
disp('error call')
return %isempty(v2) 表示都是整数 2002-12-7非常重要
end
%下面分别根据返回标志exitflag确定最终的最优解
%case 1
if exitflag1>0 & exitflag2<=0 %【左分支有】最优解，右分支无最优解
x = x1;
fval = fval1;
exitflag = exitflag1;
output = output1;
lambda = lambda1;
%case 2
elseif exitflag2>0 & exitflag1<=0 %【右分支有】最优解，左分支无最优解
x = x2;
fval = fval2;
exitflag = exitflag2;
output = output2;
lambda = lambda2;
%case 3
elseif exitflag1>0 & exitflag2>0 %【左、右分支均有】最优解，则比较选优
if fval1 fval = fval2;
exitflag = exitflag2;
output = output2;
lambda = lambda2;
end
%fval1如果x(i)是整数，则返回r(i)=1；否则返回r(i)=0
%输入参数：x 向量
%输出参数：r 向量
for i=1:length(x),
if min(abs(x(i)-floor(x(i))),abs(x(i)-ceil(x(i))))<1e-03
%这里用于判定是否为0的参数可以调整，如改为1e-6
r(i)=1;
else
r(i)=0;
end
end

function r=ifrowinmat(arow,amat)
%输入参数：
% arow 向量，
% amat 矩阵
%
%设计：如果 arow与矩阵amat中的某一行相等，则返回1，如果找不到现等的一行，则返回0
r = 0;
rows = size(amat,1);
for i=1:rows,
temp = (amat(i,:)==arow);%利用 Matlab的==操作，如果相等，则为1向量；
if length(find(temp==0))==0,%没有为0的，即temp元素全部是1
r=1;
return
end%
end %for