Allocation optimization problem using python

Question

I have a list of sectors each with certain capacity. For instance:

Sectors = { 
       "1A": 80,
       "1B": 20, 
       "2A": 10, 
       "3A": 50,
       "3B": 20,
       "3C": 110
     }

I have also a list of teams each with certain amount of members, eg.:

Teams = {
   "TeamA":20, 
   "TeamB":15, 
   "TeamC":100, 
   "TeamD":20,
   "TeamF":85, 
   "TeamG":35,
   "TeamH":25,
   "TeamI":7,
 } 

Notice that there are more Teams than Sectors and total number of team members in each team is higher than total number of sector capacity so there will be teams without a sector.

Now I need to match each of the teams to sector, so that:

1. Sum of team members in teams that have allocated sectors is maximized

For the constraints there are:

1. Team size cannot surpass allocated sector size (example: TeamC can only fit into sector 3C as a whole unless its split into two which is described in constraint 2)

Now this is harder one:

2. One Team can be allocated to more than one sector, but the sector need to start with the same number (example: its possible to allocate TeamC to 1A and 1B but not to 1A and 3B)
3. One sector can only have one team allocated to it.
4. Teams need to be either fully accommodated or fully left out.

Second constraint makes makes this problem something like multiple knapsack problem if Im thinking right, but Im not sure if a team is a knapsack and the sectors are the items that we pack into the sack. Because the item size needs to be equal or more than the sack capacity in this case.

Because of this constraint I’m quite at loss what kind of algorithm I should be using to solve this problem or if this problem is actually solvable in this form. Does anyone know as simple as possible python algorithm for solving this problem?

EDIT: ADDED CODE:

Currently I am usint ORTools library to solve this problem: Variable:

x = {}
for i in data['sector_number']:
    for j in data['teams']:
        x[(i, j)] = solver.IntVar(0, 1, 'x_%i_%i' % (i, j))

# y[j] = 1 if team j has allocated sector.
y = {}
for j in data['teams']:
    y[j] = solver.IntVar(0, 1, 'y[%i]' % j)

Here are my constraints:

# Each sector can be allocated to at most one team.
for i in data['sector_number']:
    solver.Add(sum(x[i, j] for j in data['teams']) <= 1)
# The sector seats amount allocated to each team needs to be minimum this team capacity
for j in data['teams']:
    solver.Add(sum(x[(i, j)] * data['sector_capacity'][i] for i in data['sector_number']) >= y[j] * data['team_members'][j])

And finally objective:

# Objective
objective = solver.Objective()

solver.Maximize(solver.Sum([y[j] for j in data['teams']]))

So actually the only missing constraint is the one that takes into account that if multiple sectors are allocated to one team then only sectors with same number at the beginning can be allocated there.

Here is data input for this:

{'sector_capacity': [80, 20, 10, 50, 20, 110],
 'sector_number': [0, 1, 2, 3, 4, 5],
 'sector_names': ['1A', '1B', '2A', '3A', '3B', '3C']
 'num_sectors': 6,
 'teams': [0, 1, 2, 3, 4, 5, 6, 7],
 'team_names': ['TeamA',
  'TeamB',
  'TeamC',
  'TeamD',
  'TeamE',
  'TeamF',
  'TeamG',
  'TeamH',
  'TeamI'],
 'team_members': [20, 15, 100, 20, 85, 35, 25, 7]}

Answer 1

A quick hack based on or-tools (like you tried).

Code

from collections import defaultdict
import numpy as np
from ortools.sat.python import cp_model

# INPUT
# -----
capacities = np.array([80, 20, 10, 50, 20, 110])
sector_partition = np.array([0, 0, 1, 2, 2, 2])  # 11 2 333
members = np.array([20, 15, 100, 20, 85, 35, 25, 7])

# PREPROC
# -------
n_teams = len(members)
n_sectors = len(capacities)

incompatible_sectors = defaultdict(list)  # list of incompatible sectors for each sector
for sLeft in range(n_sectors):
  for sRight in range(n_sectors):
    if sLeft < sRight:                    # symmetry
      if sector_partition[sLeft] != sector_partition[sRight]:
        incompatible_sectors[sLeft].append(sRight)
        incompatible_sectors[sRight].append(sLeft)

# MODEL
# -----
model = cp_model.CpModel()

var_assign_int = np.empty((n_teams, n_sectors), dtype=object)
var_assign_nnz_ind = np.empty((n_teams, n_sectors), dtype=object)

for t in range(n_teams):
  for s in range(n_sectors):
    # VAR assignment-matrix: member <->- team mapping : integer / amount
    var_assign_int[t, s] = model.NewIntVar(0, min(capacities[s].item(), members[t].item()), '')  # .item() because of types!
    # VAR assignment-matrix: boolean indicator of above being strictly positiv
    var_assign_nnz_ind[t, s] = model.NewBoolVar('')

for t in range(n_teams):
  for s in range(n_sectors):
    # CONSTRAINT: sector not used / not strictly positive -> no amount allowed to distribute
    model.Add(var_assign_int[t, s] == 0).OnlyEnforceIf(var_assign_nnz_ind[t, s].Not())
    # CONSTRAINT: sector used / strictly positive -> amount distributed over sectors must equal member cardinality
    model.Add(var_assign_int[t, :].sum() == members[t].item()).OnlyEnforceIf(var_assign_nnz_ind[t, s])

for s in range(n_sectors):
  # CONSTRAINT: only one team maps to this sector
  model.Add(var_assign_nnz_ind[:, s].sum() <= 1)

  for t in range(n_teams):
    # CONSTRAINT: sector used -> every other incompat sector not used!
    for incompat_s in incompatible_sectors[s]:
      model.AddImplication(var_assign_nnz_ind[t, s], var_assign_nnz_ind[t, incompat_s].Not())

# OBJECTIVE
model.Maximize(var_assign_int.sum())

# SOLVE
# -----
solver = cp_model.CpSolver()
solver.parameters.log_search_progress = True
status = solver.Solve(model)

# PRINT SOL
# ---------
vfunc = np.vectorize(lambda x: solver.Value(x), otypes=[int])
print(vfunc(var_assign_int))
print(vfunc(var_assign_nnz_ind))

Output

...
CpSolverResponse summary:
status: OPTIMAL
objective: 247
best_bound: 247
booleans: 603
conflicts: 3612
branches: 3938
propagations: 557068
integer_propagations: 70215
restarts: 424
lp_iterations: 0
walltime: 0.292556
usertime: 0.292556
deterministic_time: 0.117155
primal_integral: 0.7551

[[ 0  0  0  0 20  0]
 [ 0  0  0  0  0  0]
 [80 20  0  0  0  0]
 [ 0  0  0  0  0  0]
 [ 0  0  0  0  0 85]
 [ 0  0  0 35  0  0]
 [ 0  0  0  0  0  0]
 [ 0  0  7  0  0  0]]

 [[0 0 0 0 1 0]
  [0 0 0 0 0 0]
  [1 1 0 0 0 0]
  [0 0 0 0 0 0]
  [0 0 0 0 0 1]
  [0 0 0 1 0 0]
  [0 0 0 0 0 0]
  [0 0 1 0 0 0]]

Remarks

Correctness

Who knows... ;-)

But (after recognizing, that Magnus used different data) this code generates the same solution as the one presented by Magnus (objective 212) when his data is used.

Approach

There are tons of approaches and it's hard to select the best without a lot of information.

• This one is rather simple and somewhat dense: it's aggressively using SAT-based clauses (which or-tools solver likes)
  • We don't model sector-hierarchies in a sparse/compact form but just forbid all a-priori calculated conflicts
    • This can get very costly for bigger instances!

Answer 2

Here's another hack at this in pyomo. I peeled the sector apart into sectors and zones (the first digit) to help clarify the constraints.

# Team Assignment

import pyomo.environ as pyo

### DATA
sectors = { 
       "1A": 80,
       "1B": 20, 
       "2A": 10, 
       "3A": 50,
       "3B": 20,
       "3C": 110
     }

teams = {
   "TeamA":20, 
   "TeamB":15, 
   "TeamC":100, 
   "TeamD":20,
   "TeamF":85, 
   "TeamG":35,
   "TeamH":25,
   "TeamI":7,
 } 

sector_zones = { s:s[0] for s in sectors} # a pairing of the sector with the zone (first digit)

### SETS

m = pyo.ConcreteModel("Team Assignment")

m.S = pyo.Set(initialize=sectors.keys())
m.T = pyo.Set(initialize=teams.keys())
m.Z = pyo.Set(initialize=set(sector_zones.values()))
m.SZ = pyo.Set(within=m.S*m.Z,initialize=[(s,z) for (s,z) in sector_zones.items()])

### PARAMS

m.sector_cap = pyo.Param(m.S, initialize=sectors)
m.team_size  = pyo.Param(m.T, initialize=teams)

### VARS

# assign X players from team T to sector S in zone Z
m.X = pyo.Var(m.T, m.SZ, domain=pyo.NonNegativeIntegers)
# include team T at sector S
m.team_sector = pyo.Var(m.T, m.S, domain=pyo.Binary)
# include team T in zone Z
m.team_zone   = pyo.Var(m.T, m.Z, domain=pyo.Binary)

### OBJ

m.obj = pyo.Objective(expr=sum(m.X[t,s,z] for t in m.T for s,z in m.SZ), sense=pyo.maximize)

### CONSTRAINTS

# All-or-nothing in a particular zone
def all_or_not(m, t):
    return sum(m.X[t,s,z] for s,z in m.SZ) == sum(m.team_zone[t,z] * m.team_size[t] for z in m.Z)
m.C1 = pyo.Constraint(m.T, rule=all_or_not)

# Don't bust sector limit
def sector_lim(m, t, s, z):
    return m.X[t,s,z]  <= m.team_sector[t,s] * m.sector_cap[s]
m.C2 = pyo.Constraint(m.T, m.SZ, rule=sector_lim)

# 1 team per sector max
def one_per_sector(m, s):
    return sum(m.team_sector[t,s] for t in m.T) <= 1
m.C3 = pyo.Constraint(m.S, rule=one_per_sector)

# link sector assignment to zone
def sector_zone(m, t, z):
    return sum(m.team_sector[t,s] for s in m.S if (s,z) in m.SZ) <= \
           m.team_zone[t, z] * len(m.S)
m.C4 = pyo.Constraint(m.T, m.Z, rule=sector_zone)

# 1 zone assignment per team
def one_zone(m, t):
    return sum(m.team_zone[t,z] for z in m.Z) <= 1
m.C5 = pyo.Constraint(m.T, rule=one_zone)


### Solve
solver = pyo.SolverFactory('cbc')
res = solver.solve(m)
print(res)

#m.X.display()

assigned = 0
for idx in m.X.keys():
   val = m.X[idx].value
   if val:
      print(f'assign {val} from {idx[0]} to {idx[1]}')
      assigned += val

print(f'total assigned: {assigned}')

Yields:

assign 20.0 from TeamA to 3B
assign 80.0 from TeamC to 1A
assign 20.0 from TeamC to 1B
assign 85.0 from TeamF to 3C
assign 35.0 from TeamG to 3A
assign 7.0 from TeamI to 2A
total assigned: 247.0

Answer 3

Here is an attempt using a MILP-model implemented in MiniZinc:

int: nSectors = 5;
int: nSectorGroups = 3;
int: nTeams = 8;

set of int: SECTOR = 1..nSectors;
set of int: SECTOR_GROUP = 1..nSectorGroups;
set of int: TEAM = 1..nTeams;

array[SECTOR] of int: sectorSize = [80, 20, 10, 20, 110];
array[SECTOR] of int: sectorGroup = [1, 1, 2, 3, 3];
array[TEAM] of int: teamSize = [20, 15, 100, 20, 85, 35, 25, 7];

array[SECTOR, TEAM] of var 0..1: z; % sector is used by team? 
array[SECTOR, TEAM] of var int: x; % amount of team in sector
array[TEAM] of var 0..1: w; % team is distributed on any sector?
array[SECTOR_GROUP, TEAM] of var 0..1: g; % team is distributed on a sector group?

% maximum one team per sector
constraint forall(s in SECTOR)
    (sum(t in TEAM)(z[s,t]) <= 1);

% a team can not be split over different sector groups
constraint forall(t in TEAM)
    (sum(sg in SECTOR_GROUP)(g[sg,t]) <= 1);

% connect g and z
constraint forall(t in TEAM, s in SECTOR)
    (g[sectorGroup[s],t] >= z[s,t]);

% connect x and z
constraint forall(s in SECTOR, t in TEAM)
    (x[s,t] <= max(teamSize)*z[s,t] /\ 
     x[s,t] >= z[s,t]);

% connect w and z
constraint forall(t in TEAM)
    (nTeams*w[t] >= sum(s in SECTOR)(z[s,t]) /\ 
     w[t] <= nTeams*sum(s in SECTOR)(z[s,t]));

% team size constraint
constraint forall(t in TEAM)
    (sum(s in SECTOR)(x[s,t]) = teamSize[t]*w[t]);

% sector size constraint
constraint forall(s in SECTOR)
    (sum(t in TEAM)(x[s,t]) <= sectorSize[s]);

var int: obj = sum(x);

solve
maximize obj;

output ["obj=\(obj)\n"] ++
["w="] ++ [show(w)] ++ ["\n"] ++
["z=\n"] ++ [show2d(z)] ++
["x=\n"] ++ [show2d(x)] ++
["g=\n"] ++ [show2d(g)];

Running gives:

obj=212
w=[0, 0, 1, 1, 1, 0, 0, 1]
z=
[| 0, 0, 0, 0, 1, 0, 0, 0 |
   0, 0, 0, 0, 1, 0, 0, 0 |
   0, 0, 0, 0, 0, 0, 0, 1 |
   0, 0, 0, 1, 0, 0, 0, 0 |
   0, 0, 1, 0, 0, 0, 0, 0 |]
x=
[|   0,   0,   0,   0,  65,   0,   0,   0 |
     0,   0,   0,   0,  20,   0,   0,   0 |
     0,   0,   0,   0,   0,   0,   0,   7 |
     0,   0,   0,  20,   0,   0,   0,   0 |
     0,   0, 100,   0,   0,   0,   0,   0 |]
g=
[| 0, 0, 0, 0, 1, 0, 0, 0 |
   0, 0, 0, 0, 0, 0, 0, 1 |
   0, 0, 1, 1, 0, 0, 0, 0 |]