作者lifeisshine (芋头)
看板MATLAB
标题Re: [运算] 最佳"整数"解
时间Mon Sep 2 22:41:47 2013
※ 引述《lifeisshine (芋头)》之铭言：
: 请教各位先进们，
: 现有二元一次方程式:
: f(x,y) = 4055x + 387y
: 小弟要如何下指令计算
: 当x和y分别等於多少时(x和y都要是"整数")
: f(x,y)小於等於且最接近9000
: 这题答案是当 x=2, y=2 时，
: f(x,y) = 8884 最接近9000，
: 换句话说，9000-8884=116 为最小的残余量!
: ----------
: 目前使用左除法，以矩阵的形式解答，
: 但得到的解当然是非整数 x = 2.2195, y = 0
: 恳请各位先进帮解惑Orz
後来发现这样的一个数学问题，
可以用整数线性规划来求得，
说穿了就是最佳化(优化)，
Matlab 线性规划的指令为 linprog，
也在板上爬过文，不过遗憾的是，
没办法求"整数"线性规划，
因此我找到许多网友利用 linprog，
并加写其他 code 改成"整数"线性规划，
目前也有 Copy 一些来测试，不过尚未成功...
以下 code 我也是Copy来的，不过input多了一个n，
让我百思不得其解阿...
希望有接触过这方面数学演算的前辈们，
能提点小弟如何运用此 Code
-------------------------------------------------
实际应用题:
每个桌子:需木工 4时，漆工 2时，利润 5
每个椅子:需木工 1时，漆工 2时，利润 4
木工上限为 60 时，漆工上限为 48 时，
令 x1 为桌子的生产个数，令 x2 为椅子的生产个数
-------------------------------------------------
数学模型:
f ＝ 5x1＋ 4x2 ……(利润最大化)
s.t. 4x1 ＋ 1x2 ≦ 60 …….(木工时数限制)
2x1 ＋ 2x2 ≦ 48 …….(漆工时数限制)
x1 ≧ 0， x2 ≧ 0
-------------------------------------------------
将以上数学式转换成矩阵形式，
原 Matlab 内建线性规划函数→linprog(f,a,b,aeq,beq,lb,ub);
因为linprog是求函数最小值，所以当要求最大值时需将函数乘上-1 (you know..)
f = [-5 ; -4];
b = [60 ; 48];
a = [4 1; 2 2];
aeq = [];
beq = [];
lb = [];
ub = [];
答案是 x1=12, x2=12
但是被网友改写後，多了一个Input参数"n"，
小弟就是卡在这Orz
-------------------------------------------------
网路上所 Copy 程式码:
function [x,val]=fzdj(n,f,a,b,aeq,beq,lb,ub)
x=zeros(n,1);
x1=zeros(n,1);
m1=2;
m2=1;
[x1,val1]=linprog(f,a,b,aeq,beq,lb,ub);
if (x1==0)
x=x1;
val=val1;
else
e1={0,a,b,aeq,beq,lb,ub,x1,val1};
e(1,1)={e1}; zl=0; zu=-val1;
while (zu~=zl)
for c=1:1:m2
if (m1~=2)
if (cell2mat(e{m1-1,c}(1))==1)
e1={1,[],[],[],[],[],[],[],0};
e(m1,c*2-1)={e1};
e(m1,c*2)={e1};
continue;
end;
end;
x1=cell2mat(e{m1-1,c}(8));
x2=zeros(n,1);
s=0;
s1=1;
s2=1;
lb1=cell2mat(e{m1-1,c}(6));
ub1=cell2mat(e{m1-1,c}(7));
lb2=cell2mat(e{m1-1,c}(6));
ub2=cell2mat(e{m1-1,c}(7));
for d=1:1:n
if (abs((round(x1(d))-x1(d)))>0.0001)&(s==0)
s=1;
lb1(d)=fix(x1(d))+1;
if (a*lb1<=b)
s1=0;
end;
ub2(d)=fix(x1(d));
if (a*lb2<=b)
s2=0;
end;
end;
end;
e1={s1,a,b,aeq,beq,lb1,ub1,[],0};
e2={s2,a,b,aeq,beq,lb2,ub2,[],0};
e(m1,c*2-1)={e1}; e(m1,c*2)={e2};
end;
m1=m1+1;
m2=m2*2;
for c=1:1:m2
if (cell2mat(e{m1-1,c}(1))==0)
[x1,val1]=linprog(f,cell2mat(e{m1-1,c}(
2)),cell2mat(e{m1-1,c}(3)),cell2mat(e{m1-1,c}(4)),cell2mat(e{m1-1,c}(5)),cell2mat(e{m1-1,c}(6)),cell2mat(e{m1-1,c}(7)));
e1={cell2mat(e{m1-1,c}(1)),cell2mat(e{m1-1,c}(2)),cell2mat(e{m1-1,c}(3)),cell2mat(e{m1-1,c}(4)),cell2mat(e{m1-1,c}(5)),cell2mat(e{m1-1,c}(6)),cell2mat(e{m1-1,c}(7)),x1,val1};
e(m1-1,c)={e1};
end;
z=val1;
if ((-z)<(-zl))
e1={1,[],[],[],[],[],[],[],0};
e(m1-1,c)={e1};
elseif (abs(round(x1)-x1)<=0.0001)
zl=z;
end;
end;
for c=1:1:m2
if (cell2mat(e{m1-1,c}(1))==0)
zu=cell2mat(e{m1-1,c}(9));
end;
end;
for c=1:1:m2
if (-cell2mat(e{m1-1,c}(9))>(-zu))
zu=cell2mat(e{m1-1,c}(9));
end;
end;
end;
for c=1:1:m2
if (cell2mat(e{m1-1,c}(1))==0)&(cell2mat(e{m1-1,c}(9))==zu)
x=cell2mat(e{m1-1,c}(8));
end;
end;
val=zu;
end;
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 115.80.134.155