使用 python 的分配优化问题

Allocation optimization problem using python

我有一个扇区列表,每个扇区都有一定的容量。例如:

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,
 } 

请注意,团队数量多于扇区,并且每个团队中的团队成员总数高于扇区容量总数,因此会有没有扇区的团队。

现在我需要将每个团队与扇区相匹配,以便:

  1. 已分配扇区的团队成员总数最大化

对于约束有:

  1. 团队规模不能超过分配的扇区大小(例如:TeamC 只能作为一个整体放入扇区 3C,除非它分成两个,这在约束 2 中有描述)

现在这个更难了:

  1. 一个团队可以分配到多个扇区,但扇区需要以相同的编号开始(例如:可以将TeamC分配到1A和1B,但不能分配到1A和3B)
  2. 一个扇区只能分配一个团队。
  3. 团队需要完全适应或完全被排除在外。

第二个约束使得这个问题有点像多背包问题,如果我的想法是正确的,但我不确定一个团队是否是一个背包,并且扇区是我们装进背包的物品。因为在这种情况下,物品尺寸需要等于或大于麻袋容量。

由于这个限制,我很困惑我应该使用什么样的算法来解决这个问题,或者这个问题是否真的可以用这种形式解决。 有谁知道解决这个问题的尽可能简单的python算法?

编辑:添加代码:

目前我正在使用ORTools库来解决这个问题: 变量:

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)

这是我的约束条件:

# 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])

最后 objective:

# Objective
objective = solver.Objective()

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

所以实际上唯一缺少的约束是考虑到如果多个扇区被分配给一个团队,那么只有开始时编号相同的扇区才能被分配到那里。

这是为此输入的数据:

{'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]}

这是使用 MiniZinc 中实现的 MILP 模型的尝试:

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)];

运行 给出:

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 |]

基于 or-tools 的快速破解(如您所试)。

代码

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))

输出

...
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]]

备注

正确性

谁知道呢……;-)

但是(在认识到 Magnus 使用了不同的数据之后)当他的数据是用过。

方法

有很多方法,如果没有大量信息,很难 select 找到最好的方法。

  • 这个相当简单而且有点密集:它积极使用基于 SAT 的子句(or-tools 求解器喜欢)
    • 我们不以 sparse/compact 形式对扇区层次结构建模,只是禁止所有先验计算的冲突
      • 对于更大的实例,这可能会变得非常昂贵!

这是 pyomo 中的另一个技巧。我将扇区拆分为扇区和区域(第一个数字)以帮助阐明约束。

# 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}')

产量:

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