title: How to make linear programming optimization faster
body: <p>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:</p>
  
  <p>P1=[x1,x2,x3] , P2=[x4,x5,x6], P3=[x7,x8,x9]</p>
  
  <p>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).</p>
  
  <p>The formulation for this problem is therefore:
  <img src="https://imgur.com/y3rg8ND.png" alt="formulation"></p>
  
  <p>where yij=1 if we allocate customer Ci to adviser Pj, 0 otherwise</p>
  
  <p>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.</p>
  
  <p>I have a working solution using the optimization package: <code>PuLP</code> that finds the optimal allocation.</p>
  
  <pre class="lang-py prettyprint-override"><code>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] &lt;= z_max
      model += y_sum[j] &gt;= 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))
  </code></pre>
  
  <p>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</p>
  
  <p><strong>The problem</strong></p>
  
  <p>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?</p>