myConstrains[constrains_] := 
 Block[{num, s, mycon}, num = Length[constrains];
  mycon = 
   Table[If[Head[constrains[[i]]] === LessEqual, 
     constrains[[i]][[1]] + s[i] == constrains[[i]][[2]], 
     If[Head[constrains[[i]]] === Equal, 
      constrains[[i]][[1]] + R[i] == constrains[[i]][[2]], 
      constrains[[i]][[1]] - s[i] + R[i] == 
       constrains[[i]][[2]]]], {i, num}];
  mycon /. x__[y_] :> ToExpression[ToString[x] <> ToString[y]]]

LPstd[obj_, mycons_] := 
 Block[{myvar, basicvars, nonbasicvars, AB, BM, c, A, B, rowindex1, 
   rowindex2}, myvar = Variables[mycons[[All, 1]]];
  myvar = 
   Flatten[{Sort@
      Select[myvar, Not@StringStartsQ[SymbolName[#], "s" | "R"] &], 
     Sort@Select[myvar, StringStartsQ[SymbolName[#], "s"] &], 
     Sort@Select[myvar, StringStartsQ[SymbolName[#], "R"] &]}];
  AB = {Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, Not@StringStartsQ[SymbolName[#], "s" | "R"] &], 
    Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, StringStartsQ[SymbolName[#], "s"] &], 
    Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, StringStartsQ[SymbolName[#], "R"] &]};
  AB = Transpose@Flatten[AB, 1];
  AB = Join[{Flatten[
      Table[Select[myvar, 
        StringStartsQ[SymbolName[#], i] &], {i, {"x", "s", "R"}}]]}, 
    AB];
  c = Flatten[{Table[
      Coefficient[obj, #] & /@ 
       Select[myvar, StringStartsQ[SymbolName[#], i] &], {i, {"x", 
        "s"}}], Table[-BM, {Length[
        Select[myvar, StringStartsQ[SymbolName[#], "R"] &]]}], 0}];
  basicvars = 
   Select[AB[[All, #]] & /@ (Flatten /@ 
        Map[Position[myvar, #] &, 
         Tuples[myvar, {Length[mycons]}], {2}]), #[[2 ;; -1]] === 
       IdentityMatrix[Length[mycons]] &][[1, 1]];
  nonbasicvars = 
   Flatten[Table[
     Select[Complement[myvar, basicvars], 
      StringStartsQ[SymbolName[#], i] &], {i, {"x", "s", "R"}}]];
  c = Flatten@{c[[Flatten[Position[myvar, #] & /@ nonbasicvars]]], 
     c[[Flatten[Position[myvar, #] & /@ basicvars]]], 0};
  AB = Join[{Flatten@{nonbasicvars, basicvars}}, 
    Flatten /@ 
     Transpose@{AB[[2 ;; -1]][[All, 
         Sort[Flatten[Position[myvar, #] & /@ nonbasicvars]]]], 
       IdentityMatrix[Length[mycons]]}];
  rowindex1 = c[[-Length[basicvars] - 1 ;; -2]];
  rowindex2 = Range[Length[mycons]];
  Flatten[{c - 
       Total[rowindex1[[#]]*
           Flatten[{AB[[2 ;; -1]][[rowindex2[[#]]]], 
             mycons[[rowindex2[[#]], 2]]}] & /@ 
         Range[Length@rowindex1]], 
      AB[[2 ;; -1]][[All, 
        Sort[Flatten[Position[myvar, #] & /@ nonbasicvars]]]], 
      IdentityMatrix[Length[mycons]], mycons[[All, 2]], 
      Flatten@{Sort@
         Select[nonbasicvars, StringStartsQ[SymbolName[#], "x"] &], 
        Sort@Select[nonbasicvars, 
          StringStartsQ[SymbolName[#], "s"] &], 
        Sort@Select[nonbasicvars, 
          StringStartsQ[SymbolName[#], "R"] &]}, basicvars} & /@ 
    Range[2], 1]]

simplexmethod[c_, A_, B_, b_, xnb_, xb_, cp_, Ap_, Bp_, bp_, xnbp_, 
  xbp_] := 
 Block[{cB, nA, nB, mysymbol, pivotCol, AB, leaveVarpos, Bnew, xbnew, 
   xnbnew, Anew, cnew, znew, z, BM, ABp},
   (*初始化變數*)
   nA = Length[xnb];
   nB = Length[xb];
   
  (*建立符號變數*)
  mysymbol = Flatten@{xnbp, xbp};

  (*找出進入變數（最負的檢驗數）*)
  pivotCol = PositionSmallest[-c[[1 ;; -2]] /. BM -> 10^6][[1]];
  
  (*合併A和B矩陣*)
  AB = Flatten /@ Transpose[{A, B}];
  Print@MatrixForm@AB;
  ABp = Flatten /@ Transpose[{Ap, Bp}];
  leaveVarpos = 
   PositionSmallest[
     b/(If[# <= 0, 10^-6, #] & /@ AB[[All, pivotCol]] /. BM -> 10^6) /. 
      x_ /; x <= 0 -> 10^6][[1]];
  
  (*更新基變數和非基變數*)
  xbnew = xb /. xb[[leaveVarpos]] -> mysymbol[[pivotCol]];
  xnbnew = xnb /. mysymbol[[pivotCol]] -> xb[[leaveVarpos]];
  Bnew = 
   ABp[[All, Flatten[Position[mysymbol, #] & /@ Flatten@{xbnew}]]];
  cB = cp[[Flatten[Position[mysymbol, #] & /@ xbnew]]];
  
  (*計算新的檢驗數和目標函數值*)
  znew = Flatten[{cB . Inverse[Bnew] . Ap - cp[[1 ;; nA]], 
      cB . Inverse[Bnew], cB . Inverse[Bnew] . bp - cp[[-1]]}] // 
    FullSimplify;
  
  (*顯示單形法表格*)
  Print[Grid[
    Insert[Flatten /@ 
      Transpose[{Flatten[{z, xbnew}], 
        Insert[Map[Flatten, 
          Transpose[{Inverse[Bnew] . Ap, Inverse[Bnew], 
            Inverse[Bnew] . bp}]], znew, 1]}], 
     Flatten[{Item["BasicVariable", 
        Background -> Lighter[Yellow, 0.9]], mysymbol, "Solution"}], 
     1], Frame -> All, 
    Background -> {{Lighter[Gray, 0.8]}, {Lighter[Yellow, 0.9], 
       Lighter[Pink, 0.9]}}, ItemStyle -> Directive[FontSize -> 16], 
    Spacings -> {Table[5, {Length[mysymbol] + 3}], 2}]];
  
  (*返回更新後的參數*)
  {-Flatten@{cB . Inverse[Bnew] . Ap - cp[[1 ;; nA]], 
       cB . Inverse[Bnew]}, Inverse[Bnew] . Ap, Inverse[Bnew], 
    Inverse[Bnew] . bp, xnbnew, xbnew, cp, Ap, Bp, bp, xnbp, xbp} // 
   FullSimplify]

Nest[simplexmethod @@ # &, 
 LPstd[-15 x1 - 20 x2, 
  myConstrains[{x1 + 2 x2 >= 10, 2 x1 - 3 x2 <= 6, x1 + x2 >= 6}]], 2]

繼續閱讀全文　Mathematica教學：Simplex method in tableau form

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

You can download searoute from https://github.com/eurostat/searoute .

port = {
 {"Shanghai, China", 31.2304, 121.4737},
 {"Guangzhou, China", 23.1291, 113.2644},
 {"Rotterdam, Netherlands", 51.9225, 4.4792},
 {"Antwerp, Belgium", 51.2194, 4.4025},
 {"Valencia, Spain", 39.4699, -0.3763},
 {"Los Angeles, USA", 33.7405, -118.2719},
 {"New York/New Jersey, USA", 40.6687, -74.0451},
 {"Savannah, USA", 32.0809, -81.0912},
 {"Santos, Brazil", -23.9608, -46.3336},
 {"Buenos Aires, Argentina", -34.6037, -58.3816},
 {"Itajai, Brazil", -26.9101, -48.6703},
 {"Tangier Med, Morocco", 35.8844, -5.5036},
 {"Port Said, Egypt", 31.2653, 32.3019},
 {"Durban, South Africa", -29.8587, 31.0218},
 {"Melbourne, Australia", -37.8136, 144.9631},
 {"Auckland, New Zealand", -36.8485, 174.7633},
 {"Brisbane, Australia", -27.4698, 153.0251},
 {"Singapore, Singapore", 1.3521, 103.8198},
 {"Kaohsiung, Taiwan", 22.6273, 120.3014},
 {"Busan, South Korea", 35.1796, 129.0756},
 {"Hong Kong, China", 22.3193, 114.1694},
 {"Tanjung Pelepas, Malaysia", 1.3673, 103.5386},
 {"Jebel Ali, UAE", 25.015, 55.0633},
 {"Colombo, Sri Lanka", 6.9271, 79.8612},
 {"Felixstowe, UK", 51.9556, 1.3111},
 {"Algeciras, Spain", 36.1408, -5.4562},
 {"Vancouver, Canada", 49.2827, -123.1207},
 {"Port Klang, Malaysia", 3.0089, 101.3678},
 {"Laem Chabang, Thailand", 13.0846, 100.8956},
 {"Tokyo, Japan", 35.6828, 139.7595},
 {"Yokohama, Japan", 35.4437, 139.6380}
};

Export["/Users/yourname/Desktop/searoute/text.csv",
Insert[
Flatten@{#[[1,1]]<>" to "<>#[[2,1]],
Round[#[[1,3]],0.1],Round[#[[1,2]],0.1],
Round[#[[2,3]],0.1],Round[#[[2,2]],0.1]}&/@
Select[
Flatten[Table[{port[[i]],port[[j]]},
{i,Length@port},{j,Length@port}],1],
#[[1]]!=#[[2]]&
],
{"routename","olon","olat","dlon","dlat"},1]
];

RunProcess[{"java","-jar",
"/Users/yourname/Desktop/searoute/searoute.jar",
"/Users/yourname/Desktop/searoute/text.csv",
"-res","5","-panama","1","-suez","1"}]

route=Import["/Users/yourname/Desktop/searoute/out.geojson",
"Data"][[2,2]];

mypath=Flatten@{#[[2,2]],StringSplit["routename"/.#[[3,2]]," to "],
"distKM"/.#[[3,2]]}&/@route;

myrule=Rule@@@Transpose@{Tally[mypath[[All,2]]][[All,1]],
Range[Length@port]};

distMatrix=ParallelTable[If[i==j,Infinity,
ToExpression@Select[mypath/.myrule,
#[[2]]==i&&#[[3]]==j&][[1,-1]]],
{i,Length@port},{j,Length@port}];

tour=FindShortestTour[Range[Length@port],
DistanceFunction->(distMatrix[[#1,#2]]&)];

myline=ParallelTable[
If[i==j,{},Select[mypath/.myrule,#[[2]]==i&&#[[3]]==j&][[1,1]]],
{i,Length@port},{j,Length@port}];

myline=myline[[#[[1]],#[[2]]]]&/@Partition[tour[[2]],2,1];

GeoGraphics[{
Red,PointSize[Large],
Point[GeoPosition[port[[All,{2,3}]]]],
Blue,Thick,myline},
GeoRange->"World",
GeoCenter->GeoPosition[{20,-150}],
GeoRange->Quantity[180,"AngularDegrees"]
]

繼續閱讀全文　Mathematica 教學：Searoute for Traveling Salesman Problem

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

Concise and elegant!

g=Graph[
{1\[DirectedEdge]2,1\[DirectedEdge]3,2\[DirectedEdge]3,
2\[DirectedEdge]4,3\[DirectedEdge]5,4\[DirectedEdge]3,
4\[DirectedEdge]5,4\[DirectedEdge]6,5\[DirectedEdge]6},
EdgeWeight->{1,12,9,3,5,3,13,15,4},
EdgeLabels->Rule@@@Transpose@{
{1\[DirectedEdge]2,1\[DirectedEdge]3,2\[DirectedEdge]3,
2\[DirectedEdge]4,3\[DirectedEdge]5,4\[DirectedEdge]3,
4\[DirectedEdge]5,4\[DirectedEdge]6,5\[DirectedEdge]6},
{1,12,9,3,5,3,13,15,4}},
VertexStyle->Yellow,VertexSize->0.5,
VertexShapeFunction->{{Yellow,Disk[#1,0.1],Black,Text[#2,#1]}&}]

route={{1,2,1},{1,3,12},{2,3,9},{2,4,3},{3,5,5},{4,3,4},{4,5,13},{5,6,4},{4,6,15}};

myroute=Select[Table[Select[route,#[[2]]==z&],{z,1,6}],#=!={}&]

Table[Set[v[temp], 
SortBy[If[#[[1]]==1,{v[temp],#[[3]]},{v[temp],v[#[[1]]]+#[[3]]}]&/@
Select[Flatten[myroute,1],#[[2]]==temp&],Last][[1,2]]],{temp,2,6}]&/@Range[2]

繼續閱讀全文　Mathematica 教學：Dijkstra Algorithm

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

A Sustainable Dynamic Optimization Model of Pricing and Advertising in the Presence of Green Innovation Investment

Chung-Yuan Dye and Tsu-Pang Hsieh

Finally, the acceptance from EJOR has arrived. The Julia code for the numerical examples in this paper is provided below. Alternatively, you can also implement it directly in Julia. Please leave a message if you have any problems with the implementation of the code.

julia=StartExternalSession["Julia"]

ExternalEvaluate[julia,"
using JuMP,Ipopt
dp(x)=(140-3*x);
K1(y)=(y/5-y^2/1000);
f1(x)=100*exp(-0.03*x);
f2(x)=2*exp(-0.045*x);
v1=25;
v2=10;
omega=300;
c=16.5;
h=0;
z=0.8;
"];

dynCP11=ExternalEvaluate[julia, "
	function dynCP11(G0,gamma1,gamma2,delta,k1coef,k2coef,e,r,varpi,n)
    m1 = Model(Ipopt.Optimizer) 
    set_silent(m1)
    #set_optimizer_attributes(m1, \"tol\" => 1e-4, \"max_iter\" => 100)
    register(m1, :dp, 1, dp; autodiff = true)
    register(m1, :K1, 1, K1; autodiff = true)
    register(m1, :f1, 1, f1; autodiff = true)
    register(m1, :f2, 1, f2; autodiff = true)
    @variable(m1, p[1:n], start = 30)
    @variable(m1, w[1:n], start = 10)
    @variable(m1, u[1:n], start = 10)

    @NLexpression(m1,newG[i=1:n],(1-delta)^(-1+i)*((1-delta)*G0
	+sum((gamma1*u[j+1]^k1coef+gamma2*w[j+1]^k2coef)/((1-delta)^j) 
    for j = 0:i-1)))
    
    @NLexpression(m1,newCE[i=1:n],f1(w[i])+f2(w[i])*dp(p[i])*K1(newG[i])-varpi)

    @NLobjective(
	m1,
	Max,
	sum(z^(i-1)*((p[i]-c-h/2)*dp(p[i])*K1(newG[i])
    -v1*u[i]^2/2-v2*w[i]^2/2-e*newCE[i])
                for i = 1:n))
	for i = 1:n
        @constraint(m1, p[i] >= c + h / 2)
	end
	for i = 1:n
        @constraint(m1, p[i] <= 140 / 3)
	end
	for i = 1:n
        @constraint(m1, w[i] >= 0)
	end
	for i = 1:n
        @constraint(m1, u[i] >= 0)
	end
	for i = 1:n
        @NLconstraint(m1, newCE[i] >= 0)
	end

	JuMP.optimize!(m1)

	p_val = JuMP.value.(p)
	w_val = JuMP.value.(w)
	u_val = JuMP.value.(u)

	G_val=(1-delta)*G0+gamma1*u_val[1]^k1coef+gamma2*w_val[1]^k2coef

	CE=f1(w_val[1])+f2(w_val[1])*dp(p_val[1])*K1(G_val)

	profit=(p_val[1]-c-h/2)*dp(p_val[1])*K1(G_val)
    -v1*u_val[1]^2/2-v2*w_val[1]^2/2-e*(CE-varpi)

	return[p_val,w_val,u_val,G_val,CE,profit,objective_value(m1)]

end"];

dynCP22=ExternalEvaluate[julia,"
	function dynCP22(G0,gamma1,gamma2,delta,k1coef,k2coef,e,r,varpi,n)
    m1 = Model(Ipopt.Optimizer) 
    set_silent(m1)
    #set_optimizer_attributes(m1, \"tol\" => 1e-4, \"max_iter\" => 100)
    register(m1, :dp, 1, dp; autodiff = true)
    register(m1, :K1, 1, K1; autodiff = true)
    register(m1, :f1, 1, f1; autodiff = true)
    register(m1, :f2, 1, f2; autodiff = true)
    @variable(m1, p[1:n], start = 30)
    @variable(m1, w[1:n], start = 10)
    @variable(m1, u[1:n], start = 10)

    @NLexpression(m1,newG[i=1:n],(1-delta)^(-1+i)*((1-delta)*G0
	+sum((gamma1*u[j+1]^k1coef+gamma2*w[j+1]^k2coef)/((1-delta)^j) 
    for j = 0:i-1)))
    
    @NLexpression(m1,newCE[i=1:n],f1(w[i])+f2(w[i])*dp(p[i])*K1(newG[i])-varpi)

    @NLobjective(
	m1,
	Max,
	sum(z^(i-1)*((p[i]-c-h/2)*dp(p[i])*K1(newG[i])
    -v1*u[i]^2/2-v2*w[i]^2/2-r*newCE[i])
                for i = 1:n))
	for i = 1:n
        @constraint(m1, p[i] >= c + h / 2)
	end
	for i = 1:n
        @constraint(m1, p[i] <= 140 / 3)
	end
	for i = 1:n
        @constraint(m1, w[i] >= 0)
	end
	for i = 1:n
        @constraint(m1, u[i] >= 0)
	end
	for i = 1:n
        @NLconstraint(m1, newCE[i] <= 0)
	end

	JuMP.optimize!(m1)

	p_val = JuMP.value.(p)
	w_val = JuMP.value.(w)
	u_val = JuMP.value.(u)

	G_val=(1-delta)*G0+gamma1*u_val[1]^k1coef+gamma2*w_val[1]^k2coef

	CE=f1(w_val[1])+f2(w_val[1])*dp(p_val[1])*K1(G_val)

	profit=(p_val[1]-c-h/2)*dp(p_val[1])*K1(G_val)
    -v1*u_val[1]^2/2-v2*w_val[1]^2/2-r*(CE-varpi)

	return[p_val,w_val,u_val,G_val,CE,profit,objective_value(m1)]

end"];

mysol13[g0_,\[Gamma]1_,\[Gamma]2_,\[Delta]_,k1_,k2_,e_,r_,cap_]:=
Block[{n=40,numiter=100,iter=0,ans1,ans2},
	Print[{Dynamic[iter],Dynamic[ans2]}];
	NestList[(
	iter=iter+1;
	ans1=
	SortBy[{
		dynCP11[#[[4]],\[Gamma]1,\[Gamma]2,\[Delta],k1,k2,e,r,cap,n],
		dynCP22[#[[4]],\[Gamma]1,\[Gamma]2,\[Delta],k1,k2,e,r,cap,n]},
          Last][[-1]];
	ans2={ans1[[1,1]],ans1[[2,1]],ans1[[3,1]],ans1[[4]],
          ans1[[5]],#[[-1]]+0.8^(iter-1)*ans1[[6]]})&,
	{0,0,0,g0,0,0},numiter][[2;;-1]]
]

myexample=mysol13[#[[1]],0.5,0.8,0.2,0.5,0.8,8,2,#[[2]]]&/@
Flatten[Table[{i,j},{j,{250,300,400}},{i,{10,40}}],1]

ListLinePlot[myexample[[All,All,#]],PlotRange->All,
Frame->True]&/@{1,2,3,4,5}//Column

繼續閱讀全文　Finally, the first acceptance of 2024 comes.

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

這是我這幾年來數量方法這一門課所自編的 Mathematica 入門教學講義。Mathematica 相較於其他軟體的特點在於強大的符號運算及動態幾何繪圖的科學運算能力。此外，Mathematica 尚包含數值運算、圖表可視化、動態操作介面等。由於擁有許多的內建程式語言及內建系統，故不需要繁雜且複雜的編寫，有能夠輕易上手的便利性。本教學講義相關檔案、資料檔及未來更新維護皆可由下列網址下載： 

https://github.com/chungyuandye/Mathematica

繼續閱讀全文　Mathematica 教學講義

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有1 意見