Nonlinear Optimization with Many Constraints + Autodifferentiation: Which Julia Solution?

SebKrantz · February 24, 2024, 12:00am

I’m trying to rewrite a Matlab Package in Julia which solves a high-dimensional non-linear optimal transport problem in order to maximize economic welfare as predicted by a linked quantiative trade model. In particular, there is a vector of optimization variables x, say length 1000, a function objective(x) which computes a scalar welfare measure (the matlab implementation negates it so that it becomes a minimization problem), and a function constraints(x) which computes around 500 constraints, which should be satisfied with inequality i.e. all values should be <= 0 in the optimum. In addition, there are (optional) upper and lower bounds on the variables x, and an option to provide initial values x0.

The MatLab implementation uses the Ipopt optimizer and the ADiGator autodifferentiation package. The basic code looks like this:

n = size(x0, 1);
% Setup structure for adigatorGenFiles4Ipopt
setup.numvar=n;
setup.objective='objective'; % objective function as string
setup.order = 2;
setup.constraint='constraints'; % constraints function as string

% adigatorGenFiles4Ipopt generates everything required by ipopt (including derivatives)
funcs = adigatorGenFiles4Ipopt(setup);

% Call ipopt 
options.ipopt.tol = sqrt(eps);
options.lb = -Inf*ones(n,1); % lower bounds
options.ub = Inf*ones(n,1); % upper bounds
[z, info] = ipopt(x0, funcs, options);

Having tried this with JuMP, I am encountering errors, particularly regarding the constraints function, which is vector-valued. My basic setup would be

n = length(x0)
model = Model(Ipopt.Optimizer)
@variable(model, x[1:n])
JuMP.register(model, :objective, n, objective, autodiff=true)
JuMP.register(model, :constraints, n, constraints, autodiff=true)
# Set up the objective function to be minimized
@NLobjective(model, Min, objective(x...))

I already get an error in JuMP.register that has something to do with ForwardDiff. Then, GPT4 suggested that the correct way to specify the constraint is

for i in 1:length(constraints(x))
    @NLconstraint(model, constraints(x...)[i] <= 0)
end

There were also 2 sggestions for the bounds, either

@variable(model, x_lower[i] <= x[i=1:n] <= x_upper[i])

or

set_optimizer_attribute(model, "lb", x_lower)
set_optimizer_attribute(model, "ub", x_lower)

It seems the latter attempts to pass values directly to the optimizer. In any case, I didn’t manage to register the functions, and the loop-like way of dealing with the constraints seemed highly suboptimal. I wondered if there is a more direct way of calling Ipopt (the documentation suggests that is possible via the C API). I also wonder if using a package like Enzyme.jl might not be substantially better than ForwardDiff for the differentiation part. Anway, I’m pretty new to Julia and would be keen to know how a versed user would implement this problem.

mohamed82008 · February 24, 2024, 5:58am

JuMP may not be the most appropriate package for large vector-valued constraints although you can probably get it to work with help from JuMP developers here. You may also want to consider Nonconvex.jl (disclaimer: I am the main developer of Nonconvex), Optimization.jl or ADNLPModels.jl all of which allow the use of automatic differentiation with Ipopt (beside other solvers).

SebKrantz · February 24, 2024, 7:06am

Thanks, I figured this would be the way to go then?

  model = Nonconvex.Model(objective)
  Nonconvex.addvar!(model, x_lower, x_upper, init = x0, integer = repeat([false], length(x0))
  Nonconvex.add_ineq_constraint!(model, constraints)
  alg = IpoptAlg()
  opts = IpoptOptions()
  Nonconvex.optimize(model, alg, options = opts)

I get

 Mutating arrays is not supported -- called setindex!(Matrix{Float64}, ...)
This error occurs when you ask Zygote to differentiate operations that change
the elements of arrays in place (e.g. setting values with x .= ...)

I don’t see a mutating operation at first sight and the constraint function is complex, but I can compute gradient, jacobian and hessian using ForwardDiff. Is there a way to use ForwardDiff with Nonconvex?

That being said, I probably won’t be using Nonconvex for my production use because it has many dependencies. Currently looking into using the Ipopt C API wrapper directly. But your implementation of NonconvexIpopt is interesting to look at, in particular regarding determining the sparsity structure…

mohamed82008 · February 24, 2024, 7:13am

For other AD packages, see Using other AD packages · Nonconvex.jl, in particular the forwarddiffy(model) function maybe the most convenient. By default, Nonconvex uses Zygote which is picky about which functions it can differentiate.

mohamed82008 · February 24, 2024, 7:15am

If you want to hack around the Ipopt C API yourself, see Ipopt.jl and its Nonconvex wrapper for ideas.

gdalle · February 24, 2024, 8:27am

It seems that GPT4 is out of date: the current version of JuMP has very much improved and simplified nonlinear modeling, removing the need for @NL macros or registering.

In particular, check out the section on autodiff.

This is a compromise between speed and ease of use. If you have many input variables and one output variable, forward mode autodiff (as implemented in ForwardDiff.jl) will be very slow to compute gradients. Reverse mode autodiff (as implemented in Enzyme.jl or Zygote.jl) will be much faster. However, I don’t think JuMP supports anything beyond ForwardDiff.jl at the moment, so maybe Optimization.jl would be more suited?

In any case, I predict a complete and relevant answer from Oscar Dowson in about 3… 2… 1…

SebKrantz · February 24, 2024, 8:59am

Thanks! I’m still not managing to pass my vector-valued constraints function to the @constraints macro though. Or at least I continue to get errors. The docs as far as I can see continue to only cover cases where the constraints are explicitly laid out.

gdalle · February 24, 2024, 9:24am

Can you give a complete runnable example, including the definition of the variables, objective and constraints?

SebKrantz · February 24, 2024, 7:34pm

using JuMP, Ipopt, LinearAlgebra

n = 1000
x0 = abs.(rand(n))
auxdata = abs.(rand(Int(n / 2), n))
x_lower = repeat([-1e5], n)
x_upper = repeat([1e5], n)

function objective(x)
    return -sum(x)
end

function constraints(x)
    return auxdata * x
end

model = Model(Ipopt.Optimizer)
@variable(model, x_lower[i] <= x[i=1:n] <= x_upper[i])
@objective(model, Min, objective(x))
@constraint(model, constraints(x[1:n]) .<= 0)
optimize!(model)
objective_value(model)
value.(x)

This example works, but it seems it is hard-coding the values of auxdata in the constraints. My real constraints function is significantly more complex and there this gives an error. In any case, not really concinved that the interface is suitable to put large scale tasks into production.

gdalle · February 24, 2024, 9:02pm

I think you could directly do

@constraint(model, auxdata * x[1:n] .<= 0)

to avoid closures over global variables.

One of the things I like about JuMP is that it lets you add tons of small constraints instead of a single big one. So maybe if your function is too complex you can decompose it?

What kind of error? Can you give the stack trace and the line where it occurs? It is hard to help you with this level of detail.

I have personally helped put large tasks into production using JuMP, so I respectfully disagree Happy to help you debug yours as best I can though!

odow · February 24, 2024, 9:16pm

Hi @SebKrantz,

My real constraints function is significantly more complex and there this gives an error.

There are a lot of replies in this thread that are talking in generalities. It will be much easier if you can provide a reproducible example of your true problem.

The official suggestions are probably:

If you can write out your expressions algebraically. Do that instead.
If you must use the functional form, use Nonlinear Modeling · JuMP like you have tried:

model = Model(Ipopt.Optimizer)
@variable(model, x[1:n])
@objective(model, Min, objective(x...))
@constraint(model, constraints(x...) .<= 0)
optimize!(model)

But this will only work if you can trace your functions, which is seems like you cannot.

If neither of those things are appropriate. JuMP is not the right tool for the job.

In any case, not really concinved that the interface is suitable to put large scale tasks into production.

JuMP can and is used to solve very large scale nonlinear programs. See, e.g., AC Optimal Power Flow in Various Nonlinear Optimization Frameworks - #81 by ccoffrin.

However, one downside to Julia (and JuMP) is that a direct translation from MATLAB is often the wrong approach. If you can provide a reproducible example, people may have more concrete suggestions for improvements.

SebKrantz · February 24, 2024, 9:29pm

Ok, thanks! Sure, I can give you some details. In one of the simpler model cases the problem look like this. auxdata contains information pertaining to the graph structure and the economic model.


function objective(x)
    return -x[1]
end

function constraints_mobility(x, auxdata)

    # Extract parameters
    param = dict_to_namedtuple(auxdata[:param])
    graph = auxdata[:graph]
    kappa_ex = auxdata[:kappa_ex]
    A = auxdata[:A]

    # Extract optimization variables
    u = x[1]
    Cjn = reshape(x[2:graph.J*param.N+1], graph.J, param.N)
    Qin = reshape(x[graph.J*param.N+2:graph.J*param.N+graph.ndeg*param.N+1], graph.ndeg, param.N)
    Lj = x[graph.J*param.N+graph.ndeg*param.N+2:graph.J*param.N+graph.ndeg*param.N+graph.J+1]
    Cj = (sum(Cjn.^((param.sigma-1)/param.sigma), dims=2)).^(param.sigma/(param.sigma-1))
    Ljn = reshape(x[graph.J*param.N+graph.ndeg*param.N+graph.J+2:end], graph.J, param.N)
    Yjn = param.Zjn .* Ljn.^param.a

    # Utility constraint (Lj*u <= ... )
    cons_u = Lj .* u - (Cj / param.alpha).^param.alpha .* (param.Hj / (1 - param.alpha)).^(1 - param.alpha)

    # balanced flow constraints
    cons_Q = zeros(eltype(x), graph.J, param.N)
    for n in 1:param.N
        M = max.(A .* (ones(graph.J, 1) * sign.(Qin[:, n]')), 0)
        cons_Q[:, n] = Cjn[:, n] + A * Qin[:, n] - Yjn[:, n] + M * (abs.(Qin[:, n]).^(1 + param.beta) ./ kappa_ex)
    end

    # labor resource constraint
    cons_L = sum(Lj) - 1

    # Local labor availability constraints ( sum Ljn <= Lj )
    cons_Ljn = sum(Ljn, dims=2) .- Lj

    # return whole vector of constraints
    cons = vec(vcat(cons_u[:], cons_Q[:], cons_L, cons_Ljn))
    return cons
end

constraints = (x) -> constraints_mobility(x, auxdata)

The error I am getting is

ERROR: MethodError: Cannot `convert` an object of type AffExpr to an object of type VariableRef
Closest candidates are:
  convert(::Type{T}, ::T) where T at Base.jl:61
  GenericVariableRef{T}(::Any, ::Any) where T at ~/.julia/packages/JuMP/kSaGf/src/variables.jl:247
Stacktrace:
  [1] fill!(dest::Matrix{VariableRef}, x::AffExpr)
    @ Base ./array.jl:351
  [2] zeros(#unused#::Type{VariableRef}, dims::Tuple{Int64, Int64})
    @ Base ./array.jl:589
  [3] zeros(::Type{VariableRef}, ::Int64, ::Int64)
    @ Base ./array.jl:584

So in other words, it it not possible to generate a zeros() matrix filled with variable type references. I managed to somehow circumvent this problem but then get the next one in sign.() and abs.(). In general, my remark about this being impractical regards the fact that JuMP seems not able to simply treat constraints() as a black box julia function, but needs to digest it somehow to create an array of constraints. In the MatLab formulation this is not a problem. That’s why I am presently inclined to try and write a wrapper around the Ipopt C API myself to be able to seamlessly deal with the multiple versions of this problem (labour mobility yes, no, or only within regions, cross-good-congestion yes, no, convex or non-convex solution, etc.) in a straightforward manner.

odow · February 24, 2024, 9:45pm

The “JuMP” approach would be something like this (I can’t run, so there might be typos, etc):

function build_model(auxdata)
    param = dict_to_namedtuple(auxdata[:param])
    graph = auxdata[:graph]
    kappa_ex = auxdata[:kappa_ex]
    A = auxdata[:A]
    model = Model(Ipopt.Optimizer)
    @variable(model, u)
    @variable(model, Cjn[1:graph.J, 1:param.N])
    @variable(model, Qin[1:graph.ndeg, 1:param.N])
    @variable(model, Lj[1:graph.J])
    @objective(model, Max u)
    # Utility constraint (Lj*u <= ... )
    Cj = (sum(Cjn.^((param.sigma-1)/param.sigma), dims=2)).^(param.sigma/(param.sigma-1))
    @constraint(model, Lj .* u .<= (Cj / param.alpha).^param.alpha .* (param.Hj / (1 - param.alpha)).^(1 - param.alpha))
    Yjn = param.Zjn .* Ljn.^param.a
    # balanced flow constraints
    for n in 1:param.N
        M = max.(A .* (ones(graph.J, 1) * sign.(Qin[:, n]')), 0)
        @constraint(model, Cjn[:, n] + A * Qin[:, n] .<= Yjn[:, n] + M * (abs.(Qin[:, n]).^(1 + param.beta) ./ kappa_ex))
    end
    # labor resource constraint
    @constraint(model, sum(Lj) <= 1)
    # Local labor availability constraints ( sum Ljn <= Lj )
    @constraint(model, sum(Ljn, dims=2) .<= Lj)
    return model
end

I think this ends up being much more readable compared with the MATLAB approach of a single vector x and a f(x) <= 0 constraint.

Your error message was because
cons_Q = zeros(eltype(x), graph.J, param.N) create a matrix of variables, and you cannot store an expression in a matrix of variables. You could try instead cons_Q = Matrix{Any}(undef, graph.J, param.N).

You might want to read these parts of the JuMP documentation:

Design patterns for larger models · JuMP

and

SebKrantz · February 24, 2024, 10:26pm

Thanks, this is certainly very illuminating! I still have a problem with sign.() though.

ERROR: MethodError: no method matching sign(::VariableRef)
Closest candidates are:
  sign(::Unsigned) at number.jl:163
  sign(::Rational) at rational.jl:256
  sign(::Dates.Period) at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Dates/src/periods.jl:106
  ...

And I am wondering on how to add the bounds to the vector constraints. I presumes just adding the respective numbers to the left and right side of each constraint using .<=?

    options = Dict(
        # Bounds on x
        "lb" => vcat(-Inf, 1e-8 * ones(graph.J * param[:N]), -Inf * ones(graph.ndeg * param[:N]), 1e-8 * ones(graph.J), 1e-8 * ones(graph.J * param[:N])),
        "ub" => vcat(Inf, Inf * ones(graph.J * param[:N]), Inf * ones(graph.ndeg * param[:N]), ones(graph.J), Inf * ones(graph.J * param[:N])),
        # On the constraints
        "cl" => vcat(-Inf * ones(graph.J), -Inf * ones(graph.J * param[:N]), -1e-8, -1e-8 * ones(graph.J)),
        "cu" => vcat(-1e-8 * ones(graph.J), -1e-8 * ones(graph.J * param[:N]), 1e-8, 1e-8 * ones(graph.J))
    )

odow · February 24, 2024, 10:53pm

Ah, I didn’t see sign. This is not a function supported by JuMP.

Try something like:
M = @expression(model, max.(A .* (ones(graph.J, 1) * ifelse.(Qin[:, n]') .> 0, 1, -1), 0))

And I am wondering on how to add the bounds to the vector constraints.

These should all be supported:

@constraint(model, l .<= f(x) .<= u)
@constraint(model, f(x) .<= u)
@constraint(model, f(x) .>= l)

I think in general, your code could be much simpler if you tried to avoid MATLAB-esque vectorization of everything.

So instead of

    Cj = (sum(Cjn.^((param.sigma-1)/param.sigma), dims=2)).^(param.sigma/(param.sigma-1))
    @constraint(model, Lj .* u .<= (Cj / param.alpha).^param.alpha .* (param.Hj / (1 - param.alpha)).^(1 - param.alpha))

do

    psigma = (param.sigma - 1) / param.sigma
    for j in 1:graph.J
        Cj = sum(Cjn[j, n]^psigma for n in 1:param.N)^(1 / psigma)
        @constraint(model, Lj[j] * u <= (Cj / param.alpha)^param.alpha * (param.Hj / (1 - param.alpha))^(1 - param.alpha))
    end

I don’t really understand and the balanced flow constraint without data. It’s probably too easy to make a mistake with shapes etc.

Is it something like this?

@constraint(
    model, 
    [n in 1:param.N, j in 1:param.J],
    Cjn[j, n] + sum(A[j, i] * Qin[i, n] for i in 1:graph.ndeg) <=
    Yjn[j, n] + sum(
        max(ifelse(Qin[i, n] > 0, A[i, j], -A[i, j]), 0) *
        abs(Qin[i, n])^(1 + param.beta) / kappa_ex
        for i in 1:graph.ndeg
    )
)

Note that you do not need to use vectorization for performance in JuMP or Julia.

For variable bounds, do:

@variable(model, Cjn[1:graph.J, 1:param.N] >= 1e-8)
@variable(model, 1e-8 <= Lj[1:graph.J] <= 1)

I guess all these points get back to my comment that

However, one downside to Julia (and JuMP) is that a direct translation from MATLAB is often the wrong approach.

You can write nice JuMP code that scales to large problems, but you can’t write it like you do in MATLAB. There are pros and cons to that.

We should add some documentation for users transitioning from other modeling languages: Suggestions for documentation improvements · Issue #2348 · jump-dev/JuMP.jl · GitHub

I’ve also opened an issue to support sign(x::VariableRef): Add support for nonlinear sign(x) · Issue #3682 · jump-dev/JuMP.jl · GitHub

odow · February 24, 2024, 11:13pm

Also, it looks like you want to have constraints -1e-8 <= f(x) <= 1e-8? Do not do this. Use f(x) == 0 instead. Solvers like Ipopt will automatically add tolerances to equality constraints.

gdalle · February 25, 2024, 10:19am

@odow while you’re in the neighborhood, did I read correctly that no other autodiff backend can be used except ForwardDiff?

odow · February 25, 2024, 9:33pm

For algebraic expressions, JuMP does not use ForwardDiff, but a custom sparse reverse mode AD.

You can change the AD backend by setting the MOI.AutomaticDifferentiationBackend attribute, but currently the only other implementation is GitHub - odow/MathOptSymbolicAD.jl.

For user-defined operators, JuMP uses ForwardDiff by default. To use a different AD system, you can provide your own gradients and Hessians: Nonlinear Modeling · JuMP

SebKrantz · February 26, 2024, 3:48am

Hi @odow, thanks so much again for your help. After some trial and error I finally got things to work out, and numerically practically identical results. A couple of remaining questions. The first relates to the question of @gdalle about autodifferentiation. I find that the JuMP solution does not seem to detect the sparsity of the hessian (like ADiGator in MatLab). The Matlab output is:

This is Ipopt version 3.11.8, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:    51788
Number of nonzeros in Lagrangian Hessian.............:      783

Total number of variables............................:      784
                     variables with only lower bounds:      242
                variables with lower and upper bounds:      121
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:      364
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:      122
        inequality constraints with only upper bounds:      242

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -0.0000000e+00 2.10e-01 1.61e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
Reallocating memory for MA57: lfact (245982)
Reallocating memory for MA57: lfact (307722)
Reallocating memory for MA57: lfact (468079)
Reallocating memory for MA57: lfact (575242)
   1 -1.9344740e-05 2.10e-01 7.81e+00  -1.0 9.67e+01  -4.0 3.79e-01 2.00e-07f  2
   2 -1.8633038e+01 9.08e-01 5.49e+00  -1.7 3.75e+01    -  6.43e-01 4.97e-01f  1
   3 -1.7498744e+01 5.32e-01 4.99e+00  -2.5 2.12e+00    -  8.71e-01 5.34e-01h  1
   4 -1.4461303e+01 1.73e-01 2.50e+00  -2.5 4.13e+00    -  1.00e+00 7.35e-01h  1
   5 -1.1547139e+01 6.69e-02 2.38e+00  -2.5 2.91e+00    -  1.00e+00 1.00e+00h  1
   6 -1.1768472e+01 2.06e-02 7.31e-01  -3.8 2.70e-01    -  6.47e-01 8.20e-01f  1
   7 -1.1793453e+01 5.08e-04 8.70e-02  -3.8 2.50e-02    -  1.00e+00 1.00e+00h  1
   8 -1.1791451e+01 0.00e+00 2.14e-04  -3.8 2.00e-03    -  1.00e+00 1.00e+00h  1
   9 -1.1827654e+01 3.04e-07 2.34e-04  -8.6 3.63e-02    -  9.96e-01 9.98e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -1.1827780e+01 0.00e+00 8.75e-05  -8.6 1.26e-04    -  7.76e-01 1.00e+00h  1
  11 -1.1827784e+01 0.00e+00 2.65e-05  -8.6 3.57e-06    -  7.15e-01 1.00e+00h  1
  12 -1.1827789e+01 0.00e+00 6.76e-06  -8.6 5.34e-06    -  8.07e-01 1.00e+00f  1
  13 -1.1827790e+01 0.00e+00 1.18e-13  -8.6 5.97e-07    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 13

                                   (scaled)                 (unscaled)
Objective...............:  -1.1827789684933572e+01   -1.1827789684933572e+01
Dual infeasibility......:   1.1752412313047493e-13    1.1752412313047493e-13
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.6791510920999567e-09    2.6791510920999567e-09
Overall NLP error.......:   2.6791510920999567e-09    2.6791510920999567e-09


Number of objective function evaluations             = 15
Number of objective gradient evaluations             = 14
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 15
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 14
Number of Lagrangian Hessian evaluations             = 13
Total CPU secs in IPOPT (w/o function evaluations)   =      0.271
Total CPU secs in NLP function evaluations           =      0.863

EXIT: Optimal Solution Found.

The Julia Output is:

This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:    51788
Number of nonzeros in Lagrangian Hessian.............:    51425

Total number of variables............................:      784
                     variables with only lower bounds:      242
                variables with lower and upper bounds:      121
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:      364
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:      122
        inequality constraints with only upper bounds:      242

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 7.49e-03 1.38e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.8355390e+01 8.31e-01 8.65e+01  -1.0 7.90e+03    -  9.90e-01 1.24e-02f  1
   2  3.6539782e+01 1.89e-01 5.74e+01  -1.0 8.53e+01    -  9.90e-01 7.24e-01h  1
   3  1.4808214e+01 2.37e-02 2.42e+01  -1.7 2.64e+01    -  8.92e-01 8.25e-01h  1
   4  7.9511084e+00 1.61e-03 1.30e+01  -1.7 6.86e+00    -  7.53e-01 1.00e+00h  1
   5  8.0558347e+00 1.55e-03 1.88e+01  -2.5 3.55e+00    -  8.41e-01 2.95e-02f  1
   6  1.0874882e+01 1.95e-02 1.40e+00  -2.5 2.82e+00    -  9.97e-01 1.00e+00f  1
   7  1.1565232e+01 8.65e-03 5.65e-01  -3.8 9.02e-01    -  7.44e-01 7.66e-01f  1
   8  1.1790340e+01 2.14e-04 2.33e-02  -3.8 2.25e-01    -  1.00e+00 1.00e+00f  1
   9  1.1791296e+01 0.00e+00 4.08e-05  -3.8 9.56e-04    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.1827237e+01 0.00e+00 1.33e-03  -8.6 3.64e-02    -  9.96e-01 9.87e-01f  1
  11  1.1827780e+01 0.00e+00 1.79e-04  -8.6 5.43e-04    -  7.75e-01 1.00e+00h  1
  12  1.1827784e+01 0.00e+00 5.21e-05  -8.6 3.50e-06    -  7.80e-01 1.00e+00h  1
  13  1.1827790e+01 0.00e+00 3.19e-06  -8.6 5.90e-06    -  9.46e-01 1.00e+00f  1
  14  1.1827790e+01 0.00e+00 3.67e-14  -8.6 2.77e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 14

                                   (scaled)                 (unscaled)
Objective...............:  -1.1827789680100521e+01    1.1827789680100521e+01
Dual infeasibility......:   3.6653988763450662e-14    3.6653988763450662e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5669821264002667e-09    2.5669821264002667e-09
Overall NLP error.......:   2.5669821264002667e-09    2.5669821264002667e-09


Number of objective function evaluations             = 15
Number of objective gradient evaluations             = 15
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 15
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 15
Number of Lagrangian Hessian evaluations             = 14
Total seconds in IPOPT                               = 0.723

EXIT: Optimal Solution Found.

Here I’ve set up the problem exactly the way it is in Matlab. with upper and lower bounds (starting values are not set in either Matlab or Julia in this case):

function build_model(auxdata)
    param = dict_to_namedtuple(auxdata[:param])
    graph = auxdata[:graph]
    kappa_ex = auxdata[:kappa_ex]
    A = auxdata[:A]
    psigma = (param.sigma - 1) / param.sigma
    Hj = param.Hj
    # Model
    model = Model(Ipopt.Optimizer)
    # Variables + Bounds
    @variable(model, u)                                    # Overall utility
    # set_start_value(u, 0.0)
    @variable(model, Cjn[1:graph.J, 1:param.N] >= 1e-8)    # Good specific consumption
    # set_start_value.(Cjn, 1.0e-6)
    @variable(model, Qin[1:graph.ndeg, 1:param.N])         # Good specific flow
    # set_start_value.(Qin, 0.0)
    @variable(model, 1e-8 <= Lj[1:graph.J] <= 1)           # Total Labour
    # set_start_value.(Lj, 1 / graph.J)
    @variable(model, Ljn[1:graph.J] >= 1e-8)               # Good specific labour
    # set_start_value.(Ljn, 1 / (graph.J * param.N))
    @objective(model, Max, u)
    # Utility constraint (Lj*u <= ... )
    for j in 1:graph.J
        Cj = sum(Cjn[j, n]^psigma for n in 1:param.N)^(1 / psigma)
        @constraint(model, Lj[j] * u - (Cj / param.alpha)^param.alpha * (Hj[j] / (1 - param.alpha))^(1 - param.alpha) <= -1e-8)
    end
    Yjn = param.Zjn .* Ljn.^param.a
    # balanced flow constraints
    @constraint(
        model, 
        [n in 1:param.N, j in 1:param.J],
        Cjn[j, n] + sum(A[j, i] * Qin[i, n] for i in 1:graph.ndeg) -
        Yjn[j, n] + sum(
            max(ifelse(Qin[i, n] > 0, A[j, i], -A[j, i]), 0) *
            abs(Qin[i, n])^(1 + param.beta) / kappa_ex[i]
            for i in 1:graph.ndeg
        ) <= -1e-8
    )
    # labor resource constraint
    @constraint(model, -1e-8 <= sum(Lj) - 1 <= 1e8)
    # Local labor availability constraints ( sum Ljn <= Lj )
    @constraint(model, -1e-8 .<= sum(Ljn, dims=2) .- Lj .<= 1e-8)
    return model
end

I have experimented around with using equality constraints or relinguishing the bounds, but then the solution takes much longer and the results are not the same, so I guess the authors of the model knew what they were doing when specifying bounds on everything. But yeah, for larger problems a way of autodifferentiation that detects sparseness in the Hessian would be great.

I was also not sure if this is the recommended way to do so, but I am currently accessing the results using

results = Dict(k => value.(v) for (k, v) in model.obj_dict)

Also: the library actually iterates over these solutions for different values of auxdata (the task is actually finding the economically optimal transport network, so each iteration does some changes to the weights attached to graph edges: kappa_ex represents the cost of moving goods along each edge which depends on road network - one reason I am very interested in computational efficiency in derivative evaluation). I thus also wondered if there is a more efficient way of setting this up than calling build_model(auxdata) in each iteration (i.e. building the model once and just changing its parameters in each iteration)? Many thanks again for all your help!

odow · February 26, 2024, 4:22am

I find that the JuMP solution does not seem to detect the sparsity of the hessian

JuMP does detect and exploit sparsity in the Hessian.

The difference is probably that ADiGator sums elements in the Hessian and so passes one element or each unique non-zero, whereas JuMP passes them as a list with duplicates that Ipopt sums internally. They are probably equally efficient.

(You can also see that the dense Hessian would have O(784^2 / 2) elements, so JuMP is detecting sparsity.

I was also not sure if this is the recommended way to do so

That’s okay, but note that if you name constraints or have expressions, these will also be evaluated.

There is no “standard” way to collate all results because models can vary so much between users. We recommend that people think about their data structures and return what is most suited to the particular model.

so each iteration does some changes to the weights attached to graph edges representing infrastructure, kappa_ex represents the cost of moving goods along each edge which depends on road network

Use a parameter for kappa_ex: