..:: 科学软件专区-规划,例子,整数-一个整数线性规划的例子(分枝定界法)[有代码]，{Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛-(手机验证注册)}

qibbxxt 发表于 2010-8-31 13:40:54

一个整数线性规划的例子(分枝定界法)[有代码]

本帖最后由 messenger 于 2010-9-6 16:14 编辑

整数现行规划用lingo等很多软件都可以解，Matlab没有现成的函数（有0-1规划的函数bintprog）,一般的整数规划是用分支定界法来解，有个好处就是可以求出所有解，而lingo只能求出一个解，下面是程序，其中分支定界的程序源于网络（具体的地址忘了）

这是原问题的模型，matlab代码如下：clear;clc;close all
N=7;
pos0=;
p=ones(1,N);
p(pos0)=0;
c=hankel(1:N,0:N-1);
a=fliplr(p(c));
b=-;
t=0.4;
h = text;
h1 = text;
set(h,'Interpreter','latex');
set(h,'string','$$\overbrace{}$$','Rotation',90,'position',,'fontsize',90);
set(h1,'Interpreter','latex','string','min=$$\sum_{i=1}^{7} x_{i}$$','position',,'fontsize',20,'Rotation',0);
for i=1:N
ii=find(a(i,:));
S=sprintf('\\fontsize{20}x_{%d}+x_{%d}+x_{%d}+x_{%d}+x_{%d}>=%d',ii(1),ii(2),ii(3),ii(4),ii(5),-b(i));
text(0.5,t,S);
t=t-0.1;
end
S='\fontsize{20}x_{i}>=0';
text(0.5,t,S);
axis()
axis off以下是求解的程序代码 clear;clc;close all
N = 7;
f = ones(N,1);
N=7;
M=;
a=diag(-ones(1,N-M(1)+1),M(1)-1)+diag(-ones(1,N-M(2)+1),M(2)-1)+...
diag(-ones(1,M(3)),M(3)-N)+diag(-ones(1,M(4)),M(4)-N)+ones(N);
A=-a;
b=-;
lb = zeros(N,1);
% = linprog(f,A,b,[],[],lb);
Aeq=[];beq=[];
LB=lb;UB=[];
=ILp(f,A,b,Aeq,beq,LB,UB,[],1,[]) ;
xn=round(xn);
function =ILp(f,G,h,Geq,heq,lb,ub,x,id,options)
%整数线性规划分支定界法，可求解纯整数规划和混合整数规划。
%y=minf’*x s.t. G*x<=h Geq*x=heq x为全整数或混合整数列向量
%用法
%=ILp(f,G,h,Geq,heq,lb,ub,x,id,options)
%参数说明
%lb:解的下界列向量（Default:-int）
%ub:解的上界列向量（Default:int）
%x:迭代初值列向量
%id：整数变量指标列向量，1-整数，0-实数（Default:1）
global upper opt c x0 A b Aeq beq ID options;
warning off all;
if nargin<10
options=optimset({});
options.Display='off';
options.LargeScale='off';
end
if nargin<9,id=ones(size(f));end
if nargin<8,x=[];end
if nargin<7 ||isempty(ub),ub=inf*ones(size(f));end
if nargin<6 ||isempty(lb),lb=zeros(size(f));end
if nargin<5,heq=[];end
if nargin<4,Geq=[];end
upper=inf;c=f;x0=x;A=G;b=h;Aeq=Geq;beq=heq;ID=id;
ftemp=ILP(lb(:),ub(:));
x=opt;y=upper;
%下面是子函数
function ftemp=ILP(vlb,vub)
warning off all;
global upper opt c x0 A b Aeq beq ID options;
=linprog(c,A,b,Aeq,beq,vlb,vub,x0,options);
if how <=0
return;
end;
if ftemp-upper>0.00005 %in order to avoid error
return;
end;
if max(abs(x.*ID-round(x.*ID)))<0.00005
if upper-ftemp>0.00005 %in order to avoid error
   opt=x';upper=ftemp;
   return;
else
   opt=;
   return;
end;
end;
notintx=find(abs(x-round(x))>=0.00005); %in order to avoid error
intx=fix(x);tempvlb=vlb;tempvub=vub;
if vub(notintx(1,1),1)>=intx(notintx(1,1),1)+1;
tempvlb(notintx(1,1),1)=intx(notintx(1,1),1)+1;
ftemp=ILP(tempvlb,vub);
end;
if vlb(notintx(1,1),1)<=intx(notintx(1,1),1)
tempvub(notintx(1,1),1)=intx(notintx(1,1),1);
ftemp=ILP(vlb,tempvub);
end

shamohu 发表于 2010-8-31 18:08:32

和1stOpt对比下：
代码：
MultiRun = 15;
IntParameter x(7);
Constant c=;
MaxFunction Sum(i=1:7,x)(x);
For(i=1:7)(Sum(j=i:i+4)(x)<=c);
用线性算法只有一个结果，用全局算法，自动运行15次得14组结果：

1: 1,4,12,0,7,5,2; 目标函数 = 31
2: 1,4,12,0,9,3,2; 目标函数 = 31
3: 1,8,8,0,10,2,2; 目标函数 = 31
4: 3,4,12,0,10,2,0; 目标函数 = 31
5: 1,4,12,0,8,4,2; 目标函数 = 31
6: 1,8,8,0,6,6,2; 目标函数 = 31
7: 2,6,10,0,6,6,1; 目标函数 = 31
8: 2,8,8,0,6,6,1; 目标函数 = 31
9: 2,3,13,0,6,6,1; 目标函数 = 31
10: 2,7,9,0,11,1,1; 目标函数 = 31
11: 1,8,8,0,7,5,2; 目标函数 = 31
12: 1,6,10,0,8,4,2; 目标函数 = 31
13: 2,8,8,0,10,2,1; 目标函数 = 31
14: 3,8,8,0,11,1,0; 目标函数 = 31

qibbxxt 发表于 2010-8-31 19:58:19

本帖最后由 qibbxxt 于 2010-8-31 20:00 编辑

2# shamohu
不好意思，那会太着急，把模型写错了，现在修改了一下，我计算的结果如下：
xn =
8 0 12 0 11 5 0
7 1 12 0 11 5 0
7 1 12 0 11 4 1
6 2 12 0 11 5 0
6 2 12 0 11 4 1
6 2 12 0 11 3 2
5 3 12 0 11 5 0
5 3 12 0 11 4 1
5 3 12 0 11 3 2
5 3 12 0 11 2 3
4 4 12 0 11 5 0
4 4 12 0 11 4 1
4 4 12 0 11 3 2
4 4 12 0 11 2 3
4 4 12 0 11 1 4
3 5 12 0 11 5 0
3 5 12 0 11 4 1
3 5 12 0 11 3 2
3 5 12 0 11 2 3
3 5 12 0 11 1 4
3 5 12 0 11 0 5
2 6 12 0 11 5 0
2 6 12 0 11 4 1
2 6 12 0 11 3 2
2 6 12 0 11 2 3
2 6 12 0 11 1 4
2 6 12 0 11 0 5
1 7 12 0 11 5 0
1 7 12 0 11 4 1
1 7 12 0 11 3 2
1 7 12 0 11 2 3
1 7 12 0 11 1 4
1 7 12 0 11 0 5
0 8 12 0 11 5 0
0 8 12 0 11 4 1
0 8 12 0 11 3 2
0 8 12 0 11 2 3
0 8 12 0 11 1 4
0 8 12 0 11 0 5

最小值均为36
另外有空我也学学1stOpt,看来挺好用的

lin2009 发表于 2010-9-1 08:32:08

3# qibbxxt
2＃给出的最小值为31，而不是36。
另外，用1stopt 求解，可能会遗漏正确的解。

shamohu 发表于 2010-9-1 08:58:11

一楼好像把模型改动了，原先目标函数是求最大，约束都是小于等于。因此2楼求的是最大，而3楼是最小，约束也不同。现修改与目前的一楼一样：
MultiRun = 100;
IntParameter x(7);
Constant c=;
MinFunction Sum(i=1:7,x)(x);
For(i=1:7)(Sum(j=i:i+4)(x)>=c);
自动运行100次，得到37组不同结果：
1: 5,3,12,0,11,4,1; 目标函数 = 36
2: 1,7,12,0,11,0,5; 目标函数 = 36
3: 1,7,12,0,11,2,3; 目标函数 = 36
4: 3,5,12,0,11,0,5; 目标函数 = 36
5: 0,8,12,0,11,4,1; 目标函数 = 36
6: 5,3,12,0,11,3,2; 目标函数 = 36
7: 2,6,12,0,11,1,4; 目标函数 = 36
8: 6,2,12,0,11,4,1; 目标函数 = 36
9: 6,2,12,0,11,5,0; 目标函数 = 36
10: 5,3,12,0,11,5,0; 目标函数 = 36
11: 1,7,12,0,11,1,4; 目标函数 = 36
12: 7,1,12,0,11,4,1; 目标函数 = 36
13: 4,4,12,0,11,5,0; 目标函数 = 36
14: 1,7,12,0,11,3,2; 目标函数 = 36
15: 2,6,12,0,11,4,1; 目标函数 = 36
16: 4,4,12,0,11,4,1; 目标函数 = 36
17: 3,5,12,0,11,5,0; 目标函数 = 36
18: 3,5,12,0,11,1,4; 目标函数 = 36
19: 2,6,12,0,11,0,5; 目标函数 = 36
20: 6,2,12,0,11,3,2; 目标函数 = 36
21: 3,5,12,0,11,4,1; 目标函数 = 36
22: 0,8,12,0,11,2,3; 目标函数 = 36
23: 4,4,12,0,11,1,4; 目标函数 = 36
24: 2,6,12,0,11,5,0; 目标函数 = 36
25: 1,7,12,0,11,5,0; 目标函数 = 36
26: 3,5,12,0,11,2,3; 目标函数 = 36
27: 4,4,12,0,11,2,3; 目标函数 = 36
28: 7,1,12,0,11,5,0; 目标函数 = 36
29: 8,0,12,0,11,5,0; 目标函数 = 36
30: 0,8,12,0,11,0,5; 目标函数 = 36
31: 1,7,12,0,11,4,1; 目标函数 = 36
32: 0,8,12,0,11,3,2; 目标函数 = 36
33: 2,6,12,0,11,2,3; 目标函数 = 36
34: 0,8,12,0,11,1,4; 目标函数 = 36
35: 0,8,12,0,11,5,0; 目标函数 = 36
36: 3,5,12,0,11,3,2; 目标函数 = 36
37: 2,6,12,0,11,3,2; 目标函数 = 36

qibbxxt 发表于 2010-9-1 09:42:55

4# lin2009
对比了一下，1stOpt比matlab少两组解

5 3 12 0 11 2 3
4 4 12 0 11 3 2

对比的代码如下：

clear;clc;close all
stOpt_re=[5,3,12,0,11,4,1
1,7,12,0,11,0,5
1,7,12,0,11,2,3
3,5,12,0,11,0,5
0,8,12,0,11,4,1
5,3,12,0,11,3,2
2,6,12,0,11,1,4
6,2,12,0,11,4,1
6,2,12,0,11,5,0
5,3,12,0,11,5,0
1,7,12,0,11,1,4
7,1,12,0,11,4,1
4,4,12,0,11,5,0
1,7,12,0,11,3,2
2,6,12,0,11,4,1
4,4,12,0,11,4,1
3,5,12,0,11,5,0
3,5,12,0,11,1,4
2,6,12,0,11,0,5
6,2,12,0,11,3,2
3,5,12,0,11,4,1
0,8,12,0,11,2,3
4,4,12,0,11,1,4
2,6,12,0,11,5,0
1,7,12,0,11,5,0
3,5,12,0,11,2,3
4,4,12,0,11,2,3
7,1,12,0,11,5,0
8,0,12,0,11,5,0
0,8,12,0,11,0,5
1,7,12,0,11,4,1
0,8,12,0,11,3,2
2,6,12,0,11,2,3
0,8,12,0,11,1,4
0,8,12,0,11,5,0
3,5,12,0,11,3,2
2,6,12,0,11,3,2];
matlab_re=[8 0 12 0 11 5 0
7 1 12 0 11 5 0
7 1 12 0 11 4 1
6 2 12 0 11 5 0
6 2 12 0 11 4 1
6 2 12 0 11 3 2
5 3 12 0 11 5 0
5 3 12 0 11 4 1
5 3 12 0 11 3 2
5 3 12 0 11 2 3
4 4 12 0 11 5 0
4 4 12 0 11 4 1
4 4 12 0 11 3 2
4 4 12 0 11 2 3
4 4 12 0 11 1 4
3 5 12 0 11 5 0
3 5 12 0 11 4 1
3 5 12 0 11 3 2
3 5 12 0 11 2 3
3 5 12 0 11 1 4
3 5 12 0 11 0 5
2 6 12 0 11 5 0
2 6 12 0 11 4 1
2 6 12 0 11 3 2
2 6 12 0 11 2 3
2 6 12 0 11 1 4
2 6 12 0 11 0 5
1 7 12 0 11 5 0
1 7 12 0 11 4 1
1 7 12 0 11 3 2
1 7 12 0 11 2 3
1 7 12 0 11 1 4
1 7 12 0 11 0 5
0 8 12 0 11 5 0
0 8 12 0 11 4 1
0 8 12 0 11 3 2
0 8 12 0 11 2 3
0 8 12 0 11 1 4
0 8 12 0 11 0 5];
matlabIsIn1stOPt=ismember(matlab_re,stOpt_re,'rows');
stOPtIsInmatlab=ismember(stOpt_re,matlab_re,'rows');
re=matlab_re(not(matlabIsIn1stOPt),:);

lin2009 发表于 2010-9-1 10:04:58

6# qibbxxt
"对比的代码"可以直接用集合差运算：

setdiff(matlab_re,stOpt_re,'rows')
ans =
4 4 12 0 11 3 2
5 3 12 0 11 2 3

shamohu 发表于 2010-9-1 10:10:16

1stOpt的确不能保证给出全部最优解。估计它用的全局最优算法，初值都是随机赋的。

1stOpt的代码是很精短且易懂的。

qibbxxt 发表于 2010-9-1 10:40:23

7# lin2009
恩，互相作差，这样就全面一些

qibbxxt 发表于 2010-9-1 10:43:23

8# shamohu
恩，matlab由于基于矩阵操作，所以用起来稍麻烦一些
1stOpt在做规划问题时代码简洁易懂，就和lingo一样，但是就是很难保证找到所有解

upc1984 发表于 2010-9-1 10:56:35

本帖最后由 upc1984 于 2010-9-1 10:59 编辑

1# qibbxxt

我新学的也编了一个
function y=myfun0901(x)
y=x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7);

A=[1 0 0 1 1 1 1;...
   1 1 0 0 1 1 1;...
   1 1 1 0 0 1 1;...
   1 1 1 1 0 0 1;...
   1 1 1 1 1 0 0;...
   0 1 1 1 1 1 1;...
   0 0 1 1 1 1 1];
A=-A;
b=;
b=-b';
lb=;
lu=;
x0=;
x=fmincon(@myfun0901,x0,A,b,[],[],lb,lu)
计算结果
Active inequalities (to within options.TolCon = 1e-006):
lower upper ineqlin ineqnonlin
4 2
3
5
7

x =

4.0022 3.9978 12.0000 0 11.0000 2.8835 2.1165

y =

36.0000

不过有小数，不好办~

qibbxxt 发表于 2010-9-1 11:08:17

11# upc1984
1.这个问题是整数线性规划，你用的是非线性规划
2.matlab没有自带的求解整数规划的函数（0-1规划除外），所以整数线性规划用分支定界法
3.建议你好好看下运筹学关于规划部分

bainhome 发表于 2010-9-4 01:34:03

本帖最后由 bainhome 于 2010-9-4 01:37 编辑

嗯，整数线规用非线性方法做是个糟糕的误区，以前初学时我也犯过这种错误，这是应该去避免的。
这个分枝定界法的程序看起很不错的样子，以前找过整数线规的程序都不是很好用，意外看到这个，谢谢qibbxxt。同时建议加技术分：整数线规的MATLAB程序在版内好像没有出现过，加上应用的讨论则肯定没有，尤其在特殊的个案中和1stopt比较时能不落下风，难能可贵。

qibbxxt 发表于 2010-9-4 10:37:18

13# bainhome
谢谢bainhome兄，这个程序是不好找，我以前也找过很多，都用不了

页: [1]

Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛---(手机验证注册)'s Archiver

一个整数线性规划的例子(分枝定界法)[有代码]