如何使用 Pyomo 制作指标函数？答案

【问题标题】：How can I make an indicator function with Pyomo?如何使用 Pyomo 制作指标函数？
【发布时间】：2018-08-23 23:17:32
【问题描述】：

我希望在 Pyomo 中创建一个简单的指标变量。假设我有一个变量 x，如果 x > 0，此指标函数将取值 1，否则取值为 0。

这是我尝试过的方法：

model = ConcreteModel()  
model.A = Set(initialize=[1,2,3])
model.B = Set(initialize=['J', 'K'])

model.x = Var(model.A, model.B, domain = NonNegativeIntegers)
model.ix = Var(model.A, model.B, domain = Binary)

def ix_indicator_rule(model, a, b):
    return model.ix[a, b] == int(model.x[a, b] > 0)

model.ix_constraint = Constraint(model.A, model.B,
                             rule = ix_indicator_rule)

我收到的错误消息类似于Avoid this error by using Pyomo-provided math functions，根据this link 可以在pyomo.environ 找到...但我不知道该怎么做。我试过使用validate_PositiveValues()，像这样：

def ix_indicator_rule(model, a, b):
    return model.ix[a, b] == validate_PositiveValues(model.x[a, b])

model.ix_constraint = Constraint(model.A, model.B,
                             rule = ix_indicator_rule)

没有运气。任何帮助表示赞赏！

【问题讨论】：

标签： pyomo

【解决方案1】：

您可以使用“big-M”约束来实现这一点，如下所示：

model = ConcreteModel()  
model.A = Set(initialize=[1, 2, 3])
model.B = Set(initialize=['J', 'K'])

# m must be larger than largest allowed value of x, but it should
# also be as small as possible to improve numerical stability
model.m = Param(initialize=1e9)

model.x = Var(model.A, model.B, domain=NonNegativeIntegers)
model.ix = Var(model.A, model.B, domain=Binary)

# force model.ix to be 1 if model.x > 0
def ix_indicator_rule(model, a, b):
    return model.x <= model.ix[a, b] * model.m

model.ix_constraint = Constraint(
    model.A, model.B, rule=ix_indicator_rule
)

但请注意，大 M 约束是片面的。在这个例子中，它在model.x > 0 时强制model.ix 打开，但在model.x == 0 时不强制它关闭。您可以通过将不等式翻转为model.x >= model.ix[a, b] * model.m 来实现后者（但不是前者）。但是你不能在同一个模型中同时做这两个。通常，您只需选择适合您模型的版本，例如，如果将 model.ix 设置为 1 会使您的目标函数变差，那么您将选择上面显示的版本，求解器将负责将 model.ix 设置为 @987654333尽可能@。

Pyomo 还提供可满足您需求的分离式编程功能（请参阅 here 和 here）。 cplex求解器提供indicator constraints，但不知道Pyomo是否支持。

【讨论】：

【解决方案2】：

我最终使用了 Piecewise 函数并做了这样的事情：

DOMAIN_PTS = [0,0,1,1000000000]
RANGE_PTS  = [0,0,1,1]
model.ix_constraint = Piecewise(
     model.A, model.B,
     model.ix, model.x,
     pw_pts=DOMAIN_PTS,
     pw_repn='INC',
     pw_constr_type = 'EQ',
     f_rule = RANGE_PTS,
     unbounded_domain_var = True)

def objective_rule(model):
    return sum(model.ix[a,b] for a in model.A for b in model.B)

model.objective = Objective(rule = objective_rule, sense=minimize)

好像没问题。

【讨论】：