MIP is already exact. Let there be - $l = 26$ letters; - $m = 5$ maximum group items; - $n = 7$ as an example item count, also equal to the maximum number of groups; - $G_{i,j} \in \lbrace 0,1 \rbrace$ binary group-item assignments, $0 \le i < n$, $0 \le j < n$; - $I_{j,k} \in \lbrace 0,1 \rbrace$ binary item-letter assignments, $0 \le j < n$, $0 \le k < l$; - $L_{i,k}\in \lbrace 0,1 \rbrace$ binary group-letter assignments, $0 \le i < n$, $0 \le k < l$. Minimize $$\sum_i \sum_k L_{i,k}$$ subject to $$1=\sum_i G_{i,j} \;\forall j$$ $$m \ge \sum_j G_{i,j} \;\forall i$$ $$n L_{i,k} \ge \sum_j I_{j,k} G_{i,j} \;\forall i,k$$ Works fine: ```python from string import ascii_uppercase import scipy.sparse as sp import numpy as np from scipy.optimize import milp, Bounds, LinearConstraint items = ( 'ABC', 'CDF', 'EHI', 'AHP', 'EPZ', 'CIZ', 'ITW', ) L = len(ascii_uppercase) M = 5 # max items per group N = len(items) # Variables: NN item-group assignments; NL letter-group assignments cost = np.concatenate(( np.zeros(N*N), # group-item assignments non-optimized np.ones(N*L), # group-letter assignments minimized )) # Every item must be assigned to exactly one group item_excl = LinearConstraint( A=sp.hstack( (sp.eye(N),)*N + (sp.csc_matrix((N, N*L)),), format='csc', ), lb=np.ones(N), ub=np.ones(N), ) # Every group can have at most M items (Kronecker) group_capacity = LinearConstraint( A=sp.hstack(( sp.kron(sp.eye(N), np.ones(N)), sp.csc_matrix((N, N*L)), ), format='csc'), ub=np.full(shape=N, fill_value=M), ) # For each letter of each item, that item must have at least one group assignment that in turn has # a compatible letter assignment item_block = sp.dok_array((L, N)) letter_block = sp.dok_array((L, L)) n_letters = 0 for letter in ascii_uppercase: items_with_letter = [i for i, word in enumerate(items) if letter in word] if len(items_with_letter) > 0: item_block[n_letters, items_with_letter] = -1 letter_block[n_letters, ord(letter) - ord('A')] = N n_letters += 1 item_block = sp.block_diag((item_block[:n_letters, :],)*N, format='csc') letter_block = sp.block_diag((letter_block[:n_letters, :],)*N, format='csc') letter_assigns = LinearConstraint( A=sp.hstack((item_block, letter_block), format='csc'), lb=np.zeros(n_letters*N), ) result = milp( c=cost, integrality=np.ones(N*N + N*L), # all variables integral bounds=Bounds(lb=0, ub=1), # all variables binary constraints=(item_excl, group_capacity, letter_assigns), ) assert result.success, result.message item_assign, letter_assign = np.split(result.x, (N*N,)) item_assign = item_assign.reshape((N, N)).round().astype(int) letter_assign = letter_assign.reshape((N, L)).round().astype(int) print('Item assignments:') print(item_assign, end='\n\n') print('Letter assignments:') print(letter_assign, end='\n\n') for i_group in range(N): group_assigns = item_assign[i_group] if group_assigns.any(): group_letters = ''.join(ascii_uppercase[i] for i in letter_assign[i_group].nonzero()[0]) print(f'Group {i_group} has letters {group_letters}, items', ', '.join(items[i] for i in group_assigns.nonzero()[0])) ``` ```none Item assignments: [[1 0 1 1 1 1 0] [0 0 0 0 0 0 1] [0 0 0 0 0 0 0] [0 0 0 0 0 0 0] [0 0 0 0 0 0 0] [0 0 0 0 0 0 0] [0 1 0 0 0 0 0]] Letter assignments: [[1 1 1 0 1 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1] [0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]] Group 0 has letters ABCEHIPZ, items ABC, EHI, AHP, EPZ, CIZ Group 1 has letters ITW, items ITW Group 6 has letters CDF, items CDF ```