Whether support formulating contraints with notation commonly used in mathematical optimization

[Zhanwei-Liu]

Dear PyOptInterface's Development Team,

Thank you for providing such a valuable open-source package.

I often need to handle models with tens of millions of variables and constraints and am currently struggling with the slow performance of Pyomo in large-scale models. Although I know that JuMP also is excellent, I prefer not to mix Julia code in my Python programme. Therefore, I am very interested in PyOptInterface.

I would like to know if PyOptInterface supports formulating optimization problems in Python in a manner similar to the notation commonly used in mathematical optimization. I am aware that Pyomo and JuMP have this feature, and I am wondering if PyOptInterface offers similar functionality.

Thank you!

[metab0t]

Thanks for your kind words!

Firstly, the development of PyOptInterface just originates from our need to construct large optimization model in power system, and we are not satisfied with other solutions in Python, so you are in the right place!

I am not quite sure what "formulating optimization problems in Python in a manner similar to the notation commonly used in mathematical optimization" you mean. Do you mean the expression-based way to construct constraint or multidimensional constraint construction?

If you have some specific examples in Pyomo or JuMP to illustrate the features you need, I can write corresponding code in PyOptInterface to show the usage.

[metab0t]

If you mean adding constraint like model.add_linear_constraint(x<=y), then we choose not to implement it deliberately.

Although it can be implemented by overloading the __le__ function or <= operator of the expression object, I haven't figured out a good way to avoid copying the expression when constructing a temporary constraint object, which might be expensive when the expression is huge (each expression is a C++ object).

Therefore, we adopt a more explicit API like model.add_linear_constraint(x - y, poi.Leq, 0.0). It is not the perfect solution, but we think the efficiency to add new constraint is more important.

JuMP.jl implements this by using macro to rewrite x<=y to x-y<=0.0 then calls add_constraint(x-y, LessThan(0.0)) internally without loss of runtime efficiency. It is a clever way to combine user-friendly expression formula and efficient low-level API, but we do not have similar macro system in Python.

[Zhanwei-Liu]

Thank you for your prompt and thorough response. It has been incredibly helpful in clarifying my understanding. In our model, we anticipate the frequent use of a large number of sparse sets. However, I've found that make_tupledict is not as convenient for building variables, as it requires the use of rules. I'm wondering if it would be possible to create a wrapper for make_tupledict similar to the Var and Constraint functions in Pyomo. Additionally, I believe such a function should be able to take the model class as an input.

Another aspect I have tried to incorporate variables and constraints as attributes of the model class, rather than creating separate a container to store them.

I’m not sure if these practices are slowing down my modeling. I have following examples for your reference. Thank you very much!

from pyomo.environ import *
import pyoptinterface as poi
from pyoptinterface import gurobi 

############################ Define parameters ################################
foods = [
    "Cheeseburger", "Ham Sandwich", "Hamburger", "Fish Sandwich",
    "Chicken Sandwich", "Fries", "Sausage Biscuit", "Lowfat Milk",
    "Orange Juice"
]
nutrients = [
    "Cal", "Carbo", "Protein", "VitA", "VitC", "Calc", "Iron"
]

c = {k: v for k, v in zip(foods, 
    [1.84, 2.19, 1.84, 1.44, 2.29,0.77, 1.29, 0.60, 0.72])}
V = {k: v for k, v in zip(foods, 
    [4.0, 7.5, 3.5, 5.0, 7.3, 2.6, 4.1, 8.0, 12.0])}
Nmin = {k: v for k, v in zip(nutrients, 
    [2000, 350, 55, 100, 100, 100, 100])}
Nmax = {k: v for k, v in zip(nutrients, 
    [
        float("inf"), 375, float("inf"), 
        float("inf"), float("inf"), float("inf"), float("inf")
    ])}
Vmax= 75.0
a_values = [
        (510, 34, 28, 15, 6, 30, 20),
        (370, 35, 24, 15, 10, 20, 20),
        (500, 42, 25, 6, 2, 25, 20),
        (370, 38, 14, 2, 0, 15, 10),
        (400, 42, 31, 8, 15, 15, 8),
        (220, 26, 3, 0, 15, 0, 2),
        (345, 27, 15, 4, 0, 20, 15),
        (110, 12, 9, 10, 4, 30, 0),
        (80, 20, 1, 2, 120, 2, 2)
    ]
a = {
    (food, nutrient): a_values[ix][iy]
    for ix,food  in enumerate(foods) for iy, nutrient in enumerate(nutrients) 
}

################################### Pyomo #####################################
model =  ConcreteModel()
model.x = Var(foods, within=NonNegativeIntegers)
def cost_rule(model):
    return sum(c[i]*model.x[i] for i in foods)
model.cost = Objective(rule=cost_rule)
def nutrient_rule(model, j):
    value = sum(a[i,j]*model.x[i] for i in foods)
    return inequality(Nmin[j], value, Nmax[j])
model.nutrient_limit = Constraint(nutrients, rule=nutrient_rule)
def volume_rule(model):
    return sum(V[i]*model.x[i] for i in foods) <= Vmax
Volume = Constraint(rule=volume_rule)

solver = SolverFactory('gurobi')
solver.solve(model)
print("x = ", [value(model.x[i]) for i in foods])
print("Cost = ", model.cost())

############################### pyoptinterface ################################
# pyoptinterface
model = gurobi.Model()
model.set_model_attribute(poi.ModelAttribute.Silent, True)
import itertools

def Var(*indexes, **kwargs):
    if len(indexes) == 0:
        return model.add_variable(**kwargs)
    elif len(indexes) == 1:
        return {index[0]: model.add_variable(**kwargs) 
            for index in itertools.product(*indexes)}
    else:
        return {index:model.add_variable(**kwargs) 
            for index in itertools.product(*indexes)}

def Constraint(*args, rule):
    if len(args) == 0:
        return rule()
    return poi.make_tupledict(*args, rule=rule)

model.x = Var(foods, lb=0, domain=poi.VariableDomain.Integer)
def cost_rule(model):
    return sum(c[i]*model.x[i] for i in foods)
model.cost = cost_rule(model)
model.set_objective(model.cost, poi.ObjectiveSense.Minimize)

def nutrient_upper_rule(j):
    value = sum(a[i,j]*model.x[i] for i in foods)
    return model.add_linear_constraint(value, poi.Geq, Nmin[j])
def nutrient_lower_rule(j):
    value = sum(a[i,j]*model.x[i] for i in foods)
    return model.add_linear_constraint(value, poi.Leq, Nmax[j])
def volume_rule():
    value = sum(V[i]*model.x[i] for i in foods)
    return model.add_linear_constraint(value, poi.Leq, Vmax)

model.nutrient_upper_limit = Constraint(nutrients, rule=nutrient_upper_rule)
model.nutrient_lower_limit = Constraint(nutrients, rule=nutrient_lower_rule)
model.volume = Constraint(rule=volume_rule)
model.optimize()
cost = model.get_value(model.cost)
print("x = ", [model.get_value(model.x[i]) for i in foods])
print("Cost = ", cost)

[metab0t]

For the first issue, your Var and Constraint wrappers are unnecessary except for len(indexes)=0, but I think if there are no indexes to traverse, why not just call the rule function? Anyway, the return value is just a single element.

In PyOptInterface, make_tupledict is no magic, it just constructs a dict like what you do in Var and Constraint wrapper.

For the second issue, writing model.x=... does not harm your performance, because model is a plain Python object and do whatever you like such as reading/writing/adding its attributes.

This is my version of code for the PyOptInterface part:

############################### pyoptinterface ################################
# pyoptinterface
model = gurobi.Model()
model.set_model_attribute(poi.ModelAttribute.Silent, True)

model.x = model.add_variables(foods, lb=0, domain=poi.VariableDomain.Integer)
model.cost = poi.quicksum(c[i]*model.x[i] for i in foods)
model.set_objective(model.cost, poi.ObjectiveSense.Minimize)

def nutrient_upper_rule(j):
    value = poi.quicksum(a[i,j]*model.x[i] for i in foods)
    return model.add_linear_constraint(value, poi.Geq, Nmin[j])
def nutrient_lower_rule(j):
    value = poi.quicksum(a[i,j]*model.x[i] for i in foods)
    return model.add_linear_constraint(value, poi.Leq, Nmax[j])
def volume_rule():
    value = poi.quicksum(V[i]*model.x[i] for i in foods)
    return model.add_linear_constraint(value, poi.Leq, Vmax)

model.nutrient_upper_limit = poi.make_tupledict(nutrients, rule=nutrient_upper_rule)
model.nutrient_lower_limit = poi.make_tupledict(nutrients, rule=nutrient_lower_rule)
model.volume = volume_rule()
model.optimize()
cost = model.get_value(model.cost)
print("x = ", [model.get_value(model.x[i]) for i in foods])
print("Cost = ", cost)

[Zhanwei-Liu]

Great! Your modified version is perfect for me. I will use pyoptinterface to complete my upcoming tasks. Since I am not very familiar with pyoptinterface at the moment, I will thoroughly review the documentation. By the way, I am also troubled by Pyomo's significant memory usage and am curious to see how pyoptinterface performs in this regard. Thank you!

[metab0t]

You are welcome! If you have any questions, bug reports or feature requests when using PyOptInterface, do not hesitate to share it with us.

By the way, the memory consumption of PyOptInterface is specially optimized as described in our paper, so it is expected to reduce the memory usage for you.

It would be great if you find the performance of PyOptInterface satisfying and share us with empirical data on how it accelerates your model construction compared with Pyomo.