Halcom 发表于 2019-4-9 00:08:38

分枝(支)定界法

分枝(支)定界法:百度网盘 链接:https://pan.baidu.com/s/1gtX4sNDCzvTMYZnSIVHGWg 提取码:xq06
主程序如下:
%   This program solves mixed integer problems with a branch and bound method
clc,clear,close all
warning off
%Small test problem, optimal solution should be -21
c = -;
A = ;
b = 14;
Aeq = [];
beq = [];
lb = ';
ub = ';
yidx = true(4,1);
= miprog(c,A,b,Aeq,beq,lb,ub,yidx);
disp(['X for Objective function value: ',num2str(x_best')]);
disp(['Objective function value: ',num2str(f_best)]);lp-relaxation problem程序如下:
function = lpr(c,A,b,Aeq,beq,e,s,d,yidx,sol,opts)
%LPR
% Function to solve the lp-relaxation problem. Used by branch and bound algorithm.
%Calculate the extra constraints imposed by s
%if s is all NaN then no new constraints
sidx = ~isnan(s);
ydiag = double(yidx);
if any(sidx)
    ydiag(yidx) = d;
    l = length(yidx);
    Y = spdiags(ydiag,0,l,l);
    %Remove zero entries
    Y = Y(any(Y,2),:);
end

%linprog
if any(sidx)
    A = ;
    b = ;
end
= linprog(c,A,b,Aeq,beq,[],[],[],opts);分支定界法程序如下:function = miprog(c,A,b,Aeq,beq,lb,ub,yidx,o)
%MIPROG
% Branch and bound algorithm for linear mixed integer programs
%
%   min c'*x s.t. A*x <= b Aeq*x == beq lb <= x <= ub x(yidx) in
%
%   where yidx is a logical index vector such that
%   x(yidx) are the integer variables

o.display = 'off';
o.iterplot = false;
o.solver = 'linprog';
%Branch and bound specific options
o.Delta = 1e-8; %Stopping tolerance of the gap (f_integer-f_lp)/(f_integer+f_lp)
o.maxNodes = 1e5; %Maximum number of nodes in the branch and bound tree to visit
o.intTol = 1e-6; %Integer tolerance

%% Init
%add upper and lower bounds to A
Ab = speye(length(yidx));
if ~isempty(lb)
    A = ;
    b = ;
end
if ~isempty(ub)
    A = ;
    b = ;
end
%delete all-zero rows
b = b(any(A,2),:);
A = A(any(A,2),:);

%linprog
e=[];
opts = optimset('Display','off');
sol=2;

%Assume no initial best integer solution
%Add your own heuristic here to find a good incumbent solution, store it in
%f_best,y_best,x_best
f_best = inf;
y_best = [];
x_best = [];

%Variable for holding the objective function variables of the lp-relaxation
%problems
f = inf(o.maxNodes,1);
f(1) = 0;
fs = inf;
numIntSol = double(~isempty(y_best));

%Set of problems
S = nan(sum(yidx),1);
D = zeros(sum(yidx),1);

%The priority in which the problems shall be solved
priority = ;
%The indices of the problems that have been visited
visited = nan(o.maxNodes,1);

i=0;
%% Branch and bound loop
while i==0 || isinf(f_best) || (~isempty(priority) &&((f_best-min(fs(priority)))/abs(f_best+min(fs(priority))) > o.Delta) &&i<o.maxNodes)
    %Is the parent node less expensive than the current best
    if i==0 || fs(priority(1))<f_best
      %Solve the LP-relaxation problem
      i=i+1;
      s = S(:,priority(1));
      d = D(:,priority(1));
       = lpr(c,A,b,Aeq,beq,e,s,d,yidx,sol,opts);
      
      %Visit this node
      visited(i) = priority(1);
      priority(1) = [];
      if flag==-2 || flag==-5
            %infeasible, dont branch
            if i==1
                error('MIPROG: Infeasible initial LP problem. Try another LP solver.')
            end
            f(i) = inf;
      elseif flag==1
            %convergence
            f(i) = this_f;
            if this_f<f_best
                y = x(yidx);
                %fractional part of the integer variables -> diff
                diff = abs(round(y)-y);
                if all(diff<o.intTol)
                  %all fractions less than the integer tolerance
                  %we have integer solution
                  numIntSol = numIntSol+1;
                  f_best = this_f;
                  y_best = round(x(yidx));
                  x_best = x;
                else
                  %branch on the most fractional variable
                   = max(diff,[],1);
                  
                  %Branch into two subproblems
                  s1 = s;
                  s2 = s;
                  d1 = d;
                  d2 = d;
                  s1(branch_idx)=ceil(y(branch_idx));
                  d1(branch_idx)=-1; %direction of bound is GE
                  s2(branch_idx)=floor(y(branch_idx));
                  d2(branch_idx)=1; %direction of bound is LE
                  nsold = size(S,2);
                  
                  % add subproblems to the problem tree
                  S = ;
                  D = ;
                  fs = ;
                  
                  nsnew = nsold+2;
                  
                  %branch on the best lp solution
                  priority = ;
                   = sort(fs(priority));
                  priority=priority(pidx);
                  
                end
            end
            
      else
            error('MIPROG: Problem neither infeasible nor solved, try another solver or reformulate problem!')
      end

    else %parent node is more expensive than current f-best -> don't evaluate this node
      priority(1) = [];
    end
end
求解结果如下:
X for Objective function value: -8.1911e-11         1         1         1
Objective function value: -21




















页: [1]
查看完整版本: 分枝(支)定界法