Too many nonzeroes for NLP in JuMP and Examodels

Creating a topic as suggested by @odow and mentioned in topic.

The issue is that JuMP creates Jacobian and Hessian matrices with nonzeroes for derivatives that are fixed zero.

Consider the following two example constraints:
@constraint(m,P1, P == Vn^2 - Vn*Vm*(cos(fin-fim) + sin(fin-fim)) )
@constraint(m,P2, P == Vn^2 - Vn*Vm*(fin-fim) )

If only single constraint at the time is run with Ipopt (without objective function), one can find that JuMP makes a model with 5 nonzeroes in Jacobian and 10 nonzeroes in Hessian for both of them. For the same test AMPL makes model (even with “option presolve 0”):

  1. constraint - 5 nonzeroes in Jacobian and 9 nonzeroes in Hessian;
  2. constraint - 5 nonzeroes in Jacobian and 6 nonzeroes in Hessian.

I observed the same with Examodels native api.

In addition, to my memory, Examodels generated nonzeroes for the RHS of the following constraint despite multiplication by zero.
P == 0*Vn^2

Overall I observed worse convergence with Julia NLP optimization than with AMPL for difficult cases (OPF cases that need soft constraints to have a solution). My assumption is that excessive number of nonzeroes could impact convergence (due to solver scaling factors). How can we fix this?

Can you provide a fully reproducible example of the convergence difference for the problem with soft constraints, including the input data?

Note that the number of non-zeros can be misleading. JuMP has a conservative algorithm for sparsity detection that generates false positives. In reality, some of those will be zero. But this shouldn’t really impact the solve that much? Are the variables and constraints in the same order in JuMP and AMPL? Are the solver options the same?

Also note that you can use AMPL’s automatic differentiation instead of JuMP’s using:

using JuMP, AmplNLWriter, Ipopt_jll
model = Model(() -> AmplNLWriter.Optimizer(Ipopt_jll.amplexe))
1 Like

This is the same issue as [Nonlinear] sparsity pattern of Hessian with :(x * y) · Issue #2527 · jump-dev/MathOptInterface.jl · GitHub

julia> using JuMP

julia> begin
           model = Model()
           @variable(model, P)
           @variable(model, Vn)
           @variable(model, Vm)
           @variable(model, fin)
           @variable(model, fim)
           @constraint(model, P == Vn^2 - Vn*Vm*(fin-fim))
           nothing
       end

julia> begin
           rows = Any[]
           nlp = MOI.Nonlinear.Model()
           for (F, S) in list_of_constraint_types(model)
               for ci in all_constraints(model, F, S)
                   push!(rows, ci)
                   object = constraint_object(ci)
                   MOI.Nonlinear.add_constraint(nlp, object.func, object.set)
               end
           end
           MOI.Nonlinear.set_objective(nlp, objective_function(model))
           x = all_variables(model)
           backend = MOI.Nonlinear.SparseReverseMode()
           evaluator = MOI.Nonlinear.Evaluator(nlp, backend, index.(x))
           MOI.initialize(evaluator, [:Hess])
           hessian_sparsity = MOI.hessian_lagrangian_structure(evaluator)
           I = [i for (i, _) in hessian_sparsity]
           J = [j for (_, j) in hessian_sparsity]
           V = ones(length(hessian_sparsity))
           MOI.eval_hessian_lagrangian(evaluator, V, ones(length(x)), 1.0, ones(length(rows)))
           H = SparseArrays.sparse(I, J, V, length(x), length(x))
       end
5Ă—5 SparseMatrixCSC{Float64, Int64} with 10 stored entries:
  â‹…     â‹…     â‹…    â‹…    â‹… 
  â‹…   -2.0    â‹…    â‹…    â‹… 
  â‹…    0.0   0.0   â‹…    â‹… 
  â‹…    1.0   1.0  0.0   â‹… 
  â‹…   -1.0  -1.0  0.0  0.0

Maybe there is an argument for special casing sparsity of *. But I still don’t think this will be the cause of the difference with AMPL. But we’d need a reproducible example to confirm.

I am aware that MWE is needed to help with this. However, I am not able to provide it. Convergence issues that I had were with Examodels actually (which has the same described nonzero issue). The model was not “minimal” and I did not use it in considerable time. It would take long for me to reproduce it again 1:1 with AMPL.

JuMP overall at default is generally considerably slower than AMPL (if not set to use symbolic derivatives in JuMP at least). Examodels on CPU is slightly faster then AMPL and that was the reason why I tried it. AmplNLWriter uses AMPL’s AD, but also requires AMPL license so I get better possibilities if I write directly in AMPL.

This is incorrect. We use AMPLs open source ASL library which does not require an AMPL license.

If ExaModels.jl does not converge im sure they would like a bug report.

And yes, JuMP is slower, which is why we wrote MathOptSymbolicAD.jl but the choice of AD backend should not change convergence.

You’ve probably already seen it, but see

1 Like

Thank you. I did not know ASL was open source. I think I used Knitro when I tried AmplNLWriter and in that case, I think, it needs the Knitro license from AMPL. This is valuable to me.

I have seen very similar powerpoint presentation, but not the video. A video has some more information.

When I will have clear examole that does not converge I will post an issue. I am active in this regard. I don’t let many things “just pass”.

I think I used Knitro

KNITRO is commercial software. But you don’t need a license for JuMP+AmplNLWriter+Ipopt_jll