Hi guys, in last questions (link is shown below), I asked how to supply gradient of objective and jacobian of constraints in nonlinear optimization. But @odow told me I dont have to do that because JuMP.jl will automatically calculate them and feed them to Ipopt.jl. However, my model still runs slowly. So I show my code to find any opportunities to speed up.

Thanks for your help.

```
# ******************************************************** #
# MATHEMATICAL PROGRAM WITH EQUILIBRIUM CONSTRAINTS #
# ******************************************************** #
using Base.Threads
using Random, Distributions, Statistics
using LinearAlgebra, LogExpFunctions
using JuMP, Ipopt, ForwardDiff
const seed = Random.seed!(2023 - 11 - 24)
# ******************************************************** #
# DATA GENERATION #
# ******************************************************** #
Ps = 25 # number of products
Mkts = 40 # number of markets
Inds = 100 # number of individuals in each market
vcovXc = [
[1, -0.8, 0.3] [-0.8, 1, 0.3] [0.3, 0.3, 1]
]
Eπ½ = [-1.0, 1.5, 1.5, 0.5, -3.0] # π1
ππ½ = sqrt.([0.5, 0.5, 0.5, 0.5, 0.2]) # π2
π½i = rand(seed, MvNormal(ππ½), Inds)
struct Market
J::Int64 # number of products
π::AbstractVector{<:Real} # price vector, [J, 1]
Xc::AbstractMatrix{<:Real} # observed product characteristics, [J, B]
share::AbstractVector{<:Real} # market share exluding outshare, [J, 1]
out_share::Float64
v::AbstractMatrix{<:Real} # simulated individuals, [B, Inds]
iv::AbstractMatrix{<:Real} # instrumental variables
end
d = Vector{Market}(undef, Mkts);
Xc = rand(seed, MvNormal(vcovXc), Ps)'
let
for mkt β eachindex(d)
π, e = randn(seed, Ps), randn(seed, Ps)
price = abs.(0.5 * π + e + 1.1 * sum(Xc; dims = 2)[:])
Z = rand(seed, Ps, 6) .+ 0.25 * (e + 1.1 * sum(Xc; dims = 2)[:])
IV = hcat(ones(Ps), Z, Xc, Z .^ 2, Z .^ 3, Xc .^ 2, Xc .^ 3, sum(Z; dims = 2), sum(Xc; dims = 2), kron(Xc[:, 1], ones(size(Z, 1))') * Z, kron(Xc[:, 2], ones(size(Z, 1))') * Z)
X = hcat(ones(Ps), Xc, price)
U = X * Eπ½ .+ X * π½i
denom = 1 .+ sum(exp.(U); dims = 1)
pchoice = exp.(U) ./ denom
share = mean(pchoice; dims = 2)[:]
out_share = 1 - sum(share)
v = rand(seed, MvNormal(ones(size(X, 2))), Inds)
d[mkt] = Market(Ps, price, Xc, share, out_share, v, IV)
end
end
# ======================================================== #
# ESTIMATION #
# ======================================================== #
model = Model(Ipopt.Optimizer)
@variable(model, theta2[1:length(ππ½)], start = 1.0)
@variable(model, theta1[1:length(Eπ½)])
π = Vector{Vector{VariableRef}}(undef, Mkts)
for mkt in 1:Mkts
π[mkt] = @variable(model, [1:Ps], base_name = "π$(mkt)")
end
pXc = Vector{Matrix{Float64}}(undef, Mkts)
for m in 1:Mkts
pXc[m] = hcat(ones(d[m].J), d[m].Xc, d[m].π)
end
π = @expression(model, [m = 1:Mkts], pXc[m] * Diagonal(theta2) * d[m].v);
πΏ = @expression(model, [m = 1:Mkts], pXc[m] * theta1 .+ π[m]);
π = @expression(model, [m = 1:Mkts], πΏ[m] .+ π[m]);
expπ = @expression(model, [m = 1:Mkts], exp.(π[m]));
denom = @expression(model, [m = 1:Mkts], mapreduce(x -> 1 + sum(x), hcat, eachcol(expπ[m])));
pchoice = @expression(model, [m = 1:Mkts], expπ[m] ./ denom[m]);
share_eq = @constraint(model, [m = 1:Mkts, p = 1:Ps], mean(pchoice[m][p, :]) == d[m].share[p]);
π_vec = @expression(model, reduce(vcat, π))
IV = reduce(vcat, getproperty.(d, :iv))
π· = IV' * IV
EM = @expression(model, mean(IV .* π_vec; dims = 1));
obj = @expression(model, EM * π· * EM');
@objective(model, Min, obj[1, 1]);
optimize!(model)
value.(theta1)
value.(theta2)
```