search

jump-dev / Ipopt.jl

A Julia interface to the Ipopt nonlinear solver
https://github.com/coin-or/ipopt
Other
152 stars 58 forks source link

Incorrect number of Hessian structure (nonzero entries) #353

Closed marco83pt closed 1 year ago

marco83pt commented 1 year ago

Dear All, I'm experiencing an issue using Julia/JuMP with Ipopt. In particular, the number of non-zero entries (also reported at the beginning on the solver printout) are 10x more than expected (evaluated analytically). While using *.nl file gives the righ problem structure.

I did a small investigation (by checking the Hessian structure MOI.hessian_lagrangian_structure) and I realized that duplicated nonzero entries (same matrix position) are not summed up.

The solver behaves appropriately, but the solution is not scalable (while, again, proper problem formulation is created by using *.nl file with Ipopt).

Thanks for your support. Marco

odow commented 1 year ago

I realized that duplicated nonzero entries (same matrix position) are not summed up.

This is a known feature and intentional on our part.

It arises because of how MathOptInterface computes the Hessian using automatic differentiation. Can you share a reproducible example that demonstrates the issue? There might be somethings we can do to improve.

We have some work underway to allow switching the AD backend for some problems. See:

At the moment, you can use AmplNLWriter to switch to the NL formulation:

using JuMP, Ipopt
model = Model(Ipopt.Optimizer)

becomes

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

If the problem is nicely structured, you could also try https://github.com/odow/MathOptSymbolicAD.jl

marco83pt commented 1 year ago

Dear Odow, Thanks a lot for your support. It is really appreciated. Sure, here a replicable test case.

First piece of code:

m = read_from_file("demos\\problem.nl") set_optimizer(m, Ipopt.Optimizer) @time jump_solve_time = optimize!(m)

Returns: 12.683538 seconds (6.26 M allocations: 542.467 MiB, 3.35% gc time) Number of nonzeros in equality constraint Jacobian...: 65471 Number of nonzeros in inequality constraint Jacobian.: 13976 Number of nonzeros in Lagrangian Hessian.............: 314291 Total seconds in IPOPT = 11.674 Objective...............: 2.3836722884185951e+03 2.3836722884185951e+03 Number of Iterations....: 41

Second piece of code:

m = read_from_file("demos\\problem.nl") set_optimizer(m, () -> AmplNLWriter.Optimizer(Ipopt_jll.amplexe)) @time ampl_solve_time = optimize!(m)

Returns: 8.320853 seconds (5.85 M allocations: 402.056 MiB, 1.75% gc time) Number of nonzeros in equality constraint Jacobian...: 65471 Number of nonzeros in inequality constraint Jacobian.: 13976 Number of nonzeros in Lagrangian Hessian.............: 36439 Total seconds in IPOPT = 6.914 Objective...............: 2.3836722884185870e+03 2.3836722884185870e+03 Number of Iterations....: 41

Considerations:

I hope this data set can help you understand this behaviour better. With bigger cases, this slowdown is also amplified.

problem.zip

odow commented 1 year ago

Do you have the JuMP script that builds the model?

odow commented 1 year ago

Here's what I get. It's hard to say more without the JuMP script. The NL file doesn't preserve all of the structure. For example, it treats linear constraints as nonlinear, where as JuMP+Ipopt.jl directly has special detection and support for linear and quadratic functions.

But I'd mark this as expected behavior. The choice of AD system can have a large impact on the solve time, and no single AD system is best for all models. So it's nice that we can easily try out multiple different ones:

julia> using JuMP

julia> import AmplNLWriter

julia> import Ipopt

julia> import Ipopt_jll

julia> import MathOptSymbolicAD

julia> model = read_from_file("/Users/Oscar/Downloads/problem.nl")
A JuMP Model
Minimization problem with:
Variables: 8744
Objective function type: Nonlinear
`VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 4724 constraints
`VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 4724 constraints
Nonlinear: 15029 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> set_optimizer(model, Ipopt.Optimizer)

julia> @time optimize!(model)
This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.

Number of nonzeros in equality constraint Jacobian...:    65471
Number of nonzeros in inequality constraint Jacobian.:    13976
Number of nonzeros in Lagrangian Hessian.............:   314291

Total number of variables............................:     8718
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     4698
                     variables with only upper bounds:        0
Total number of equality constraints.................:     8041
Total number of inequality constraints...............:     6988
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:     6988
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.2061341e+03 1.31e+01 2.52e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.3795678e+03 1.14e+01 1.33e+01  -1.0 6.12e+00    -  1.04e-02 1.31e-01h  1
   2  1.3651384e+03 1.11e+01 1.22e+01  -1.0 1.35e+01    -  6.53e-02 2.22e-02f  1
   3  1.4170592e+03 1.03e+01 8.58e+00  -1.0 1.24e+01    -  4.36e-01 7.11e-02h  1
   4  1.8085839e+03 4.98e+00 5.47e+00  -1.0 9.05e+00    -  1.90e-01 4.87e-01h  1
   5  2.0721391e+03 2.29e+00 2.62e+00  -1.0 3.45e+00    -  5.16e-01 5.30e-01h  1
   6  2.1845834e+03 1.34e+00 1.79e+00  -1.0 4.89e+00    -  3.31e-01 4.09e-01h  1
   7  2.2408122e+03 9.41e-01 1.67e+00  -1.0 2.76e+00    -  1.48e-01 2.98e-01h  1
   8  2.2606240e+03 8.02e-01 1.18e+00  -1.0 1.90e+00    -  4.36e-01 1.48e-01h  1
   9  2.2878726e+03 6.16e-01 1.06e+00  -1.0 1.50e+00    -  3.65e-01 2.32e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.3122630e+03 4.42e-01 7.08e-01  -1.0 1.09e+00    -  1.50e-01 2.82e-01h  1
  11  2.3328322e+03 2.91e-01 6.89e-01  -1.7 1.22e+00    -  2.85e-01 3.42e-01h  1
  12  2.3392720e+03 2.53e-01 4.92e-01  -1.7 2.49e+00    -  1.72e-01 1.29e-01h  1
  13  2.3434938e+03 2.28e-01 3.69e-01  -1.7 1.38e+00    -  1.90e-01 1.01e-01h  1
  14  2.3477968e+03 2.02e-01 3.51e-01  -1.7 1.07e+00    -  1.18e-01 1.12e-01h  1
  15  2.3502357e+03 1.86e-01 6.20e+00  -1.7 9.47e-01    -  7.57e-01 8.01e-02h  1
  16  2.3702492e+03 1.05e-01 1.09e+00  -1.7 1.97e+00    -  1.08e-02 4.49e-01h  1
  17  2.3784950e+03 7.88e-02 8.06e-01  -1.7 1.59e+00    -  1.00e+00 2.53e-01h  1
  18  2.4020780e+03 9.54e-03 1.33e+00  -1.7 7.65e-01    -  1.30e-01 1.00e+00h  1
  19  2.4029110e+03 9.74e-05 1.30e-02  -1.7 9.84e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  2.3928236e+03 1.22e-02 6.46e-02  -3.8 8.91e-01    -  8.21e-01 9.62e-01f  1
  21  2.3897847e+03 1.41e-02 1.57e+00  -3.8 1.56e+00    -  7.10e-01 2.85e-01f  1
  22  2.3875809e+03 1.30e-02 4.98e+00  -3.8 1.97e+00    -  8.70e-01 2.19e-01f  1
  23  2.3852591e+03 1.43e-02 8.37e+00  -3.8 1.17e+00    -  1.00e+00 3.18e-01f  1
  24  2.3847405e+03 1.24e-02 6.89e+00  -3.8 2.72e+00    -  1.00e+00 1.55e-01f  1
  25  2.3841607e+03 9.14e-03 8.39e+00  -3.8 2.76e+00    -  5.17e-01 3.20e-01f  1
  26  2.3838886e+03 6.33e-03 1.35e+01  -3.8 2.34e+00    -  5.83e-01 3.34e-01h  1
  27  2.3838154e+03 4.33e-03 1.17e+01  -3.8 1.99e+00    -  5.73e-01 3.32e-01h  1
  28  2.3837640e+03 2.58e-03 9.41e+00  -3.8 1.00e+00    -  9.75e-01 4.24e-01h  1
  29  2.3837400e+03 7.95e-04 2.74e+00  -3.8 6.09e-01    -  7.86e-01 7.07e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  2.3837407e+03 5.22e-05 2.00e-04  -3.8 1.71e-01    -  1.00e+00 1.00e+00h  1
  31  2.3837030e+03 9.63e-05 2.18e+00  -5.7 9.49e-01    -  7.83e-01 5.35e-01f  1
  32  2.3836811e+03 3.37e-04 5.59e-01  -5.7 2.00e+00    -  6.92e-01 7.10e-01f  1
  33  2.3836736e+03 2.13e-04 2.58e-01  -5.7 1.02e+00    -  7.61e-01 9.29e-01f  1
  34  2.3836732e+03 3.40e-05 3.14e-02  -5.7 2.27e-01    -  1.00e+00 8.54e-01f  1
  35  2.3836731e+03 1.70e-05 1.66e-02  -5.7 5.13e-02    -  1.00e+00 5.00e-01f  2
  36  2.3836731e+03 2.20e-07 2.82e-07  -5.7 2.91e-02    -  1.00e+00 1.00e+00h  1
  37  2.3836725e+03 7.11e-07 1.67e-03  -8.6 8.63e-02    -  7.90e-01 7.81e-01f  1
  38  2.3836723e+03 1.59e-06 3.15e-03  -8.6 2.05e-01    -  8.92e-01 9.72e-01f  1
  39  2.3836723e+03 2.43e-08 6.03e-08  -8.6 5.45e-03    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  2.3836723e+03 4.49e-10 1.03e-08  -8.6 1.32e-01    -  1.00e+00 1.00e+00h  1
  41  2.3836723e+03 2.74e-11 6.32e-10  -9.0 6.90e-03    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 41

                                   (scaled)                 (unscaled)
Objective...............:   2.3836722884186038e+03    2.3836722884186038e+03
Dual infeasibility......:   6.3247362678235533e-10    6.3247362678235533e-10
Constraint violation....:   2.7384317036194261e-11    2.7384317036194261e-11
Variable bound violation:   1.2994724585269068e-07    1.2994724585269068e-07
Complementarity.........:   1.2578999124635368e-09    1.2578999124635368e-09
Overall NLP error.......:   1.2578999124635368e-09    1.2578999124635368e-09

Number of objective function evaluations             = 43
Number of objective gradient evaluations             = 42
Number of equality constraint evaluations            = 43
Number of inequality constraint evaluations          = 43
Number of equality constraint Jacobian evaluations   = 42
Number of inequality constraint Jacobian evaluations = 42
Number of Lagrangian Hessian evaluations             = 41
Total seconds in IPOPT                               = 9.758

EXIT: Optimal Solution Found.
 10.555810 seconds (10.55 M allocations: 669.945 MiB, 2.58% gc time)

julia> @time optimize!(model; _differentiation_backend = MathOptSymbolicAD.DefaultBackend())
This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.

Number of nonzeros in equality constraint Jacobian...:    65471
Number of nonzeros in inequality constraint Jacobian.:    13976
Number of nonzeros in Lagrangian Hessian.............:    36439

Total number of variables............................:     8718
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     4698
                     variables with only upper bounds:        0
Total number of equality constraints.................:     8041
Total number of inequality constraints...............:     6988
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:     6988
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.2061341e+03 1.31e+01 2.52e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.3795678e+03 1.14e+01 1.33e+01  -1.0 6.12e+00    -  1.04e-02 1.31e-01h  1
   2  1.3651384e+03 1.11e+01 1.22e+01  -1.0 1.35e+01    -  6.53e-02 2.22e-02f  1
   3  1.4170592e+03 1.03e+01 8.58e+00  -1.0 1.24e+01    -  4.36e-01 7.11e-02h  1
   4  1.8085839e+03 4.98e+00 5.47e+00  -1.0 9.05e+00    -  1.90e-01 4.87e-01h  1
   5  2.0721391e+03 2.29e+00 2.62e+00  -1.0 3.45e+00    -  5.16e-01 5.30e-01h  1
   6  2.1845834e+03 1.34e+00 1.79e+00  -1.0 4.89e+00    -  3.31e-01 4.09e-01h  1
   7  2.2408122e+03 9.41e-01 1.67e+00  -1.0 2.76e+00    -  1.48e-01 2.98e-01h  1
   8  2.2606240e+03 8.02e-01 1.18e+00  -1.0 1.90e+00    -  4.36e-01 1.48e-01h  1
   9  2.2878726e+03 6.16e-01 1.06e+00  -1.0 1.50e+00    -  3.65e-01 2.32e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.3122630e+03 4.42e-01 7.08e-01  -1.0 1.09e+00    -  1.50e-01 2.82e-01h  1
  11  2.3328322e+03 2.91e-01 6.89e-01  -1.7 1.22e+00    -  2.85e-01 3.42e-01h  1
  12  2.3392720e+03 2.53e-01 4.92e-01  -1.7 2.49e+00    -  1.72e-01 1.29e-01h  1
  13  2.3434938e+03 2.28e-01 3.69e-01  -1.7 1.38e+00    -  1.90e-01 1.01e-01h  1
  14  2.3477968e+03 2.02e-01 3.51e-01  -1.7 1.07e+00    -  1.18e-01 1.12e-01h  1
  15  2.3502357e+03 1.86e-01 6.20e+00  -1.7 9.47e-01    -  7.57e-01 8.01e-02h  1
  16  2.3702492e+03 1.05e-01 1.09e+00  -1.7 1.97e+00    -  1.08e-02 4.49e-01h  1
  17  2.3784950e+03 7.88e-02 8.06e-01  -1.7 1.59e+00    -  1.00e+00 2.53e-01h  1
  18  2.4020780e+03 9.54e-03 1.33e+00  -1.7 7.65e-01    -  1.30e-01 1.00e+00h  1
  19  2.4029110e+03 9.74e-05 1.30e-02  -1.7 9.84e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  2.3928236e+03 1.22e-02 6.46e-02  -3.8 8.91e-01    -  8.21e-01 9.62e-01f  1
  21  2.3897847e+03 1.41e-02 1.57e+00  -3.8 1.56e+00    -  7.10e-01 2.85e-01f  1
  22  2.3875809e+03 1.30e-02 4.98e+00  -3.8 1.97e+00    -  8.70e-01 2.19e-01f  1
  23  2.3852591e+03 1.43e-02 8.37e+00  -3.8 1.17e+00    -  1.00e+00 3.18e-01f  1
  24  2.3847405e+03 1.24e-02 6.89e+00  -3.8 2.72e+00    -  1.00e+00 1.55e-01f  1
  25  2.3841607e+03 9.14e-03 8.39e+00  -3.8 2.76e+00    -  5.17e-01 3.20e-01f  1
  26  2.3838886e+03 6.33e-03 1.35e+01  -3.8 2.34e+00    -  5.83e-01 3.34e-01h  1
  27  2.3838154e+03 4.33e-03 1.17e+01  -3.8 1.99e+00    -  5.73e-01 3.32e-01h  1
  28  2.3837640e+03 2.58e-03 9.41e+00  -3.8 1.00e+00    -  9.75e-01 4.24e-01h  1
  29  2.3837400e+03 7.95e-04 2.74e+00  -3.8 6.09e-01    -  7.86e-01 7.07e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  2.3837407e+03 5.22e-05 2.00e-04  -3.8 1.71e-01    -  1.00e+00 1.00e+00h  1
  31  2.3837030e+03 9.63e-05 2.18e+00  -5.7 9.49e-01    -  7.83e-01 5.35e-01f  1
  32  2.3836811e+03 3.37e-04 5.59e-01  -5.7 2.00e+00    -  6.92e-01 7.10e-01f  1
  33  2.3836736e+03 2.13e-04 2.58e-01  -5.7 1.02e+00    -  7.61e-01 9.29e-01f  1
  34  2.3836732e+03 3.40e-05 3.14e-02  -5.7 2.27e-01    -  1.00e+00 8.54e-01f  1
  35  2.3836731e+03 1.70e-05 1.66e-02  -5.7 5.13e-02    -  1.00e+00 5.00e-01f  2
  36  2.3836731e+03 2.20e-07 2.82e-07  -5.7 2.91e-02    -  1.00e+00 1.00e+00h  1
  37  2.3836725e+03 7.11e-07 1.67e-03  -8.6 8.63e-02    -  7.90e-01 7.81e-01f  1
  38  2.3836723e+03 1.59e-06 3.15e-03  -8.6 2.05e-01    -  8.92e-01 9.72e-01f  1
  39  2.3836723e+03 2.43e-08 6.03e-08  -8.6 5.45e-03    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  2.3836723e+03 4.49e-10 1.03e-08  -8.6 1.32e-01    -  1.00e+00 1.00e+00h  1
  41  2.3836723e+03 2.74e-11 6.33e-10  -9.0 6.90e-03    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 41

                                   (scaled)                 (unscaled)
Objective...............:   2.3836722884186111e+03    2.3836722884186111e+03
Dual infeasibility......:   6.3252805544783373e-10    6.3252805544783373e-10
Constraint violation....:   2.7385524403733541e-11    2.7385524403733541e-11
Variable bound violation:   1.2994724585269068e-07    1.2994724585269068e-07
Complementarity.........:   1.2578999131652345e-09    1.2578999131652345e-09
Overall NLP error.......:   1.2578999131652345e-09    1.2578999131652345e-09

Number of objective function evaluations             = 43
Number of objective gradient evaluations             = 42
Number of equality constraint evaluations            = 43
Number of inequality constraint evaluations          = 43
Number of equality constraint Jacobian evaluations   = 42
Number of inequality constraint Jacobian evaluations = 42
Number of Lagrangian Hessian evaluations             = 41
Total seconds in IPOPT                               = 4.730

EXIT: Optimal Solution Found.
 18.420762 seconds (67.16 M allocations: 2.622 GiB, 3.35% gc time)

julia> set_optimizer(model, () -> AmplNLWriter.Optimizer(Ipopt_jll.amplexe))

julia> @time optimize!(model)
Ipopt 3.14.4: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

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

Number of nonzeros in equality constraint Jacobian...:    65471
Number of nonzeros in inequality constraint Jacobian.:    13976
Number of nonzeros in Lagrangian Hessian.............:    36439

Total number of variables............................:     8718
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     4698
                     variables with only upper bounds:        0
Total number of equality constraints.................:     8041
Total number of inequality constraints...............:     6988
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:     6988
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.2061341e+03 1.31e+01 2.52e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.3795678e+03 1.14e+01 1.33e+01  -1.0 6.12e+00    -  1.04e-02 1.31e-01h  1
   2  1.3651384e+03 1.11e+01 1.22e+01  -1.0 1.35e+01    -  6.53e-02 2.22e-02f  1
   3  1.4170592e+03 1.03e+01 8.58e+00  -1.0 1.24e+01    -  4.36e-01 7.11e-02h  1
   4  1.8085839e+03 4.98e+00 5.47e+00  -1.0 9.05e+00    -  1.90e-01 4.87e-01h  1
   5  2.0721391e+03 2.29e+00 2.62e+00  -1.0 3.45e+00    -  5.16e-01 5.30e-01h  1
   6  2.1845834e+03 1.34e+00 1.79e+00  -1.0 4.89e+00    -  3.31e-01 4.09e-01h  1
   7  2.2408122e+03 9.41e-01 1.67e+00  -1.0 2.76e+00    -  1.48e-01 2.98e-01h  1
   8  2.2606240e+03 8.02e-01 1.18e+00  -1.0 1.90e+00    -  4.36e-01 1.48e-01h  1
   9  2.2878726e+03 6.16e-01 1.06e+00  -1.0 1.50e+00    -  3.65e-01 2.32e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.3122630e+03 4.42e-01 7.08e-01  -1.0 1.09e+00    -  1.50e-01 2.82e-01h  1
  11  2.3328322e+03 2.91e-01 6.89e-01  -1.7 1.22e+00    -  2.85e-01 3.42e-01h  1
  12  2.3392720e+03 2.53e-01 4.92e-01  -1.7 2.49e+00    -  1.72e-01 1.29e-01h  1
  13  2.3434938e+03 2.28e-01 3.69e-01  -1.7 1.38e+00    -  1.90e-01 1.01e-01h  1
  14  2.3477968e+03 2.02e-01 3.51e-01  -1.7 1.07e+00    -  1.18e-01 1.12e-01h  1
  15  2.3502357e+03 1.86e-01 6.20e+00  -1.7 9.47e-01    -  7.57e-01 8.01e-02h  1
  16  2.3702492e+03 1.05e-01 1.09e+00  -1.7 1.97e+00    -  1.08e-02 4.49e-01h  1
  17  2.3784950e+03 7.88e-02 8.06e-01  -1.7 1.59e+00    -  1.00e+00 2.53e-01h  1
  18  2.4020780e+03 9.54e-03 1.33e+00  -1.7 7.65e-01    -  1.30e-01 1.00e+00h  1
  19  2.4029110e+03 9.74e-05 1.30e-02  -1.7 9.84e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  2.3928236e+03 1.22e-02 6.46e-02  -3.8 8.91e-01    -  8.21e-01 9.62e-01f  1
  21  2.3897847e+03 1.41e-02 1.57e+00  -3.8 1.56e+00    -  7.10e-01 2.85e-01f  1
  22  2.3875809e+03 1.30e-02 4.98e+00  -3.8 1.97e+00    -  8.70e-01 2.19e-01f  1
  23  2.3852591e+03 1.43e-02 8.37e+00  -3.8 1.17e+00    -  1.00e+00 3.18e-01f  1
  24  2.3847405e+03 1.24e-02 6.89e+00  -3.8 2.72e+00    -  1.00e+00 1.55e-01f  1
  25  2.3841607e+03 9.14e-03 8.39e+00  -3.8 2.76e+00    -  5.17e-01 3.20e-01f  1
  26  2.3838886e+03 6.33e-03 1.35e+01  -3.8 2.34e+00    -  5.83e-01 3.34e-01h  1
  27  2.3838154e+03 4.33e-03 1.17e+01  -3.8 1.99e+00    -  5.73e-01 3.32e-01h  1
  28  2.3837640e+03 2.58e-03 9.41e+00  -3.8 1.00e+00    -  9.75e-01 4.24e-01h  1
  29  2.3837400e+03 7.95e-04 2.74e+00  -3.8 6.09e-01    -  7.86e-01 7.07e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  2.3837407e+03 5.22e-05 2.00e-04  -3.8 1.71e-01    -  1.00e+00 1.00e+00h  1
  31  2.3837030e+03 9.63e-05 2.18e+00  -5.7 9.49e-01    -  7.83e-01 5.35e-01f  1
  32  2.3836811e+03 3.37e-04 5.59e-01  -5.7 2.00e+00    -  6.92e-01 7.10e-01f  1
  33  2.3836736e+03 2.13e-04 2.58e-01  -5.7 1.02e+00    -  7.61e-01 9.29e-01f  1
  34  2.3836732e+03 3.40e-05 3.14e-02  -5.7 2.27e-01    -  1.00e+00 8.54e-01f  1
  35  2.3836731e+03 1.70e-05 1.66e-02  -5.7 5.13e-02    -  1.00e+00 5.00e-01f  2
  36  2.3836731e+03 2.20e-07 2.82e-07  -5.7 2.91e-02    -  1.00e+00 1.00e+00h  1
  37  2.3836725e+03 7.11e-07 1.67e-03  -8.6 8.63e-02    -  7.90e-01 7.81e-01f  1
  38  2.3836723e+03 1.59e-06 3.15e-03  -8.6 2.05e-01    -  8.92e-01 9.72e-01f  1
  39  2.3836723e+03 2.43e-08 6.03e-08  -8.6 5.45e-03    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  2.3836723e+03 4.49e-10 1.03e-08  -8.6 1.32e-01    -  1.00e+00 1.00e+00h  1
  41  2.3836723e+03 2.74e-11 6.32e-10  -9.0 6.90e-03    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 41

                                   (scaled)                 (unscaled)
Objective...............:   2.3836722884185924e+03    2.3836722884185924e+03
Dual infeasibility......:   6.3249250062102252e-10    6.3249250062102252e-10
Constraint violation....:   2.7393740054115767e-11    2.7393740054115767e-11
Variable bound violation:   1.2994724585269068e-07    1.2994724585269068e-07
Complementarity.........:   1.2578999130247670e-09    1.2578999130247670e-09
Overall NLP error.......:   1.2578999130247670e-09    1.2578999130247670e-09

Number of objective function evaluations             = 43
Number of objective gradient evaluations             = 42
Number of equality constraint evaluations            = 43
Number of inequality constraint evaluations          = 43
Number of equality constraint Jacobian evaluations   = 42
Number of inequality constraint Jacobian evaluations = 42
Number of Lagrangian Hessian evaluations             = 41
Total seconds in IPOPT                               = 5.219

EXIT: Optimal Solution Found.
 10.613051 seconds (14.07 M allocations: 896.936 MiB, 2.69% gc time, 1.43% compilation time)
marco83pt commented 1 year ago

Hi Oscar, Thanks a lot for your description. Unfortunately, I cannot share the JuMP script, but I can give you an overall description of the problem. It is an energy application (AKA optimal power flow) where there are two sets of variables, complex generation from one side and complex voltages from the other side. Each variable has an upper and a lower bound. The relationship between variables is driven by the power balance for each node (or bus). Power balance is a sparse set of non-linear equality constraints. The shape is the product of voltages and trigonometric functions (because they are complex numbers) equal to power injection for each bus. That's all. What are the other options that I could test for the "_differentiation_backend"?

Best Regards, Marco

odow commented 1 year ago

Unfortunately, I cannot share the JuMP script

It's hard to offer advice without the JuMP script. Can you remove any confidential data and share a minimal reproducible example?

What are the other options that I could test for the "_differentiation_backend"?

That' mostly it, at present. The idea is that someone could write a specialized AD system for their model class and easily integrate it, but I've only written the symbolic AD as a proof of concept.

marco83pt commented 1 year ago

Hi Oscar, here I summary of what has been done at the modelling level: Example of variable:

@variable(m, vmin[b+1] <= vmag[b=0:nbus-1] <= vmax[b+1], start = 1.)
@variable(m, pmim[g+1] <= genp[g=0:ngen-1] <= pmax[g+1], start = 0.)

Example of non-linear constraint:

# Active Power Balance Constraints
@NLconstraint(m, pbal[b=0:nbus-1], -vmag[b] * sum(vmag[j] * (yr[b+1][ij] * cos(vang[b]-vang[j]) + yi[b+1][ij] * sin(vang[b]-vang[j])) for (ij, j) in enumerate(yidx[b+1])) - vmag[b]^2 * ydr[b+1] + sum(genp[g] for g in gmap[b+1]) == demp[b+1])

The objective function is has only linear and quadratic terms, then I used

@objective(m, Min, ...)

That's all.

BR/ Marco

odow commented 1 year ago

Ah okay. I think this is just a bad case for JuMP's reverse AD, because there are many constraints which are quite dense. We compute the Hessian constraint-by-constraint, and so if you have the same variable pairs appearing in multiple constraints then they will appear as two non-zeros, instead of being collated into a single entry.

Have you tried using the formulations in https://github.com/lanl-ansi/PowerModels.jl?

marco83pt commented 1 year ago

Hi Oscar, Thanks a lot for your feedback. Related to the problem modelling, this is not dense because for significant instances, the nnz are <1% compared to the matrix size. However, based on your sentences: "if you have the same variable pairs appearing in multiple constraints, then they will appear as two non-zeros" I will switch to AmplNLWriter.Optimizer(Ipopt_jll.amplexe) definitively. Then, using this way to communicate to IPOPT, is it possible to initialize (warm start) the solution properly (i.e., passing x0 and also dual and Lagrange multiplier)?

BR/ Marco

odow commented 1 year ago

Related to the problem modelling, this is not dense because for significant instances, the nnz are <1% compared to the matrix size

I guess "quite" is relative. Dense enough to cause an issue for the AD system.

. is it possible to initialize (warm start) the solution properly

You can use variable primal starts with the NL interface, but not the dual Lagrange multipliers.