mip.m

运用分支定界法求解混合整数线性规划

%%%运用分枝定界法求解混合整数规划程序
%%%各输入参数的含义：
%  min f*x
% A*x<=b
% Aeq*x=beq
% iter_num_max:迭代次数的上限
% lb:各变量的下限
% ub:各变量的上限
% x0:初值
% ID:0、1变量，ID(k)=1,表示第k个变量为整数变量，ID(k)=0，表示第k个变量为非整数变量
% ftemp_up:初始上界,无法找到初始上界，默认为无穷大
% IP：初始可行解，无法找到初始可行解，默认为空
%%%各输出参数的含义：
% x:最优解
% y:最优解对应的目标函数值
% iter_num:整个分支定界过程的迭代次数
% Time:运行本程序所耗费的时间
function [x,y,iter_num]=MIP(f,A,b,Aeq,beq,iter_num_max,lb,ub,x0,ID,ftemp_up,IP,options)
if nargin<13
    options=optimset({});
    options.Display='off';
    options.LargeScale='on';
end
if nargin<12
    IP=[];
end
if nargin<11
    ftemp_up=inf;
end
if nargin<10
    ID=ones(size(f));
end
if nargin<9
    x0=[];
end
if nargin<8
    % lisempty(ub),
    ub=inf*ones(size(f));
end
if nargin<7
    % lisempty(lb),
    lb=zeros(size(f));
end
if nargin<6
    iter_num_max=100;
end
if nargin<5
    beq=[];
end
if nargin<4
    Aeq=[];
end
iter_num=0;
[x,ftemp,how]=linprog(f,A,b,Aeq,beq,lb,ub,x0,options);
iter_num=iter_num+1;
if how<=0
    fprintf('该优化模型无可行解\n');
    y=ftemp;
    Time=0;
    return;
end
if max(abs(x'.*ID-round(x'.*ID)))<0.00005
    y=ftemp;
    return; 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%找出需要分支的变量
notintx=[];
notintx_flag=1;
bin_var=find(ID==1);
for k=1:length(bin_var)
    if abs(x(bin_var(k))-round(x(bin_var(k))))>=0.00005
        notintx(notintx_flag)=bin_var(k);
        notintx_flag=notintx_flag+1;
    end
end

for k=1:length(notintx)
    xx(k)=x(notintx(k));
    mean_temp=[lb(notintx(k)),ub(notintx(k))];
    mean_value(k)=mean(mean_temp);
end  

branch_var=find(abs(xx-mean_value)==min(abs(xx-mean_value)));
branch_var_index=notintx(branch_var(1));        %%待分支的变量序号
mean_temp=[];
mean_value=[];
xx=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%初始化需要存放非整数解信息的结构体和需要存放整数解信息的矩阵
RP=struct('x',[],'ftemp',[],'lb',[],'ub',[]);   %%存放非整数解信息
RP(1)=[];
RP_flag=1;
if isempty(IP)
%IP=[];                                          %%存放整数解信息
   IP_flag=1;
else
    IP_flag=2;
    IP(1,:)=IP;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%对需要分支的变量进行第一次分支，即改变该变量的下界再一次进行优化
x_temp=x;
vlb=lb;
vlb(branch_var_index)=fix(x_temp(branch_var_index))+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[x,ftemp,how]=linprog(f,A,b,Aeq,beq,vlb,ub,x0,options);
iter_num=iter_num+1;
if  how>0
    if ftemp<ftemp_up
       if max(abs(x'.*ID-round(x'.*ID)))<0.00005
           IP(IP_flag,1:length(f))=x';                   %%存放整数解信息
           IP(IP_flag,length(f)+1)=ftemp;
           IP_flag=IP_flag+1;
       else
           RP(RP_flag)=struct('x',x','ftemp',ftemp,'lb',vlb,'ub',ub);   %%存放非整数解信息
           RP_flag=RP_flag+1;
       end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%对需要分支的变量进行第二次分支，即改变该变量的上界再一次进行优化
vub=ub;
vub(branch_var_index)=fix(x_temp(branch_var_index));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[x,ftemp,how]=linprog(f,A,b,Aeq,beq,lb,vub,x0,options);
iter_num=iter_num+1;
if  how>0
    if ftemp<ftemp_up
       if max(abs(x'.*ID-round(x'.*ID)))<0.00005
           IP(IP_flag,1:length(f))=x';
           IP(IP_flag,length(f)+1)=ftemp;
           IP_flag=IP_flag+1;
       else
           RP(RP_flag)=struct('x',x','ftemp',ftemp,'lb',lb,'ub',vub);
           RP_flag=RP_flag+1;
       end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

while ~isempty(RP) && iter_num<iter_num_max
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%从结构体RP中找到目标函数值最小的那个记录作为分支的对象
    temp=[];
    for k=1:length(RP)
        temp(k)=RP(k).ftemp;
    end
    temp_index=find(temp==min(temp));
    temp_branch=RP(temp_index(1));       %%需要分支的记录
    RP(temp_index)=[];
    RP_flag=length(RP)+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%找出需要分支的变量
    notintx=[];
    notintx_flag=1;
    for k=1:length(bin_var)         
      if  abs(temp_branch.x(bin_var(k))-round(temp_branch.x(bin_var(k))))>=0.00005
          notintx(notintx_flag)=bin_var(k);
          notintx_flag=notintx_flag+1;
      end
    end
    for k=1:length(notintx)
        xx(k)=temp_branch.x(notintx(k));
        mean_temp=[temp_branch.lb(notintx(k)),temp_branch.ub(notintx(k))];
        mean_value(k)=mean(mean_temp);
    end  

       branch_var=find(abs(xx-mean_value)==min(abs(xx-mean_value)));
       branch_var_index=notintx(branch_var(1));        %%待分支的变量序号
       mean_temp=[];
       mean_value=[];
       xx=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%对需要分支的变量进行第一次分支，即改变该变量的下界再一次进行优化
     vlb=temp_branch.lb;
     vub=temp_branch.ub;
     vlb(branch_var_index)=fix(temp_branch.x(branch_var_index))+1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   [x,ftemp,how]=linprog(f,A,b,Aeq,beq,vlb,vub,x0,options);
   iter_num=iter_num+1;
   if   how>0
       if ftemp<ftemp_up
         if max(abs(x'.*ID-round(x'.*ID)))<0.00005
             IP(IP_flag,1:length(f))=x';                   %%存放整数解信息
             IP(IP_flag,length(f)+1)=ftemp;
             IP_flag=IP_flag+1;
         else
            RP(RP_flag)=struct('x',x','ftemp',ftemp,'lb',vlb,'ub',vub);   %%存放非整数解信息
            RP_flag=RP_flag+1;
         end
       end
   end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%对需要分支的变量进行第二次分支，即改变该变量的上界再一次进行优化
      vlb=temp_branch.lb;
      vub=temp_branch.ub;
      vub(branch_var_index)=fix(temp_branch.x(branch_var_index));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     [x,ftemp,how]=linprog(f,A,b,Aeq,beq,vlb,vub,x0,options);
     iter_num=iter_num+1;
    if  how>0
        if ftemp<ftemp_up
          if max(abs(x'.*ID-round(x'.*ID)))<0.00005
             IP(IP_flag,1:length(f))=x';
             IP(IP_flag,length(f)+1)=ftemp;
             IP_flag=IP_flag+1;
          else
            RP(RP_flag)=struct('x',x','ftemp',ftemp,'lb',vlb,'ub',vub);
            RP_flag=RP_flag+1;
          end
       end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    if(~isempty(IP))
         ftemp_up_flag=find(IP(:,length(f)+1)==min(IP(:,length(f)+1)));
         ftemp_up=IP(ftemp_up_flag,length(f)+1);
    end
    iter_num=iter_num+1;
end
if isempty(IP)==0
     opt_flag=find(IP(:,length(f)+1)==min(IP(:,length(f)+1)));
     x=IP(opt_flag,1:size(IP,2)-1);              %%返回最优解
     y=IP(opt_flag,size(IP,2));                  %%返回最优解对应的目标函数值
elseif iter_num<iter_num_max
    x=[];
    y=[];
    fprintf('无可行解\n');
else
    x=[];
    y=[];
    fprintf('达到迭代次数上限还未找到可行解\n');
end