Run with open-energy-modeling-benchmarks

[odow]

Our claim, and that of @Opt-Mucca, is that the energy modeling benchmarks are "bad" (numerically).

https://github.com/jump-dev/open-energy-modeling-benchmarks/tree/main/instances

We should quantify this.

The resolution for this issue is a report that summaries the set of issues that MathOptAnalyzer can find in these instances.

[Opt-Mucca]

Here's a list of ways an instance is marked with the numerics flag in MIPLIB (may be helpful):

• A single row of the matrix has max / min coefficient ratio of 10**6
• A single column of the matrix has max / min coefficient ratio of 10**6
• The matrix has max / min coefficient ratio of 10**8
• Two solvers disagree on the optimal solution value by more than some tolerance

[odow]

With no data other than general vibes: I'd expect most energy models to fail most of those checks...

I think our pitch to BT is. Why yes. These models suck. But if we provide better tooling, modelers may be able to write better models, and therefore the models may be easier to solve. So in addition to improving HiGHS, we can help make all models better and bring modelers closer to the solver.

[rschwarz]

Here's a list of ways an instance is marked with the numerics flag in MIPLIB (may be helpful):

Sorry for hijacking this discussion, but I have a question here: is this analysis done before or after presolving is applied? I remember that MIPLIB does perform some limited preprocessing on all submitted, right?

The question seems relevant because sometimes a modeler can "cheat" by, eg, splitting long power lines into smaller segments, that each have smaller coefficients in power loss, but that a solver might aggregate again by eliminating variables and equations. So, coefficient ranges before and after presolving can differ dramatically.

[Opt-Mucca]

Sorry for hijacking this discussion, but I have a question here: is this analysis done before or after presolving is applied? I remember that MIPLIB does perform some limited preprocessing on all submitted, right?

You're remembering correctly (there's some limited presolving applied). Such presolve should be applied, although I think it will not do much for the energy instances if the generator is decent. It's more there to catch modelling mistakes and truly basic potential simplifications (it will not catch many "cheats").
There definitely are modelling techniques to disguise a troublesome instance essentially (even if not explicitly intended by the modeller themself). My recommendation is not focus on that too much because that problem will always exist to some degree and assuming bad faith on the modellers part is pointless if you're building a benchmark from those instances.

[rschwarz]

Thanks for the explanation.

My recommendation is not focus on that too much because that problem will always exist to some degree and assuming bad faith on the modellers part is pointless if you're building a benchmark from those instances.

Yes, I was not thinking of modellers trying to cheat a benchmark, but rather cheating themselves, because they think (according to the metrics above) they are improving the formulation numerically by applying these changes (been there).

[odow]

I think a related key point is: we're building a benchmark from these because they're the models people write. Not because they're good benchmark instances from a MIPLIB perspective.

People do trick themselves. GenX want to apply their own prescaling to make the warnings go away and make the solve go faster. Even though with post scaling the results might violate their expected tolerances. In some cases I don't know what the right fix is. People except more precision from a solver than it can provide.

[odow]

and truly basic potential simplifications

@jd-lara and I talked about some of these today. Solar power plants being added to a model at midnight with everything fixed to 0. People write these models, and I don't know we should be benchmarking them after presolve. Or at least. It's a different motivation. Both are useful.

[Opt-Mucca]

I think a related key point is: we're building a benchmark from these because they're the models people write. Not because they're good benchmark instances from a MIPLIB perspective.

Strong agree.

Another point to add: The initial presolve for MIPLIB is not important for improving the numerics, but for creating a feature presentation of the instance (which is then used for numeric flags).

Another point: I think people can hack their instance to be numerically friendly from the solver's POV (like discretising pipe segments). It's not the solver's problem so to speak that the returned solution is not valid to some larger framework it's not aware of.

[jd-lara]

There are some "realities" in energy modeling that make it harder to write good models out of the box. Ee incur in some of these in Sienna too because they are hard to catch before build times and are very system and scenario dependent. I'll listen them for the sake of making the point that the effort to improve solve times in MILP solvers as it applies to energy systems has to inescapably consider the realities of energy modeling.

1. Setting penalties on slacks. There isn't a clear way to define a slack penalty before the model is built, probably the only way is to build the cost function without the slack, determine the range and then add the penalty. Making this a behind the scenes model build feature is not trivial at all. See this improve the addition of slack penalties NREL-Sienna/PowerSimulations.jl#1336.

2. Big M constraints. These are used a lot even when modelers don't know they are using them explicitly. I think there is a HUGE literature (a lot of it unhelpful) on how set big M's. We have some of these since for most solvers is not possible to fix the binaries within an SOS constraint to get an LP relaxed model.

3. Difficult to catch scaled data issues for the b vector. This is particularly bad when you combine regionally diverse models and have situation when a time series is close to 0.0 but not 0.0. For instance 0.1 MW solar power output (in the morning or tail end of the day) in 100 MVA Base per unit is 1e-4 and a huge aggregate plant 20 GW plant that turns into 2e5 on the same base. I have tried my best to address this and it is extremely hard. Something similar happens for lines when they can span subsystems at 69 kV with 2-3 MW capacity and high voltage 750 KV subsystems with 1 GW capacities.

This can partially be solved with configurable implementations of ParameterOptInterface.jl if we can define a threshold for 0.0 to make a parameter on the RHS 0.0 in parametric optimization cases for renewable but it isn't a bulletproof solution.

[jd-lara]

Here's a list of ways an instance is marked with the numerics flag in MIPLIB (may be helpful):

• A single row of the matrix has max / min coefficient ratio of 10**6
• A single column of the matrix has max / min coefficient ratio of 10**6
• The matrix has max / min coefficient ratio of 10**8
• Two solvers disagree on the optimal solution value by more than some tolerance

some additional ones that can reflect the quality of the model could be

1. Constraint matrix sparsity/density.

I'd suggest adding some other metric that can only be calculated post-solve

1. Integrality gap: difference between LP relaxation and optimal MILP solution.
2. Tightness of constraints. If we can estimate how often all the constraints in the model are binding it will reflect that the models can also be poorly defined and ineffective.

[odow]

I'll listen them for the sake of making the point that the effort to improve solve times in MILP solvers as it applies to energy systems has to inescapably consider the realities of energy modeling.

There's an inescapable tension. The question is where on the trade-off curve do we want to sit.

We have some of these since for most solvers is not possible to fix the binaries within an SOS constraint to get an LP relaxed model.

Okaaaay. So this is an interesting one. We could add this to JuMP's relax_integrality.

What if you added an SOS constraint but used a solver like HiGHS which does not support them natively?

[joaquimg]

Here's a list of ways an instance is marked with the numerics flag in MIPLIB

Thanks @Opt-Mucca !

Is this list exhaustive? I will plan to add these features.

Question:
Is there a version of such a list for quadratic objectives or constraints?

[joaquimg]

About @jd-lara points,

Setting penalties on slacks

Extremely complex problem. Our focus (at least for now) should be on pointing when they are bad and educating people through tutorials and blog posts. Proposing a solution to the user automatically is way to complex.

Big M constraints

SOS is fun. However, if users can provide good big-M's, the solve times tend to be much much better. The main issue is SOS does help with root node. Thats why people in bilevel have been using: SOS + valid constraints: https://link.springer.com/article/10.1007/s10287-023-00435-5 Note that the issue with bilevel is that big M's can be as hard to find as solving the problem. This is not the case for many energy systems, where food big Ms are known.

Difficult to catch scaled data issues for the b vector

Yeah, this is important. Zeros are less of an issue from what I have seen. However, people use random number gerators that return 1e-6 for some solar and wind plant. Setting those to zero usually improves solve time, solution quality and do not change the objective. I would not want JuMP to do it automatically, though. Again we should report, and point to tutorials etc.

Constraint matrix sparsity/density

We have this per constraint in this package already. I would like to hear @Opt-Mucca take on us presenting this info also and what would a reasonable threshold for a warning be.

Integrality gap

I am not sure how to use that info to guide a user. Yes, small gap is great, but hard to tell the user how to "fix the issue". Also, this could be more of a recommendation as this analysis would require solving the problem.
From the perspective of presenting an analysis about instances, this is great, but for the automatic analysers in this package, I don't know.

Tightness of constraints

This could go into the Feasibility analyzer. However, I am not sure how to differentiate issues from expected stuff and report suggestions to the user. A simple example: bounds on a variable. On of the bounds will typically be tight while the other will be veeeeery loose.

[joaquimg]

So this is an interesting one. We could add this to JuMP's relax_integrality.

I think he meant he needed https://jump.dev/JuMP.jl/stable/api/JuMP/#fix_discrete_variables

[Opt-Mucca]

Is this list exhaustive? I will plan to add these features.

The objective dynamism should also be no larger than 10**6
We also would flag anything with a coefficient larger than 10**10 (My personal opinion is that this is not needed if you have the matrix dynamism though)
There's also the basis condition number of the coefficient matrix, but we did not incorporate this as not all solvers support such functionality or agree that it's an important flag (Personal experience points to it being a good indicator for numerical issues). See kappa in Xpress.
The solvers disagreeing is going to be the hardest to compute, and is straight up incorrect if the solver simply has a bug. Row / col / matrix / obj dynamism together will eliminate the majority of bad actors

Constraint matrix sparsity/density

That's a more difficult one..... I don't think the density is the inherent reason behind the numeric difficulties, e.g. dense 0-1 problems of normal sizes shouldn't have issues. Rather dense rows / columns exacerbate numeric issues, and especially so for larger problems, where I think density can get more dangerous. I would not recommend using density by itself as a flag for numerical issues (and I'd have no clue how to make it a conditional flag). I would maybe recommend adding a "consider another formulation flag" though as dense rows / columns also can often slow everything down even if they're not causing numerical problems.

Big M constraints. These are used a lot even when modelers don't know they are using them explicitly. I think there is a HUGE literature (a lot of it unhelpful) on how set big M's. We have some of these since for most solvers is not possible to fix the binaries within an SOS constraint to get an LP relaxed model.

What does this mean exactly? I imagine if your fixes violate the SOS constraint then the MILP solver would return infeasible instantly? The "MILP" solve after fixing should also be as quick as an LP solve (although I guess that limits your range of solvers?)