【数学建模入门】整数规划求最优解｜分枝定界法求最优解

整数规划

intprog()

function [x,fval,status] = intprog(f,A,B,I,Aeq,Beq,lb,ub,e)

其中f为目标函数向量
A和B为不等式约束Aeq和Beq为等式约束
I为整数约束
lb与ub分别为变量下界与上界
x为最优解，fval为最优值

---

例子：

maximiaze 20*x1+10*x2
S.T.
		5 x1 + 4 x2 <= 24
		2 x1 + 5 x2 <=13
			x1, x2 >=0
			x1,x2是整数

分枝定届法

intprog

function [x,fval,status] = intprog(f,A,B,I,Aeq,Beq,lb,ub,e)
%整数规划求解函数 intprog()
%     其中 f为目标函数向量
%     A和B为不等式约束 Aeq与Beq为等式约束
%     I为整数约束
%     lb与ub分别为变量下界与上界
%     x为最优解，fval为最优值
%例子:
%        maximize 20 x1 + 10 x2 
%        S.T.
% 	             5 x1 + 4 x2 <=24
%                2 x1 + 5 x2 <=13
%                   x1, x2 >=0 
%                   x1, x2是整数
% f=[-20, -10];
% A=[ 5  4; 2 5];
% B=[24; 13];
% lb=[0 0];
% ub=[inf inf];
% I=[1,2];
% e=0.000001;
% [x v s]= IP(f,A,B,I,[],[],lb,ub,,e)
% x = 4     1  v = -90.0000   s = 1

% 控制输入参数
if nargin < 9, e = 0.00001;
   if nargin < 8, ub = []; 
      if nargin < 7, lb = []; 
         if nargin < 6, Beq = []; 
            if nargin < 5, Aeq = [];
               if nargin < 4, I = [1:length(f)];
               end, end, end, end, end, end
  
%求解整数规划对应的线性规划,判断是否有解
options = optimset('display','off');
[x0,fval0,exitflag] = linprog(f,A,B,Aeq,Beq,lb,ub,[],options);
if exitflag < 0
    disp('没有合适整数解');
    x = x0;
    fval = fval0;
    status = exitflag;
    return;
else
    %采用分支定界法求解
    bound = inf;
    [x,fval,status] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
end

branchbound

function [newx,newfval,status,newbound] = branchbound(f,A,B,I,x,fval,bound,Aeq,Beq,lb,ub,e)

% 分支定界法求解整数规划
% f,A,B,Aeq,Beq,lb,ub与线性规划相同
% I为整数限制变量的向量
% x为初始解，fval为初始值

options = optimset('display','off');
[x0,fval0,status0]=linprog(f,A,B,Aeq,Beq,lb,ub,[],options);

%递归中的最终退出条件
%无解或者解比现有上界大则返回原解
if status0 <= 0 || fval0 >= bound
    newx = x;
    newfval = fval;
    newbound = bound;
    status = status0;
    return;
end

%是否为整数解,如果是整数解则返回
intindex = find(abs(x0(I) - round(x0(I))) > e);
if isempty(intindex) %判断是否为空值
    newx(I) = round(x0(I));
    newfval = fval0;
    newbound = fval0;
    status = 1;
    return;
end

%当有非整可行解时，则进行分支求解
%此时必定会有整数解或空解
%找到第一个不满足整数要求的变量
n = I(intindex(1));
addA = zeros(1,length(f));
addA(n) = 1;
%构造第一个分支 x<=floor(x(n))
A = [A;addA];
B = [B,floor(x(n))];%向下取整
[x1,fval1,status1,bound1] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
A(end,:) = [];
B(:,end) = [];
%解得第一个分支，若为更优解则替换，若不是则保持原状

status = status1;
if status1 > 0 && bound1 < bound
    newx = x1;
    newfval = fval1;
    bound = fval1;
    newbound = bound1;
else
    newx = x0;
    newfval = fval0;
    newbound = bound;
end

%构造第二分支
A = [A;-addA];
B = [B,-ceil(x(n))];%向上取整
[x2,fval2,status2,bound2] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);    
A(end,:) = [];
B(:,end) = [];

%解得第二分支，并与第一分支做比较，如果更优则替换
if status2 > 0 && bound2 < bound
    status = status2;
    newx = x2;
    newfval = fval2;
    newbound = bound2;
end