DeveloperError: Internal Pyomo implementation error

[abb-omidi]

Dear support team,

I have tried to run the following code. The first part is compiled without any issues, but the last shows the following note:

DeveloperError: Internal Pyomo implementation error:
    "sympy expression type '<class 'sympy.core.relational.LessThan'>' not found in the operator map"
    Please report this to the Pyomo Developers.

The code:

ml = ConcreteModel()
ml.I = RangeSet(3)
ml.M = RangeSet(2)
ml.PAIRS = Set(initialize = ml.I * ml.I, dimen=2, filter=lambda ml, i, j : i < j)

ml.x = BooleanVar(ml.I, ml.I)
ml.start = Var(ml.I, bounds=(0, 10))
ml.finish = Var(ml.I, bounds=(0, 10))

ml.c1 = LogicalConstraint(ml.PAIRS, ml.M, rule=lambda ml, i, j, m: land(ml.x[i,j], ml.x[j,m]).implies(lor(ml.finish[j]<=ml.start[i], ml.finish[i]<=ml.start[j])))
ml.c1.pprint()

# where the issue is occurred: 
TransformationFactory('core.logical_to_linear').apply_to(ml)
ml.logic_to_linear.transformed_constraints.pprint()

Would you please, how can I fix this? Regards

[mrmundt]

Thanks for the report, @abb-omidi . Please provide the following information:

Pyomo version: Python version: Operating system: How Pyomo was installed (PyPI, conda, source): Solver (if applicable):

For future reference, please use the Bug report or Feature/enhancement request templates we supply. They ask all of the information we need to get started on debugging or prioritizing requests.

[abb-omidi]

Thanks @mrmundt for your attention. Actually, I am working with Pyomo with Google-Colab. But, as far as I can see:

Pyomo version: 6.5.0 Python version: 3.9 Operating system: Google-colab How Pyomo was installed (PyPI, conda, source): !pip install Pyomo Solver (if applicable): -

I hope this would be helpful.

Regards Abbas

[jsiirola]

In this case, it is correct that you are getting an error, but we could be generating a better (different) message. I'm working on a PR to improve the error message.

The crux of the problem is that you are mixing Logical and Relational expressions. That is, you are putting a relational expression as a node in a logical expression. The core.logical_to_linear transformation was designed for convert purely logical expressions (i.e., combinations of BooleanVars with and, or, xor, implies, equivalent, atleast, atmost, exactly expressions). We are in the process of moving away from core.logical_to_linear and to a new transformation (contrib.logical_to_disjunctive) that is significantly faster and doesn't rely on sympy. However, that transformation does not yet support "mixed" logical expressions, either (although it is on the relatively short list to add). The challenge here is how to convert the relational expression to a boolean variable. The simple approach is to add a disjunction for each term. That is, convert:

(ml.x[i,j] & ml.x[j,m]) --> (ml.finish[j]<=ml.start[i] | ml.finish[i]<=ml.start[j])`

to (pardon the abuse of notation)

[ Y1[i,j]; ml.finish[j]<=ml.start[i] ] V [ ~Y1[i,j] ]
[ Y2[i,j]; ml.finish[i]<=ml.start[j] ] V [ ~Y2[i,j] ]
(ml.x[i,j] & ml.x[j,m]) --> (Y1[i,j] | Y2[i,j])

...unfortunately, while easy, that is a pretty poor relaxation. A better approach is to convert the inner lor to a disjunction:

(ml.x[i,j] & ml.x[j,m]) --> ([ ml.finish[j]<=ml.start[i] ] V [ ml.finish[j]<=ml.start[i] ])

The relaxation here is better, but the transformation is more complicated (you have to recognize the inner disjunction as a nested disjunction and relax the implicit "or constraint' by a new Boolean variable. That is, something like:

[ Y1; ml.finish[j]<=ml.start[i] ] V [ Y2; ml.finish[j]<=ml.start[i] ]
Y3 = Y1 | Y2
(ml.x[i,j] & ml.x[j,m]) --> Y3

OK, so how to work around this right now? The "easiest" thing is to explicitly write the nested disjunction:

[Y1; [ ml.finish[j]<=ml.start[i] ] V [ ml.finish[j]<=ml.start[i] ]] V [ ~Y1 ]
(ml.x[i,j] & ml.x[j,m]) --> Y1

In Pyomo, we could do it with something like:

from pyomo.gdp import Disjunction

@ml.Disjunct(ml.PAIRS)
def enforce_no_overlap(b, i, j):
    ml = b.model()
    @b.Disjunction()
    def no_overlap(b):
        return [ml.finish[i]<=ml.start[j], ml.finish[j]<=ml.start[i]]

@ml.Disjunction(ml.PAIRS)
def disj(ml, i, j):
    return [ml.enforce_no_overlap[i,j], []]

@ml.LogicalConstraint(ml.PAIRS, ml.M)
def c1(ml, i, j, m):
    return land(ml.x[i,j], ml.x[j,m]).implies(ml.enforce_no_overlap[i,j].indicator_var)

To solve the problem, you would need to convert the whole thing to a MIP, e.g., using a BigM relaxation (which will call contrib.logical_to_disjunctive for you):

TransformationFactory('gdp.bigm').apply_to(ml)

[abb-omidi]

Dear @jsiirola,

Many thanks for your detailed explanation. I confirm that the above form works, at least without throwing any error, but I need a hard time debugging what really is happening as I am pretty new to using Pyomo. Let's say, I have worked on some variants of the parallel machine scheduling problem. The mentioned logical expression, as you correctly pointed out, came back to well-known machine_no_overlap_constraints. This can be written as:

comp[j] <= start[i] + M*(1-x[i,m]) + M*(1-x[j,m]) + M*(alpha[i,j]),    forall m, i, j : i < j;
comp[i] <= start[j] + M*(1-x[i,m]) + M*(1-x[j,m]) + M*(1- alpha[i,j]),    forall m, i, j : i < j;

It is also, very similar to the implementation of Pyomo scheduling books. And it can be easily transformed into its equivalent GDF. But, it produces a weak relaxation. Now, I am seeking a tighter form of the above transformation. I have tried the following derivations:

Iff (x[i,m] and x[j,m]) => (c[j]<=s[i] or c[i]<=s[j]) => 
Iff (x[i,m] and x[j,m]) => (w_1[i,j] <=> c[j]<=s[i] or w_2[i,j] <=>c[i]<=s[j]) =>
1) (1-x[i,m]) + (1-x[j,m]) + w_1[i,j] + w_2[i,j] >= 1;
2) c[j] <= s[i] + M*(1-w_1[i,j]);
3) c[j] >= s[i] -M*w_1[i,j];
4) c[i] <= s[j] + M*(1-w_2[i,j]);
5) c[i] >= s[j] -M*w_2[i,j];

After solving the model with both approaches, I found that the second one can be solved so faster that the first one. The first transformation takes around 28 sec, while the second just takes around 2 sec. In both cases, the optimal solution is the same and the results are reasonable. 1) At the moment what I am willing to know is, does my second transformation and your last formulation implies the same thing? 2) Would you pleaes, how can I write my second derivations as a GDF? I tried something like this, but it cannot implies no_overlap_constraint:

m.z = Var(m.J, m.M, domain=Binary)
m.y = Var(m.PAIRS, domain=Binary)    
BigM = max([JOBS[j]['release'] for j in m.J]) + sum([JOBS[j]['duration'] for j in m.J])
m.w_1 = Var(m.PAIRS, domain=Binary)
m.w_2 = Var(m.PAIRS, domain=Binary)

m.c7 = LogicalConstraint(m.PAIRS, m.M, rule=lambda m, i, j, mach: land(m.z[i,mach], m.z[j,mach]).implies(lor(m.w_1[i,j], m.w_2[i,j])))
m.c8 = Disjunction(m.PAIRS, rule=lambda m, j, k: [ m.start[j] + JOBS[j]['duration'] <= m.start[k] + BigM*(1- m.w_1[j,k]), m.start[j] + JOBS[j]['duration'] >= m.start[k] - (BigM)*(m.w_1[j,k]) ])
m.c9 = Disjunction(m.PAIRS, rule=lambda m, j, k: [ m.start[k] + JOBS[k]['duration'] <= m.start[j] + BigM*(1- m.w_2[j,k]), m.start[k] + JOBS[k]['duration'] >= m.start[j] - (BigM)*(m.w_1[j,k]) ])

Best regards Abbas