0

A while ago I started trying to make a fair performance comparison of Pyomo and JuMP. I created a GitHub repository to share the code I use and I am happy if anyone wants to contribute. As of now, I am trying to find a more efficient implementation for the Pyomo models than the intuitive one i came up so far:

IJKLM

def pyomo(I, IJK, JKL, KLM, solve):
    model = pyo.ConcreteModel()

    model.I = pyo.Set(initialize=I)

    x_list = [
        (i, j, k, l, m) for (i, j, k) in IJK for l in JKL[j, k] for m in KLM[k, l]
    ]

    constraint_dict_i = {
        ii: ((i, j, k, l, m) for (i, j, k, l, m) in x_list if i == ii) for ii in I
    }

    model.x_list = pyo.Set(initialize=x_list)
    model.c_dict_i = pyo.Set(model.I, initialize=constraint_dict_i)

    model.z = pyo.Param(default=1)

    model.x = pyo.Var(model.x_list, domain=pyo.NonNegativeReals)

    model.OBJ = pyo.Objective(expr=model.z)

    model.ei = pyo.Constraint(model.I, rule=ei_rule)

    if solve:
        opt = pyo.SolverFactory("gurobi")
        opt.solve(model)


def ei_rule(model, i):
    return sum(model.x[i, j, k, l, m] for i, j, k, l, m in model.c_dict_i[i]) >= 0

Supply Chain

def intuitive_pyomo(I, L, M, IJ, JK, IK, KL, LM, D, solve):
    model = pyo.ConcreteModel()

    model.I = pyo.Set(initialize=I)
    model.L = pyo.Set(initialize=L)
    model.M = pyo.Set(initialize=M)
    model.IJ = pyo.Set(initialize=IJ)
    model.JK = pyo.Set(initialize=JK)
    model.IK = pyo.Set(initialize=IK)
    model.KL = pyo.Set(initialize=KL)
    model.LM = pyo.Set(initialize=LM)

    model.f = pyo.Param(default=1)
    model.d = pyo.Param(model.I, model.M, initialize=D)

    model.x = pyo.Var(
        [
            (i, j, k)
            for (i, j) in model.IJ
            for (jj, k) in model.JK
            if jj == j
            for (ii, kk) in model.IK
            if (ii == i) and (kk == k)
        ],
        domain=pyo.NonNegativeReals,
    )

    model.y = pyo.Var(
        [(i, k, l) for i in model.I for (k, l) in model.KL], domain=pyo.NonNegativeReals
    )

    model.z = pyo.Var(
        [(i, l, m) for i in model.I for (l, m) in model.LM], domain=pyo.NonNegativeReals
    )

    model.OBJ = pyo.Objective(expr=model.f)

    model.production = pyo.Constraint(model.IK, rule=intuitive_production_rule)
    model.transport = pyo.Constraint(model.I, model.L, rule=intuitive_transport_rule)
    model.demand = pyo.Constraint(model.I, model.M, rule=intuitive_demand_rule)

    # model.write("int.lp")

    if solve:
        opt = pyo.SolverFactory("gurobi")
        opt.solve(model)


def intuitive_production_rule(model, i, k):
    lhs = [
        model.x[i, j, k]
        for (ii, j) in model.IJ
        if ii == i
        for (jj, kk) in model.JK
        if (jj == j) and (kk == k)
    ]
    rhs = [model.y[i, k, l] for (kk, l) in model.KL if kk == k]

    if lhs or rhs:
        return sum(lhs) >= sum(rhs)
    else:
        return pyo.Constraint.Skip


def intuitive_transport_rule(model, i, l):
    lhs = [model.y[i, k, l] for (k, ll) in model.KL if ll == l]
    rhs = [model.z[i, l, m] for (lll, m) in model.LM if lll == l]

    if lhs or rhs:
        return sum(lhs) >= sum(rhs)
    else:
        return pyo.Constraint.Skip


def intuitive_demand_rule(model, i, m):
    return sum(model.z[i, l, m] for (l, mm) in model.LM if mm == m) >= model.d[i, m]

I measure performance in model generation time for two exemplary models I refer to as IJKLM and Supply Chain.

Models:

IJKLM IJKLM model

Supply Chain Supply Chain model

These are the results for increasing instance sizes:

IJKLM Model generation IJKLM

Supply Chain Model generation supply chain

Can anyone help to improve Pyomos performance?

5
  • There needs to be code here. Commented May 24, 2023 at 14:15
  • I tinkered with this a bit.... There might be a few things to implement on the pyomo side, but one of the driving issues here is that your data is somewhat nonsensical. At least on the Supply Chain Model. (I didn't look at the other). If you inspect the sizes of some of the sets you are generating, when the size of I is 4600, the domain of x is only ~230. Meaning less than 5% of your items are even capable of being manufactured/delivered. I would expect the domain of x to me many times larger than |I|, somewhere closer to |I| x |J| x |K|. Commented May 24, 2023 at 17:09
  • Due to the fact that your problem scale is minuscule compared to what I think the data should represent, it is difficult to gage any accelerants in construction because you are only generating maybe 2% (?) of the expected constraints. For same reason, your solve times are negligible where they should be driving the discussion. Commented May 24, 2023 at 17:12
  • Could you also include the code to generate your test data sets? Commented May 24, 2023 at 18:58
  • Data and code is here: github.com/justine18/performance_experiment/blob/master/… python main_supply_chain.py. I had to mkdir supply_chain/data before running though. Commented May 24, 2023 at 20:05

3 Answers 3

3

Here are a couple mods/notions for the Supply Chain Pyomo model for you to consider aside from the data issues that you have (under-representative scope and infeasibilities).

  1. Making the lhs set inside the production constraint is expensive, and likely much more expensive if the data were representative. And in the current construct, you are executing that expensive operation inside of a function that is going to execute |model.IK| times. You can wrangle the data separately from the elements that you have to pre-compute the valid indices for model.x as shown with an indexed set in the model, or just a dictionary if you don't want to make the indexed set in the model. This same construct can be used to generate the domain for model.x as shown, for a tiny speedup (because the domain for x is only computed once) for free.

  2. You can reduce the number of production constraints to only what is required by the demand (rhs in your construct). Presumably, in the production model, you would minimize x, so the only constraints needed are where there is a rhs. Note that I dropped in a little print statement to show infeasibilities where there is demand, but no means to produce (defect in the data).

With those 2 tweaks, I'm getting runs of |model.I| at 7900 in 3.8 seconds build time. I think that will improve (relative to original) with better data.

Your code w/ mods:

import pyomo.environ as pyo
import logging
import timeit
import pandas as pd
import numpy as np
from collections import defaultdict
from itertools import chain

logging.getLogger("pyomo.core").setLevel(logging.ERROR)


########## Intuitive Pyomo ##########
def run_intuitive_pyomo(I, L, M, IJ, JK, IK, KL, LM, D, solve, repeats, number):
    setup = {
        "I": I,
        "L": L,
        "M": M,
        "IJ": IJ,
        "JK": JK,
        "IK": IK,
        "KL": KL,
        "LM": LM,
        "D": D,
        "solve": solve,
        "model_function": intuitive_pyomo,
    }
    r = timeit.repeat(
        "model_function(I, L, M, IJ, JK, IK, KL, LM, D, solve)",
        repeat=repeats,
        number=number,
        globals=setup,
    )

    result = pd.DataFrame(
        {
            "I": [len(I)],
            "Language": ["Intuitive Pyomo"],
            "MinTime": [np.min(r)],
            "MeanTime": [np.mean(r)],
            "MedianTime": [np.median(r)],
        }
    )
    return result


def intuitive_pyomo(I, L, M, IJ, JK, IK, KL, LM, D, solve):

    # some data wrangling
    IJ_dict = defaultdict(set)
    for i, j in IJ:  IJ_dict[i].add(j)
    KJ_dict = defaultdict(set)
    for j, k in JK:  KJ_dict[k].add(j)
    # make a dictionary of (i, k) : {(i, j, k) tuples}
    IK_IJK = {(i, k): {(i, j, k) for j in IJ_dict.get(i, set()) & KJ_dict.get(k, set())}
             for (i, k) in IK}

    model = pyo.ConcreteModel()

    model.I = pyo.Set(initialize=I)
    model.L = pyo.Set(initialize=L)
    model.M = pyo.Set(initialize=M)
    model.IJ = pyo.Set(initialize=IJ)
    model.JK = pyo.Set(initialize=JK)
    model.IK = pyo.Set(initialize=IK)
    model.KL = pyo.Set(initialize=KL)
    model.LM = pyo.Set(initialize=LM)

    model.IK_IJK = pyo.Set(IK_IJK.keys(), initialize=IK_IJK)

    model.f = pyo.Param(default=1)
    model.d = pyo.Param(model.I, model.M, initialize=D)

    # x_idx = [
    #         (i, j, k)
    #         for (i, j) in model.IJ
    #         for (jj, k) in model.JK
    #         if jj == j
    #         for (ii, kk) in model.IK
    #         if (ii == i) and (kk == k)
    #     ]
    x_idx_quick = list(chain(*IK_IJK.values()))
    # assert set(x_idx) == set(x_idx_quick)   # sanity check.  Make sure it is same...
    model.x = pyo.Var(x_idx_quick,
        domain=pyo.NonNegativeReals,
    )
    print(f'length of model.I: {len(I)}')
    print(f'length of modek.IK: {len(IK)}')
    print(f'size of model.x: {len(x_idx_quick)}')

    model.y = pyo.Var(
        [(i, k, l) for i in model.I for (k, l) in model.KL], domain=pyo.NonNegativeReals
    )

    model.z = pyo.Var(
        [(i, l, m) for i in model.I for (l, m) in model.LM], domain=pyo.NonNegativeReals
    )

    model.OBJ = pyo.Objective(expr=model.f)


    # model.write("int.lp")

    if solve:
        opt = pyo.SolverFactory("cbc")
        opt.solve(model)


    def intuitive_production_rule(model, i, k):
        lhs = [model.x[i, j, k] for i, j, k in model.IK_IJK[i, k]]

        rhs = [model.y[i, k, l] for (kk, l) in model.KL if kk == k]

        # show where the data is infeasible...
        # if rhs and not lhs:
        #     print(f'infeasible for (i, k): {i}, {k}')

        if lhs and rhs:
            return sum(lhs) >= sum(rhs)
        else:
            return pyo.Constraint.Skip


    def intuitive_transport_rule(model, i, l):
        lhs = [model.y[i, k, l] for (k, ll) in model.KL if ll == l]
        rhs = [model.z[i, l, m] for (lll, m) in model.LM if lll == l]

        if lhs or rhs:
            return sum(lhs) >= sum(rhs)
        else:
            return pyo.Constraint.Skip


    def intuitive_demand_rule(model, i, m):
        return sum(model.z[i, l, m] for (l, mm) in model.LM if mm == m) >= model.d[i, m]

    model.production = pyo.Constraint(IK_IJK.keys(), rule=intuitive_production_rule)
    print(f'created {len(model.production)} production constraints')
    model.transport = pyo.Constraint(model.I, model.L, rule=intuitive_transport_rule)
    model.demand = pyo.Constraint(model.I, model.M, rule=intuitive_demand_rule)
Sign up to request clarification or add additional context in comments.

Comments

3

IJKLM

The first thing to mention is to update Pyomo. Running on Pyomo 6.6.0 under Python 3.11 gives a slightly different result for the IJKLM model (showing the model generation time; JuMP 1.11.1 on Julia 1.9.0; RHEL7):

enter image description here

Now, to add to @AirSquid's answer, The IJKLM model scales linearly with |I|, and you should expect that the Pyomo runtime should also scale linearly as well. The fact that it is scaling quadratically indicates a modeling problem. The first thing to do is find out where. Pyomo has a built-in tool to report the timing information for common operations, including the time it takes to construct each component. To enable that, call the following before creating a model:

from pyomo.common.timing import report_timing
report_timing()

Rerunning python main_IJKLM.py gives:

           0 seconds to construct Block ConcreteModel; 1 index total
           0 seconds to construct Set I; 1 index total
        0.10 seconds to construct Set x_list; 1 index total
        5.62 seconds to construct Set c_dict_i; 3700 indices total
           0 seconds to construct Set Any; 1 index total
           0 seconds to construct Param z; 1 index total
        0.03 seconds to construct Var x; 69382 indices total
           0 seconds to construct Objective OBJ; 1 index total
        0.05 seconds to construct Constraint ei; 3700 indices total
           0 seconds to construct Block ConcreteModel; 1 index total
           0 seconds to construct Set I; 1 index total
        0.10 seconds to construct Set x_list; 1 index total
        5.71 seconds to construct Set c_dict_i; 3700 indices total
           0 seconds to construct Set Any; 1 index total
           0 seconds to construct Param z; 1 index total
        0.03 seconds to construct Var x; 69382 indices total
           0 seconds to construct Objective OBJ; 1 index total
        0.05 seconds to construct Constraint ei; 3700 indices total
Pyomo             done   3700 in 5.82s

The bulk of the time is actually creating the c_dict_i IndexedSet. If we look at it, this starts to make sense:

constraint_dict_i = {
    ii: ((i, j, k, l, m) for (i, j, k, l, m) in x_list if i == ii) for ii in I
}
model.c_dict_i = pyo.Set(model.I, initialize=constraint_dict_i)

You are initializing the IndexedSet with a dictionary of generators, and each generator has to walk through the entire x_list to extract the relevant subset (giving the observed quadratic behavior). We would be better off to explicitly create the dictionary of lists in a single pass with something like:

constraint_dict_i = {i: [] for i in I}
for idx in x_list:
    constraint_dict_i[idx[0]].append(idx)

Finally, your model generates an error because some of the lists in c_dict_i are actually empty (resulting in a trivial constraint). Your "Intuitive Pyomo" model skips those constraints, and modifying the ei_rule in the "Pyomo" model will resolve that error:

def ei_rule(model, i):
    if not model.c_dict_i[i]:
        return pyo.Constraint.Skip
    return sum(model.x[idx] for idx in model.c_dict_i[i]) >= 0

(As a sidebar, also note that you don't have to unpack the indicies if you are not going to manipulate them)

With these changes, we now see the linear scaling that we should expect:

enter image description here

As a final comment, the JuMP model would also benefit from a similar change to only iterate through the indexing set once.


Supply Chain

We can apply the same trick to the supply chain model. Starting from @AirSquid's answer, we see

enter image description here

From report_timing(), we see that 40% of the time is in the transport constraint, 30% in the demand constraint, 20% in the d Param, and 10% everywhere else. We can leverage IndexedSets to avoid repeated iteration by defining 3 new indexed sets:

L_M = {l: [] for l in model.L}
M_L = {m: [] for m in model.M}
for l, m in model.LM:
    L_M[l].append(m)
    M_L[m].append(l)
L_K = {l: [] for l in model.L}
for k, l in model.KL:
    L_K[l].append(k)

model.L_M = pyo.Set(model.L, initialize=L_M)
model.M_L = pyo.Set(model.M, initialize=M_L)
model.L_K = pyo.Set(model.L, initialize=L_K)

and rewriting the constraint rules:

def intuitive_transport_rule(model, i, l):
    if len(model.L_K[l]) or len(model.L_M[l]):
        return (
            sum(model.y[i, k, l] for k in model.L_K[l]) >= 
            sum(model.z[i, l, m] for m in model.L_M[l])
        )
    return pyo.Constraint.Skip

def intuitive_demand_rule(model, i, m):
    return sum(model.z[i, l, m] for l in model.M_L[m]) >= model.d[i, m]

While not as dramatic as with the IJLKM model, this still gives a ~15% improvement:

enter image description here

1 Comment

0

Thanks for all those helpful tips and tricks. I followed the advice given by @jsiirola on the pyomo side and @Oscar Dowson on the jump side and want to give an update on the new results for anyone interested.

IJKLM Model generation IJKLM

Supply Chain Model generation supply chain

@AirSquid: Thanks for pointing this out. Data generation for the supply chain model is an artifact from the IJKLM model and needs to be adapted. I will implement a new data generator hopefully soonish.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.