(BranchBound)- 线性整数规划

moutqvw

1. function [x,y]=lpint(f,G,h,lb,ub,x,n,id)
   
2. %整数线性规划分枝定界法，可求解线性全整数或线性混合整数规划。
   
3. %         y = min f'x    subject to:   Gx <= h  x为整
   
4. %              x
   
5. %用法
   
6. %    [x,y]=lpint(f,G,h)
   
7. %    [x,y]=lpint(f,G,h,lb,ub)
   
8. %    [x,y]=lpint(f,G,h,lb,ub,x)
   
9. %    [x,y]=lpint(f,G,h,lb,ub,x,n)
   
10. %    [x,y]=lpint(f,G,h,lb,ub,x,n,id)
    
11. %参数说明
    
12. %    x: 最优解列向量
    
13. %    y: 目标函数最小值
    
14. %    f: 目标函数系数列向量
    
15. %    G: 约束条件系数矩阵
    
16. %    h: 约束条件右端列向量
    
17. %    lb: 解的的下界列向量(Default: -inf)
    
18. %    ub: 解的的上界列向量(Default: inf)
    
19. %    x: 迭代初值列向量
    
20. %    n: 等式约束数(Default: 0)
    
21. %    id: 整数变量指标列向量。1-整数，0-实数(Default: 1)
    
22. %例: min Z=x1+4x2
    
23. % s.t.  2x1+x2<=8
    
24. %      x1+2x2>=6
    
25. %      x1, x2>=0且为整数
    
26. %先将x1+2x2>=6化为 - x1 - 2x2<= -6
    
27. %[x,y]=lpint([1;4],[2 1;-1 -2],[8;-6],[0;0])
    
28. 
    
29. % Y. MA & L.J. HU 1999
    
30. 
    
31. global upper opt c N x0 A b ID;
    
32. if nargin<8, id=ones(size(f));end
    
33. if nargin<7|isempty(n), n=0;end
    
34. if nargin<6, x=[];end
    
35. if nargin<5|isempty(ub), ub=inf*ones(size(f));end
    
36. if nargin<4|isempty(lb), lb=zeros(size(f));end
    
37. 
    
38. upper=inf;
    
39. c=f;N=n;x0=x;A=G;b=h;ID=id;
    
40. temp=ILP(lb(:),ub(:));
    
41. x=opt;y=upper;
    
42. 
    
43. %以下子函数
    
44. function y=ILP(vlb,vub)
    
45. global upper opt c N x0 A b ID;
    
46.   warning off;
    
47.   [x,temp,how]=lp(c,A,b,vlb,vub,x0,N,-1);
    
48.   if strcmp(how,'ok')~=1
    
49.      return;
    
50.   end;
    
51.   if c'*x-upper>0.00005  %in order to avoid error
    
52.      return;
    
53.   end;
    
54.    
    
55.   if max(abs(x.*ID-round(x.*ID)))<0.00005
    
56.     if upper-c'*x>0.00005 %in order to avoid error
    
57.       opt=x';
    
58.       upper=c'*x;
    
59.       return;
    
60.    else
    
61.       opt=[opt;x'];
    
62.       return;
    
63.     end;
    
64.   end;
    
65.   notintx=find(abs(x-round(x))>=0.00005); %in order to avoid error
    
66.   intx=fix(x);
    
67.   tempvlb=vlb;
    
68.   tempvub=vub;
    
69.   if vub(notintx(1,1),1)>=intx(notintx(1,1),1)+1
    
70.      tempvlb(notintx(1,1),1)=intx(notintx(1,1),1)+1;
    
71.      temp=ILP(tempvlb,vub);
    
72.   end;
    
73.   if vlb(notintx(1,1),1)<=intx(notintx(1,1),1)
    
74.     tempvub(notintx(1,1),1)=intx(notintx(1,1),1);
    
75.      temp=ILP(vlb,tempvub);
    
76.   end;

复制代码