標題: Re: [求救]Mathematica求有限制式的極大值指令問題
時間: Fri Oct 23 19:54:05 2009
※ 引述《sbtpc (翻身)》之銘言：
: 這題我已經卡很多天了
: 想請高手幫個忙
: 謝謝~~
: 題目我用圖片檔
: http://yfrog.com/0j74082627j
: 我的想法是用以下這個指令
: Maximize[{目標函數,限制式},{變數}]
: 但都說找不到符合的解
: 麻煩高手了...
Maximize比較適合用在你的解是有close form的時候！這邊比較適合用NMaximze 或是
FindMaximum
myrule = x[i_, j_] :> ToExpression["x" <> ToString[i] <> ToString[j]];
obj = Plus @@
Flatten[Table[-Log[
E^-x[i, j] Sqrt[2 \[Pi]] x[i, j]^(0.5 + x[i, j])], {i, 5}, {j,
5}], 1] //. myrule;
vars = Flatten[Table[ x[i, j], {i, 5}, {j, 5}], 1] //. myrule;
cons = Flatten[{Thread[
Table[Plus @@ Table[x[i, j], {j, 5}], {i, 5}] <= {5.78, 2.61,
1.71, 0.12, 0.27}],
Thread[Table[Plus @@ Table[x[i, j], {i, 5}], {j, 5}] >= {5.94,
1.55, 2.54, 0.23, 0.23}],
Flatten[Table[x[i, j], {i, 5}, {j, 5}], 1].{14.7, 16.44, 32.39,
58.68, 51.62, 20.47, 13.88, 21.15, 57.25, 52.78, 33.69, 19.14,
11.8, 51.11, 50.25, 42.45, 34.69, 41.82, 15.04, 46.05, 42.69,
28.8, 40.9, 42.8, 14.43} <= 177.42,
Thread[Flatten[Table[x[i, j], {i, 5}, {j, 5}], 1] >= 10^(-10)]},
1] //. myrule;
FindMinimum[{obj, cons}, vars]
{10.966, {x11 -> 5.54243, x12 -> 0.0864126, x13 -> 0.0679551,
x14 -> 0.0427504, x15 -> 0.0404569, x21 -> 0.262848, x22 -> 1.30528,
x23 -> 0.9365, x24 -> 0.0559572, x25 -> 0.049412, x31 -> 0.0670608,
x32 -> 0.0844314, x33 -> 1.46502, x34 -> 0.0507783,
x35 -> 0.0427054, x41 -> 0.0225331, x42 -> 0.022662,
x43 -> 0.0228854, x44 -> 0.0309131, x45 -> 0.0210063,
x51 -> 0.0451326, x52 -> 0.0512112, x53 -> 0.0476358,
x54 -> 0.049601, x55 -> 0.0764195}}
NMinimize[{obj, cons}, vars,
Method -> {"NelderMead",
"PenaltyFunction" -> Function[{d, i}, 10^i*d^2],
"InitialPoints" ->
Table[Flatten[Table[ Random[], {i, 5}, {j, 5}]], {26}]}]
{10.966, {x11 -> 5.54243, x12 -> 0.0864126, x13 -> 0.0679551,
x14 -> 0.0427504, x15 -> 0.0404569, x21 -> 0.262848, x22 -> 1.30528,
x23 -> 0.9365, x24 -> 0.0559572, x25 -> 0.049412, x31 -> 0.0670608,
x32 -> 0.0844314, x33 -> 1.46502, x34 -> 0.0507783,
x35 -> 0.0427054, x41 -> 0.0225331, x42 -> 0.022662,
x43 -> 0.0228854, x44 -> 0.0309131, x45 -> 0.0210063,
x51 -> 0.0451326, x52 -> 0.0512112, x53 -> 0.0476358,
x54 -> 0.049601, x55 -> 0.0764194}}
最小值應該是確定可以找到！
最大值就有點問題！
NMaximize[{obj, cons}, vars, MaxIterations -> 2000,
Method -> {"NelderMead", "PostProcess" -> True,
"PenaltyFunction" -> Function[{d, i}, 10 i*d^2],
"InitialPoints" ->
Table[Flatten[Table[ Random[], {i, 5}, {j, 5}]], {26}]}]
{44.8349, {x11 -> 4.20224, x12 -> 1.2057, x13 -> 0.313754,
x14 -> 0.0569259, x15 -> 0.00138137, x21 -> 1.73577,
x22 -> 0.339281, x23 -> 0.45701, x24 -> 0.0486815, x25 -> 0.0292547,
x31 -> 0.00103741, x32 -> 0.00170413, x33 -> 1.70649,
x34 -> 0.000700329, x35 -> 0.000064355, x41 -> 0.000847514,
x42 -> 2.03391*10^-10, x43 -> 0.00169768, x44 -> 0.117455,
x45 -> 2.78825*10^-8, x51 -> 0.0000991133, x52 -> 0.00331898,
x53 -> 0.0610449, x54 -> 0.00623745, x55 -> 0.1993}}
這邊要多試幾次∼∼
--
養花種魚數月亮賞星星
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.130.237
我想請問
回覆刪除我這樣到底哪裡有問題
為什麼一直跑不出來
SR = 5; SM1 = 2; SM2 = 2; c = 1; h = 0.05; H = 15; \[Lambda] = -0.1; \
Subscript[M, 1] = 0.7; Subscript[M, 2] = 0.3; \[Beta] = 0.25; \
\[Theta] = 0.05; \[CapitalDelta] = 0.05; \[Delta] = 0.2;
Subscript[c, 1] = c;
Subscript[c, 2] = c*\[Delta];
piB1 = \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(j\), \(k\)]\(\((
\*SubscriptBox[\(p\), \(2\)] -
\*SubscriptBox[\(w\), \(2\)]*Exp[\[Theta]*\((t - j)\)] - h*\(
\*SubsuperscriptBox[\(\[Integral]\), \(j\), \(t\)]Exp[\[Theta]*\((t -
u)\)] \[DifferentialD]u\))\)*\((
\*SubscriptBox[\(M\), \(2\)] -
\*SubscriptBox[\(p\), \(2\)] + \[Beta]*
\*SubscriptBox[\(p\), \(1\)])\)*
Exp[\[Lambda]*t] \[DifferentialD]t\)\);(*零售商再製品利潤*)
piB2 = \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(j\), \(k\)]\(\((
\*SubscriptBox[\(p\), \(1\)] -
\*SubscriptBox[\(w\), \(1\)]*Exp[\[Theta]*\((t - j)\)] - h*\(
\*SubsuperscriptBox[\(\[Integral]\), \(j\), \(t\)]Exp[\[Theta]*\((t -
u)\)] \[DifferentialD]u\))\)*\((
\*SubscriptBox[\(M\), \(1\)] -
\*SubscriptBox[\(p\), \(1\)] + \[Beta]*
\*SubscriptBox[\(p\), \(2\)])\)*
Exp[\[Lambda]*t] \[DifferentialD]t\)\);(*零售商新品利潤*)
piB = piB1 + piB2 - SR;(*零售商總利潤，減固定成本*)
funpidpB1 = D[piB, Subscript[p, 1]];(*對p1作一階*)
funpidpB2 = D[piB, Subscript[p, 2]];(*對p2作一階*)
solpidpB1 =
Solve[{funpidpB1 == 0, funpidpB2 == 0}, {Subscript[p, 1], Subscript[
p, 2]}];(*求p1和p2*)
ansp1 = Subscript[p, 1] /. solpidpB1;
ansp2 = Subscript[p, 2] /. solpidpB1;
piV1 = \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(j\), \(k\)]\(\((
\*SubscriptBox[\(w\), \(1\)] -
\*SubscriptBox[\(c\), \(1\)])\)*\((
\*SubscriptBox[\(M\), \(1\)] -
\*SubscriptBox[\(p\), \(1\)] + \[Beta]*
\*SubscriptBox[\(p\), \(2\)])\)*
Exp[\[Lambda]*t + \[Theta]*\((t - j)\)] \[DifferentialD]t\)\) -
SM1;(*製造商新品利潤*)
piV2 = \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(j\), \(k\)]\(\((
\*SubscriptBox[\(w\), \(2\)] -
\*SubscriptBox[\(c\), \(2\)])\)*\((
\*SubscriptBox[\(M\), \(2\)] -
\*SubscriptBox[\(p\), \(2\)] + \[Beta]*
\*SubscriptBox[\(p\), \(1\)])\)*
Exp[\[Lambda]*t + \[Theta]*\((t - j)\)] \[DifferentialD]t\)\) -
SM2;(*製造商再製品利潤*)
funwp1 = piV2 /. {Subscript[p, 2] -> ansp2};(*把p2代入*)
funwp2 = funwp1 /. {Subscript[p, 1] -> ansp1};(*把p1代入*)
pidwpV1 = D[funwp2, Subscript[w, 2]];(*對w2做一階*)
solpidwV1 = Solve[pidwpV1 == 0, Subscript[w, 2]];(*求w2*)
ansp3 = Subscript[w, 2] /. solpidwV1;(*將w2替代為anps3*)
funwp3 = piV1 /. {Subscript[p, 2] -> ansp2};(*把p2代入*)
funwp4 = funwp3 /. {Subscript[p, 1] -> ansp1};(*把p1代入*)
funwp5 = funwp4 /. {Subscript[w, 2] -> ansp3};(*將w2代回利潤式*)
pidwpV2 = D[funwp5, Subscript[w, 1]];(*對w1做一階*)
solpidwV2 = Solve[pidwpV2 == 0, Subscript[w, 1]];(*求出w1*)
funwp6 = funwp2 /. solpidwV2;(*帶回利潤式*)
pidwpV3 = D[funwp6, Subscript[w, 2]];(*對w2做一階*)
solpidwV3 = Solve[pidwpV3 == 0, Subscript[w, 2]];(*重新求w2，進行比對驗證之用*)
j = 0;
For[k = 1, k <= H, k = k + 1,
optw1 = Maximize[{funwp5, 0 <= Subscript[w, 1] <= ansp2}, Subscript[
w, 1]];
Subscript[WholeSalePrice1, j, k] = Subscript[w, 1] /. Last[optw1];