I am trying to allocate customers Ci to financial advisers Pj. Each customer has a policy value xi. I'm assuming that the number of customers (n) allocated to each adviser is the same, and that the same customer cannot be assigned to multiple advisers. Therefore each partner will have an allocation of policy values like so:
P1=[x1,x2,x3] , P2=[x4,x5,x6], P3=[x7,x8,x9]
I am trying to find the optimal allocation to minimise dispersion in fund value between the advisers. I am defining dispersion as the difference between the adviser with the highest fund value (z_max) and the lowest fund value (z_min).
The formulation for this problem is therefore:

where yij=1 if we allocate customer Ci to adviser Pj, 0 otherwise
The first constraint says that zmax has to be greater than or equal to each policy value; since the objective function encourages smaller values of zmax, this means that zmax will equal the largest policy value. Similarly, the second constraint sets zmin equal to the smallest policy value. The third constraint says that each customer must be assigned to exactly one adviser. The fourth says that each adviser must have n customers assigned to him/her.
I have a working solution using the optimization package: PuLP that finds the optimal allocation.
import random
import pulp
import time
# DATA
n = 5 # number of customers for each financial adviser
c = 25 # number of customers
p = 5 # number of financial adviser
policy_values = random.sample(range(1, 1000000), c) # random generated policy values
# INDEXES
set_I = range(c)
set_J = range(p)
set_N = range(n)
x = {i: policy_values[i] for i in set_I} #customer policy values
y = {(i,j): random.randint(0, 1) for i in set_I for j in set_J} # allocation dummies
# DECISION VARIABLES
model = pulp.LpProblem("Allocation Model", pulp.LpMinimize)
y_sum = {}
y_vars = pulp.LpVariable.dicts('y_vars',((i,j) for i in set_I for j in set_J), lowBound=0, upBound = 1, cat=pulp.LpInteger)
z_max = pulp.LpVariable("Max Policy Value")
z_min = pulp.LpVariable("Min Policy Value")
for j in set_J:
y_sum[j] = pulp.lpSum([y_vars[i,j] * x[i] for i in set_I])
# OBJECTIVE FUNCTION
model += z_max - z_min
# CONSTRAINTS
for j in set_J:
model += pulp.lpSum([y_vars[i,j] for i in set_I]) == n
model += y_sum[j] <= z_max
model += y_sum[j] >= z_min
for i in set_I:
model += pulp.lpSum([y_vars[i,j] for j in set_J]) == 1
# SOLVE MODEL
start = time.clock()
model.solve()
print('Optimised model status: '+str(pulp.LpStatus[model.status]))
print('Time elapsed: '+str(time.clock() - start))
Note that I have implemented constraints 1 and 2 slightly differently by including an additional variable y_sum to prevent duplicating an expression with a large number of nonzero elements
The problem
The issue is that for larger values of n,p and c the model takes far too long to optimise. Is it possible to make any changes to how I've implemented the objective function/constraints to make the solution faster?
w; minimizewand usew >= 0.