整数规划模型与求解
1.整数规划问题分类
（1）纯整数规划：决策变量要求取整，但松弛变量、剩余变量不需取整。
（2）全整数规划：决策变量、松弛变量、剩余变量都要取整。
（3）混合整数变量：部分决策变量取整。
（4）0-1整数规划：指派问题，决策变量只取0或1.
若线性规划无可行解，则相应的整数变量也无可行解。
2.整数规划模型
我们在这里只讨论线性整数规划，其标准模型为标准模型：
3.算法与matlab求解
(1)分枝定界法
a.算法思路：
不考虑整数限制，先求相应松弛问题的最优解；
若松弛问题无可行解，则整数规划问题无可行解；
若松弛问题最优解符合整数要求，则该解也是整数规划的最优解；
若不满足整数条件
则任选一个不满足整数条件的变量
x
i
0
x_i^0
xi0来构造一个新的约束，添加到松弛问题中形成两个子问题
x
i
≤
[
x
i
0
]
;
x
i
≥
[
x
i
0
]
+
1
x_i\le[x_i^0] \ ;\quad x_i\ge[x_i^0]+1
xi≤[xi0] ;xi≥[xi0]+1
依次缩小可行域，求解新构造的线性规划最优解，并重复上述过程，直到问题误解或有整数最优解（被查清）
若最后还未成功，则考虑其他方法，分枝定界法失效。
分支定界法可用于解纯整数和混合整数问题
b.matlab求解
主函数：
%例子
f = [-20 -10];
A = [5 4;2 5];
B = [24 13];
lb = [0 0];
[x,fval,status] = intprog(f,A,B,[1 2],[],[],lb)
二叉树函数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
整数规划求解函数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
[x,fval,status] = intprog(f,A,B,I,Aeq,Beq,lb,ub,e)
整数规划求解函数 intprog
其中 f为目标函数向量
A和B为不等式约束 Aeq与Beq为等式约束
I为整数约束
lb与ub分别为变量下界与上界
x为最优解，fval为最优值
（2）割平面法
a.算法思路
若松弛算法（
P
0
P_0
P0）无解，则（
P
P
P)无解；
若松弛问题最优解符合整数要求，则该解也是（
P
P
P）最优解；
若不能满足整数条件，则对（
P
0
P_0
P0）增加割平面条件：即对(
P
0
P_0
P0)增加一个线性约束，将(
P
0
P_0
P0)的可行区域割掉一块，但又没有割掉原问题的可行解，得到问题(
P
1
P_1
P1)，然后重复上述过程。
切割方程构建方法：
1.将方程化作：一个决策变量+多个松弛变量的形式：
x
i
+
∑
a
i
k
x
k
=
b
i
x_i+\sum a_{ik}x_{k}=b_i
xi+∑aikxk=bi
2.把每个松弛变量写成整数+小数的形式,把
b
i
b_i
bi也做同样处理：
a
i
k
=
1
+
(
a
i
k
−
1
)
=
[
a
i
k
]
+
f
i
k
a_{ik}=1+(a_{ik}-1)=[a_{ik}]+f_{ik}
aik=1+(aik−1)=[aik]+fik
b
i
=
1
+
(
b
i
k
−
1
)
=
[
b
i
k
]
+
f
i
b_{i}=1+(b_{ik}-1)=[b_{ik}]+f_i
bi=1+(bik−1)=[bik]+fi
3.把整数都放到左边，小数都放到右边，得到下式
x
i
+
∑
[
a
i
k
]
−
∑
[
b
i
k
]
=
f
i
−
∑
[
f
i
k
]
x_i+\sum [a_{ik}]-\sum [b_{ik}]=fi-\sum [f_{ik}]
xi+∑[aik]−∑[bik]=fi−∑[fik]
易知：
f
i
−
∑
[
f
i
k
]
≤
0
fi-\sum [f_{ik}]\le0
fi−∑[fik]≤0
这个就是我们用割平面构建的线性条件。
b.matlab求解
创建：[intx,intf] = DividePlane(A,c,b,baseVector)函数来解决
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
注意，该函数的A是加上松弛变量后的系数矩阵
应用举例
例：
设A、B、C、D、E、F六个车间这个季度各生产量为
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
x_1,x_2,x_3,x_4,x_5,x_6
x1,x2,x3,x4,x5,x6建立数学模型如下：
写成含有松弛条件的形式：
求解代码如下（要先创建DividePlane(函数）：
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];
c=[-20;-14;-16;-36;-32;-30];
b=[850;700;100;900];
baseVector=[7,8,9,10];
[intx,intf]=DividePlane(A,c,b,baseVector)
(3)0-1规划问题与匈牙利算法
a.常见0-1规划问题
1.投资问题：一般有
2.互斥约束问题及其推广：
3.固定费用问题（生产线安排）
4.指派问题
b.指派问题的标准形式
某单位有项任务，正好需n个人去完成，由于每项任务的性质和每个人的能力和专长的不同，假设分配每个人仅完成一项任务。设 c i j c_{ij} cij表示分配第 i i i个人去完成第 j j j项任务的费用（时间等），问应如何指派，完成任务的总费用最小？
建立数学模型
现将工作代价列成矩阵
C
C
C
所以可以抽象出以下问题：
在x取值限制为0，1的情况下
求解：
我们称矩阵
X
X
X：
为指派问题的解矩阵。
c.非标准形式的指派问题
1.最大化指派问题： 找到
C
=
(
c
i
j
)
n
×
n
C=(c_{ij})_{n \times n}
C=(cij)n×n中最大元素m，令矩阵
B
=
(
b
i
j
)
n
×
n
=
(
m
−
c
i
j
)
n
×
n
B=(b_{ij})_{n\times n}=(m-c_{ij})_{n\times n}
B=(bij)n×n=(m−cij)n×n ,这时候把
B
B
B当作系数矩阵来求解。*
2.人数和工作数不等：
如果人少工作多：添加虚拟的人，代价都为0
如果人多工作少：添加虚拟的工作，代价都为0
3.一个人做多个工作：
该人可化作多个相同的人，对虚拟人不限制0，1
4.某工作一定不能由某人做：
该人做该工作的相应代价取足够大的M,这样可以保证
x
i
j
≤
0
x_{ij}\le0
xij≤0。
d.匈牙利算法
第一步：变换指派问题的系数(也称效率)矩阵(cij) 为(bij)，使在 (bij) 的各行各列中都出现0元素，
即：(1) 从 (cij) 的每行元素都减去该行的最小元素；
(2) 再从所得新系数矩阵的每列元素中减去该列的最小元素
第二步：进行试指派，以寻求最优解。 在(bij)中找尽可能多的独立0元素，若能找出n个独立0元素，就以这n 个独立0元素对应解矩阵(xij)中的元素为1，其余为0，这就得到最优解。找 独立0元素，常用的步骤为： (1)从只有一个0元素的行(列)开始，给这个0元素加圈，记作◎ 。然后 划去◎ 所在列(行)的其它0元素，记作Ø ；这表示这列所代表的任务已指派 完，不必再考虑别人了。 (2)给只有一个0元素的列(行)中的0元素加圈，记作◎；然后划去◎ 所 在行的0元素，记作Ø ． (3)反复进行(1)，(2)两步，直到尽可能多的0元素都被圈出和划掉为止。
(4)若仍有没有划圈的0元素，且同行(列)的0元素至少有两个， 则从剩有0元素最少的行(列)开始，比较这行各0元素所在列中0元 素的数目，选择0元素少的那列的这个0元素加圈(表示选择性多的 要“礼让”选择性少的)。然后划掉同行同列的其它0元素。可反 复进行，直到所有0元素都已圈出和划掉为止。 （5）若◎ 元素的数目m 等于矩阵的阶数n，那么这指派问题的 最优解已得到。若m < n, 则转入下一步。 第三步：作最少的直线覆盖所有0元素。 (1)对没有◎的行打√号； (2)对已打√号的行中所有含Ø元素的列打√号； (3)再对打有√号的列中含◎ 元素的行打√号；
(4)重复(2)，(3)直到得不出新的打√号的行、列为止； (5)对没有打√号的行画横线，有打√号的列画纵线，这就得到覆盖 所有0元素的最少直线数 l 。l 应等于m，若不相等，说明试指派过 程有误，回到第二步(4)，另行试指派；若 l＝m < n，须再变换当前 的系数矩阵，以找到n个独立的0元素，为此转第四步。 第四步：变换矩阵(bij)以增加0元素。 在没有被直线覆盖的所有元素中找出最小元素，然后打√各行都减去 这最小元素；打√各列都加上这最小元素（以保证系数矩阵中不出现 负元素）。新系数矩阵的最优解和原问题仍相同。转回第二步
匈牙利算法用于求解标准指派问题
c.matlab求解
• 例: 求解下列指派问题，已知指派矩阵为
• 解：编写Matlab程序如下：
c=[3 8 2 10 3;8 7 2 9 7;6 4 2 7 5; 8 4 2 3 5;9 10 6 9 10];
c=c(:); %把矩阵c转化为向量
a=zeros(10,25);
for i=1:5 %实现循环运算
a(i,(i-1)*5+1:5*i)=1;
a(5+i,i:5:25)=1;
end %此循环把指派问题转化为线性规划问题 b=ones(10,1);
[x,y]=linprog(c,[],[],a,b,zeros(25,1),ones(25,1));
X=reshape(x,5,5)
opt=y
求得最优指派方案是 x 15 = x 23 = x 32 = x 51 = 1 x_{15}=x_{23}=x_{32}=x_{51}=1 x15=x23=x32=x51=1,最优解为21.
4.旅行商问题
本问题可以用图论中的赋权图来解，这里用纯粹的0-1整数规划来解
对每一对城市
v
i
,
v
j
v_i,v_j
vi,vj,定义变量
x
i
j
x_{ij}
xij来表示是否要从
v
i
v_i
vi出发访问
v
j
v_j
vj：
其中
i
,
j
=
1
,
2
,
…
，
n
i,j=1,2,…，n
i,j=1,2,…，n
目标函数是总旅费最小
M的选定是为了防止从自己本来所在地走到本来所在地（不让从自己走到自己）
约束条件
为了防止遍历过程出现一个子回路，即回不到出发点，附加一个条件：
综上，得到0-1规划模型如下
第三个约束条件的证明可以用反证法证明，其中
u
u
u可以不是整数
5.固定费用问题
对于每种方式的成本可以写作：
则目标函数可以写作：
引入0-1变量
y
i
y_i
yi
则约束条件为：
ε
\varepsilon
ε是充分小实数，
M
M
M是无穷大实数
6.蒙特卡洛法求解非线性整数规划
蒙特卡洛法代码如下：
clc, clear
%rng('shuffle') %根据当前时间为随机数生成器提供种子
rng(0) %进行一致性比较，每次产生的随机数是一样的
p0=0; n=10^6; tic %计时开始
for i=1:n
x=randi([0,99],1,5); %产生一行五列的区间[0,99]上的随机整数
[f,g]=mengte(x);
if all(g<=0)
if p0<f
x0=x; p0=f; %记录下当前较好的解
end
end
end
x0, p0, toc %计时结束
function [f,g]=mengte(x); %定义目标函数和约束条件
f=x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)-8*x(1)-2*x(2)-3*x(3)-...
x(4)-2*x(5);
g=[sum(x)-400
x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800
2*x(1)+x(2)+6*x(3)-200
x(3)+x(4)+5*x(5)-200];
end
用Lingo软件可以求得最优解，程序如下
model:
sets:
row/1..4/:b;
col/1..5/:c1,c2c,x;
link(row,col):a;
endsets
data:
c1=1,1,3,4,2;
c2=-8,-2,-3,-1,-2;
a=1 1 1 1 1
1 2 2 1 6
2 1 6 0 0
0 0 1 1 5;
b=400,800,200,200;
enddata
max=@sum(col:c1*x^2+c2*x);
@for(row(i):@sum(col(j):a(i,j)*x(j))<b(i));
@for(col:@gin(x));
@for(col:@bnd(0,x,99));
end
求得全局最优解 x 1 = 50 , x 2 = 99 , x 3 = 0 , x 4 = 99 , x 5 = 20 x_1=50,x_2=99,x_3=0,x_4=99,x_5=20 x1=50,x2=99,x3=0,x4=99,x5=20,最优值 z = 51568 z=51568 z=51568
410
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
