# How to formulate: each pair of elements in $A$ has one common unit in $B$

We have two sets, $$A$$ and $$B$$. Some elements of $$A$$ must be connected to some elements of $$B$$, but no element of a given set is connected to another element of the same set. (Think of a bipartite graph, and we are trying to establish the edges between the two sets.)

I have a constraint that says:

For every pair of elements of $$A$$, there is (at least/exactly) one element of $$B$$ that both elements of $$A$$ are connected to.

For now, I have formulated this constraint by introducing a binary variable $$x_{a_1,a_2,b}$$ that equals 1 if $$a_1,a_2\in A$$ are both connected to $$b\in B$$, and then adding a constraint that says:

$$\sum_{b\in B} x_{a_1,a_2,b} \ge 1 \qquad \forall a_1,a_2\in A$$

(Or replace $$\ge$$ with $$=$$ for the "exactly one" version.)

But I'd really like to avoid having these triply-indexed variables. Can anyone think of a way to enforce this constraint without the $$x_{a_1,a_2,b}$$ variables?

If you don’t require linear constraints, you can introduce (or reuse) binary variables $$y_{a.b}$$ and quadratic constraints $$\sum_{b\in B} y_{a_1,b} y_{a_2,b} \ge 1 \qquad \forall a_1,a_2\in A,$$ for which your triply-indexed formulation is a linearization. If you do need linear and are worried about the problem size, you can generate the violated constraints dynamically. Depending on your other constraints, you might also be able to use compact linearization a la Liberti.

• Would a commercial MIP solver be able to handle these constraints OK? – LarrySnyder610 Nov 4 '19 at 1:57
• Some solvers support them, and some don’t. – Rob Pratt Nov 4 '19 at 12:36

Can you consider an extended formulation? In this case you could consider variables $$x_{A',b} \in \{0,1\}$$ with value 1 iff

\begin{align} \{a,b\} \in E & \quad & \forall a \in A' \\ \{a,b\} \not\in E & \quad & \forall a \not\in A' \end{align}

(where $$E$$ is the edge set of the bipartite graph). In other words it determines which subset $$A' \subseteq A$$ is connected with vertex $$b$$.

Then, I believe, your constraint becomes: $$\sum_{b \in B} \sum_{A' \supseteq \{a_1, a_2\}} x_{A',b} = 1 \quad \forall a_1,a_2 \in A \;\; (a_1 \neq a_2)$$

and you need to ensure that

$$\sum_{A' \subseteq A} x_{A',b} = 1 \quad \forall b \in B$$

Of course enumerating the variables is not feasible, and you should resort to some column generation approach.

• Hmm. OK, I think this is possible, but I don't want to use column generation. I was hoping for a simpler set of logical constraints that I could just send to a commercial MIP solver. – LarrySnyder610 Nov 4 '19 at 2:04

For the at-least variant, it seems you can write down this problem as propositional logic and convert your constraints into Conjunctive Normal Form, which gives you the constraints you need only using the edge-variables.

Given sets $$A$$ and $$B$$ and variables $$\forall i\in A, j \in B: x_{i,j}$$, the constraints can be written as $$\bigwedge_{i,j \in A} \bigvee_{k \in B} (x_{ik} \wedge x_{jk})$$ which are not the constraints of the language of Integer Programming. However, these constraints can be converted into the Conjunctive Normal Form to obtain constraints that are easily mapped to IP constraints.

I made a short Python script that uses the boolean algebra module of Sympy to automate this process.

from sympy.core import Symbol
from sympy.logic.boolalg import *

# Define the sizes of sets A and B
a_size = 3
b_size = 2

# Define sets A and B
A = [i for i in range(a_size)]
B = [i for i in range(b_size)]

# Build a dictionary for the variables
X = {i : {j : Symbol('x{}_{}'.format(i,j)) for j in B} for i in A}

# Construct a logic expression capturing all constraints
expr = True
for i in A:
for j in A:
if i != j:
subExpr = False
for k in B:
subSubExpr = And(X[i][k],X[j][k])
subExpr = Or(subExpr, subSubExpr)
expr = And(expr, subExpr)

print('The logic expression is:')
print(expr)

# Convert the logic expression to CNF
cnf_expr = to_cnf(expr)

print('In CNF the logic expression is:')
print(cnf_expr)

# Iterate over the 'And' terms of the CNF expression
for conjunctive_term in cnf_expr.args:
# Extract variable names of the terms
varnames = [str(a) for a in conjunctive_term.args]
# Print the IP-constraint corresponding to the conjunctive term
lhs = ' + '.join(varnames)
print(lhs + ' >= 1')


Applied to a case where $$|A| = 3$$ and $$|B| = 2$$, this gives us the logical constraint $$((x_{0,0} \wedge x_{1,0}) \vee (x_{0,1} \wedge x_{1,1})) \wedge ((x_{0,0} \wedge x_{2,0}) \vee (x_{0,1} \wedge x_{2,1})) \wedge ((x_{1,0} \wedge x_{2,0}) \vee (x_{1,1} \wedge x_{2,1}))$$

converted to CNF this constraint becomes

$$(x_{0,0} \vee x_{0,1}) \wedge (x_{0,0} \vee x_{1,1}) \wedge (x_{0,0} \vee x_{2,1}) \wedge (x_{0,1} \vee x_{1,0}) \wedge (x_{0,1} \vee x_{2,0}) \wedge (x_{1,0} \vee x_{1,1}) \wedge (x_{1,0} \vee x_{2,1}) \wedge (x_{1,1} \vee x_{2,0}) \wedge (x_{2,0} \vee x_{2,1})$$

And converted to linear constraints, this gives us \begin{align*} x_{0,0} + x_{0,1} &>= 1 \\ x_{0,0} + x_{1,1} &>= 1 \\ x_{0,0} + x_{2,1} &>= 1 \\ x_{0,1} + x_{1,0} &>= 1 \\ x_{0,1} + x_{2,0} &>= 1 \\ x_{1,0} + x_{1,1} &>= 1 \\ x_{1,0} + x_{2,1} &>= 1 \\ x_{1,1} + x_{2,0} &>= 1 \\ x_{2,0} + x_{2,1} &>= 1 \end{align*}

Note that the number of constraints generated by this approach grows quite rapidly: for $$|A|=3$$ and $$|B|$$ is $$3, 4, 5, 6$$ the numbers of constraints generated by this script seem to be $$21, 45, 93, 189$$. The Online Encyclopedia of Integer Sequences seems to contain some information on this sequence which suggest some binomial-related growth, but I did not analyze this in more detail. It may thus be necessary to figure out the actual structure of the generated constraints and use a cutting plane method to solve larger instances in practice.

However, I do think does avoid the need to generate the triplet variables for the at-least variant.

• Seems like the sequence is $3+3\sum\limits_{i=1}^{|B|-1}2^i=3(2^{|B|}-1)$. It might be interesting to test this for different $|A|,|B|$ and see if there's a more general pattern. – TheSimpliFire Nov 5 '19 at 15:50
• @TheSimpliFire I ran some more checks. For $|A|=3$ and $|B|=5,6,7,8$ the sequence is $93, 189, 381, 765$. For $|A|=4$ and $|B|=5,6,7,8$ the sequence is $184, 376, 760, 1528$. For $|A|=5$ and $|B|=5,6,7,8$ the sequence is $305, 625, 1265, 2545$. For $|A|=6$ and $|B|=5,6,7,8$ the sequence is $456, 936, 1896, 3816$. So substituting the $3$ in that formula for $|A|$ is not good enough, it seems. – Paul Bouman Nov 5 '19 at 17:16
• Thanks, this is an interesting approach. I'm not sure I will wind up implementing this approach, but it was interesting to learn about it from your answer! – LarrySnyder610 Nov 5 '19 at 17:36
• I have found the formula. The number of constraints $N(|A|,|B|)$ is given by $$N=|A|+|A|(|A|-1)\left(2^{|B|-1}-1\right).$$ – TheSimpliFire Nov 5 '19 at 19:25
• I tried out some small examples which seem to suggest that these constraints can also be found by applying Fourier–Motzkin elimination to the original system. – Rolf van Lieshout Nov 6 '19 at 15:03

How about we write your constraint, just replacing the tuple of $$(a_1, a_2)$$ with, say, $$a'$$. So, you'll have two indexed variables of $$x_{a',b}$$.

What I mean is, first enumerate all the pairs of $$A$$, and then, do your sum over those (and I think you'll have a faster write-time of your constraint in python, for instance, too)

• This doesn't really fix anything, unfortunately. My objection isn't to the fact that there are 3 indices (although that's what I said). It's to the fact that the $x$ variable has $O(|A|^2|B|)$ terms, which is true whether we formulate it as $x_{a_1,a_2,b}$ or $x_{a',b}$ with $a' \in A \times A$. – LarrySnyder610 Nov 4 '19 at 2:01
• @LarrySnyder610 Then, if you don't mind too many constraints, I go with what Maarten suggested. Maybe linearize it by introducing a new binary variable $y_i$ for every pair in $A$. Then two constraint for every pair: i) $x_{a_1,b} + x_{a_2,b} \le 2y_i+1$ and ii) $x_{a_1,b} + x_{a_2,b} > 2y_i -1$ – EhsanK Nov 4 '19 at 3:47
• What is $y_i$? It equals 1 if ... ? – LarrySnyder610 Nov 4 '19 at 14:52
• It equals 1 if $x _{a1,b}+x_{a2,b} = 2$ and 0 otherwise. As I said, it's essentially Maarten constraint. – EhsanK Nov 4 '19 at 15:02
• But then it's really $y_{a1,a2,b}$, so there are still $O(|A|^2|B|)$ terms. :( – LarrySnyder610 Nov 4 '19 at 15:25

If $$B_i$$ is the set whose elements are connected by $$a_i$$, then for the "at least" condition we want $$\sum_{b_k\in B_i\cap B_j\\\quad i\ne j}|b_k|>0\qquad\forall a_i,a_j\in A$$ where $$k=1,\cdots,K$$ with $$K=|B_i\cap B_j|$$.

For the "exactly one" condition we want $$-b_1+\sum_{b_k\in B_i\cap B_j\\\quad i\ne j}|b_k|=0\qquad\forall a_i,a_j\in A.$$

• It looks like you are using $b_k$ as a set, as an element of $B$, and as a binary variable. – Rob Pratt Nov 3 '19 at 16:26
• @RobPratt $b_k$ is an element of $B_i\cap B_j$, both of which are sets. That is, $B_i\cap B_j=\{b_1,b_2,\cdots,b_K\}$. Can you elaborate where I have defined $b_k$ as the other two? – TheSimpliFire Nov 3 '19 at 16:29
• You have $|b_k|$, which I interpret as cardinality of a set. And $b_k$ looks like a (binary?) variable because otherwise the constraint has no variables. Or is $B_i$ a set-valued variable? And what does $b_1$ mean, that is, what is special about $k=1$? – Rob Pratt Nov 3 '19 at 16:35
• @RobPratt Ah, the ambiguity lies in my defining $|b_k|$ as the absolute value of $b_k$ (which disappears if $b_k$ is positive). For your second question, $B_i$ contains all the elements of $B$ that are connected to $a_i$ and is thus a set. Finally, the choice of $b_1$ is arbitrary (to an extent) -- I could have written something like $\max b_k$. Example: $B=\{1,2,3,4,5\}$, $B_1=\{1,3,4\}$, $B_2=\{2,3,4\}$ so $B_1\cap B_2=\{3,4\}$ with $b_1=3$ and $b_2=4$. – TheSimpliFire Nov 3 '19 at 16:39
• Hmm, seems like a complicated way to express $|B_i \cap B_j|>0$. – Rob Pratt Nov 3 '19 at 16:47

(edit: now I read more carefully, I see Rob Pratt already alluded to this approach, but might as well expand on it)

Another option is the lazy constraint method: solve the problem without bothering to implement the constraint, check whether the constraint is violated, and if it is add a new constraint that prevents those particular violations and repeat.

For instance, suppose I get a solution where $$a_i$$ and $$a_j$$ have no common connection. Let $$S$$ be the set of $$k$$ such that $$a_i$$ connected to $$b_k$$ in this particular (invalid) solution.

For that particular $$i$$ and $$j$$, I can then add a constraint: if $$x_{i,k}=1 \forall {k\in S}$$ and $$x_{i,k}=0 \forall {k\notin S}$$ then $$\sum_{k\in S} x_{j,k} \ge 1$$. This is very easily linearised.

Iterating on this will eventually give a solution that satisfies the "$$\ge 1$$" condition; a "$$\le 1$$" condition can be handled in a similar way if desired.

Whether this is an efficient way to do it depends on other aspects of the problem. If the other constraints and objective/s are such that solutions will mostly satisfy the condition even without being specifically constrained to do so, with only a few violations, then this might be a relatively painless way to deal with those exceptions. It might be useful since it avoids needing to create $$O(n^3)$$ extra DVs.

OTOH, if there are a lot of violations of this rule, then you may end up iterating for a long time. Especially if solution time per iteration is long, lazy constraints might not be your best bet.