为最大连接集数创建 Pyomo 约束

Question

我做了什么

首先，这可能不是最好的论坛，如果是这样，请道歉。我正在创建一个 Pyomo 模型，我在其中创建了一个二进制矩阵，如下所示：

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

我的模型求解这个矩阵，典型输出如下：

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]]

结果被解释为 1 的坐标，即 (1,2),(1,4),(2,1),(2,3),(3,2),(3,6),本例中为 (4,1),(4,5),(5,4),(5,6),(6,3),(6,5)。

然后根据连接元素的组来考虑。在这种情况下，只有 1 个唯一组：(1,2,3,4,5,6)。

我需要什么

我希望通过引用 model.binMat 中的值来创建一个新的约束，以仅允许 2 个大小相同的唯一组。

这些最终组的示例如下：(1,5,6) 和 (2,3,4)。对应的坐标可以是：(1,5),(1,6),(2,3),(2,4),(3,2),(3,4),(4,2), (4,3),(5,1),(5,6),(6,1),(6,5)

我目前正在尝试使用 Pyomo 集来解决这个问题，但由于这些对我来说是新的，所以我没有运气。

编辑

对于那些对相同问题的替代方法感兴趣的人，我也发布了这个here

Answer 1

可能有一种更简单的方法，但我能想到的最好方法是添加二进制约束来检查每个可能的此类集合，并强制选择其中一组大小相同的唯一组件。请注意，这种方法会导致成倍增加的约束，因此它不是解决较大问题的好方法。

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