как моделировать логическое ИЛИ или кванторы в ограничениях Pyomo

Ошибка возникает из-за того, что вы пытаетесь указать неалгебраическое ограничение. Концептуально следующее будет определять логическую дизъюнкцию:

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

Для конкретных синтаксических проблем:

• | - это двоичное ИЛИ. Он связывает более жестко, чем реляционные операторы, и поэтому не будет делать то, что вы хотите.
• || недопустимый синтаксис Python
• концептуально вам нужно логическое или, которое в Python реализовано с помощью or. Было бы неплохо поддерживать такой синтаксис, однако в настоящее время он не поддерживается Pyomo.

У вас есть два варианта указания таких ограничений: либо (1) указать его как дизъюнкцию с использованием расширения pyomo.gdp, а затем использовать преобразования в pyomo.gdp, чтобы ослабить дизъюнктивную программу обратно в MIP, либо (2) явно ослабить дизъюнкцию, используя, например, релаксация Big-M. Чтобы сделать первое, вам нужно будет определить два дизъюнкта, а затем дизъюнкцию:

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)

Затем вам нужно вызвать преобразование, чтобы преобразовать дизъюнкции обратно в алгебраические ограничения. Примечание: преобразования либо требуют, чтобы все ваши переменные Pyomo имели нижнюю и верхнюю границы, либо вам необходимо указать допустимые значения «Big-M» через суффикс BigM в вашей модели. Вы также можете:

• укажите преобразование в командной строке Pyomo (например, --transform=gdp.bigm или --transform=gdp.chull)

• укажите преобразование в BuildAction

  def xfrm(m):
      TransformationFactory('gdp.bigm').apply_to(m)
  model.xfrm = BuildAction(rule=xfrm)
  

• явным образом вызовите преобразование как часть сценария настраиваемого драйвера.

Альтернативой pyomo.gdp является явное выполнение расслабления самостоятельно. Вам нужно будет добавить двоичную переменную (назовем ее y), которая указывает, какая сторона дизъюнкции должна быть True, а затем явно ослабить оба дизъюнкта, используя эту двоичную переменную. Концептуально вы превратите

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

в

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

Обратите внимание, что в зависимости от границ для expr и aff некоторые из этих ограничений могут быть избыточными. Кроме того, четыре «больших М» (большие константы) не обязательно должны быть разными, хотя проблема будет решена лучше, если вы сможете определить наименьшее допустимое значение для каждого из М.

спасибо @jsiirola.

Я реализовал ограничение с помощью Big-M-Relaxation. оно работает.

Кроме того, я упростил выражения expr- и return-statement.

для каждого f (бывшего f1, f2) сумма по S и M из (model.b [s, f, 1, m] * model.x [s, f, 1, m]) должна равняться 0 OR model.af [f] (бывший model.aff [f1, f2]).

model.af [f] также используется как граница M.

## Big-M Relaxation
def P1_constraint_rule (model, f):
    expr = 0
    for (s,m) in model.index_sm:
        expr += model.b[s,f,1,m] * model.x[s,f,1,m]

    return expr <= model.af[f] * model.xf[f]

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


def P2_constraint_rule (model, f):
    expr = 0
    for (s,m) in model.index_sm:
        expr += model.b[s,f,1,m] * model.x[s,f,1,m]

    return expr >= model.af[f] * model.xf[f] 

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


def P3_constraint_rule (model, f):
    expr = 0
    for (s,m) in model.index_sm:
        expr += model.b[s,f,1,m] * model.x[s,f,1,m]

    return expr - model.af[f] <= model.af[f]  * (1- model.xf[f])

model.P3Constraint = Constraint(model.F, rule=P3_constraint_rule)


def P4_constraint_rule (model, f):
    expr = 0
    for (s,m) in model.index_sm:
        expr += model.b[s,f,1,m] * model.x[s,f,1,m]

    return model.af[f] - expr <= model.af[f] * (1- model.xf[f])

model.P4Constraint = Constraint(model.F, rule=P4_constraint_rule)