How to conditionally add variables to a Pyomo constraint based on predetermined data settings

I'm new to Pyomo.

I would like to know if there is an elegant way to code a constraint that contains variables that we may or may not want to include. The option to include them will be known at the time the model is solved and will be based on settings it reads in from the database.

A good example might be having slacks in the constraint, sometimes we want them and sometimes we don't.

I have tapped out a small demo below for a warehouse location example. In equation buildLimit I have added a slack to allow the number of warehouses to exceed the build limit [the code might contain some syntax errors I haven't run it]

# Import pyomo
from pyomo.environ import *
from pyomo.opt import SolverStatus, TerminationCondition

N = ['Harlingen', 'Memphis', 'Ashland']
M = ['NYC', 'LA', 'Chicago', 'Houston']

d = {('Harlingen', 'NYC'): 1956, \
     ('Harlingen', 'LA'): 1606, \
     ('Harlingen', 'Chicago'): 1410, \
     ('Harlingen', 'Houston'): 330, \
     ('Memphis', 'NYC'): 1096, \
     ('Memphis', 'LA'): 1792, \
     ('Memphis', 'Chicago'): 531, \
     ('Memphis', 'Houston'): 567, \
     ('Ashland', 'NYC'): 485, \
     ('Ashland', 'LA'): 2322, \
     ('Ashland', 'Chicago'): 324, \
     ('Ashland', 'Houston'): 1236 }

P = 2

model = ConcreteModel("warehouse location problem")

model.N = Set(dimen=1, initialize=N)
model.M = Set(dimen=1, initialize=M)
model.d = Param(model.N, model.M, within=PositiveIntegers, initialize=d)
model.P = Param(initialize=P)
model.y = Var(model.N, within=Binary)
model.x = Var(model.N, model.M, bounds=(0,1))
##########################
model.buildLimitSlack = Var(within=NonNegativeIntegers)
model.useSlacks = Param() # assume some data read will populate this at some stage before the solve
##########################

# Objective, minimise delivery costs
def obj_rule(model):
    return sum(model.d[n,m] * model.x[n,m] for n in model.N for m in model.M) + 99*model.buildLimitSlack
model.obj = Objective(rule=obj_rule)
# All customer demand must be met
def demand_rule(model, m):
    return sum(model.x[n,m] for n in model.N) == 1
model.demand = Constraint(model.M, rule=demand_rule)
# Can only ship from a warehouse if that warehouse is built
def supplyOnlyIfBuilt_rule(model, m, n):
    return model.x[n,m] <= model.y[n]
model.supplyOnlyIfBuilt = Constraint(model.M, model.N, rule=supplyOnlyIfBuilt_rule)
##############################
#### WE WANT THE SLACK IN THIS EQUATION TO BE OPTIONAL BASED ON DATA SETTINGS
def buildLimit_rule(model):
    return sum(model.y[n] for n in model.N) <= model.P + model.buildLimitSlack
model.buildLimit = Constraint(rule=buildLimit_rule)
##############################

I imagine we could have an if statement in the constraint, something like the following. But we don't want that as our model equations will likely have many such optional variables in the same constraint and I don't want to have tons of nested if statements [unless there is a nice way to do this?].

def buildLimit_rule(model):
  if model.useSlacks:    
    return sum(model.y[n] for n in model.N) <= model.P + model.buildLimitSlack
  else:
    return sum(model.y[n] for n in model.N) <= model.P
model.buildLimit = Constraint(rule=buildLimit_rule)

Any advice?

Thanks in advance

Solution

(The question is quite general, so I'll propose some approaches but it's up to you to decide which ones work best in which situation.)

EDIT 1 - sparse components: for sparsity-related issues, see this question: Create a variable with sparse index in pyomo Essentially, you should initialize your sets so that they only contain the "valid" values, and that will automatically ensure that you are only initializing the required components.

EDIT 2 - optional components: there are several ways to add optional components to the model. A simple if statement is a strategy, or you could look into BuildActions. Another option is to create Expressions (combinations of parameters and variables), and "optionally" change those instead of the constraint definition itself.

If all your constraints have that form (i.e. if 0 is the "neutral" state for your variables), what you can do is re-write them as:

def buildLimit_rule(m):  
    return sum(m.y[n] for n in m.N) <= m.P + m.buildLimitSlack * m.useSlacks

If you have more than one slack variable with independent on/off parameters:

def buildLimit_rule(m):  
    return sum(m.y[n] for n in m.N) <= m.P + m.v1 * m.use_v1 + m.v2 * m.use_v2 + ...

The most general (but less readable) approach is to index all the potential slack variables:

model.SLACKS = Set(...) # set for slack variables
model.slack = Var(model.SLACKS)
model.use_slack = Param(model.SLACKS, within=Binary)

def buildLimit_rule(m):  
    return sum(m.y[n] for n in m.N) <= m.P + sum(m.slack[i] * m.use_slack[i] for i in m.SLACKS)

One other way is to fix the variables you don't need:

if not m.useSlacks:
    m.Slack.fix(0)
    m.OtherSlack.fix(0)