This problem is trivial and $Q_i = 20 \, \forall i$ if $C_i \ge 0 \, \forall i$. If $C$ may be negative, then keep your $y_{ij}$ binary assignment variables, add three group value variables $G_j$, and only enforce that the group value variables are constrained to integers - the original $Q_{i}$ are easier solved as bounded continuous. Then add constraints like the following:
Error: $$\epsilon_{ij} := Q_i - G_j $$
Tolerance: $$\tau_{ij} := 20 \left( 1 - y_{ij} \right) $$
$$\epsilon_{ij} \le \tau_{ij}$$
$$-\epsilon_{ij} \le \tau_{ij}$$
import pulp
from numpy.random import default_rng
n = 50
m = 3
group_bound = 5
Qbound = 20
rand = default_rng(seed=399)
C = rand.uniform(low=-1, high=1, size=n)
Q = pulp.LpVariable.matrix(
name='Q', indices=range(n), cat=pulp.LpContinuous, lowBound=0, upBound=Qbound)
Qgroup = pulp.LpVariable.matrix(
name='Qgroup', indices=range(m), cat=pulp.LpInteger, lowBound=0, upBound=Qbound)
y = pulp.LpVariable.matrix(
name='y', indices=(range(n), range(m)), cat=pulp.LpBinary)
prob = pulp.LpProblem(name='var_grouping', sense=pulp.LpMaximize)
prob.objective = pulp.lpDot(Q, C)
for i in range(n):
prob.addConstraint(
name=f'groupexcl{i}',
constraint=pulp.lpSum(y[i]) == 1)
for j in range(m):
prob.addConstraint(
name=f'groupmin{j}',
constraint=pulp.lpSum([
y[i][j] for i in range(n)
]) >= group_bound)
for i in range(n):
for j in range(m):
# if y == 0, Q == 0
# if y == 1, Q == Qgroup
error = Q[i] - Qgroup[j]
tol = Qbound*(1 - y[i][j])
prob.addConstraint(name=f'assign_{i}_{j}_hi', constraint=error <= tol)
# Only needed if C may be negative
prob.addConstraint(name=f'assign_{i}_{j}_lo', constraint=-error <= tol)
print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal
print('Qgroup =')
print([int(Qgroup[j].value()) for j in range(m)])
print('y =')
for j in range(m):
print([int(y[i][j].value()) for i in range(n)])
print()
print('Q =')
print([int(Q[i].value()) for i in range(n)])
print()
var_grouping:
MAXIMIZE
0.3626220607951214*Q_0 + 0.42196536951115493*Q_1 + ...
SUBJECT TO
groupexcl0: y_0_0 + y_0_1 + y_0_2 = 1
groupexcl1: y_1_0 + y_1_1 + y_1_2 = 1
groupexcl2: y_2_0 + y_2_1 + y_2_2 = 1
groupexcl3: y_3_0 + y_3_1 + y_3_2 = 1
...
groupmin0: y_0_0 + y_10_0 + y_11_0 + y_12_0 + y_13_0 + y_14_0 + y_15_0
+ y_16_0 + y_17_0 + y_18_0 + y_19_0 + y_1_0 + y_20_0 + y_21_0 + y_22_0
+ y_23_0 + y_24_0 + y_25_0 + y_26_0 + y_27_0 + y_28_0 + y_29_0 + y_2_0
+ y_30_0 + y_31_0 + y_32_0 + y_33_0 + y_34_0 + y_35_0 + y_36_0 + y_37_0
+ y_38_0 + y_39_0 + y_3_0 + y_40_0 + y_41_0 + y_42_0 + y_43_0 + y_44_0
+ y_45_0 + y_46_0 + y_47_0 + y_48_0 + y_49_0 + y_4_0 + y_5_0 + y_6_0 + y_7_0
+ y_8_0 + y_9_0 >= 5
groupmin1: y_0_1 + y_10_1 + y_11_1 + y_12_1 + y_13_1 + y_14_1 + y_15_1
+ y_16_1 + y_17_1 + y_18_1 + y_19_1 + y_1_1 + y_20_1 + y_21_1 + y_22_1
+ y_23_1 + y_24_1 + y_25_1 + y_26_1 + y_27_1 + y_28_1 + y_29_1 + y_2_1
+ y_30_1 + y_31_1 + y_32_1 + y_33_1 + y_34_1 + y_35_1 + y_36_1 + y_37_1
+ y_38_1 + y_39_1 + y_3_1 + y_40_1 + y_41_1 + y_42_1 + y_43_1 + y_44_1
+ y_45_1 + y_46_1 + y_47_1 + y_48_1 + y_49_1 + y_4_1 + y_5_1 + y_6_1 + y_7_1
+ y_8_1 + y_9_1 >= 5
groupmin2: y_0_2 + y_10_2 + y_11_2 + y_12_2 + y_13_2 + y_14_2 + y_15_2
+ y_16_2 + y_17_2 + y_18_2 + y_19_2 + y_1_2 + y_20_2 + y_21_2 + y_22_2
+ y_23_2 + y_24_2 + y_25_2 + y_26_2 + y_27_2 + y_28_2 + y_29_2 + y_2_2
+ y_30_2 + y_31_2 + y_32_2 + y_33_2 + y_34_2 + y_35_2 + y_36_2 + y_37_2
+ y_38_2 + y_39_2 + y_3_2 + y_40_2 + y_41_2 + y_42_2 + y_43_2 + y_44_2
+ y_45_2 + y_46_2 + y_47_2 + y_48_2 + y_49_2 + y_4_2 + y_5_2 + y_6_2 + y_7_2
+ y_8_2 + y_9_2 >= 5
assign_0_0_hi: Q_0 - Qgroup_0 + 20 y_0_0 <= 20
assign_0_0_lo: - Q_0 + Qgroup_0 + 20 y_0_0 <= 20
assign_0_1_hi: Q_0 - Qgroup_1 + 20 y_0_1 <= 20
assign_0_1_lo: - Q_0 + Qgroup_1 + 20 y_0_1 <= 20
assign_0_2_hi: Q_0 - Qgroup_2 + 20 y_0_2 <= 20
assign_0_2_lo: - Q_0 + Qgroup_2 + 20 y_0_2 <= 20
...
VARIABLES
Q_0 <= 20 Continuous
Q_1 <= 20 Continuous
Q_10 <= 20 Continuous
Q_11 <= 20 Continuous
...
0 <= Qgroup_0 <= 20 Integer
0 <= Qgroup_1 <= 20 Integer
0 <= Qgroup_2 <= 20 Integer
0 <= y_0_0 <= 1 Integer
0 <= y_0_1 <= 1 Integer
0 <= y_0_2 <= 1 Integer
0 <= y_10_0 <= 1 Integer
0 <= y_10_1 <= 1 Integer
0 <= y_10_2 <= 1 Integer
...
At line 2 NAME MODEL
At line 3 ROWS
At line 358 COLUMNS
At line 1915 RHS
At line 2269 BOUNDS
At line 2473 ENDATA
Problem MODEL has 353 rows, 203 columns and 1200 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 207.128 - 0.01 seconds
Cgl0004I processed model has 353 rows, 203 columns (153 integer (150 of which binary)) and 1200 elements
Cbc0038I Initial state - 100 integers unsatisfied sum - 4.63636
Cbc0038I Pass 1: suminf. 0.00000 (0) obj. -207.128 iterations 129
Cbc0038I Solution found of -207.128
Cbc0038I Relaxing continuous gives -207.128
Cbc0038I Cleaned solution of -207.128
Cbc0038I Before mini branch and bound, 52 integers at bound fixed and 50 continuous
Cbc0038I Mini branch and bound did not improve solution (0.04 seconds)
Cbc0038I After 0.04 seconds - Feasibility pump exiting with objective of -207.128 - took 0.01 seconds
Cbc0012I Integer solution of -207.12836 found by feasibility pump after 0 iterations and 0 nodes (0.04 seconds)
Cbc0001I Search completed - best objective -207.1283647272854, took 0 iterations and 0 nodes (0.04 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from -207.128 to -207.128
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Result - Optimal solution found
Objective value: 207.12836473
Enumerated nodes: 0
Total iterations: 0
Time (CPU seconds): 0.05
Time (Wallclock seconds): 0.04
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.06 (Wallclock seconds): 0.06
Qgroup =
[0, 20, 20]
y =
[0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1]
[1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0]
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 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]
Q =
[20, 20, 20, 0, 20, 20, 0, 20, 20, 0, 0, 0, 20, 20, 0, 20, 20, 0, 0, 20, 0, 0, 20, 0, 20, 20, 20, 0, 20, 0, 0, 0, 0, 20, 20, 20, 0, 20, 0, 0, 0, 0, 0, 0, 20, 20, 20, 0, 0, 0]