Create Pyomo constraint for maximum number of connected sets

Question

What I have done

Firstly, this may not be the best forum, so apologies if that's the case. I am creating a Pyomo model, where I have created a binary matrix as follows:

model.binMat = Var(range(6),range(6),domain=Binary)

My model solves for this matrix, with a typical output like this:

binaryMatrix =  [[0 1 0 1 0 0]
                 [1 0 1 0 0 0]
                 [0 1 0 0 0 1]
                 [1 0 0 0 1 0]
                 [0 0 0 1 0 1]
                 [0 0 1 0 1 0]]

Results are interpreted as the coordinates of the 1's, i.e. (1,2),(1,4),(2,1),(2,3),(3,2),(3,6),(4,1),(4,5),(5,4),(5,6),(6,3),(6,5) in this example.

This is then thought of in terms of groups of connected elements. In this case, there would only be 1 unique group: (1,2,3,4,5,6).

What I need

I would like help to create a new constraint to only allow 2 unique groups that are equally sized by referencing the values in model.binMat.

An example of what these final groups could look like is: (1,5,6) and (2,3,4). The corresponding coordinates for this could be: (1,5),(1,6),(2,3),(2,4),(3,2),(3,4),(4,2),(4,3),(5,1),(5,6),(6,1),(6,5)

I am currently attempting to solve this using Pyomo sets, but as these are new to me, I haven't had any luck.

Edit

For those interested in alternative approaches to the same problem, I also posted this here

Answer 1

There may be a simpler way, but the best way I could think of is to add binary constraints to check each possible such set and force one of those sets of equally sized unique components to be chosen. Note, this approach results in an exponential number of constraints so it's not a good solution for larger problems.

import pyomo.environ as pyo
import itertools

nodes = set(range(6))
# the possible sets of components of length 3
full_comp_list = [(set(i),nodes-set(i)) for i in itertools.combinations(nodes,3)]
# only take the first half because it's symmetric with six nodes and equal size
comp_list = full_comp_list[:int(len(full_comp_list)/2)]

num_poss_component_sets = len(comp_list)

#%% Build model
model = pyo.ConcreteModel()
model.binMat = pyo.Var(nodes,nodes,domain=pyo.Binary)

#%% Additional Variables
# binaries to track if each component connected
model.comp1_connected= pyo.Var(range(num_poss_component_sets),within=pyo.Binary)
model.comp2_connected= pyo.Var(range(num_poss_component_sets),within=pyo.Binary)
# tracks if the two components are disjoint
model.comps_disjoint = pyo.Var(range(num_poss_component_sets),within=pyo.Binary)
# tracks if the criteria met for any set of components
model.meet_criteria = pyo.Var(range(num_poss_component_sets),within=pyo.Binary)

#%% Additional constraints
def is_comp1_connected_rule(model,comp_num):
    ''' The component is complete iff the number of (directed) edges between ==6 (all three undirected edges selected)'''
    return(sum(model.binMat[i,j] for i,j in itertools.combinations(comp_list[comp_num][0],2))
    >=3*model.comp1_connected[comp_num])
   
model.comp1_connected_constraint = pyo.Constraint(range(num_poss_component_sets),
                                                  rule=is_comp1_connected_rule)

# Check if each component set is a complete graph
def is_comp2_connected_rule(model,comp_num):
    ''' The component is complete iff the number of (directed) edges between == 6 (all three undirected edges selected)'''
    return(sum(model.binMat[i,j] for i,j in itertools.combinations(comp_list[comp_num][1],2))
    >= 3*model.comp2_connected[comp_num])
   
model.comp2_connected_constraint = pyo.Constraint(range(num_poss_component_sets),
                                                  rule=is_comp2_connected_rule)

# Check if components are separate from each other (no edges between)
def are_both_disjoint_rule(model,comp_num):
    '''Disjoint if no edges between any nodes in different component
    If there are ANY edges between, then not disjoint (model.both_comps_connected must be 0)
    '''
    return(sum([model.binMat[i,j] for i in comp_list[comp_num][0] for j in comp_list[comp_num][1]])
    <= 9 * (1-model.comps_disjoint[comp_num]))
   
model.comps_disjoint_constraint = pyo.Constraint(range(num_poss_component_sets),
                                                      rule=are_both_disjoint_rule)

# Determines if a given set of components meet the rule
def meet_criteria_rule(model,comp_num):
    '''Rule met if both components are connected and separate from each other'''
    return(model.comp1_connected[comp_num] + model.comp2_connected[comp_num]
    + model.comps_disjoint[comp_num] >= 3 * model.meet_criteria[comp_num])
   
model.comp_meets_criteria_constraint = pyo.Constraint(range(num_poss_component_sets),
                                                rule=meet_criteria_rule)

# at least one component must meet rule that theyre separate connected components
model.must_meet_criteria_constraint = pyo.Constraint(expr = sum(model.meet_criteria[comp_num]
for comp_num in range(num_poss_component_sets)) >= 1)

### New constraint to make adjacency matrix symmetric (binMat_{i,j} == binMat_{j,i})
def edges_symmetric_rule(model,node1,node2):
    '''Rule requiring both directions for edges to be the same'''
    return(model.binMat[node1,node2] == model.binMat[node2,node1])
model.edges_symmetric_constraint = pyo.Constraint(nodes,nodes,rule=edges_symmetric_rule)

#%% Add objective and solve
des_edges = [(4,0),(1,2),(1,3),(2,1),(2,3),(3,1),(3,2)]
pos_c_dict = {e:1 for e in des_edges}
c = [[pos_c_dict.get((i,j),-1) for i in nodes] for j in nodes]
model.obj = pyo.Objective(expr = sum([c[i][j]*model.binMat[i,j] for i in nodes for j in nodes]),
                          sense=pyo.maximize)

solver = pyo.SolverFactory('glpk')
res = solver.solve(model)

# get the components and the index for what's chosen
[comp_list[i] for i in range(len(comp_list)) if pyo.value(model.meet_criteria[i])]
# [({0, 4, 5}, {1, 2, 3})]
[i for i in range(len(comp_list)) if pyo.value(model.meet_criteria[i])]
# 9

# View the final binMat
final_binMat = pd.DataFrame({'source':list(nodes)*len(nodes),
                             'target':[j for i in nodes for j in [i]*len(nodes)]})
final_binMat['value'] = [pyo.value(model.binMat[i,j]) for i,j in final_binMat.values]
final_binMat['cost'] = [c[i][j] for i,j in final_binMat[['source','target']].values]
final_binMat_wide = pd.pivot(data=final_binMat,index='source',columns='target',values='value')

# target    0    1    2    3    4    5
# source                              
# 0       0.0  0.0  0.0  0.0  1.0  1.0
# 1       0.0  0.0  1.0  1.0  0.0  0.0
# 2       0.0  1.0  0.0  1.0  0.0  0.0
# 3       0.0  1.0  1.0  0.0  0.0  0.0
# 4       1.0  0.0  0.0  0.0  0.0  1.0
# 5       1.0  0.0  0.0  0.0  1.0  0.0