Search code examples
pythonconstraintslogical-operatorspyomo

how to model logical-OR or quantifiers in Pyomo-constraints


I am using PYOMO and I want to implement a logical-OR within the "P_constraint_rule", but I cannot make it.

there are some parts of my model:

model.s = Param(within=NonNegativeIntegers)
model.S = RangeSet(1, model.s)
model.f = Param(within=NonNegativeIntegers)
model.F = RangeSet(1, model.f)
model.p = Param(within=NonNegativeIntegers)
model.P = RangeSet(0, model.p)
model.m = Param(within=NonNegativeIntegers)
model.M = RangeSet(1, model.m)
model.g = Param(model.M, default=5.0)
model.b2 = Param(model.F,model.F, within=Binary, default=0)
model.xf1f2 = Var(model.F, model.F, within=Binary)
model.y = Var(model.M, within=Binary)
model.xf = Var(model.F, within=Binary)
model.aff = Param(model.F,model.F, within=NonNegativeIntegers, default=0)

...

model.x = Var(model.index_sfpm, within=Binary)
model.b1 = Param(model.index_sfpm, within=Binary, default=0)


def Obj_rule(model):
    expr = 0.0
    for (s,f,p,m) in model.index_sfpm:
        expr += model.g[m] * model.xf[f] * model.b1[s,f,p,m] * model.x[s,f,p,m]    
    for m in model.M:
        expr += model.g[m] * model.y[m]

    return expr

model.OBJ = Objective(rule=Obj_rule, sense=maximize)


def P_constraint_rule (model, f1, f2):
    expr = 0.0
    for (s,m) in model.index_sm:
        expr += model.b2[f1,f2] * model.xf[f1] * model.b1[s,f1,1,m] * model.x[s,f1,1,m]
        expr += model.xf[f2] * model.b1[s,f2,1,m] * model.x[s,f2,1,m]

e.g. in my .dat : param aff:= 1 7 10

    return expr == model.aff[f1,f2] | expr == 0

model.PConstraint = Constraint(model.index_f1f2, rule=P_constraint_rule)

When I use" | ", I get following error:

ERROR: Rule failed when generating expression for constraint PConstraint with index (7, 3):
        TypeError: unsupported operand type(s) for |: 'int' and '_SumExpression' ERROR: Constructing component 'PConstraint' from data=None failed:
        TypeError: unsupported operand type(s) for |: 'int' and '_SumExpression' [    1.72] Pyomo Finished ERROR: Unexpected exception while running model:
        unsupported operand type(s) for |: 'int' and '_SumExpression'

and when I use " || "

ERROR: Unexpected exception while loading model:
        invalid syntax 

When that constraint is commented, the model and gurobi run fine.

Could anybody assist me with these errors?

Is there another possibility to use a quantifier? The inequation P1.constraint should be valid for model.index_f1f2 The equation P2Constraint should be valid for 2 elements of model.F or 1 element of model.index_f1f2 something like this:

def P1_constraint_rule (model, f1, f2):
expr = 0.0
for (s,m) in model.index_sm:
    expr += model.b2[f1,f2] * model.xf[f1] * model.b1[s,f1,1,m] * model.x[s,f1,1,m]
    expr += model.xf[f2] * model.b1[s,f2,1,m] * model.x[s,f2,1,m]
return expr <= model.aff[f1,f2]

model.P1Constraint = Constraint(model.index_f1f2, rule=P1_constraint_rule)


def P2_constraint_rule (model, f1, f2):
    expr = 0.0
    for (s,m) in model.index_sm:
        expr += model.b2[f1,f2] * model.xf[f1] * model.b1[s,f1,1,m] * model.x[s,f1,1,m]
        expr += model.xf[f2] * model.b1[s,f2,1,m] * model.x[s,f2,1,m]
        #this equation should be valid for 2 elements of model.F or 1 element of model.index_f1f2
     return expr == model.aff[f1,f2]

model.P2Constraint = Constraint(model.index_f1f2, rule=P2_constraint_rule)

thank you in advance, Lala


Solution

  • The error is because you are attempting to specify a non-algebraic constraint. Conceptually, the following would define a logical disjunction:

    expr == model.aff[f1,f2] | expr == 0
    

    For specific syntactic issues:

    • | is binary OR. It binds tighter than relational operators and thus would not do what you want.
    • || is not valid Python syntax
    • conceptually what you would want is logical or, which in Python is implemented with or. That would be nice syntax to support - however, it is not currently supported by Pyomo.

    You have two options for specifying constraints like this: either (1) specify it as a disjunction using the pyomo.gdp extension and then leverage the transformations in pyomo.gdp to relax the disjunctive program back to a MIP, or (2) explicitly relax the disjunction using, for example, a Big-M relaxation. To do the former, you would need to define the two disjuncts and then the disjunction:

    from pyomo.gdp import *
    
    def P_disjunct_rule (b, f1, f2, i):
        model = b.model()
        expr = 0.0
        for (s,m) in model.index_sm:
            expr += model.b2[f1,f2] * model.xf[f1] * model.b1[s,f1,1,m] * model.x[s,f1,1,m]
            expr += model.xf[f2] * model.b1[s,f2,1,m] * model.x[s,f2,1,m]
        if i:
            return expr == model.aff[f1,f2]
        else:
            return expr == 0
    model.PDisjunct = Disjunct(model.index_f1f2, [0,1], rule=P_constraint_rule)
    
    def P_disjunction_rule(m,f1,f2):
        return [ m.PDisjunct[f1,f2,i] for i in [0,1] ]
    model.PDisjunction = Disjunction(model.index_f1f2, rule=P_Disjunction_rule)
    

    You then need to invoke a transformation to convert the disjunctions back to algebraic constraints. Note: the transformations either need your Pyomo variables to all have lower and upper bounds, OR you need to specify valid "Big-M" values through a BigM suffix on your model. You can either:

    • specify the transformation on the Pyomo command line (e.g., --transform=gdp.bigm or --transform=gdp.chull)
    • specify the transformation in a BuildAction

      def xfrm(m):
          TransformationFactory('gdp.bigm').apply_to(m)
      model.xfrm = BuildAction(rule=xfrm)
      
    • explicitly call the transformation as part of your custom driver script.

    The alternative to pyomo.gdp is to explicitly implement the relaxation yourself. You would need to add a binary variable (let's call it y) that indicates which side of the disjunction must be True and then explicitly relax both disjuncts using that binary variable. Conceptually you would turn

    expr == model.aff[f1,f2] | expr == 0
    

    into

    expr - model.aff[f1.f2] <= M1 * y
    model.aff[f1.f2] - expr <= M2 * y
    expr <= M3 * (1-y)
    expr >- M4 * (1-y)
    

    Note that depending on the bounds for expr and aff, some of those constraints may be redundant. Also, the four "Big-M's" (large constants) may not necessarily need to be different, although the problem will solve better if you can determine the smallest valid value for each of the M's.