Avoid re-building models

[klamike]

In sampler.jl, we build a new JuMP model for each formulation for each seed. Instead, we should build one model per formulation and use POI to change the loads for each seed. This should reduce build times substantially after the first one.

In the package code, this will also require updating all of the formulations to make the loads parameters, wrapping the optimizer with POI.Optimizer, and adding ParametricOptInterface to the dependencies. Since we add the dependency, the code for changing the parameter values can go in the package as well.

In the SLURM pipeline, this would still build one model per formulation per minibatch. There may be a way to create each formulation's model once in the make_sysimage or ref job, then reuse it in each minibatch.

[mtanneau]

I fully agree that we should improve the current data-generation pipeline, but I would like to have a deeper discussion about this before we merge such changes.

A few notes:

• JuMP currently supports parameters (see, e.g., https://jump.dev/JuMP.jl/stable/manual/variables/#Parameters). I am not 100% sure how that differs from POI

• I would also like to understand how things work under the hood, namely, what problem ends up being presented to the solver. Specifically, I'm referring to this part of the JuMP docs:

  If the solver supports the MOI.Parameter set, it may decide to replace all instances of the parameter variable by the associated constant. If the solver does not support parameters, it will add the parameter as a decision variable with fixed bounds.

  This is particularly relevant for Ipopt, because it does not have a presolve. Therefore, if parameters are passed as fixed variables (i.e., same lower and upper bound), Ipopt may struggle (because there's no strictly feasible solution anymore).

• In a similar vein (same section of the docs):

  JuMP treats a parameter multiplied by a decision variable as a quadratic expression, even though it is equivalent to a linear expression.

  Not an immediate issue if we only modify right-hand sides, but it might be more problematic if/when we start considering changes in topology (namely removing branches). Granted, it's far-sighting, but I want to make sure we are as future-proof on that front.

[mtanneau]

An intermediate step (which would probably need to happen anyway) would be to define an interface for updating an existing OPF model, given new input data. We can have an agnostic interface like, e.g., update!(opf::OPFModel, data), the internals of which we can always modify in the future

[klamike]

Agreed. I'll see if I can get an update! working without POI...

[mtanneau]

For the loads, we should be able to use something like set_normalized_rhs

[klamike]

With the normalization, I think it would be cleaner to delete+unregister the old constraints and add new ones

(though the docs do explicitly recommend against doing this...)

[mtanneau]

Normalization?

Is that related to the fact that, the active/reactive power balance constraints will have -pd / -qd on the right-hand side?

For instance, on the 14-ieee system, bus #2 has a load with +0.217p.u. but... in the ACPPowerModel problem the power balance constraint at that node will have -0.217 as RHS

julia> using PGLib, Ipopt, PowerModels, JuMP
julia> data = make_basic_network(pglib("14_ieee"));
julia> data["load"]["1"]["load_bus"]
2
julia> data["load"]["1"]["pd"]
0.217
julia> opf = OPFGenerator.build_opf(ACPPowerModel, data, Ipopt.Optimizer);
julia> opf.model[:kirchhoff_active][2];
kirchhoff_active[2] : -pg[2] + pf[(5, 2, 5)] + pf[(3, 2, 3)] + pf[(4, 2, 4)] + pf[(1, 2, 1)] = -0.217
julia> normalized_rhs(opf.model[:kirchhoff_active][2])
-0.217

This is because the original code for active/reactive power balance https://github.com/AI4OPT/OPFGenerator/blob/f7364bc04d2d075ac4999ff1025b9dc4a38a885f/src/opf/acp.jl#L57-L71 has variables on both sides of the equation. JuMP automatically moves all the variables to the same side and all the constants to the other side... but we don't fully control which side is left and which is right.

We can fix this by writing all the variables on the same side of the equation, when declaring those constraints. I've been hesitant to do that, because it would effectively negate the corresponding dual variable, and I didn't want to mess up our existing datasets.

BTW, we have the same issue with DC and SOC. Easy fix, I just want to make sure our code / formulation / data / etc. are all in sync.