Implementation of the COPS benchmark in JuMP

Specific Domains Optimization (Mathematical)
1

Together with @tmigot , we have re-implemented the COPS benchmark in pure JuMP: GitHub - MadNLP/COPSBenchmark.jl: Implementation of the COPS benchmark with JuMP.

The COPS benchmark is a collection of challenging nonlinear programs. The dimension of each instance can be parameterized, meaning we can generate very large-scale nonlinear instances.

On the contrary to the OPF benchmark used in rosetta opf, the instances are coming from various domains: PDE-constrained optimization, optimal control, identification of parameters. As a result, they are a good way to test the robustness of a nonlinear modeler or a nonlinear solver.

You can find attached a preliminary benchmark comparing the performance we obtain with AMPL (using AmplNLWriter.jl) and with bare JuMP. The instances are also used in the AMPL-NLP benchmark generated by Hans Mittelmann.

In both cases, we use the solver Ipopt with HSL MA27 for the linear solver. On this benchmark, most of the running time is spent in the linear solver.

Instance #vars #cons Ampl (total) JuMP (total) JuMP (AD)
bearing_160000 161604 1608 12.0 6.2 -
camshape_6400 6400 12803 1.6 0.8 -
dirichlet_120 64783 11143 499.6 508.0 10.1
elec_400 1200 400 98.5 1811.4 1712.1
gasoil_3200 83203 83200 6.2 6.5 0.9
henon_120 48557 16397 1354.7 1143.5 26.8
lane_120 66491 9011 556.2 490.8 11.4
marine_1600 51215 51192 1.0 1.0 -
pinene_3200 160005 160000 6.0 6.6 -
robot_1600 14410 9612 0.7 0.9 0.3
rocket_12800 51205 38404 13.9 8.5 1.6
steering_12800 64006 51208 1.7 1.7 0.6

A few notes:

  • The performance of AMPL and JuMP are relatively similar, except on elec_400, a dense instance.
  • Despite solving the same instances, we can obtain different convergence patterns between AMPL and JuMP. Indeed, we observe very large primal-dual regularization within Ipopt (lg(rg)), leading to difference in the floating point arithmetic inside the linear solver (here HSL MA57).
  • MathOptSymbolicAD crashes with a segfault on the PDE-constrained instances (dirichlet, henon, lane). With the large number of nonlinear terms, I think we are pushing Julia’s compiler to its limit.
  • This package is provided just to help people developing sparse AD backend and optimization solvers in Julia. We do not want to use it as a new benchmark for different optimization solvers. We believe the Mittelmann benchmark already fulfills that purpose :slight_smile:

We hope this benchmark would be useful for the community!

9 Likes
2

That is awesome, thank you!

As you may know, I’m currently trying to develop a generic sparse AD framework in Julia, together with @hill. Our collection of packages SparseConnectivityTracer.jl + SparseMatrixColorings.jl + DifferentiationInterface.jl is starting to look very promising, and we would like to challenge it.

Suppose I am given a JuMP Model like those in the COPS suite, and I want to benchmark my methods against JuMP in the fairest way possible:

  • How can I convert the JuMP Model to pure Julia functions that return the objective value and the constraint values, such that the resulting functions are as efficient as possible?
  • To compute the Hessian of the Lagrangian (or just its sparsity pattern), is the method outlined here the fastest, with MOI.Nonlinear.Evaluator?

Since it’s a rather generic JuMP question, perhaps @odow can help

2 Likes
3

For @galle’s questions:

  1. You cannot
  2. Yes
1 Like
4

Okay so if I have a sparse AD approach that is not based on an algebraic modeling language, the only option I have to compare it with JuMP is to have 2 separate versions of the benchmark problems, a pure Julia one and a JuMP model?

5

In the Hessian tutorial, you can evaluate the objective and constraint with:

MOI.eval_objective(evaluator, x)
MOI.eval_constraint(evaluator, g, x)

But these functions are interpreted from the tape inside evaluator. They aren’t regular Julia functions like you may expect.

1 Like
6

@gdalle These problems will also be in OptimizationProblems.jl (50% already are) with both JuMP models and Julia function models (in ADNLPModel format, but it is trivial to get the functions back).

2 Likes
7

How do Examodels.jl perform on these problems compared with JuMP and AMPL (on CPU and GPU)?

8

How do they differ from a “regular” Julia function?

9

There are a bunch of calls like

We don’t compile a Julia function from the expression.

10

@metab0t We used COPSBenchmark.jl in a recent paper: [2405.14236] Condensed-space methods for nonlinear programming on GPUs (Figure 5 and Table 4, page 28).

We had to use a subset of the COPS benchmark, as ExaModels is not able to parse the large-scale PDE-constrained problems. It suffers from the same caveat as MathOptSymbolicAD.jl: if the problem has long nonlinear expressions, we put too much burden on Julia’s compiler.

2 Likes
11

@frapac
Thanks! I have read the inspiring preprint. Is the code to replicate result of Examodels.jl in this paper open source?

12

@metab0t Thank you for the kind words! The code is open-source, indeed. A subset of the COPS instances have been converted to ExaModels here.

In the paper, we used the JuMP instances directly, and converted them to ExaModels using e.g.,

model = COPSBenchmark.rocket_model(12800)
nlp = ExaModels.ExaModel(model; backend=CUDABackend())

to instantiate the model on the GPU, or

nlp = ExaModels.ExaModel(model)

to instantiate the model on the CPU with ExaModels.

Compared to the raw ExaModels instances, there is a slight overhead when we convert the JuMP model to ExaModels.

1 Like
13

Great! I will look into these examples.

1 Like