|
一、整数规划基本原理
数学规划中的变量（部分或全部）限制为整数时，称为整数规划。若在线性规划模型中，变量限制为整数，则称为整数线性规划。目前所流行的求解整数规划的方法，往往只适用于整数线性规划。目前还没有一种方法能有效地求解一切整数规划。
1.1 整数规划的分类
(1)变量全限制为整数时，称为纯（完全）整数规划。
(2)变量部分限制为整数的，称为混合整数规划。
(3)变量取值要么为0要么为1，称为0-1规划。
1.2 整数规划的特点
1. 原线性规划有最优解，当自变量限制为整数后，其整数规划解出现下述情况：
(1)原线性规划最优解全是整数，则整数规划最优解与线性规划最优解一致。
(2)整数规划可能不存在可行解。
(3)有可行解(当然就存在最优解)，但最优解值变差。
2. 整数规划最优解不能按照实数最优解简单取整而获得
1.3 案例
合理下料问题：
设用某型号的圆钢下零件A1,A2...,Am的毛坯。在一根圆钢上下料的方式有B1,B2,... Bn种，每种下料方式可以得到各种零件的毛坯数以及每种零件的需要量，如表所示。问怎样安排下料方式，使得即满足需要，所用的原材料又最少?
如表所示：
模型：
其中，xj表示用Bj种方式下料根数，aij为每种下料方式可以得到各种零件的毛坯数，bi每种零件的需要量。
1.4 整数规划的数学模型一般形式
另外补充： 整数规划与松弛的线性规划之间的关系。
不难看出两者之间的关系。
二、整数线性规划的求解方法
从数学模型上看整数规划似乎是线性规划的一种特殊形式，求解只需在线性规划的基础上，通过舍入取整，寻求满足整数要求的解即可。但实际上两者却有很大的不同，通过舍入得到的解（整数）也不一定就是最优解，有时甚至不能保证所得到的解是整数可行解。
2.1 分枝定界法
分枝定界法（branch and bound）是一种求解整数规划问题的最常用算法。分支定界法是一种搜索与迭代的方法，选择不同的分支变量和子问题进行分支。
2.1.1 分枝定界法的求解过程
1. 先求出相应松弛问题的最优解。
2. 若松弛问题无可行解，则整数线性规划问题无可行解。
3. 若求得的松弛问题最优解符合整数要求，则是整数线性规划问题的最优解。若不满足整数条件，则任选一个不满足整数条件的变量 来构造新的约束，分别添加到松弛问题中形成两个子问题：
和
依次在缩小的可行域中求解新构造的线性规划的最优解，并重复上述过程，直到子问题无解或有整数最优解（被查清）。
Q：为什么在步骤3中不满足整数条件时要添加上述两个约束条件？
A：因为变量是一个不满足整数条件的变量，它注定无法构成整数规划的最优解，因此我们要向变量的两边去寻找合适的整数最优解。
注意：构造新的约束条件是分别到整数规划问题对应的松弛问题中，独立构成两个不同的子问题再求解。
2.1.2 案例
分枝定界法总体上有点像二叉树，能看懂下面的案例基本就能理解分枝定界法。
2.1.3 代码实现
数学模型如下：
1. 先创建intprog.m文件：
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
2. 再创建branchbound.m文件
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
3. 测试
function[x. fval.stetus] = intprog(f,A,B,I,Aeq,Beq,lb,ub,e)
f = [-20, -10]
A = [5,4;2,5]
B = [24,13]
lb = [0,0]
[x, fval, status] = intprog(f, A, B, [1,2],[],[],lb)
2.2 割平面法
2.2.1 割平面法的基本思想
1. 如果松弛问题(P0)无解，则(P)无解;
2. 如果(P0)的最优解为整数向量，则也是(P)的最优解;
3. 如果(P0)的解含有非整数分量，则对(P0)增加割平面条件：即对(P0)增加一个线性约束，将(P0)的可行区域割掉一块，使得非整数解恰好在割掉的一块中，但又没有割掉原问题(P)的可行解，得到问题(P1)，重复上述的过程。
2.2.2 割平面法的基本步骤
注：这里的第二、三步是推理过程，简略来看，只需要看第一步、第二步中1和第四步即可。
第一步：求解整数规划对应松弛问题是否有整数最优解。若有整数最优解，则它也是整数规划的最优解，否则，转到第二步。
第二步：1. 引入松弛变量转化为等式约束：
其中，xi为决策变量，xk为松弛变量，aik为松弛变量的系数，bi与原约束条件一致。
然后，将松弛变量的系数分别划分为整数部分和小数部分，如下：
其中，[aik]表示向下取整，Laik也表示aik向下取整，fik表示小数部分（aik减去它整数部分得到的小数，小于0）。
同样，将求和结果划分为整数部分和小数部分，如下：
其中，[bi]表示向下取整，Lbi也表示aik向下取整，fi表示小数部分（同上，小于0）。
第三步：将划分后的aik和bi代入所构建的方程中，可得：
我们将所有的整数部分放入公式左侧，所有的小数部分放入公式右侧：
因为 <0，所以 <0，则可推出：
这样就对决策变量进行了一次割平面（增加约束条件）。
第四步：将约束条件添加到原数学模型中，重复第一步。
2.2.3 案例和代码实现
已知AM工厂是一个拥有四个车间的玩具生产厂商，该厂商今年新设计出A、B、C、D、E、F六种玩具模型，根据以前的生产情况及市场调查预测，得知生产每单位产品所需的工时、每个车间在一季度的工时上限以及产品的预测价格，如下表所示。问:每种设计产品在这个季度各应生产多少，才能使AM工厂这个季度的生产总值达到最大?
1. 建立模型
其中，x1，x2，x3，x4，x5，x6 是A、B、C、D、E、F六种玩具模型的生产数量。注意，是车间之间是合作生产，而不是独立生产。
2. 引入松弛变量转化为等式约束
3. 代码实现
DividePlane.m
function [intx,intf] = DividePlane(A,c,b,baseVector)
%功能：用割平面法求解整数规划
%调用格式：[intx,intf]=DividePlane(A,c,b,baseVector)
%其中, A：约束矩阵；
% c：目标函数系数向量；
% b：约束右端向量；
% baseVector：初始基向量；
% intx：目标函数取最值时的自变量值；
% intf：目标函数的最值；
sz = size(A);
nVia = sz(2);%获取有多少决策变量
n = sz(1);%获取有多少约束条件
xx = 1:nVia;
if length(baseVector) ~= n
disp('基变量的个数要与约束矩阵的行数相等！');
mx = NaN;
mf = NaN;
return;
end
M = 0;
sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
xb = b;
%首先用单纯形法求出最优解
while 1
[maxs,ind] = max(sigma);
%--------------------用单纯形法求最优解--------------------------------------
if maxs <= 0 %当检验数均小于0时，求得最优解。
vr = find(c~=0 ,1,'last');
for l=1:vr
ele = find(baseVector == l,1);
if(isempty(ele))
mx(l) = 0;
else
mx(l)=xb(ele);
end
end
if max(abs(round(mx) - mx))<1.0e-7 %判断最优解是否为整数解，如果是整数解。
intx = mx;
intf = mx*c;
return;
else %如果最优解不是整数解时，构建切割方程
sz = size(A);
sr = sz(1);
sc = sz(2);
[max_x, index_x] = max(abs(round(mx) - mx));
[isB, num] = find(index_x == baseVector);
fi = xb(num) - floor(xb(num));
for i=1:(index_x-1)
Atmp(1,i) = A(num,i) - floor(A(num,i));
end
for i=(index_x+1):sc
Atmp(1,i) = A(num,i) - floor(A(num,i));
end
Atmp(1,index_x) = 0; %构建对偶单纯形法的初始表格
A = [A zeros(sr,1);-Atmp(1,:) 1];
xb = [xb;-fi];
baseVector = [baseVector sc+1];
sigma = [sigma 0];
%-------------------对偶单纯形法的迭代过程----------------------
while 1
%----------------------------------------------------------
if xb >= 0 %判断如果右端向量均大于0，求得最优解
if max(abs(round(xb) - xb))<1.0e-7 %如果用对偶单纯形法求得了整数解，则返回最优整数解
vr = find(c~=0 ,1,'last');
for l=1:vr
ele = find(baseVector == l,1);
if(isempty(ele))
mx_1(l) = 0;
else
mx_1(l)=xb(ele);
end
end
intx = mx_1;
intf = mx_1*c;
return;
else %如果对偶单纯形法求得的最优解不是整数解，继续添加切割方程
sz = size(A);
sr = sz(1);
sc = sz(2);
[max_x, index_x] = max(abs(round(mx_1) - mx_1));
[isB, num] = find(index_x == baseVector);
fi = xb(num) - floor(xb(num));
for i=1:(index_x-1)
Atmp(1,i) = A(num,i) - floor(A(num,i));
end
for i=(index_x+1):sc
Atmp(1,i) = A(num,i) - floor(A(num,i));
end
Atmp(1,index_x) = 0; %下一次对偶单纯形迭代的初始表格
A = [A zeros(sr,1);-Atmp(1,:) 1];
xb = [xb;-fi];
baseVector = [baseVector sc+1];
sigma = [sigma 0];
continue;
end
else %如果右端向量不全大于0，则进行对偶单纯形法的换基变量过程
minb_1 = inf;
chagB_1 = inf;
sA = size(A);
[br,idb] = min(xb);
for j=1:sA(2)
if A(idb,j)<0
bm = sigma(j)/A(idb,j);
if bm<minb_1
minb_1 = bm;
chagB_1 = j;
end
end
end
sigma = sigma -A(idb,:)*minb_1;
xb(idb) = xb(idb)/A(idb,chagB_1);
A(idb,:) = A(idb,:)/A(idb,chagB_1);
for i =1:sA(1)
if i ~= idb
xb(i) = xb(i)-A(i,chagB_1)*xb(idb);
A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);
end
end
baseVector(idb) = chagB_1;
end
%------------------------------------------------------------
end
%--------------------对偶单纯形法的迭代过程---------------------
end
else %如果检验数有不小于0的，则进行单纯形算法的迭代过程
minb = inf;
chagB = inf;
for j=1:n
if A(j,ind)>0
bz = xb(j)/A(j,ind);
if bz<minb
minb = bz;
chagB = j;
end
end
end
sigma = sigma -A(chagB,:)*maxs/A(chagB,ind);
xb(chagB) = xb(chagB)/A(chagB,ind);
A(chagB,:) = A(chagB,:)/A(chagB,ind);
for i =1:n
if i ~= chagB
xb(i) = xb(i)-A(i,ind)*xb(chagB);
A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);
end
end
baseVector(chagB) = ind;
end
M = M + 1;
if (M == 1000000)
disp('找不到最优解！');
mx = NaN;
minf = NaN;
return;
end
end
test.m
c = [-20 ; -14 ; -16 ; -36 ; -32 ; -30]; % 目标函数系数向量
A = [0.01 0.01 0.01 0.03 0.03 0.03 1 0 0 0;
0.02 0 0 0.05 0 0 0 1 0 0;
0 0.02 0 0 0.05 0 0 0 1 0;
0 0 0.03 0 0 0.08 0 0 0 1]; % 加上松弛变量后约束条件的系数矩阵
b = [850; 700; 100; 900;]; % 约束右端向量
[intx , intf] = DividePlane(A,c,b,[7 8 9 10]); % 初始基向量的下标 7 8 9 10
intx
intf = -intf % 取反求最大
|
