Introduction of general constraint in Gurobipy breaks Cobrapy assumption [BUG]

[Palaract]

Is there an existing issue for this?

• [x] I have searched the existing issues

Problem description

I tried to incorporate a piecewise linear constraint into my model. Because Gurobipy offers a function for that, I used the Gurobipy API directly to add it by accessing model.solver.problem. After calling model.solver.problem.update() the change is in the model and it can be optimized. But if I try to call model.optimize(), Gurobipy throws an error: gurobipy.GurobiError: Unable to retrieve attribute 'RC'.

In search of fixing the error I found the culprit. I also prepared a "Proof of Concept" code snippet which shows the problem directly.

Please bear with me as I explain the problem: During normal operation, model.optimize() is eventually called and therefore the get_solution() function in Cobrapy core gets called: https://github.com/opencobra/cobrapy/blob/79a16d445635103343915a0ed2c2a33074eafd27/src/cobra/core/solution.py#L177-L196 There are two cases: The first case is, that the model is an integer model which is just checking the condition if Gurobi sets the NumIntVars higher than zero. The second case is everything else. Cobrapy assumes that every non-integer programming model, like mixed integer programming (MIP) is always a convex, continous model and therefore always tries to get the reduced costs of the model. In the case of adding general constraints, expecially piecewise linear constraints, Gurobi states:

Note that adding any of these constraints to an otherwise continuous model will transform it into a MIP

So Gurobi converts it to a non-continuous, MIP model for which reduced costs are not available.

As I don't know how to get specific variables from Gurobi through optlang, I can't check the NumGenConstrs attribute from Gurobi in the library code, so I could just make another case in the if clause. If someone would take the time to look at this issue and maybe check in with me in how to enable this rather solver specific feature in a solver agnostic way, I would be delighted to hand in a pull request.

Code sample

Code run:

# Minimal Working Example
from cobra.io import load_json_model
from cobra import Reaction
import gurobipy

model = load_json_model('model.json')
print("----- Working Example START -----")
# Check if it is a mixed integer problem
print(f"Is this a mixed integer problem? { 'Yes' if model.solver.problem.getAttr('IsMIP') == 1 else 'No'}")

solution = model.optimize()
print(f"Objective value: {solution.objective_value}")
print("----- Working Example END -----\n")
print("----- Not Working Example START -----")
# Introducing a new general constraint, specifically a PWL constraint
model = load_json_model('model.json')
# Creating new variable before R1
reaction = Reaction('R0')
model.add_reactions([reaction])
# Building and adding the PWL constraint
xvar = model.solver.problem.getVarByName('R0')
yvar = model.solver.problem.getVarByName('R1')
xpts = [0, 1, 2, 3, 4]
ypts = [0, 1, 2, 3, 4]
name = "PWLconstraint"
model.solver.problem.addGenConstrPWL(xvar, yvar, xpts, ypts, name)
model.solver.problem.update()

# Check if it is a mixed integer problem now
print(f"Is this a mixed integer problem? { 'Yes' if model.solver.problem.getAttr('IsMIP') == 1 else 'No'}")

# Try using the cobrapy interface
try:
    solution = model.optimize()
    print(f"Objective value: {solution.objective_value}")
except gurobipy.GurobiError as e:
    print("The Error gets thrown because it is assumed, that the model is always a convex, continous model")
    print(f"Error: {e}")

# Try using the gurobipy interface
model.solver.problem.optimize()
print("Using the Gurobi API directly, doesn't give an error and works")
print(f"Objective value: {model.solver.problem.getAttr('ObjVal')}")

print("----- Not Working Example END -----")

Environment

### Package Information | Package | Version | |:--------|--------:| | cobra | 0.29.0 | ### Dependency Information | Package | Version | |:--------------------|------------:| | appdirs | 1.4.4 | | black | 23.11.0 | | bumpversion | **missing** | | depinfo | 2.2.0 | | diskcache | 5.6.3 | | future | 0.18.3 | | httpx | 0.25.1 | | importlib-resources | 6.1.1 | | isort | **missing** | | numpy | 1.24.2 | | optlang | 1.8.1 | | pandas | 2.1.3 | | pydantic | 2.5.1 | | python-libsbml | 5.20.2 | | rich | 13.7.0 | | ruamel.yaml | 0.18.5 | | scipy | 1.11.4 | | swiglpk | 5.0.10 | | tox | **missing** | ### Build Tools Information | Package | Version | |:-----------|--------:| | pip | 23.0.1 | | setuptools | 65.5.0 | ### Platform Information | | | |:--------|-------------:| | Linux | 6.6.3-x86_64 | | CPython | 3.10.12 |

Anything else?

Attached is the json for the model.json for the script that shows the problem:

{
    "metabolites": [
        {
            "id": "Glu",
            "name": "Glucose",
            "compartment": ""
        },
        {
            "id": "ATP",
            "name": "Adenosine Tri-Phosphate",
            "compartment": ""
        },
        {
            "id": "Pre",
            "name": "Precursor",
            "compartment": ""
        },
        {
            "id": "Bio",
            "name": "Biomass",
            "compartment": ""
        }
    ],
    "reactions": [
        {
            "id": "R1",
            "name": "Glucose Uptake",
            "metabolites": {
                "Glu": 1.0
            },
            "lower_bound": 0,
            "upper_bound": 1000,
            "gene_reaction_rule": ""
        },
        {
            "id": "R2",
            "name": "Glucose to ATP",
            "metabolites": {
                "ATP": 1.0,
                "Glu": -1.0
            },
            "lower_bound": 0,
            "upper_bound": 1000,
            "gene_reaction_rule": ""
        },
        {
            "id": "R3",
            "name": "Glucose to Precursor",
            "metabolites": {
                "Glu": -1.0,
                "Pre": 1.0
            },
            "lower_bound": 0,
            "upper_bound": 1000,
            "gene_reaction_rule": ""
        },
        {
            "id": "R4",
            "name": "ATP to Biomass",
            "metabolites": {
                "ATP": -1.0,
                "Bio": 1.0
            },
            "lower_bound": 0,
            "upper_bound": 1000,
            "gene_reaction_rule": ""
        },
        {
            "id": "R5",
            "name": "Precursor to Biomass",
            "metabolites": {
                "Bio": 1.0,
                "Pre": -1.0
            },
            "lower_bound": 0,
            "upper_bound": 1000,
            "gene_reaction_rule": ""
        },
        {
            "id": "R6",
            "name": "Biomass Product",
            "metabolites": {
                "Bio": -1.0
            },
            "lower_bound": 0,
            "upper_bound": 1000,
            "gene_reaction_rule": "",
            "objective_coefficient": 1.0
        }
    ],
    "genes": [],
    "id": "PoCModel",
    "compartments": {},
    "version": "1"
}

Also an image showing the example network used in a graphical way:

[cdiener]

Hi, this is more related to optlang which only supports linear constraints. The main reason is that for COBRApy we need a unified interface to support different solver backends. This makes it hard to implement features that are only supported by one or a few solvers. For instance, what should happen to the constraint if you change the solver with model.solver = "glpk"?

That's why you need to use the Gurobi API here as well.

I do like the idea to fix the mip detection in the optlang gurobi interface though. Do you know for which Gurobi versions Model.getAttr('IsMIP') is supported? This wouldn't support the full API but as long as we don't expose those constraints in optlang it would be fine. Let me know if you want to take a stab at it or would like me to.

[Palaract]

I can easily implement an optlang function which returns if the model is a MIP or not. I'll do that gladly. The IsMIP attribute is supported since at least version 8.0, which should be backwards compatible enough.

Would you think that it would be better to make the gurobi optlang interface more foolproof by checking if reduced costs and shadow prices are actually available and returning None instead of throwing an error?

In the specific case I provided, this would also work with cobrapy because the lists reduced and shadow are then only passed into a pandas series which I think doesn't care if the data attribute is None.

If you think that this is a good solution, I will build that out tomorrow and make an issue over at optlang and also provide the necessary PR for it.

[cdiener]

For the integer check you would only have to adjust the already existing one in the Gurobi interface. And if you think we should adjust the primal/dual getters as well, feel free to do so in the functions below here.

Regarding get_solution that is a good question. Would depend on whether other solvers will give sensible results there when the problem is integer. Though I would prefer optlang throwing an error if you can't/shouldn't get them and then have get_solution just deal with that explicitly. If you do want to make PRs against get_solution it would be great if you can add a small benchmark against the previous implementation because it is one of our bottlenecks so we try to keep it fast.

But yeah, sounds good. Definitely start with the PR in optlang and we can follow that with one here. Thanks!

[Palaract]

Yeah, that sounds good! I think I will add a function like is_integer called is_mip and use that in the reduced costs and shadow prices getter functions to raise an ValueError like they did for the integer problem. After that, I will try to build an elegant error handling patch for the get_solution function which doesn't impact performance. Of course I'll follow the contribution guidelines and hand in a benchmark according to that in the PR.

Thank you for taking the time to answer my questions and giving me proper directions to where to look and planning things out with me :)