### 分枝（支）定界法

%   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