请选择 进入手机版 | 继续访问电脑版

Hello Mat

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 216|回复: 0

分枝(支)定界法

[复制链接]

687

主题

808

帖子

2万

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
29423
发表于 2019-4-9 00:08:38 | 显示全部楼层 |阅读模式
分枝(支)定界法:百度网盘 链接: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. %linprog
  16. if any(sidx)
  17.     A = [A; Y];
  18.     b = [b; s(sidx).*d(sidx)];
  19. end
  20. [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. o.display = 'off';
  10. o.iterplot = false;
  11. o.solver = 'linprog';
  12. %Branch and bound specific options
  13. o.Delta = 1e-8; %Stopping tolerance of the gap (f_integer-f_lp)/(f_integer+f_lp)
  14. o.maxNodes = 1e5; %Maximum number of nodes in the branch and bound tree to visit
  15. o.intTol = 1e-6; %Integer tolerance

  16. %% Init
  17. %add upper and lower bounds to A
  18. Ab = speye(length(yidx));
  19. if ~isempty(lb)
  20.     A = [A; -Ab];
  21.     b = [b; -lb];
  22. end
  23. if ~isempty(ub)
  24.     A = [A; Ab];
  25.     b = [b; ub];
  26. end
  27. %delete all-zero rows
  28. b = b(any(A,2),:);
  29. A = A(any(A,2),:);

  30. %linprog
  31. e=[];
  32. opts = optimset('Display','off');
  33. sol=2;

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

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

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

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

  53. i=0;
  54. %% Branch and bound loop
  55. while i==0 || isinf(f_best) || (~isempty(priority) &&  ((f_best-min(fs(priority)))/abs(f_best+min(fs(priority))) > o.Delta) &&  i<o.maxNodes)
  56.     %Is the parent node less expensive than the current best
  57.     if i==0 || fs(priority(1))<f_best
  58.         %Solve the LP-relaxation problem
  59.         i=i+1;
  60.         s = S(:,priority(1));
  61.         d = D(:,priority(1));
  62.         [x,this_f,flag] = lpr(c,A,b,Aeq,beq,e,s,d,yidx,sol,opts);
  63.         
  64.         %Visit this node
  65.         visited(i) = priority(1);
  66.         priority(1) = [];
  67.         if flag==-2 || flag==-5
  68.             %infeasible, dont branch
  69.             if i==1
  70.                 error('MIPROG: Infeasible initial LP problem. Try another LP solver.')
  71.             end
  72.             f(i) = inf;
  73.         elseif flag==1
  74.             %convergence
  75.             f(i) = this_f;
  76.             if this_f<f_best
  77.                 y = x(yidx);
  78.                 %fractional part of the integer variables -> diff
  79.                 diff = abs(round(y)-y);
  80.                 if all(diff<o.intTol)
  81.                     %all fractions less than the integer tolerance
  82.                     %we have integer solution
  83.                     numIntSol = numIntSol+1;
  84.                     f_best = this_f;
  85.                     y_best = round(x(yidx));
  86.                     x_best = x;
  87.                 else
  88.                     %branch on the most fractional variable
  89.                     [maxdiff,branch_idx] = max(diff,[],1);
  90.                     
  91.                     %Branch into two subproblems
  92.                     s1 = s;
  93.                     s2 = s;
  94.                     d1 = d;
  95.                     d2 = d;
  96.                     s1(branch_idx)=ceil(y(branch_idx));
  97.                     d1(branch_idx)=-1; %direction of bound is GE
  98.                     s2(branch_idx)=floor(y(branch_idx));
  99.                     d2(branch_idx)=1; %direction of bound is LE
  100.                     nsold = size(S,2);
  101.                     
  102.                     % add subproblems to the problem tree
  103.                     S = [S s1 s2];
  104.                     D = [D d1 d2];
  105.                     fs = [fs f(i) f(i)];
  106.                     
  107.                     nsnew = nsold+2;
  108.                     
  109.                     %branch on the best lp solution
  110.                     priority = [nsold+1:nsnew priority];
  111.                     [dum,pidx] = sort(fs(priority));
  112.                     priority=priority(pidx);
  113.                     
  114.                 end
  115.             end
  116.             
  117.         else
  118.             error('MIPROG: Problem neither infeasible nor solved, try another solver or reformulate problem!')
  119.         end

  120.     else %parent node is more expensive than current f-best -> don't evaluate this node
  121.         priority(1) = [];
  122.     end
  123. end
复制代码

求解结果如下:
  1. X for Objective function value: -8.1911e-11           1           1           1
  2. Objective function value: -21
复制代码





















程序:算法QQ  3283892722
群智能算法视频:X元,链接http://halcom.cn/forum.php?mod=forumdisplay&fid=73
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则


Python|Opencv|MATLAB|Halcom.cn  

GMT+8, 2019-4-22 16:20 , Processed in 0.124445 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表