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
