How to speed up model generation in JuMP?

Specific Domains Optimization (Mathematical)
,
1

I have a large-scale scenarios-based optimization model that takes a long time to solve on JuMP and Julia using the Gurobi solver. While solving the constructed model using Gurobi (as it uses all available cores) is quick, the constraints definition part (done by JuMP) takes all the time. Can we parallelize the model definition part on JuMP? I searched online and found that the model definition part is inherently done in a single process and can not be parallelized. I was just wondering if you have encountered the same problem and have a way around it.

** My understanding of solving an optimization model using JuMP/Julia and Gurobi solver is that JuMP defines and constructs the model and passes it to the Gurobi solver to solve it. I want to leverage a parallel computing mechanism for the model generation part (done by JuMP). Apparently, this is not possible. Also, is “model generation” the right term to refer to the part done performed by JuMP?
Thank you!

2

Hi @SantoshSharma, welcome to the forum.

Can you post a reproducible example of your code? There are almost certainly some small tweaks that we can make that will make it run faster.

Have you tried to profile your code to see which bit takes the most time?

You cannot (and should not need to) parallelize the problem creation.

1 Like
3

Hi Odow,

Thanks for the comment. The following two constraints are taking the majority of the time:

for k in 1:size(y, 1)
    for i in 1:NumBus
        for j in 1:NumBus
            if i != j
                @constraint(model, -100 <= B[i, j] * (Angles[i] + sum(B_u[i, m] * (y[k, i] - Lambda[m]) for m in 1:NumBus) - Angles[j] - sum(B_u[j, n] * (y[k, j] - Lambda[n]) for n in 1:NumBus)) <= 100)
            end
        end
    end
end

for k in 1:size(y, 1) 
    for i in 1:NumBus
        for j in 1:NumBus
            if i != j
                @constraint(model,
                    -100 <= B[i, j] * (
                        Angles[i] +
                        sum(B_u[i, m] * (y[k, i] - Lambda[m]) for m in 1:NumBus) -
                        Angles[j] -
                        sum(B_u[j, n] * (y[k, j] - Lambda[n]) for n in 1:NumBus)
                    ) <= 100
                )
            end
        end
    end
end

Here, size(y)=10,000 and NumBus=118. I have 10,000 scenarios and those two set of constraints (118*118) has to be replicated 10,000 times. So basically, will have 10k * 118 * 118 constraints. Total computation time looks like this: roughly 3% for data pre-processing, 90% for model build time, and 7% for solving the model. Since it is a scenario-based model, I have to replicate the same constraint with different parameters many times, which takes most of the time. Since those constraints are independent, I wondered if we could utilize the parallel model construction with JuMP. Thank you!

1 Like
4

A few comments:

You can format your code by putting it in triple backticks

```julia
1 + 1
```

I don’t know that I understand your first @constraint block. Is that valid JuMP syntax? Is it meant to be different to the second block?

10k * 118^2 is a 139,240,000. That is a very big model. It is made worse because Gurobi.jl does not support a <= f(x) <= b constraints, so JuMP rewrites as a <= f(x), and f(x) <= b, so you’re really adding 2 * 139,240,000 constraints.

You could try this reformulation:

@variable(model, -100 <= tmp[1:size(y, 1), 1:NumBus, 1:NumBus] <= 100)
for k in 1:size(y, 1), i in 1:NumBus, j in 1:NumBus
    if i == j
        continue
    end
    @constraint(
        model,
        B[i, j] * (
        Angles[i] +
        sum(B_u[i, m] * (y[k, i] - Lambda[m]) for m in 1:NumBus) -
        Angles[j] -
        sum(B_u[j, n] * (y[k, j] - Lambda[n]) for n in 1:NumBus)
        ) == tmp[k, i, j]
    )
end

What are the variables and what are the data? Is the problem symmetric? Do you need all i and j or is it sufficient to just have the upper or lower triangle?

5

Hi Oscar, I will get back to you with those information shortly. Is “model generation” the correct term to refer to this part where JuMP prepares a model to pass it to the solver?
Thank you!

6

Is “model generation” the correct term

Yes.

You could also try:

using JuMP, Gurobi
model = direct_model(Gurobi.Optimizer())
set_string_names_on_creation(model, false)
# ... model here

Note that this will throw a warning if you try to use interval constraints.

7

Hi Oscar, Thank you so much for the response. While I work on a new implementation of constraints that you mentioned, I wanted to ask you a question. Is there any research or project you know of for improving the speed of model generation for these modeling frameworks like JuMP or Pyomo other than what you already said? We have been running with this slow model generation problem for a while and are looking for ways to improve this step as this problem has been seen more and more for a model with a huge number of constraints. One way I am thinking of is to parallelize the workflow of model generation part, but not sure about the feasibility of this though. Thank you!

8

One way I am thinking of is to parallelize the workflow of model generation part

You cannot parallelize the workflow of model generation.

We have been running with this slow model generation problem

The problem is not that JuMP (or Pyomo) is “slow.” It is that you are building models with ~10^9 variables and constraints. These are very large. How long do they take to build? How much memory do you have available?

If the model takes longer to build than to solve, then it is highly likely that the solution is trivial, and that the problems can almost be solved independently.

What you can, and probably should consider, is some sort of decomposition instead of trying to solve the deterministic equivalent.