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 意見

session=StartExternalSession["Julia"]
test =
  "function testpi(size)
    sample = rand(size, 2)
    dist = broadcast(hypot, sample[:,1], sample[:,2])
    return sample, dist
end";
pitest=ExternalEvaluate[session,test]
size=100000
data=pitest[size];
ListPlot[GatherBy[Transpose[data[[1]]],Norm[#]<1&],Frame ->True,AspectRatio -> 1]
Length@GatherBy[Transpose[data[[1]]], Norm[#]<1&][[1]]/size//N




\documentclass[tikz]{standalone}
\usepackage{tikz}
\usepackage{ifthen}
\begin{document}
  \begin{tikzpicture}
	\draw[help lines, color=gray!30, dashed] (-1.2,1.5) grid (-1.2,1.5);
	\draw[->] (-1.2,0)--(1.1,0) node[right]{$x$};
	\draw[->] (0,-1.2)--(0,1.5) node[above]{$y$};
	\draw[domain=0:1, smooth, variable=\x, black] plot ({\x}, {sqrt(1-\x*\x)});
  \pgfmathsetseed{42}
    \foreach \i in {1,...,10000}
    {
      \pgfmathsetmacro{\x}{abs(rand)}
      \pgfmathsetmacro{\y}{abs(rand)}
      \pgfmathsetmacro{\r}{\x*\x+\y*\y}
       \pgfmathsetmacro{\col}{ifthenelse(\r<1,"blue!90","red!90")}
      \fill[\col] (\x,\y) circle (0.1pt)
      ;
    }
  \end{tikzpicture}
\end{document}

繼續閱讀全文　Mathematica 教學：Call Julia from Mathematica

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