Gurobi / gurobipy-pandas

Convenience wrapper for building optimization models from pandas data
https://gurobipy-pandas.readthedocs.io
Apache License 2.0
95 stars 15 forks source link

Support for index broadcasting in add_constrs #89

Open taka110811 opened 1 week ago

taka110811 commented 1 week ago

Description:

I am working on a problem where I need to create constraints with a right-hand side that has two indices, as shown in the following formulation:

$$ \sum_j cj y{ij} \le b_{ij} \quad \forall i \in I, j \in J $$

In the current implementation of gurobipy-pandas, I can implement a constraint with a single-indexed right-hand side using the following syntax:

constraints = gppd.add_constrs(  
    model,
    (c * y).groupby("i").sum(),  # left-hand side
    GRB.LESS_EQUAL,              # inequality (sense)
    b,                           # right-hand side
    name="constr",
)

This works for constraints of the form:

$$ \sum_j cj y{ij} \le b_i \quad \forall i \in I $$

However, I am now facing a situation where the right-hand side has two indices, b_{ij} , as shown in the first equation above.

Question: Is there a way to implement constraints with a two-indexed right-hand side b_{ij} using gurobipy-pandas? If so, could you please provide an example?

If not, would this feature be worth considering as an enhancement to gurobipy-pandas? I believe the ability to handle more complex constraints with two indices on both sides would be a valuable addition for users working on more advanced optimization problems.

Thank you for your consideration.

simonbowly commented 1 week ago

Hi @taka110811, thanks for raising this. Multi-indexes are supported. You can create a set of constraints with any arbitrary left- and right- hand side series, as long as the indices are aligned.

What's missing for your case is broadcasting, i.e. what pandas is doing to align the single-indexed c with the multiindexed y to get the result c * y. Subtraction has the same behaviour, so you can leverage pandas operators to create the constraints you want like so:

constraints = gppd.add_constrs(  
    model,
    (c * y).groupby("i").sum() - b,
    GRB.LESS_EQUAL,
    0.0,
    name="constr",
)

Does that achieve what you're looking for?

Feature request

gurobipy-pandas is intentionally strict on index alignment in add_constrs to avoid common errors. But I agree it would make sense if it were more consistent with pandas' arithmetic operators. I'll take a look.

Side note

If $b_{ij}$ is a series of constants, your constraint expression

$$ \sum_j cj y{ij} \le b_{ij} \quad \forall i \in I, j \in J $$

is really equivalent to

$$ \sum_j cj y{ij} \le \minj b{ij} \quad \forall i \in I $$

which is expressed over a single level index in gurobipy-pandas like so:

constraints = gppd.add_constrs(  
    model,
    (c * y).groupby("i").sum(),
    GRB.LESS_EQUAL,
    b.groupby("i").min(),
    name="constr",
)

This is a better approach since it will create a smaller model that still expresses the same requirement. Even better would be to do this kind of reduction of the input data in pre-processing so that the right-hand side data is already aligned on the correct index.

Of course, your request still would be a good addition in the case that the right-hand side involves variables.

taka110811 commented 1 week ago

@simonbowly Thank you for the prompt reply. Also, I would like to mention that both implementation methods you provided worked correctly.This might get a bit long, but for my own reference as well, I will document the implementation log of the sample case.

Sample Case Here’s the setup for the sample case:

sample case

c

j
0    1.0
1    2.0
2    3.0
Name: c, dtype: float64
y

i  j
0  1    <gurobi.Var y[0,1]>
   2    <gurobi.Var y[0,2]>
1  0    <gurobi.Var y[1,0]>
2  0    <gurobi.Var y[2,0]>
   1    <gurobi.Var y[2,1]>
Name: y, dtype: object
b_multi_index

i  j
0  1    1
   2    2
1  0    3
2  0    4
   1    5
dtype: int64

Constraints using pandas broadcasting

$$ \sum_j cj y{ij} \leq b_{ij} \quad \forall i \in I, j \in J $$

constraints = gppd.add_constrs(  
    model,
    (c * y).groupby("i").sum() - b_multi_index,
    GRB.LESS_EQUAL,
    0.0,
    name="constr",
)
constraints.apply(model.getRow)

i  j
0  1    2.0 y[0,1] + 3.0 y[0,2]
   2    2.0 y[0,1] + 3.0 y[0,2]
1  0                     y[1,0]
2  0        y[2,0] + 2.0 y[2,1]
   1        y[2,0] + 2.0 y[2,1]
Name: constr, dtype: object
constraints.gppd.RHS

i  j
0  1    1.0
   2    2.0
1  0    3.0
2  0    4.0
   1    5.0
Name: constr, dtype: float64

Constraints using min

$$ \sum_j cj y{ij} \leq \minj b{ij} \quad \forall i \in I $$

constraints = gppd.add_constrs(  
    model,
    (c * y).groupby("i").sum(),
    GRB.LESS_EQUAL,
    b_multi_index.groupby("i").min(),
    name="constr",
)
constraints.apply(model.getRow)

i
0    2.0 y[0,1] + 3.0 y[0,2]
1                     y[1,0]
2        y[2,0] + 2.0 y[2,1]
Name: constr, dtype: object
constraints.gppd.RHS

i
0    1.0
1    3.0
2    4.0
Name: constr, dtype: float64

Feature request

Thank you for suggesting the enhancement regarding the broadcasting functionality. If the broadcasting feature is implemented, does the following implementation match the intended usage?

constraints = gppd.add_constrs(  
    model,
    (c * y) - b,  # Here, b is a MultiIndex, but will be handled via broadcasting
    GRB.LESS_EQUAL,
    0.0,
    name="constr",
)
simonbowly commented 1 week ago

Thanks for adding the complete example! The last example is not what I intended though, sorry for the confusion. I'm proposing that this:

# Case 1
# Proposed new feature: broadcast between the left- and right- sides in add_constrs
# (In the current version, this is an error)
constraints = gppd.add_constrs(  
    model,
    (c * y).groupby("i").sum(),
    GRB.LESS_EQUAL,
    b_multi_index,
    name="constr",
)

would give the same result you can already get from this, where the broadcasting happens in the subtraction operation:

# Case 2
constraints = gppd.add_constrs(  
    model,
    (c * y).groupby("i").sum() - b_multi_index,
    GRB.LESS_EQUAL,
    0.0,
    name="constr",
)

i.e. these both create this set of constraints:

$$ \sum_j cj y{ij} \leq b_{ij} \quad \forall i \in I, j \in J. $$

By contrast, this:

# Case 3
constraints = gppd.add_constrs(  
    model,
    (c * y) - b_multi_index,
    GRB.LESS_EQUAL,
    0.0,
    name="constr",
)

and this:

# Case 4
constraints = gppd.add_constrs(  
    model,
    (c * y),
    GRB.LESS_EQUAL,
    b_multi_index,
    name="constr",
)

should already work, but they both model the following:

$$ cj y{ij} \leq b_{ij} \quad \forall i \in I, j \in J $$

which I don't think is what you are after. We won't change the behaviour for (2), (3) or (4), we would only implement the necessary changes to make (1) give the same result as (2).

Note that broadcasting is a pandas concept. In this case it refers to aligning a single-level index with a multi-index, based on the names of the series levels. The left-side entries of (c * y).groupby("i").sum(), indexed by $i$, are broadcasted onto the $(i, j)$ multi-index of the right-side b, such that the entries on the left are matched with each corresponding $b_{ij}$ having the same value of $i$. You still need to carry out the groupby+sum yourself; broadcasting only matches up entries between two series, it does not perform any reduction operations.

I hope that's clear?

Dr-Irv commented 1 week ago

There is a more fundamental issue here is that the original constraint, as mathematically written, has an unclear meaning:

$$ \sum_j cj y{ij} \le b_{ij} \quad\forall i\in I, j\in J $$

The reason this is unclear is that you are stating the constraint for each element $i\in I$ AND each element $j\in J$, but then you are summing over $j$, and it isn't clear which elements of the set $J$ you are summing over in the LHS.

When I write mathematical constraints, anything after the "forall" ($\forall$) essentially binds the values of those indices for each "iterate" of the constraint prior to the "forall". So it is inappropriate from a mathematical sense to have a sum prior to the "forall" that is summing over an index that has been bound by an expression after the "forall".

Fix your math, and then we can decide if gurobipy-pandas needs to do something.