Hi everybody, I am trying to solve a large-scale nonlinear problem that is convex and has linear inequality constraints. There are 550 unknowns and 25000 linear inequality constraints. I want to get your recommendations on which solver is better suited for solving this problem and if there are any obvious performance flaws in my current formulation.

The Mathematical Problem is the following:

\min_{x_{i}\geq 0 , y_{j}} \sum_{i=1}^{I}x_{i}C_{i} + \sum_{j=1}^{J}a_{j}{y_{j}}^{1-\sigma}

\text{s.t}\; \; y_{j} \leq \tau_{ij}\left(\frac{w_{i}}{\varepsilon_{i}}+x_{i}\right) \hspace{5.5pt}\forall i,j

with \sigma>1.

In the MWE I am providing I am using JuMP with Ipopt. It seems to me that the fact that I have the y_{j}^{1-\sigma} terms prevent me from using Convex.jl or a more powerful commercial solver like Gurobi or Mosek. One option would be to linearize this term and add an extra constraint, but that would add j additional variables and constraints and it is not clear to me that a convex solver can handle this type of concern. At this point the performance of my MWE is disappointing as I would need to improve this around 100x to make it feasible to implement.

Are there any other recommended solvers (open-source or commercial) that you would recommend?

```
using JuMP, Ipopt, LinearAlgebra, Random, Distributions, BenchmarkTools
const Ι = 500
const J = 50
const σ = 5
#Create a vector of productivities
const z=exp.(rand(Normal(0,1),Ι))
const β=exp.(rand(Normal(0,1),J))./10
C = ones(Ι).*100000000
# Create a Matrix of Distances
function distmat(J,Ι)
Distances = zeros(J,Ι)
Random.seed!(7777)
coordinates_market = 100*rand(J,2)
coordinates_plant = 100*rand(Ι,2)
for j = 1:J, l=1:Ι
Distances[j,l] = sqrt((coordinates_market[j,1]-coordinates_plant[l,1])^2+(coordinates_market[j,2]-coordinates_plant[l,2])^2)
end
return 1 .+ Distances./100
end
const τ = distmat(J,Ι)
function solving_max_pi_constraints_dual(C,β,τ,z)
(J,Ι) = size(τ)
dual_profits = Model(Ipopt.Optimizer)
#set_optimizer_attribute(dual_profits, "nlp_scaling_method", "none")
#set_silent(dual_profits)
a = ((σ-1)/σ)^(σ-1)/(σ).*β.^(σ)
#set_optimizer_attribute(primal_capacity, "nlp_scaling_method", "none")
@variable(dual_profits, x[1:Ι] >= 0)
@variable(dual_profits, y[1:J]>= 0)
@expression(dual_profits, part1, x'*C)
@NLexpression(dual_profits, part2, sum(a[j]*y[j]^(1-σ) for j=1:J) )
@NLobjective(dual_profits, Min, part1 + part2)
for i=1:Ι, j=1:J
@constraint(dual_profits, τ[j,i]/z[i] + τ[j,i]*x[i] - y[j] >= 0)
end
optimize!(dual_profits)
return objective_value(dual_profits), value.(x), value.(y)
end
(sol,xsol,ysol) = @time solving_max_pi_constraints_dual(C,β,τ,z)
Total seconds in IPOPT = 6.624
EXIT: Optimal Solution Found.
6.693756 seconds (9.70 M allocations: 235.620 MiB, 0.57% gc time)
```