# Implementation of the COPS benchmark in JuMP

**URL:** <https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555>\
**Category:** Optimization (Mathematical)\
**Created:** [June 12, 2024, 1:56pm UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555 "2024-06-12T13:56:48Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![frapac](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frapac/32/6879_2.png) [@frapac](https://discourse.julialang.org/u/frapac)\
**Post date:** [June 12, 2024, 1:56pm UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/1 "2024-06-12T13:56:48Z")

</div>

Together with @tmigot , we have re-implemented the [COPS benchmark](https://www.mcs.anl.gov/~more/cops/) in pure JuMP: [GitHub - MadNLP/COPSBenchmark.jl: Implementation of the COPS benchmark with JuMP.](https://github.com/MadNLP/COPSBenchmark.jl)

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](https://github.com/lanl-ansi/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](https://plato.asu.edu/ftp/ampl-nlp.html) 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 🙂

We hope this benchmark would be useful for the community!

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 12, 2024, 3:16pm UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/2 "2024-06-12T15:16:19Z")

</div>

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](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/querying_hessians/#Create-a-nonlinear-model) the fastest, with `MOI.Nonlinear.Evaluator`?

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

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 12, 2024, 6:56pm UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/3 "2024-06-12T18:56:52Z")

</div>

For @galle’s questions:

1. You cannot
2. Yes

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 12, 2024, 7:00pm UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/4 "2024-06-12T19:00:32Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 12, 2024, 7:05pm UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/5 "2024-06-12T19:05:04Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![tmigot](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tmigot/32/23914_2.png) [@tmigot](https://discourse.julialang.org/u/tmigot)\
**Post date:** [June 13, 2024, 12:04am UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/6 "2024-06-13T00:04:53Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![metab0t](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/metab0t/32/15551_2.png) [@metab0t](https://discourse.julialang.org/u/metab0t)\
**Post date:** [June 13, 2024, 12:10am UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/7 "2024-06-13T00:10:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 13, 2024, 6:21am UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/8 "2024-06-13T06:21:05Z")

</div>

> [@odow](#):
>
> But these functions are interpreted from the tape inside `evaluator`. They aren’t regular Julia functions like you may expect.

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

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 13, 2024, 7:15am UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/9 "2024-06-13T07:15:30Z")

</div>

There are a bunch of calls like

> <https://github.com/jump-dev/MathOptInterface.jl/blob/08d34003cd3d7500379d9440e92445bc4d9a8fa8/src/Nonlinear/operators.jl#L607-L637>

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

---

<div class="post-metadata">

**Author:** ![frapac](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frapac/32/6879_2.png) [@frapac](https://discourse.julialang.org/u/frapac)\
**Post date:** [June 13, 2024, 8:40am UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/10 "2024-06-13T08:40:31Z")

</div>

@metab0t We used COPSBenchmark.jl in a recent paper: [[2405.14236] Condensed-space methods for nonlinear programming on GPUs](https://arxiv.org/abs/2405.14236) (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.

---

<div class="post-metadata">

**Author:** ![metab0t](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/metab0t/32/15551_2.png) [@metab0t](https://discourse.julialang.org/u/metab0t)\
**Post date:** [June 13, 2024, 10:43am UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/11 "2024-06-13T10:43:48Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![frapac](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frapac/32/6879_2.png) [@frapac](https://discourse.julialang.org/u/frapac)\
**Post date:** [June 13, 2024, 11:30am UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/12 "2024-06-13T11:30:37Z")

</div>

@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](https://github.com/exanauts/ExaModelsExamples.jl/tree/main/src).

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

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

```

to instantiate the model on the GPU, or

```julia
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.

---

<div class="post-metadata">

**Author:** ![metab0t](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/metab0t/32/15551_2.png) [@metab0t](https://discourse.julialang.org/u/metab0t)\
**Post date:** [June 13, 2024, 11:55am UTC](https://discourse.julialang.org/t/implementation-of-the-cops-benchmark-in-jump/115555/13 "2024-06-13T11:55:10Z")

</div>

Great! I will look into these examples.