[Bridges] Add FixParametricVariablesBridge #2092
Conversation
|
This is mostly my fault, but the format check failed: https://github.com/jump-dev/MathOptInterface.jl/actions/runs/4136044676/jobs/7151476121. You can either follow the changes in the link, or do cd("~/.julia/dev/MathOptInterface") # or whatever the path is
using JuliaFormatter
format("src")
format("test")I don't know if I fully understand the benchmarks, or the drop in performance as T grows. I need to take a deeper look. cc @joaquimg @guilhermebodin @jd-lara who will all be interested in this. |
I'm not sure I understand this part. Why not? Was an error thrown? One thing the bridge can do is add support for MathOptInterface.jl/src/Bridges/Constraint/bridges/flip_sign.jl Lines 57 to 68 in fbf833c |
|
Formatting should be correct now.
I simplified the test model to what is actually necessary, and did additional tests. Full code at the end of this comment. New benchmark below
It "corrupted" the model resulting either in an infeasible model or one without a proper objective. I'm not sure - could be due to the way I was modifying the objective, ... . I've left relaxation out for now - will test that again + look at the BenchmarkThoughts on what to benchmark: There are two different timings we care about; that of the first iteration, and those of all subsequent iterations (the third should already be equal to the second). While it may be a relevant advantage in a real case scenario, that subsequent solves are faster when using a parametric approach (compared to throwing away the first iteration's solution; since the solver could use some information for a warm start / ...), this is not a relevant timing to be measured (since it would essentially be the same for every parametric / set_coefficient / set warm start / ... approach). Therefore I have excluded solve times from now on (by using The model looks like Here Performance seems to be: bad for the first iteration (standard and parametric seem to be
Timings (in milliseconds) are:
However, there is still a "breakpoint" iteration after which the parametric version is faster (regarding the overall time of first iteration + all subsequent iterations; except for Benchmark CodeExpandDependencies ( |
|
I have now done a "soft" test, that does not assume that each constraint contains a (bridged) parameter. The following constraints are added to the model (just to add some "arbitrary" non-parametric effort): This adds 3 new constraints (the last 2 are just a single one), which means that now 25% (1 out of 4) constraints contain a multiplicative parameter and need the bridge. The timing plot looks like:
Now, at least up to However, I'm still unsure about the first iteration performance - for small models it's probably fine, but for larger ones (and for some use cases the |
|
I think I need to profile this properly to see where the issue is. (If you haven't used it before, https://github.com/timholy/ProfileView.jl is the way to go.) There's a few potential causes:
We should also think about our aims:
|
I haven't profiled it yet - but I'll make sure to do that on the weekend!
I'll add that. On a similar note: Is there an easy way to detect if
Is there a way to prevent the other bridges, while still loading the parameter one? Using
At least for me, yes. While rebuilding the model is probably the easiest thing to do most of the time, it wastes an unnecessary amount of time. Therefore the total time over some number of iterations should be considerably better using the parametric approach (otherwise people will stick to rebuilding) - while also making sure iteration 1 is actually doable (e.g. not "exploding" the memory usage too much)
Yes. |
|
I don't really understand the benchmark results. Most of the time spent is in HiGHS calling (These numbers include a small change to cache the dictionaries so that we don't rebuild them every time.) Codeusing JuMP
import HiGHS
function build_and_pass(N::Int64, opt; manual::Bool)
model = Model(opt)
set_silent(model)
set_string_names_on_creation(model, false)
add_bridge(model, MOI.Bridges.Constraint.FixParametricVariablesBridge)
@variable(model, x[1:N])
@variable(model, p[1:N])
fix.(p, 1.5; force = true)
if manual
@constraint(model, c, x .>= 1)
set_normalized_coefficient.(c, x, 1.5)
else
@constraint(model, c, p .* x .>= 1)
end
@objective(model, Min, sum(x))
@time optimize!(model)
fix.(p, 2.0; force = true)
if manual
set_normalized_coefficient.(c, x, 2.0)
end
@time optimize!(model)
return
end
import ProfileView
ProfileView.@profview build_and_pass(50_000, HiGHS.Optimizer; manual = false)
ProfileView.@profview build_and_pass(50_000, HiGHS.Optimizer; manual = true)Outputjulia> ProfileView.@profview build_and_pass(50_000, HiGHS.Optimizer; manual = false)
1.583368 seconds (4.05 M allocations: 193.149 MiB, 5.90% gc time)
0.084514 seconds (700.01 k allocations: 10.681 MiB)
julia> ProfileView.@profview build_and_pass(50_000, HiGHS.Optimizer; manual = true)
0.540596 seconds (1.75 M allocations: 89.918 MiB, 7.51% gc time)
0.024224 seconds (2 allocations: 32 bytes)manual = false
manual = true
|
|
Ah. Once cause could be that the initial pass doesn't insert values into the constraint matrix. Then modifying the coefficient value later changes from 0 to 1, requiring a shift in the sparsity pattern of the constraint matrix. We could work around this by setting the affine part for all variables that appear. Let me have a go at a few modifications. |
|
I made some initial changes here: sstroemer#1. They haven't changed the runtime performance much, but I did catch a bug handling |
|
I swear those changes didn't make a difference. But I can back after a break, and found: using JuMP
import HiGHS
function build_and_pass(N::Int64, opt; manual::Bool)
@info "Build model"
@time begin
model = Model(opt)
set_silent(model)
set_string_names_on_creation(model, false)
add_bridge(model, MOI.Bridges.Constraint.FixParametricVariablesBridge)
@variable(model, x[1:N])
@variable(model, p[1:N])
fix.(p, 1.5; force = true)
if manual
@constraint(model, c, x .>= 1)
set_normalized_coefficient.(c, x, 1.5)
else
@constraint(model, c, x .>= 1)
end
@objective(model, Min, sum(x))
end
@info "First optimize"
@time optimize!(model)
@info "Update"
@time if manual
fix.(p, 2.0; force = true)
set_normalized_coefficient.(c, x, 2.0)
else
fix.(p, 2.0; force = true)
end
@info "Re-optimize"
@time optimize!(model)
@info "Total"
return model
end
@time build_and_pass(50_000, HiGHS.Optimizer; manual = false);
@time build_and_pass(50_000, HiGHS.Optimizer; manual = true);julia> @time build_and_pass(50_000, HiGHS.Optimizer; manual = false);
[ Info: Build model
0.135477 seconds (2.23 M allocations: 106.356 MiB, 37.81% gc time)
[ Info: First optimize
0.435786 seconds (1.75 M allocations: 89.918 MiB)
[ Info: Update
0.084833 seconds (300.00 k allocations: 5.341 MiB)
[ Info: Re-optimize
0.011754 seconds (2 allocations: 32 bytes)
[ Info: Total
0.669441 seconds (4.28 M allocations: 201.639 MiB, 7.65% gc time)
julia> @time build_and_pass(50_000, HiGHS.Optimizer; manual = true);
[ Info: Build model
0.092192 seconds (2.33 M allocations: 108.645 MiB)
[ Info: First optimize
0.696616 seconds (1.75 M allocations: 89.918 MiB, 33.52% gc time)
[ Info: Update
0.128399 seconds (500.00 k allocations: 9.155 MiB)
[ Info: Re-optimize
0.022884 seconds (2 allocations: 32 bytes)
[ Info: Total
0.941419 seconds (4.58 M allocations: 207.742 MiB, 24.81% gc time)
julia> @time build_and_pass(50_000, HiGHS.Optimizer; manual = false);
[ Info: Build model
0.275107 seconds (2.23 M allocations: 106.356 MiB, 68.17% gc time)
[ Info: First optimize
0.433582 seconds (1.75 M allocations: 89.918 MiB)
[ Info: Update
0.084287 seconds (300.00 k allocations: 5.341 MiB)
[ Info: Re-optimize
0.011777 seconds (2 allocations: 32 bytes)
[ Info: Total
0.806095 seconds (4.28 M allocations: 201.639 MiB, 23.27% gc time)
julia> @time build_and_pass(50_000, HiGHS.Optimizer; manual = true);
[ Info: Build model
0.117296 seconds (2.33 M allocations: 108.645 MiB, 27.67% gc time)
[ Info: First optimize
0.463744 seconds (1.75 M allocations: 89.918 MiB, 5.65% gc time)
[ Info: Update
0.119186 seconds (500.00 k allocations: 9.155 MiB)
[ Info: Re-optimize
0.022526 seconds (2 allocations: 32 bytes)
[ Info: Total
0.724041 seconds (4.58 M allocations: 207.742 MiB, 8.10% gc time)If I read those results correctly, there's very little difference between |
Thanks for the changes! I looked at your first profiling benchmark and the second with the timings - there are differences in the Note: I have kept the using JuMP
import HiGHS
function build_and_pass(N::Int64, opt; manual::Bool)
@info "Build model"
@time begin
model = Model(opt)
set_silent(model)
set_string_names_on_creation(model, false)
add_bridge(model, MOI.Bridges.Constraint.FixParametricVariablesBridge)
@variable(model, x[1:N])
# @variable(model, p[1:N]) # moved into the branch `manual=false`
# fix.(p, 1.5; force = true) # moved into the branch `manual=false`
if manual
@constraint(model, c, x .>= 1)
set_normalized_coefficient.(c, x, 1.5)
else
@variable(model, p[1:N]) # <--
fix.(p, 1.5; force = true) # <--
@constraint(model, c, p .* x .>= 1) # added `p .*`
end
@objective(model, Min, sum(x))
end
@info "First optimize"
@time optimize!(model)
@info "Update"
@time if manual
# fix.(p, 2.0; force = true) # no `p` in `manual=true`
set_normalized_coefficient.(c, x, 2.0)
else
fix.(p, 2.0; force = true)
end
@info "Re-optimize"
@time optimize!(model)
@info "Total"
return model
end
GC.gc()
@time build_and_pass(50_000, HiGHS.Optimizer; manual = false);
GC.gc()
@time build_and_pass(50_000, HiGHS.Optimizer; manual = true);With that, we are still left with a sizeable difference: julia> @time build_and_pass(50_000, HiGHS.Optimizer; manual = false);
[ Info: Build model
0.168728 seconds (2.05 M allocations: 163.493 MiB, 35.16% gc time)
[ Info: First optimize
0.885333 seconds (1.95 M allocations: 146.583 MiB, 4.39% gc time)
[ Info: Update
0.055777 seconds (150.00 k allocations: 3.052 MiB)
[ Info: Re-optimize
0.092662 seconds (250.01 k allocations: 3.815 MiB)
[ Info: Total
1.203645 seconds (4.40 M allocations: 316.973 MiB, 8.16% gc time)
julia> @time build_and_pass(50_000, HiGHS.Optimizer; manual = true);
[ Info: Build model
0.072086 seconds (901.52 k allocations: 68.443 MiB)
[ Info: First optimize
0.156815 seconds (400.20 k allocations: 39.966 MiB)
[ Info: Update
0.042946 seconds (150.00 k allocations: 3.052 MiB)
[ Info: Re-optimize
0.020261 seconds (2 allocations: 32 bytes)
[ Info: Total
0.293478 seconds (1.45 M allocations: 111.491 MiB) |
| bridge::FixParametricVariablesBridge{T,S}, | ||
| chg::MOI.ScalarCoefficientChange{T}, | ||
| ) where {T,S} | ||
| MOI.modify(model, bridge.affine_constraint, chg) |
There was a problem hiding this comment.
I don't this that this would work, the coefficient change should change the coefficient before we added the quadratic part. If final_touch was already called then we should add the coefficient here
There was a problem hiding this comment.
Won't final touch need to be called after this anyway?
There are two parts:
- We modify the existing
.affine_constraint. This catches adding new variables, or modifying a variable that doesn't exist in the quadratic part - We modify
bridge.f, so that the next time we callfinal_touch, it gets updated correctly.
| s::S, | ||
| ) where {T,S<:MOI.AbstractScalarSet} | ||
| affine = MOI.ScalarAffineFunction(f.affine_terms, f.constant) | ||
| ci = MOI.add_constraint(model, affine, s) |
There was a problem hiding this comment.
Wouldn't it be more efficient to not add any constraint here and add it in final_touch when it is already modified ?
There was a problem hiding this comment.
You can leave a TODO to keep this PR simple
There was a problem hiding this comment.
Then we'd need a Union{Nothing,ConstraintIndex} in the type. I guess it's a trade-off that might improve the first solve time, but it wouldn't improve the re-solve time (which is more important for this bridge).
Ah. Sure. Now I remember what I intended to do.
I think that's bad news, because it means we can't close the gap unless we find a way to not add the The good news is that I think it means this PR can be merged once we fix @blegat's comments. At a minimum, it simplifies the syntax for users, and they can evaluate whether it is worth the performance trade-off. There's also a route forward to writing a better bridge if we decide to add explicit parameters to MOI. |
I'll take care of working in the comments. But I'm still trying to dig into it a little bit more, to fully understand whats going on behind the scenes. Question 1: Why does the following example result in one call to using JuMP
import HiGHS
model = Model(HiGHS.Optimizer)
add_bridge(model, MOI.Bridges.Constraint.FixParametricVariablesBridge)
@variable(model, x)
@variable(model, p); fix(p, 1.0; force=true)
@constraint(model, p * x >= 1)
@constraint(model, x >= 0)
@objective(model, Min, x + 0)
optimize!(model)Question 2: How much does the ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}}(x.value)
MOI.get(model, MOI.ConstraintSet(), ci).valuelookup cost? Could it be benefitial to replace that by a dictionary lookup for the current parameter value? |
It should be minimal, but I guess it makes a difference. Note that we have to do this lookup at least every time we call |
odow
changed the title
The qualitative evaluation of good and bad news also depends on the problem size. What we have learned with really large models are: 1) the performance is solver dependent and now the presolve works, HiGHS isn't the best to assess the trade off. 2) with recurrent solves not rebuilding the model outweighs the additional variables except when the solver can't support more than 1e9 variables. For left hand side parameters (parameters on b) there should still be a way to avoid making the updates a la ParameterJuMP allows without additional variables. |
|
How does this compare to ParametricOptInterface performance-wise? We should at least run these benchmarks: https://github.com/jump-dev/ParametricOptInterface.jl/tree/master/benchmark If there is a selection of benchmarks that are considered better, we should certainly add to the ones in POI. |
We should benchmark thoroughly before merging. Also,
This is very relevant. |
|
|
||
| MOI.Bridges.needs_final_touch(::FixParametricVariablesBridge) = true | ||
|
|
||
| function MOI.Bridges.final_touch( |
There was a problem hiding this comment.
I think there should be some comments in this function
I'd suggest and offer to help with benchmarking different aspects of the approaches. Not only the solve time matter but as we learned with ParameterJuMP the build cost can become problematic. |
|
Many of the simple benchmarks in POI are purely build-time benchmarks. More examples there are very welcome |
Note: I am unsure if/how this compares to I have looked at function build_and_pass_poi(N::Int64, opt)
@info "Build model"
@time begin
model = direct_model(ParametricOptInterface.Optimizer(opt))
set_silent(model)
@variable(model, x[1:N])
@variable(model, p[i=1:N] in ParametricOptInterface.Parameter(1.0))
@constraint(model, c[i=1:N], (1.5 * p[i]) * x[i] >= 1)
@objective(model, Min, sum(x))
end
@info "First optimize"
@time optimize!(model)
@info "Update"
@time begin
MOI.set.(model, POI.ParameterValue(), p, 2.0)
end
@info "Re-optimize"
@time optimize!(model)
@info "Total"
return model
endNote that I have dropped using Remark: I am not sure how to use Results for POI:Results for jump-dev/JuMP.jl#3215Results for the bridge-based version (this PR)BenchmarkTools resultsFor POI: BenchmarkTools.Trial: 30 samples with 1 evaluation.
Range (min … max): 479.372 ms … 582.332 ms ┊ GC (min … max): 15.22% … 27.34%
Time (median): 535.092 ms ┊ GC (median): 23.12%
Time (mean ± σ): 534.643 ms ± 19.614 ms ┊ GC (mean ± σ): 22.56% ± 3.04%
▁ ▄ █▁▄ ▁▁
▆▁▁▁▁▆▁▁▁▁▁▁▁▁▁▁▆▁▁▁▁▁▁▆▁▁▁█▆█▆▁███▆▁▁▁██▁▆▆▁▆▆▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁
479 ms Histogram: frequency by time 582 ms <
Memory estimate: 340.02 MiB, allocs estimate: 5348999.BenchmarkTools.Trial: 30 samples with 1 evaluation.
Range (min … max): 398.501 ms … 473.766 ms ┊ GC (min … max): 15.43% … 25.29%
Time (median): 449.538 ms ┊ GC (median): 23.65%
Time (mean ± σ): 442.453 ms ± 18.666 ms ┊ GC (mean ± σ): 22.76% ± 3.24%
▃ ▃ ▃ ▃ ██ ▃ ▃█
▇▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▇▁█▇▁▁▁▁▁▁▁▇█▁▇▇▁▁██▁▁█▇██▁▁▁▇▁▁▁▁▁▁▁▇ ▁
399 ms Histogram: frequency by time 474 ms <
Memory estimate: 316.26 MiB, allocs estimate: 5097227.For the bridge-based version (this PR) ConclusionThe bridge based version and POI are pretty similar performance wise, with the bridge based version being better in terms of memory and slightly better build times. However, it is probably less easy to use and does not support everything that POI supports. The "JuMP"-based version (that basically intercepts the |
|
@sstroemer do you want to email me so we can set up a call and I can walk you through the lay of the land and the different options we've tried? |
|
POI was not optimized for memory, might be possible to improve that. |
POI was planned as a replacement for ParameterJuMP, my ideal recommendation would be to contribute to POI. Let me know if you and oscar figure something out, I can help with POI guidance if needed. |
You are correct, ParameterJuMP does not compare as it does not cover multiplicative parameters, Also, you don't want to you the benders benchmark, I meant this: or this: We should have more complete end-to-end benchmarks indeed, like the one you used here. |
|
I've renamed this to We can always remove later if it turns out not to be useful. |
|
Although we could potentially wait for #2095 and then update this bridge to work only on |
|
I've updated this to use the new |
|
Thoughts on what we should do with this? |
|
Monthly developer call suggests we leave this open for now and work on simplifying ParametricOptInterface. (It's likely that we can use the code here to simplify what's in POI.) We can always merge this later, but it shouldn't be enabled by default. |
|
I'm not sure where we landed on this, but this shouldn't be turned on by default, and it seems unlikely to get used very much if it's opt-in. Since we discussed this PR, @joaquimg put a lot of improvements into POI:
So perhaps it's best if we close and focus our efforts on https://github.com/jump-dev/ParametricOptInterface.jl. |
|
Closing as won't-merge. |
Idea
Preface: This PR was moved over from PR#3211 in JuMP.jl. It is not meant as an actual PR, but more as starting point for the discussion of a possible implementation for (partial, not generally viable) parameters in (MI)LP models. I copied the relevant paragraphs.
The main idea originated while skimming through this rework in PowerSimulations.jl: using fixed variables as parameters. While (from now on I'm only thinking about LPs) this works easily as long as the parameters are part of
b(thought as in some constraints likeAx <= b), because the solver just eliminates them during presolve. However, as soon as they are part ofAx, they actually create an expressionparam * var, and since both of them are actuallyJuMP.VariableRefthis will entail the creation of aQuadExpr(instead of anAffExpr).After asking about "multiplicative parameters", the idea (see here for the first steps) was that - since solvers like
HiGHSactually do not support quadratic terms in the constraints - adding a MathOptInterface bridge that convertsMOI.ScalarQuadraticFunctionintoMOI.ScalarAffineFunctioncouldTherefore the usage of a "ParameterRef" would then entail:
set_value(param, val) = fix(param, val; force=true)automatically takes care of updating the parameter's valueInitial prototype
I added the code that @odow proposed.
A simple test model is available in this gist. It implements a basic concept of something generating (
gen) a commodity (here: heat) that can be stored (charge/discharge) intostate(with some losses due to efficiency), while making sure that a given demand is met. All that for each timestep. The parameters are: thedemand(RHS / constant), thecost(coefficient in the objective function) of generation, and the conversion efficiencycop(coefficient in the constraints). Those are randomized. While simple, this is a pretty "hard" example, since 50% of all constraints contain a multiplicative parameter that needs to be adjusted.Benchmarks
An initial performance test compares the timing for the steps: build, solve, re-build/update, re-solve. This is done for a standard model that is rebuilt from scratch each time vs. a parametric model. Timings (all values given in milliseconds) obtained from BenchmarkTools are listed below. All values correspond to median results. "Re-***" steps are only given for the parametric model, assuming that the standard would roughly be the same time in each iteration.
Percentage time increases/decreases for the total of the first iteration (timing it. #1) as well as all subsequent iterations (timing it. #i) are given at the end ("+" means the parametric model is slower, "-" means it's faster).
Tdenotes the number of timesteps in the model; scroll to the right!.While the overhead seems okay-ish for the first iteration (basically after 2 iterations we are already faster compared to rebuilding the model), a huge drop in performance can be seen for
T=10000. I suspected some GC shenanigans, but compare:Memory usage and allocs for smaller models are basically the same ratio (~ 3:1), but I assume that some "non-model-size" related constant term seems to dominate the total time for "build+pass" there.
Notes
This is just a collection of notes/thoughts:
relax_with_penalty!(...)direct_model(...)(I'll check that later, didn't have enough time right now)