分支定界算法的MATLAB程序

Linprogdis子程序: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);ift2~=1,b=b';%将b转置为列向量end%调用线性规划求解[x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb,ub,x0,options);ifexitflag=0,%如果线性规划失败，则本求解也失败returnend%得到有整数约束的决策变量的序号v1=find(ifint==1);%整数变量的indextmp=x(v1);%【整数约束之决策变量】的当前值ifisempty(tmp),%无整数约束，则是一般的线性规划，直接返回即可returnendv2=find(checkint(tmp)==0);%寻找不是整数的indexifisempty(v2),%如果整数约束决策变量确实均为整数，则调用结束returnend%第k个决策变量还不是整数解%注意先处理第1个不满足整数约束的决策变量k=v1(v2(1));%分支1：左分支tmp1=zeros(1,length(f));%线性约束之系数向量tmp1(k)=1;low=floor(x(k));%thisA分支后实际调用线性规划的不等式约束的系数矩阵A%thisb分支后实际调用线性规划的不等式约束向量bififrowinmat([tmp1,low],[A,b])==1%如果分支的约束已经存在旧的A，b中，则不改变约束thisA=A;thisb=b;elsethisA=[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));ififrowinmat([tmp2,high],[A,b])==1thisA=A;thisb=b;elsethisA=[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);ifisempty(v2)&((exitflag10&exitflag2=0&fval=fval1)...|(exitflag20&exitflag1=0&fval=fval2)...|(exitflag10&exitflag20&fval=fval1&fval=fval2)),disp('errorcall')return%isempty(v2)表示都是整数end%下面分别根据返回标志exitflag确定最终的最优解%case1ifexitflag10&exitflag2=0%【左分支有】最优解，右分支无最优解x=x1;fval=fval1;exitflag=exitflag1;output=output1;lambda=lambda1;%case2elseifexitflag20&exitflag1=0%【右分支有】最优解，左分支无最优解x=x2;fval=fval2;exitflag=exitflag2;output=output2;lambda=lambda2;%case3elseifexitflag10&exitflag20%【左、右分支均有】最优解，则比较选优iffval1fval2,%【左】分支最优（min）x=x1;fval=fval1;exitflag=exitflag1;output=output1;lambda=lambda1;elsex=x2;,%【右】分支最优（min）fval=fval2;exitflag=exitflag2;output=output2;lambda=lambda2;end%fval1fval2EndCheckint子程序:functionr=checkint(x)%算法：如果x(i)是整数，则返回r(i)=1；否则返回r(i)=0%输入参数：x向量%输出参数：r向量fori=1:length(x),ifmin(abs(x(i)-floor(x(i))),abs(x(i)-ceil(x(i))))1e-03%这里用于判定是否为0的参数可以调整，如改为1e-6r(i)=1;elser(i)=0;endendifrowinmat子程序:functionr=ifrowinmat(arow,amat)%输入参数：%arow向量，%amat矩阵%%设计：如果arow与矩阵amat中的某一行相等，则返回1，如果找不到现等的一行，则返回0r=0;rows=size(amat,1);fori=1:rows,temp=(amat(i,:)==arow);%利用Matlab的==操作，如果相等，则为1向量；iflength(find(temp==0))==0,%没有为0的，即temp元素全部是1r=1;returnend%end%for主程序(示例):ifint=[01];f=[109];a=[10;01;-5-3];b=[810-45];[x,fval,exitflag]=linprogdis(ifint,f,a,b,[],[],[],[],[],[])