# It's running too long

I ran the following code for 7 hours without any results, how can I modify it? Thanks

``````
using Optim, LinearAlgebra, Statistics, Plots, QuantEcon, Distributions, Interpolations, NLsolve, Random, JuMP, GLPK
import Ipopt

β = 0.88
θ = 3.4
Nb = 1000
Nz = 50
v = 0.19
α = 0.33
δ = 1.5
χ = 0.05
ϵ = 0.335
γ = 0.89
ψ = 0.95
ξ = 0.8
grid_z = range(1, 10, length = 50)
grid_b = range(0, 5, length = 1000)
grid_b = collect(grid_b)
grid_z = collect(grid_z)
dist_z = Pareto(θ)
dist_b = Pareto(θ)
wage = 1
rate = 1
out = zeros(50,1000)
action = zeros(50,1000)
Tw1 = zeros(1000,1)
Tw2 = zeros(1000,1)
Tw3 = zeros(1000,1)
c_choice1 = zeros(1000,1)
k_choice1 = zeros(1000,1)
l_choice1 = zeros(1000,1)
c_choice2 = zeros(1000,1)
k_choice2 = zeros(1000,1)
l_choice2 = zeros(1000,1)
c_choice3 = zeros(1000,1)
k_choice3 = zeros(1000,1)
l_choice3 = zeros(1000,1)
#################################### bellman begin #################################

for i in 1:50

################################### choice1 ##########################
function T1(w)

# 插值函数t1
function t1(x)
w_func = LinearInterpolation(grid_b, w,extrapolation_bc=Line())
return w_func(x)
end

# 优化模型choice1
choice1 = Model(Ipopt.Optimizer)
@variable(choice1, x >= 0)
register(choice1,:w_func,1,t1; autodiff = true)
@NLparameter(choice1, b == 0)
@NLobjective(choice1, Max, (x^(1-δ))/(1-δ) + β * (w_func((b*(1+rate)+wage-x))))
for j in 1:Nb
set_value(b, grid_b[j])
optimize!(choice1)
Tw1[j,1] = objective_value(choice1)
c_choice1[j,1] = value(x)
end

return Tw1
end

################################### choice2 ##########################
function T2(w)

# 插值函数t2
function t2(x)
w_func2 = LinearInterpolation(grid_b, w, extrapolation_bc=Line())
return w_func2(x)
end

#优化模型choice2
choice2 = Model(Ipopt.Optimizer)
@variable(choice2, x >= 0)
@variable(choice2, y >= 0)
@variable(choice2, z >= 0)
register(choice2, :w_func2, 1, t2; autodiff = true)
@NLparameter(choice2, b == 1)
@NLobjective(choice2, Max, (x^(1-δ))/(1-δ) + β *  (w_func2(-x + (1+rate)*(b-y) + (1-δ)*y + grid_z[i]*((y^α)*(z^(1-α))^(1-v))-wage*z)))
@NLconstraint(choice2, y <= b)

for j in 1:Nb
set_value(b, grid_b[j])
optimize!(choice2)
Tw2[j,1] = objective_value(choice2)
c_choice2[j,1] = value(x)
k_choice2[j,1] = value(y)
l_choice2[j,1] = value(z)
end

return Tw2
end

################################### choice3 ##########################
function T3(w)

function t3(x)
w_func3 = LinearInterpolation(grid_b, w,extrapolation_bc=Line())
return w_func3(x)
end

choice3 = Model(Ipopt.Optimizer)
@variable(choice3, x >= 0)
@variable(choice3, y >= 0)
@variable(choice3, z >= 0)
register(choice3, :w_func3, 1, t3; autodiff = true)
@NLparameter(choice3, b == 1)
@NLobjective(choice3, Max, (x^(1-δ))/(1-δ) + β * (w_func3(-x - wage*z + (1-δ)*y - (1+rate+χ)*(y-b+ψ) +grid_z[i]*((y^α)*(z^(1-α))^(1-v)))))
@NLconstraint(choice3, ξ*y <= b - ψ)
@NLconstraint(choice3, y >= b)

for j in 1:Nb
set_value(b, grid_b[j])
optimize!(choice3)
Tw3[j,1] = objective_value(choice3) # 返回生产率为zi，初始财富为bj下的最优值函数
c_choice3[j,1] = value(x)
k_choice3[j,1] = value(y)
l_choice3[j,1] = value(z)
end

return Tw3
end

################################### solve ##########################
ww1 = grid_b
ww2 = grid_b
ww3 = grid_b
for mm in 1:500
Tw1 = T1(ww1)
ww1 = [Tw1[j,1] for j in 1:1000]
Tw2 = T2(ww1)
ww2 = [Tw2[j,1] for j in 1:1000]
Tw3 = T3(ww1)
ww4 = [Tw3[j,1] for j in 1:1000]
end

func1 = LinearInterpolation(grid_b, ww1, extrapolation_bc=Line())
func2 = LinearInterpolation(grid_b, ww2, extrapolation_bc=Line())
func3 = LinearInterpolation(grid_b, ww3, extrapolation_bc=Line())

for j in 1:1000

b = grid_b[j]
v1 = func1(b)
v2 = func2(b)
v3 = func3(b)

if v1 > max(v2, v3)
action = 1
elseif v2 > max(v1, v3)
action = 2
else
action = 3
end

action[i,j] = action
out[i,j] = max(v1, v2, v3)
end
end

return action, out

#################################### Bellman end #################################

``````

There are a number of performance issues here, first among them your use of global variables - see the Performance Tips section of the manual. You’ll need to either wrap the whole thing in a function (including your parameter definitions in the beginning), or pass the parameters as arguments to your functions. Parameters.jl helps with parameter-handling

When trying to tune performance, it’s best to trim down the problem as much as possible so you can iterate quickly - for example, replace `for i in 1:50` with `for i in 1:1`, and reduce the length of your grids drastically. It helps to set the grid lengths using a variable so it only needs to be modified in one place, e.g. `N = 1000; c_choice1 = zeros(N), k_choice1 = zeros(N)`, which can be quickly changed to `N = 10`. Once you have a functional ‘diet’ version of your code, you can use Julia’s profiling tools (mentioned in the Performance Tips) to find the slow spots. Pay special attention to memory allocation and type instability.

Try to avoid writing functions inside the inner loops of your code - this is sometimes necessary, but it’s often more ergonomic to move the function to the outermost level and pass any necessary data as function arguments.

You use this pattern several times:

``````function t1(x)
w_func = LinearInterpolation(grid_b, w,extrapolation_bc=Line())
return w_func(x)
end
``````

Are you trying to create an interpolation and evaluate it at one point, or create a function that can be passed to the optimizer? In either case, you should only need to create the interpolator once (not once per function call!), e.g.

``````w_func = LinearInterpolation(grid_b, w,extrapolation_bc=Line())
t1 = x -> w_func(x)
``````

You may also be able to use `Interpolations.gradient` to calculate the derivative, which will likely be faster than `:autodiff`.

6 Likes

thanks！But I need to go round the w in the linearInterpolation, how do I set it up? I’ve been reporting errors by running the following programs. The error is “Incorrect number of arguments for “tt” in nonlinear expression.”

``````grid_b =range(1,10,length = 100)

function tt(x,wx,wy)
w_func = LinearInterpolation(wx, wy,extrapolation_bc=Line())
return w_func(x)
end
choice1 = Model(Ipopt.Optimizer)
@variable(choice1, x >= 0)
register(choice1,:tt,1,tt; autodiff = true)
@NLparameter(choice1, b == 0)
@NLobjective(choice1, Max, (x^0.5)/0.5 + 0.5 * tt((0.3-x),1,1))
optimize!(choice1)
``````

The suggestion is something like

``````grid_b =range(1,10,length = 100)
w_func = LinearInterpolation(grib_b, w, extrapolation_bc=Line())
choice1 = Model(Ipopt.Optimizer)
@variable(choice1, x >= 0)
foo(x) = w_func(x)
register(choice1, :foo, 1, foo; autodiff = true)
@NLobjective(choice1, Max, (x^0.5)/0.5 + 0.5 * foo(0.3-x))
optimize!(choice1)
``````

Thank you for your answer. When I adjust the code, I can get the result quickly when I set Nb, Nz and Nm to all 2. But running for a long time at a larger Nb， Nz and Nz doesn’t work out. it says that“Excessive output truncated after 524812 bytes.”

``````
using Optim, LinearAlgebra, Statistics, Plots, QuantEcon, Distributions, Interpolations, NLsolve, Random, JuMP, GLPK
import Ipopt
function CareerWorkerProblem()
β = 0.88
θ = 3.4
Nb = 100
Nz = 4
Nm = 4
v = 0.19
α = 0.33
δ = 1.5
χ = 0.05
ϵ = 0.335
γ = 0.89
ψ = 0.95
ξ = 0.8
grid_z = collect(range(1, 10, length = Nz))
grid_b = collect(range(0, 5, length = Nb))
dist_z = Pareto(θ)
dist_b = Pareto(θ)
wage = 1
rate = 1
out = zeros(Nz,Nb)
action = zeros(Nz,Nb)
Tw1 = zeros(Nb,1)
Tw2 = zeros(Nb,1)
Tw3 = zeros(Nb,1)
c_choice1 = zeros(Nb,1)
k_choice1 = zeros(Nb,1)
l_choice1 = zeros(Nb,1)
c_choice2 = zeros(Nb,1)
k_choice2 = zeros(Nb,1)
l_choice2 = zeros(Nb,1)
c_choice3 = zeros(Nb,1)
k_choice3 = zeros(Nb,1)
l_choice3 = zeros(Nb,1)
return ( β = β, θ = θ, Nb = Nb,Nz = Nz, v = v, α = α, δ = δ, χ = χ, ϵ = ϵ, γ = γ, ψ = ψ, ξ = ξ,
grid_z = grid_z, grid_b = grid_b, dist_z = dist_z, dist_b = dist_b, wage = wage, rate = rate,
out = out, action = action, Tw1 = Tw1, Tw2 = Tw2, Tw3 = Tw3, Nm = Nm,
c_choice1 = c_choice3, k_choice1 = k_choice3, l_choice1 = l_choice3, c_choice2 = c_choice3,
k_choice2 = k_choice3, l_choice2 = l_choice3, c_choice3 = c_choice3, k_choice3 = k_choice3,
l_choice3 = l_choice3)
end

####################### Linear Interpolation ##########################

function tt(x, lx, ly)
w_func = LinearInterpolation(lx, ly,extrapolation_bc=Line())
return w_func(x)
end

################################### choice1 ##########################
function T1(cp, ty, Tw1)

# 优化模型choice1
choice1 = Model(Ipopt.Optimizer)
@variable(choice1, x >= 0)
foo(xx)=tt(xx,cp.grid_b,ty)
register(choice1,:foo,1,foo; autodiff = true)
@NLparameter(choice1, b == 0)
@NLobjective(choice1, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β * (foo(b*(1+cp.rate)+cp.wage-x)))

for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice1)
Tw1[j,1] = objective_value(choice1)
end

return Tw1
end

################################### choice2 ##########################
function T2(cp, ty, i, Tw2)

#优化模型choice2
choice2 = Model(Ipopt.Optimizer)
@variable(choice2, x >= 0)
@variable(choice2, y >= 0)
@variable(choice2, z >= 0)
foo(xx) = tt(xx,cp.grid_b,ty)
register(choice2, :foo, 1, foo; autodiff = true)
@NLparameter(choice2, b == 1)
@NLobjective(choice2, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β *  (foo(-x + (1+cp.rate)*(b-y) + (1-cp.δ)*y + cp.grid_z[i]*((y^cp.α)*(z^(1-cp.α))^(1-cp.v))-cp.wage*z)))
@NLconstraint(choice2, y <= b)

for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice2)
Tw2[j,1] = objective_value(choice2)
end

return Tw2
end

################################### choice3 ##########################
function T3(cp, ty, i, Tw3)

choice3 = Model(Ipopt.Optimizer)
@variable(choice3, x >= 0)
@variable(choice3, y >= 0)
@variable(choice3, z >= 0)
foo(xx) = tt(xx,cp.grid_b,ty)
register(choice3, :foo, 1, foo; autodiff = true)
@NLparameter(choice3, b == 1)
@NLobjective(choice3, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β * (foo(-x - cp.wage*z + (1-cp.δ)*y - (1+cp.rate+cp.χ)*(y-b+cp.ψ) +cp.grid_z[i]*((y^cp.α)*(z^(1-cp.α))^(1-cp.v)))))
@NLconstraint(choice3, cp.ξ*y <= b - cp.ψ)
@NLconstraint(choice3, y >= b)

for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice3)
Tw3[j,1] = objective_value(choice3) # 返回生产率为zi，初始财富为bj下的最优值函数
end

return Tw3
end

#################################### bellman begin #################################

function  bellman(cp)
action = cp.action
out = cp.out
for i in 1:cp.Nz

ww1 = cp.grid_b
ww2 = cp.grid_b
ww3 = cp.grid_b

for mm in 1:cp.Nm
@time Tw1 = T1(cp, ww1, cp.Tw1)
ww1 = [Tw1[j,1] for j in 1:cp.Nb]
@time Tw2 = T2(cp, ww2, mm, cp.Tw2)
ww2 = [Tw2[j,1] for j in 1:cp.Nb]
@time Tw3 = T3(cp, ww3, mm, cp.Tw3)
ww3 = [Tw3[j,1] for j in 1:cp.Nb]
end

@time func1 = LinearInterpolation(cp.grid_b, ww1, extrapolation_bc=Line())
func2 = LinearInterpolation(cp.grid_b, ww2, extrapolation_bc=Line())
func3 = LinearInterpolation(cp.grid_b, ww3, extrapolation_bc=Line())

for j in 1:cp.Nb

v1 = func1(cp.grid_b[j])
v2 = func2(cp.grid_b[j])
v3 = func3(cp.grid_b[j])

if v1 > max(v2, v3)
action[i,j] = 1
elseif v2 > max(v1, v3)
action[i,j] = 2
else
action[i,j] = 3
end

out[i,j] = max(v1, v2, v3)
end
end
return action, out
end
#################################### Bellman end #################################
wp = CareerWorkerProblem()
@time action, out = bellman(wp)
``````

Number of nonzeros in inequality constraint Jacobian.: 0.347385 seconds (67.18 k allocations: 11.137 MiB)。
thanks

It’s helpful if you can create a minimal working example that demonstrates the problem. Simplify your code as much as possible and remove any unneeded packages.

It’s also helpful if you can provide the full text of the error message.

Note that the code below is still incorrect. Why? Because every call to `foo(xx)` calls `tt(xx, cp.grib_b, ty)`, which creates a new `LinearInterpolation` object.

``````function tt(x, lx, ly)
w_func = LinearInterpolation(lx, ly,extrapolation_bc=Line())
return w_func(x)
end

function T1(cp, ty, Tw1)
choice1 = Model(Ipopt.Optimizer)
@variable(choice1, x >= 0)
foo(xx)=tt(xx,cp.grid_b,ty)
register(choice1,:foo,1,foo; autodiff = true)
``````

``````function T1(cp, ty, Tw1)
choice1 = Model(Ipopt.Optimizer)
@variable(choice1, x >= 0)
w_func = LinearInterpolation(cp.grid_b, ty, extrapolation_bc=Line())
foo(x) = w_func(x)
register(choice1, :foo, 1, foo; autodiff = true)
``````

Thanks! I have change it as belows. I need to make Nb =1000, Nz = 50, Nm =500, when I run I can’t run out of results. The following is Nb=2, Nz=2,Nm=1,I can get results instantly

``````using LinearAlgebra, Statistics, Distributions, Interpolations, JuMP, Ipopt
function CareerWorkerProblem()
β = 0.88
θ = 3.4
Nb = 2
Nz = 2
Nm = 1
v = 0.19
α = 0.33
δ = 1.5
χ = 0.05
ϵ = 0.335
γ = 0.89
ψ = 0.95
ξ = 0.8
grid_z = collect(range(1, 10, length = Nz))
grid_b = collect(range(0, 5, length = Nb))
dist_z = Pareto(θ)
dist_b = Pareto(θ)
wage = 1
rate = 1
out = zeros(Nz,Nb)
action = zeros(Nz,Nb)
Tw1 = zeros(Nb,1)
Tw2 = zeros(Nb,1)
Tw3 = zeros(Nb,1)
c_choice1 = zeros(Nb,1)
k_choice1 = zeros(Nb,1)
l_choice1 = zeros(Nb,1)
c_choice2 = zeros(Nb,1)
k_choice2 = zeros(Nb,1)
l_choice2 = zeros(Nb,1)
c_choice3 = zeros(Nb,1)
k_choice3 = zeros(Nb,1)
l_choice3 = zeros(Nb,1)
return ( β = β, θ = θ, Nb = Nb,Nz = Nz, v = v, α = α, δ = δ, χ = χ, ϵ = ϵ, γ = γ, ψ = ψ, ξ = ξ,
grid_z = grid_z, grid_b = grid_b, dist_z = dist_z, dist_b = dist_b, wage = wage, rate = rate,
out = out, action = action, Tw1 = Tw1, Tw2 = Tw2, Tw3 = Tw3, Nm = Nm,
c_choice1 = c_choice3, k_choice1 = k_choice3, l_choice1 = l_choice3, c_choice2 = c_choice3,
k_choice2 = k_choice3, l_choice2 = l_choice3, c_choice3 = c_choice3, k_choice3 = k_choice3,
l_choice3 = l_choice3)
end

################################### choice1 ##########################
function T1(cp, ty, Tw1)

choice1 = Model(Ipopt.Optimizer)
set_silent(choice1 ::Model)
@variable(choice1, x >= 0)
w_func = LinearInterpolation(cp.grid_b, ty, extrapolation_bc=Line())
foo(x) = w_func(x)
register(choice1,:foo,1,foo; autodiff = true)
@NLparameter(choice1, b == 0)
@NLobjective(choice1, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β * (foo(b*(1+cp.rate)+cp.wage-x)))

for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice1)
Tw1[j,1] = objective_value(choice1)
end

return Tw1
end

################################### choice2 ##########################
function T2(cp, ty, i, Tw2)

choice2 = Model(Ipopt.Optimizer)
set_silent(choice2 ::Model)
@variable(choice2, x >= 0)
@variable(choice2, y >= 0)
@variable(choice2, z >= 0)
w_func = LinearInterpolation(cp.grid_b, ty, extrapolation_bc=Line())
foo(x) = w_func(x)
register(choice2, :foo, 1, foo; autodiff = true)
@NLparameter(choice2, b == 1)
@NLobjective(choice2, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β *  (foo(-x + (1+cp.rate)*(b-y) + (1-cp.δ)*y + cp.grid_z[i]*((y^cp.α)*(z^(1-cp.α))^(1-cp.v))-cp.wage*z)))
@NLconstraint(choice2, y <= b)

for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice2)
Tw2[j,1] = objective_value(choice2)
end

return Tw2
end

################################### choice3 ##########################
function T3(cp, ty, i, Tw3)

choice3 = Model(Ipopt.Optimizer)
set_silent(choice3 ::Model)
@variable(choice3, x >= 0)
@variable(choice3, y >= 0)
@variable(choice3, z >= 0)
w_func = LinearInterpolation(cp.grid_b, ty, extrapolation_bc=Line())
foo(x) = w_func(x)
register(choice3, :foo, 1, foo; autodiff = true)
@NLparameter(choice3, b == 1)
@NLobjective(choice3, Max, (x^(1-cp.δ))/(1-cp.δ) + cp.β * (foo(-x - cp.wage*z + (1-cp.δ)*y - (1+cp.rate+cp.χ)*(y-b+cp.ψ) +cp.grid_z[i]*((y^cp.α)*(z^(1-cp.α))^(1-cp.v)))))
@NLconstraint(choice3, cp.ξ*y <= b - cp.ψ)
@NLconstraint(choice3, y >= b)

for j in 1:cp.Nb
set_value(b, cp.grid_b[j])
optimize!(choice3)
Tw3[j,1] = objective_value(choice3)
end

return Tw3
end

#################################### bellman begin #################################

function  bellman(cp)
action = cp.action
out = cp.out
for i in 1:cp.Nz

ww1 = cp.grid_b
ww2 = cp.grid_b
ww3 = cp.grid_b

for mm in 1:cp.Nm
Tw1 = T1(cp, ww1, cp.Tw1)
ww1 = [Tw1[j,1] for j in 1:cp.Nb]
Tw2 = T2(cp, ww2, i, cp.Tw2)
ww2 = [Tw2[j,1] for j in 1:cp.Nb]
Tw3 = T3(cp, ww3, i, cp.Tw3)
ww3 = [Tw3[j,1] for j in 1:cp.Nb]
end

w_func1=LinearInterpolation(cp.grid_b, ww1, extrapolation_bc=Line())
w_func2=LinearInterpolation(cp.grid_b, ww2, extrapolation_bc=Line())
w_func3=LinearInterpolation(cp.grid_b, ww3, extrapolation_bc=Line())

for j in 1:cp.Nb

v1 = w_func1(cp.grid_b[j])
v2 = w_func2(cp.grid_b[j])
v3 = w_func3(cp.grid_b[j])

if v1 > max(v2, v3)
action[i,j] = 1
elseif v2 > max(v1, v3)
action[i,j] = 2
else
action[i,j] = 3
end

out[i,j] = max(v1, v2, v3)
end
end
return action, out
end
#################################### Bellman end #################################
wp = CareerWorkerProblem()
action, out = bellman(wp)
``````

i have also run profile, results as below

``````Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
╎ 90  @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
╎  90  @Base/essentials.jl:706; invokelatest
╎   90  @Base/essentials.jl:708; #invokelatest#2
╎    90  .../execute_request.jl:67; execute_request(socket::ZMQ.Soc...
╎     90  .../SoftGlobalScope.jl:65; softscope_include_string(m::Mo...
1╎    ╎  90  @Base/boot.jl:360; eval
╎    ╎   16  In[1]:127; bellman(cp::NamedTuple{(:β, ...
╎    ╎    1   In[1]:53; T1(cp::NamedTuple{(:β, :θ, :...
╎    ╎     1   @JuMP/src/nlp.jl:1185; register##kw
╎    ╎    ╎ 1   @JuMP/src/nlp.jl:1191; register(m::Model, s::Sym...
╎    ╎    ╎  1   @Base/Base.jl:33; getproperty(x::JuMP._NLPDa...
╎    ╎    15  In[1]:59; T1(cp::NamedTuple{(:β, :θ, :...
╎    ╎     15  ...izer_interface.jl:106; optimize!
╎    ╎    ╎ 15  ...zer_interface.jl:106; optimize!
╎    ╎    ╎  15  ...zer_interface.jl:130; optimize!(model::Model, o...
╎    ╎    ╎   1   ...ingoptimizer.jl:248; optimize!(m::MathOptInter...
╎    ╎    ╎    1   ...ingoptimizer.jl:185; attach_optimizer(model::...
╎    ╎    ╎     1   ...ge_optimizer.jl:401; (::MathOptInterface.var"...
╎    ╎    ╎    ╎ 1   ...e_optimizer.jl:401; #copy_to#4
╎    ╎    ╎    ╎  1   ...lities/copy.jl:23; automatic_copy_to##kw
╎    ╎    ╎    ╎   1   ...lities/copy.jl:24; #automatic_copy_to#127
╎    ╎    ╎    ╎    1   ...ities/copy.jl:714; default_copy_to(dest:...
╎    ╎    ╎    ╎     1   ...ities/copy.jl:497; (::MathOptInterface.U...
╎    ╎    ╎    ╎    ╎ 1   ...ties/copy.jl:497; pass_constraints##kw
1╎    ╎    ╎    ╎    ╎  1   ...ties/copy.jl:509; pass_constraints(de...
╎    ╎    ╎   14  ...ingoptimizer.jl:252; optimize!(m::MathOptInter...
╎    ╎    ╎    14  ...ge_optimizer.jl:319; optimize!(b::MathOptInte...
╎    ╎    ╎     1   ...MOI_wrapper.jl:1282; optimize!(model::Ipopt....
1╎    ╎    ╎    ╎ 1   ...MOI_wrapper.jl:999; hessian_lagrangian_stru...
╎    ╎    ╎     13  ...MOI_wrapper.jl:1441; optimize!(model::Ipopt....
13╎    ╎    ╎    ╎ 13  ...t/src/Ipopt.jl:513; solveProblem(prob::Ipop...
╎    ╎   16  In[1]:129; bellman(cp::NamedTuple{(:β, ...
╎    ╎    1   In[1]:71; T2(cp::NamedTuple{(:β, :θ, :...
╎    ╎     1   @JuMP/src/macros.jl:91; macro expansion
╎    ╎    ╎ 1   ...src/variables.jl:801; add_variable(model::Model,...
╎    ╎    ╎  1   ...src/variables.jl:807; _moi_add_variable(backend...
1╎    ╎    ╎   1   ...rc/variables.jl:814; _moi_constrain_variable(b...
╎    ╎    15  In[1]:83; T2(cp::NamedTuple{(:β, :θ, :...
╎    ╎     15  ...izer_interface.jl:106; optimize!
╎    ╎    ╎ 15  ...zer_interface.jl:106; optimize!
╎    ╎    ╎  15  ...zer_interface.jl:130; optimize!(model::Model, o...
╎    ╎    ╎   15  ...ingoptimizer.jl:252; optimize!(m::MathOptInter...
╎    ╎    ╎    15  ...ge_optimizer.jl:319; optimize!(b::MathOptInte...
╎    ╎    ╎     2   ...MOI_wrapper.jl:1344; optimize!(model::Ipopt....
2╎    ╎    ╎    ╎ 2   ...t/src/Ipopt.jl:350; createProblem(n::Int64,...
╎    ╎    ╎     13  ...MOI_wrapper.jl:1441; optimize!(model::Ipopt....
13╎    ╎    ╎    ╎ 13  ...t/src/Ipopt.jl:513; solveProblem(prob::Ipop...
╎    ╎   57  In[1]:131; bellman(cp::NamedTuple{(:β, ...
╎    ╎    57  In[1]:108; T3(cp::NamedTuple{(:β, :θ, ...
╎    ╎     57  ...izer_interface.jl:106; optimize!
╎    ╎    ╎ 57  ...zer_interface.jl:106; optimize!
╎    ╎    ╎  57  ...zer_interface.jl:130; optimize!(model::Model, o...
╎    ╎    ╎   57  ...ingoptimizer.jl:252; optimize!(m::MathOptInter...
╎    ╎    ╎    57  ...ge_optimizer.jl:319; optimize!(b::MathOptInte...
╎    ╎    ╎     1   ...MOI_wrapper.jl:1280; optimize!(model::Ipopt....
╎    ╎    ╎    ╎ 1   @JuMP/src/nlp.jl:393; initialize(d::NLPEvalua...
╎    ╎    ╎    ╎  1   @JuMP/src/nlp.jl:280; JuMP._FunctionStorage(...
╎    ╎    ╎    ╎   1   ...s/coloring.jl:454; hessian_color_preproce...
╎    ╎    ╎    ╎    1   ...s/coloring.jl:183; acyclic_coloring(g::J...
╎    ╎    ╎    ╎     1   ...sjoint_set.jl:35; DataStructures.IntDis...
1╎    ╎    ╎    ╎    ╎ 1   ...sjoint_set.jl:27; IntDisjointSets
╎    ╎    ╎     2   ...MOI_wrapper.jl:1344; optimize!(model::Ipopt....
2╎    ╎    ╎    ╎ 2   ...t/src/Ipopt.jl:350; createProblem(n::Int64,...
╎    ╎    ╎     54  ...MOI_wrapper.jl:1441; optimize!(model::Ipopt....
50╎    ╎    ╎    ╎ 54  ...t/src/Ipopt.jl:513; solveProblem(prob::Ipop...
╎    ╎    ╎    ╎  3   ...t/src/Ipopt.jl:202; eval_g_wrapper(n::Int3...
╎    ╎    ╎    ╎   3   ...OI_wrapper.jl:1305; (::Ipopt.var"#eval_g_...
╎    ╎    ╎    ╎    3   ...OI_wrapper.jl:1113; eval_constraint(model...
╎    ╎    ╎    ╎     3   @JuMP/src/nlp.jl:569; eval_constraint(d::NL...
╎    ╎    ╎    ╎    ╎ 1   @Base/timing.jl:286; macro expansion
1╎    ╎    ╎    ╎    ╎  1   @Base/Base.jl:65; time_ns
╎    ╎    ╎    ╎    ╎ 2   @Base/timing.jl:287; macro expansion
╎    ╎    ╎    ╎    ╎  2   ...P/src/nlp.jl:571; macro expansion
╎    ╎    ╎    ╎    ╎   2   .../src/nlp.jl:498; _forward_eval_all(d...
╎    ╎    ╎    ╎    ╎    1   .../forward.jl:125; forward_eval(stora...
1╎    ╎    ╎    ╎    ╎     1   .../NaNMath.jl:23; pow
╎    ╎    ╎    ╎    ╎    1   .../forward.jl:186; forward_eval(stora...
1╎    ╎    ╎    ╎    ╎     1   .../forward.jl:5; eval_and_check_retu...
╎    ╎    ╎    ╎  1   ...t/src/Ipopt.jl:238; eval_jac_g_wrapper(n::...
╎    ╎    ╎    ╎   1   ...OI_wrapper.jl:1315; (::Ipopt.var"#eval_ja...
1╎    ╎    ╎    ╎    1   ...OI_wrapper.jl:1175; eval_constraint_jacob...
Total snapshots: 167
``````

thanks！