作者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