# Write constraints with relationship sets in Pyomo

I would like to write the following constraint in the most compact form possible using Pyomo.

$$\sum_{o\in O}y_{n,o} = \sum_{o\in O} \sum_{n'\in LN}y_{n',o}\qquad \forall n$$

Where, the important part, is that $$LN$$ is a relationship set, or multiset, such that its elements are tuples like $$(n_1,n_2), (n_2, n_4), ...$$

In GAMS, this is easy. The code for a small example would look like:

SETS
n       /n1*n6/
LN(n,n) /(n1.n2), (n2.n4)/
O       /1,2,3/
;
alias(n,nn);

VARIABLES
y(n,o)
;

EQUATIONS
EQ01;
EQ01(n)..      sum(O, y(n,o)) = sum(o, sum(nn\$LN(n,nn), y(nn,o));


For example, for $$n=n1$$, I would want the following equation. $$y_{n_1,1} + y_{n_1,2} + y_{n_1,3} = y_{n_2,1} + y_{n_2,2} + y_{n_2,3}$$

I'm not sure of how to translate the compactness of GAMS to Pyomo (if possible). I am working with some rules for a pyomo model in the form:

model    = pyo.ConcreteModel()
model.N  = pyo.Set(initialize = ["n1","n2","n3","n4","n5", "n6"]
model.LN = pyo.Set(within = model.N*model.N, initialize = [("n1", "n2"), ("n2", "n4")]
model.O  = pyo.Set(initialize = [1,2,3])
model.y  = pyo.Var(model.N, model.O)

def _rule(model,n):
return sum(model.y[n,o] for o in model.O) == sum(sum(model.y[nn,o] for nn in model.LN) for o in model.O)
model.equation = pyo.Constraint(model.N, rule = _rule)


But it is clearly not working, since in the inner sum it does not know what nn is. The set of LN also doesn't have a way of knowing that I'm referring to the outer $$n$$ in there, I think.

I think something could be done with the advantage of Pyomo's ability to use "ifs" and other flow controls, together with ConstraintSkip for the cases where there is no $$n'$$ associated to a $$n$$. However, I would like to know if there is a more elegant, compact way of writing this, such as GAMS'.

• I would just make a python dictionary, something like set_di = {'n1': ['n2'], 'n2': ['n4'] ..... Then inside your function, you should be able to define nn = set_di[n] and have an if statement to prevent a rule if there is no remaining set. Jul 4, 2022 at 17:04
• @Daniel V., actually what you mentioned as a tuple set in GAMS, LN(n), is a subset and does not perform as a nested/tuple set. Jul 5, 2022 at 7:34
• @Andy That indeed could be a solution. Also to just perform the sums in standard for loops instead of the list comprehension manner, but I was wondering if there would be a more elegant way such as GAMS' A.Omidi: Indeed it is! I have edited it from LN(n) to LN(n,n). Thank you for catching that Jul 5, 2022 at 8:09

First create a dictionary like

Allowed={"n1":"n2", "n2":"n4"}


and then

Use this

from pyomo.environ import *
Allowed={"n1":"n2", "n2":"n4"}
model    = AbstractModel()
model.N  = Set(initialize = ["n1","n2","n3","n4"])
model.O  = Set(initialize = [1,2,3])
model.y  = Var(model.N, model.O, within=Reals)
def yrule(model,n):
if n in Allowed:
return sum(model.y[n,o] for o in model.O) ==sum(model.y[nn,o] for nn in model.N for o in model.O if nn in Allowed[n])
else:
return Constraint.Skip
model.equation = Constraint(model.N, rule = yrule)

instance = model.create_instance()
instance.equation.pprint()


I will mark the answer from Optimization team, which is completely correct, as it works for the problem given, although it does not use tuple sets. I will also put here for completion how I ended up doing it, since I found that it could be made more compact with list comprehension and without Constraint.Skip.

Apparently, I was on the right track with this:

def _rule(model,n):
return sum(model.y[n,o] for o in model.O) == sum(sum(model.y[nn,o] for nn in model.LN) for o in model.O)


But I had to slightly change it to this.

def _rule(model,n):
return sum(model.y[n,o] for o in model.O) == sum(sum(model.y[nn,o] for nn in m.N if (n,nn) in model.LN) for o in model.O)


So now it knows where the nn should come from, the set N, and also knows that only do it for the pairs that are present in the tuple set.

Thank you all who helped me!

• Dear Daniel, What if for a given n you do not have any pair ? your formulation assumes that for every n in model.N there is another nn which might nobt be always correct Jul 6, 2022 at 9:15
• Apparently, the same as GAMS, pyomo is "smart" enough to then put a 0 in there, as I wanted. Wasn't sure it was going to do that, but it looks like it does, so everything works alright. For example, if I'm not mistaken, my formulation would write for nodes n3 and onwards that their summation is equal to zero, while yours would not write the constraint for those nodes. In my case, I want the first behavior, which I think is the specific one from the mathematical formulation. But I can see the usefulness of the second behavior for other problems. Jul 6, 2022 at 10:29
• I just noticed that this may seem confusing because for this case, my formulation would put everything to zero. n3 to n6 would be zero because of the condition in the right side, which will make n2 zero because n4 is zero, and n1 zero because n3 is zero. This example does not make much sense, but I just wanted to quickly put some code together that showed what I intented to do. The real case study is not equalities, but inequalities, and there are a couple more sets, but the equation I want to reproduce is analogous to this one :) Jul 6, 2022 at 10:49