LP solving time doubled by adding redundant constraints

[pierre-haessig]

Hello,
I have a LP Energy system optimization model defined from JuMP. It takes about 20s to solve with the following characteristics:

A JuMP Model
Minimization problem with:
Variables: 78849
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 43803 constraints
`AffExpr`-in-`MathOptInterface.GreaterThan{Float64}`: 8786 constraints
`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 26281 constraints
`VariableRef`-in-`MathOptInterface.EqualTo{Float64}`: 8760 constraints
`VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 43809 constraints
`VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 4 constraints
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: HiGHS

Then, in the course of model development, I just found out I had missed specifying one constraint (in fact a vector constraint on an annual time series, so it makes for 8760 scalar constraints), but fortunately I know it cannot be active (it's redundant with an existing constraint). Only, the solving time jumped from 20 s to 50 s (more than ×2).

So I made some experiments by switching solvers (with the redundant constraints):

• HiGHS with "choose", "imp" or "simplex": 50s
• HiGHS with "pdlp": terrible (1700s)
• Clp: 20s (with or without the redundant constraint)

Now, I'm really not familiar with what may influence the solving time of an LP, but I really got surprised by this result (the doubling and the discrepancy with Clp) and I have no idea on how to diagnose it. In the notebook I mentioned above, the code is not easy to reproduce, but I should be able to make a standalone version if needed.

[jajhall]

If you've solved a model then added a bunch of constraints, HiGHS will solve from an advanced basis. If you don't add many constraints this should be very much more efficient than solving from scratch, but if you add a lot, then solution time could be longer than previously - when your first model is reduced by presolve and solved from scratch .

However, if your redundant constraints can be eliminated by presolve, it might be much faster to solve from scratch. If you switched to Clp at that point, then it will be running presolve.

HiGHS with "choose" will choose simplex. "IPM" is normally faster than simplex. "Pdlp" is an experimental first order solver that needs a GPU for top performance, but currently can only use CPU in HiGHS.

If you use JuMP to generate an MPS or .LP file of your problem and share it with me, I can look further at the performance of HiGHS

[pierre-haessig]

Thanks a lot for the feedback. I wasn't specific enough about one point: there is no warm start in the results I mentioned because the JuMP model is recreated fresh in each case (I'm not actually using the notebook I referenced, but a sibling notebook, which is more complex (that's why I didn't referenced it first), but where the model creation is wrapped in functions.).

about HiGHS solvers:

• Thanks for the info on "pdlp". I though it meant "primal-dual LP"... I'll not not to use it (I'm presently on my 2017 laptop with 6th gen i5 CPU and no dedicated GPU so it's indeed not a good target).
• About IMP solver, it is indeed a bit faster than simplex on my problem, but perhaps only by 10-20%.
• By the way, the only places I found info on HiGHS solvers are in the docs List of options and in the About HiGHS section of the README. Are there other places?

Now coming back to the original problem, there may indeed something about the presolve stage. I explored a bit further the problem and I tried removing another redundant constraint and then the solving time is halved (drops to 8-9 s), for both Clp and HiGHS with its "choose" solver. Also, I suspect at least one of the two constraints I'm playing with is not strictly speaking redundant but only near redundant (parallel within a 1.05 factor)

This being said, the constraints I adding/removing are quite central to the problem since it's a about energy storage operation limits. Those make a problem genuinely intertemporal (coupling between the successive 8760 instants) so I can understand they can have a huge impact on the solving time. The only thing I don't understand is that in all the cases I mentioned, the core storage constraints are always present, so from an application point of view these are just minor variants of the same problem. How could I diagnose what makes them so different from the point of view of an LP solver?

[jajhall]

It's very hard to identify why a particular collection of constraints makes such a difference to solution time. I can think of a reason why a near parallel constraint could make a difference. If you can generate a MPS/.LP model before and after eliminating the single constraint I could have a look.

Note

• HiGHS documentation is

https://ergo-code.github.io/HiGHS/stable/

• pdlp is primal-dual LP

[pierre-haessig]

Here is zipped version of four variants, in both LP and MPS format https://filesender.renater.fr/?s=download&token=da8e926e-c7d4-4a6e-b24c-b2e78ca407b7 (valid until 2024-08-03). LP files are about 8 MB while MPS are about 15MB. Variants are:

• model_09s: with less constraints, solved in 09s
• model_20s_Psto_dis and model_20s_Psto_cha: both solved in 20 s, by adding a constraint on either Psto_cha (A) or Psto_dis (B)
• model_50s_Psto_cha_dis: solved in 50 s, adding both constraints on Psto_cha (A) and Psto_dis (B)

being a bit more specific, in the LP Energy system optimization model, those are the relevant lines. The constraints I'm switching on and off are the one marked "A" and "B"

# Charge and discharge power
@variable(model, Psto_cha[1:K] >= 0)
@variable(model, Psto_dis[1:K] >= 0)
# Energy
@variable(model, Esto[1:K+1])
@constraint(model, Esto .<= energy_rated_sto) # SoCmax = 1 implicitly
@constraint(model, Esto .>= mg.storage.SoC_min*energy_rated_sto) # SoC_min is often 0
# Power constraints
@constraint(model, Psto_cha .<= mg.storage.charge_rate    * energy_rated_sto) # A
@constraint(model, Psto_dis .<= mg.storage.discharge_rate * energy_rated_sto) # B
# Evolution of the State of Energy, with piecewise linear in P loss term (aka constant efficiency model)
a = mg.storage.loss_factor # = 0.05
@constraint(model,
    stodyn[k=1:K],
    Esto[k+1] == Esto[k] + (Psto_cha[k] - Psto_dis[k] - a*(Psto_cha[k]+Psto_dis[k]))*dt
    )

the near redundancy comes from the fact that:

• mg.storage.(dis)charge_rate is equal to dt
• and the loss factor a is << 1.

More specifically, it is the combination of:

• $0 \leq E(k) \leq E_{rated}$ for all k (in the code: Esto .<= energy_rated_sto)
• $E(k+1) =E(k)+ (P_{cha}(k) - P_{dis}(k) - a*(P_{cha}(k) +P_{dis}(k)))$

which is almost included in both $P_{cha}(k) \leq E_{rated}$ (A) and $P_{dis}(k) \leq E_{rated}$ (B). Still, those last two constraints represent 2*8760 elements. It would be equivalent in the case $a=0$

Note: my question about the doc was not about the doc url but whether inside the doc the solver choices are documented elsewhere than in the repo README and in the List of options

[pierre-haessig]

OK, another iteration: I removed the small penalty I'm adding in the objective to emulate a lexicographic biobjective problem and now the solving time drops to 16s. This penalty probably makes the problem (a bit?) numerically ill-conditioned.

For comparison, I also ran the problem without this penalty and without the constraints on Psto_cha (A) and Psto_dis (B) and it gets solved in 9 s, so in this case the penalty has no effect.

I've created a new zip file with all models (now with 6 variants...) https://filesender.renater.fr/?s=download&token=368d7f4d-d825-461c-ac8d-f6999d8630cc

[jajhall]

The small penalty will not make the problem numerically difficult to solve. However, if it's the same order of magnitude as the cost perturbation used to alleviate dual degeneracy, it may trigger some extra iterations when the perturbation is removed.

I can't say more until I run your models, and I can't do that this week as I'm travelling.

[pierre-haessig]

OK, I suggest that perhaps you don't waste of time on this, because it seems to depends on many peculiarities. E.g. this morning the running time of 50 s had dropped to 16s. After some digging, it happens that the change was a package update:

• 50 s solving time is obtained when using v1.30.0 of MathOptInterface, which was in my software setup yesterday
• 16s with MathOptInterface v1.31.0

Something must have happened in that release, although nothing struck me in the release notes.

Other relevant package versions:

• HiGHS.jl v1.9.0 or v1.9.1, corresponding to HiGHS binary v1.7.0 or v1.7.1 (no effect from this upgrade)
• JuMP v1.22.2 in all cases

Still, for the reference, I've appended the MPS and LP exports of the problems saved with MOI 1.31. Here is the link (valid until 2024-08-04). The two sibling problems are:

• model_50s_Psto_cha_dis: exported from MOI 1.30, takes 50s (model mentioned above: with power constraints A and B + a small penalty in the objective)
• model_50s-15s_Psto_cha_dis_MOI-1.31: same JuMP problem definition, exported from MOI 1.31, takes 16s

I haven't tried to load and save any of these MPS/LP files, since I'm not familiar. Only, I've checked that their checksums are different!