Я хочу создать простую индикаторную переменную в Pyomo. Предполагая, что у меня есть переменная x, эта индикаторная функция примет значение 1, если x> 0, и 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
, которые, согласно эта ссылка, находятся в 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)
без везения. Любая помощь приветствуется!
В итоге я использовал функцию 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)
Вроде нормально работает.
Этого можно добиться с помощью ограничения «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
)
Но обратите внимание, что ограничение big-M одностороннее. В этом примере он принудительно включает model.ix
, когда model.x > 0
, но не выключает его, когда model.x == 0
. Вы можете достичь последнего (но не первого), перевернув неравенство на model.x >= model.ix[a, b] * model.m
. Но вы не можете сделать и то, и другое в одной и той же модели. Как правило, вы просто выбираете версию, которая подходит вашей модели, например, если установка model.ix
на 1
ухудшает вашу целевую функцию, вы должны выбрать версию, показанную выше, и решатель позаботится о настройке model.ix
на 0
, когда это возможно.
Pyomo также предлагает функции дизъюнктивного программирования (см. здесь и здесь), которые могут удовлетворить ваши потребности. И решатель cplex предлагает индикаторные ограничения, но я не знаю, поддерживает ли их Pyomo.