分枝（支）定界法

Halcom

分枝（支）定界法：百度网盘 链接：https://pan.baidu.com/s/1gtX4sNDCzvTMYZnSIVHGWg 提取码：xq06
主程序如下：

1. %   This program solves mixed integer problems with a branch and bound method
   
2. clc,clear,close all
   
3. warning off
   
4. %Small test problem, optimal solution should be -21
   
5. c = -[8 11 6 4];
   
6. A = [5 7 4 3];
   
7. b = 14;
   
8. Aeq = [];
   
9. beq = [];
   
10. lb = [0 0 0 0]';
    
11. ub = [1 1 1 1]';
    
12. yidx = true(4,1);
    
13. [x_best,  f_best] = miprog(c,A,b,Aeq,beq,lb,ub,yidx);
    
14. disp(['X for Objective function value: ',num2str(x_best')]);
    
15. disp(['Objective function value: ',num2str(f_best)]);

复制代码

lp-relaxation problem程序如下：

1. function [x,fval,flag] = lpr(c,A,b,Aeq,beq,e,s,d,yidx,sol,opts)
   
2. %LPR
   
3. % Function to solve the lp-relaxation problem. Used by branch and bound algorithm.
   
4. %Calculate the extra constraints imposed by s
   
5. %if s is all NaN then no new constraints
   
6. sidx = ~isnan(s);
   
7. ydiag = double(yidx);
   
8. if any(sidx)
   
9.     ydiag(yidx) = d;
   
10.     l = length(yidx);  
    
11.     Y = spdiags(ydiag,0,l,l);
    
12.     %Remove zero entries
    
13.     Y = Y(any(Y,2),:);
    
14. end
    
15. 
    
16. %linprog
    
17. if any(sidx)
    
18.     A = [A; Y];
    
19.     b = [b; s(sidx).*d(sidx)];
    
20. end
    
21. [x,fval,flag] = linprog(c,A,b,Aeq,beq,[],[],[],opts);

复制代码

分支定界法程序如下：

1. function [x_best,  f_best] = miprog(c,A,b,Aeq,beq,lb,ub,yidx,o)
   
2. %MIPROG
   
3. % Branch and bound algorithm for linear mixed integer programs
   
4. %
   
5. %   min c'*x s.t. A*x <= b Aeq*x == beq lb <= x <= ub x(yidx) in [0,1,2,...]
   
6. %
   
7. %   where yidx is a logical index vector such that
   
8. %   x(yidx) are the integer variables
   
9. 
   
10. o.display = 'off';
    
11. o.iterplot = false;
    
12. o.solver = 'linprog';
    
13. %Branch and bound specific options
    
14. o.Delta = 1e-8; %Stopping tolerance of the gap (f_integer-f_lp)/(f_integer+f_lp)
    
15. o.maxNodes = 1e5; %Maximum number of nodes in the branch and bound tree to visit
    
16. o.intTol = 1e-6; %Integer tolerance
    
17. 
    
18. %% Init
    
19. %add upper and lower bounds to A
    
20. Ab = speye(length(yidx));
    
21. if ~isempty(lb)
    
22.     A = [A; -Ab];
    
23.     b = [b; -lb];
    
24. end
    
25. if ~isempty(ub)
    
26.     A = [A; Ab];
    
27.     b = [b; ub];
    
28. end
    
29. %delete all-zero rows
    
30. b = b(any(A,2),:);
    
31. A = A(any(A,2),:);
    
32. 
    
33. %linprog
    
34. e=[];
    
35. opts = optimset('Display','off');
    
36. sol=2;
    
37. 
    
38. %Assume no initial best integer solution
    
39. %Add your own heuristic here to find a good incumbent solution, store it in
    
40. %f_best,y_best,x_best
    
41. f_best = inf;
    
42. y_best = [];
    
43. x_best = [];
    
44. 
    
45. %Variable for holding the objective function variables of the lp-relaxation
    
46. %problems
    
47. f = inf(o.maxNodes,1);
    
48. f(1) = 0;
    
49. fs = inf;
    
50. numIntSol = double(~isempty(y_best));
    
51. 
    
52. %Set of problems
    
53. S = nan(sum(yidx),1);
    
54. D = zeros(sum(yidx),1);
    
55. 
    
56. %The priority in which the problems shall be solved
    
57. priority = [1];
    
58. %The indices of the problems that have been visited
    
59. visited = nan(o.maxNodes,1);
    
60. 
    
61. i=0;
    
62. %% Branch and bound loop
    
63. while i==0 || isinf(f_best) || (~isempty(priority) &&  ((f_best-min(fs(priority)))/abs(f_best+min(fs(priority))) > o.Delta) &&  i<o.maxNodes)
    
64.     %Is the parent node less expensive than the current best
    
65.     if i==0 || fs(priority(1))<f_best
    
66.         %Solve the LP-relaxation problem
    
67.         i=i+1;
    
68.         s = S(:,priority(1));
    
69.         d = D(:,priority(1));
    
70.         [x,this_f,flag] = lpr(c,A,b,Aeq,beq,e,s,d,yidx,sol,opts);
    
71.         
    
72.         %Visit this node
    
73.         visited(i) = priority(1);
    
74.         priority(1) = [];
    
75.         if flag==-2 || flag==-5
    
76.             %infeasible, dont branch
    
77.             if i==1
    
78.                 error('MIPROG: Infeasible initial LP problem. Try another LP solver.')
    
79.             end
    
80.             f(i) = inf;
    
81.         elseif flag==1
    
82.             %convergence
    
83.             f(i) = this_f;
    
84.             if this_f<f_best
    
85.                 y = x(yidx);
    
86.                 %fractional part of the integer variables -> diff
    
87.                 diff = abs(round(y)-y);
    
88.                 if all(diff<o.intTol)
    
89.                     %all fractions less than the integer tolerance
    
90.                     %we have integer solution
    
91.                     numIntSol = numIntSol+1;
    
92.                     f_best = this_f;
    
93.                     y_best = round(x(yidx));
    
94.                     x_best = x;
    
95.                 else
    
96.                     %branch on the most fractional variable
    
97.                     [maxdiff,branch_idx] = max(diff,[],1);
    
98.                     
    
99.                     %Branch into two subproblems
    
100.                     s1 = s;
     
101.                     s2 = s;
     
102.                     d1 = d;
     
103.                     d2 = d;
     
104.                     s1(branch_idx)=ceil(y(branch_idx));
     
105.                     d1(branch_idx)=-1; %direction of bound is GE
     
106.                     s2(branch_idx)=floor(y(branch_idx));
     
107.                     d2(branch_idx)=1; %direction of bound is LE
     
108.                     nsold = size(S,2);
     
109.                     
     
110.                     % add subproblems to the problem tree
     
111.                     S = [S s1 s2];
     
112.                     D = [D d1 d2];
     
113.                     fs = [fs f(i) f(i)];
     
114.                     
     
115.                     nsnew = nsold+2;
     
116.                     
     
117.                     %branch on the best lp solution
     
118.                     priority = [nsold+1:nsnew priority];
     
119.                     [dum,pidx] = sort(fs(priority));
     
120.                     priority=priority(pidx);
     
121.                     
     
122.                 end
     
123.             end
     
124.             
     
125.         else
     
126.             error('MIPROG: Problem neither infeasible nor solved, try another solver or reformulate problem!')
     
127.         end
     
128. 
     
129.     else %parent node is more expensive than current f-best -> don't evaluate this node
     
130.         priority(1) = [];
     
131.     end
     
132. end

复制代码

求解结果如下：

1. X for Objective function value: -8.1911e-11           1           1           1
   
2. Objective function value: -21

复制代码