As in the comments: as it stands, the problem is infeasible. For student 3, they could never be given a class with any teacher, because they have never had at least two classes with any teacher. You need an escape hatch for students that have a small class history.
A crude approximation to "variance" that I demonstrate below is: teachers with low load are not themselves a problem, but teachers with high above-mean load variance should have their load minimized in aggregate (this accomplishes, very roughly, the same thing).
import numpy as np
import pandas as pd
import pulp
teachers = ["t1", "t2", "t3", "t4"]
students = ["s1", "s2", "s3", "s4", "s5", "s6", "s7"]
class_history_data = {
"s1": ["t1", "t4", "t3", "t1"],
"s2": ["t3", "t2", "t4"],
"s4": ["t1", "t4", "t3", "t2", "t4"],
"s6": ["t3", "t4", "t1", "t2"],
"s7": ["t2", "t1", "t3"],
"s5": ["t1", "t3"],
"s3": ["t4"],
}
# 22 long, for every teacher-student history pair
class_history = pd.DataFrame(
[
[student, teacher]
for student, teachers in class_history_data.items()
for teacher in teachers
],
columns=['student', 'teacher'],
)
# 20-long, teacher-student history pair counts, multi-indexed
class_history['n'] = 1
class_history = (
class_history.set_index(['student', 'teacher'])
.groupby(level=['student', 'teacher'])
.sum()
)
# 28-long: for every possible teacher-student pair
combos = pd.DataFrame(
index=pd.MultiIndex.from_product((students, teachers), names=('student', 'teacher'))
)
# 28-long names by student and teacher, e.g. s1_t2
var_names = (
combos.index.get_level_values('student')
+ '_'
+ combos.index.get_level_values('teacher')
).to_series(index=combos.index)
combos['asn'] = var_names.apply(pulp.LpVariable, cat=pulp.LpBinary)
combos['history'] = class_history.n
combos['history'] = combos.history.fillna(value=0).astype(int)
# for a student, concerning their current potential teacher, and all other teachers
history_pairs = np.count_nonzero(
np.broadcast_to(
# 7x4x1
combos.history.unstack(level='teacher').values[..., np.newaxis],
# to 7x4x4 repeated
(len(students), len(teachers), len(teachers)),
)
# with self references (teacher to same teacher) removed
*(1 - np.eye(len(teachers))),
axis=1,
)
# for this student-teacher pair, the number of other teachers taught at least once
combos['history_other'] = history_pairs.ravel()
# for each student, total number of classes so far
combos['history_total'] = combos.history.groupby(level='student').transform('sum')
# A student cannot have classes with a teacher before they have at least 2 classes with other teachers
# OR IF they have not yet taken classes with at least two teachers
# Express this by row elimination and not an LP constraint
# Reduces this from 28 to 26 rows.
combos = combos[
(combos.history_other >= 2) |
(combos.history_total < 2)
]
# 4x7, with zeros where there can never be a teacher-student assignment
assignments_by_teacher = combos.asn.unstack(level='student').fillna(0)
# 4-long: deviation above mean load per teacher (deviation below mean load is ignored).
mean_load = len(students) / len(teachers)
load_dev = (
(
'load_dev_' +
assignments_by_teacher.index
.to_series(index=assignments_by_teacher.index)
)
.apply(pulp.LpVariable, cat=pulp.LpContinuous, lowBound=0)
)
# minimize total deviation above mean load
prob = pulp.LpProblem(name='class_distribution', sense=pulp.LpMinimize)
prob.objective = load_dev.sum()
# Every student must have one class
for student, group in combos.asn.groupby(level='student'):
prob.addConstraint(name=student + '_class_count', constraint=group.sum() == 1)
# The deviation above mean load has a lower bound of sum(load) - mean(load)
for teacher, load in assignments_by_teacher.sum(axis=1).items():
prob.addConstraint(
name='dev_' + teacher,
constraint=load_dev[teacher] >= load - mean_load,
)
print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal
print('Student-teacher assignments for this week:')
combos['asn'] = combos.asn.apply(pulp.LpVariable.value).astype(int)
print(combos.asn.unstack(level='teacher', fill_value=0), end='\n\n')
print('Teacher loads:')
loads = pd.DataFrame({
'load': assignments_by_teacher.sum(axis=1).apply(pulp.LpAffineExpression.value),
'deviation above mean': load_dev.apply(pulp.LpVariable.value),
})
print(loads)
class_distribution:
MINIMIZE
1*load_dev_t1 + 1*load_dev_t2 + 1*load_dev_t3 + 1*load_dev_t4 + 0
SUBJECT TO
s1_class_count: s1_t1 + s1_t2 + s1_t3 + s1_t4 = 1
s2_class_count: s2_t1 + s2_t2 + s2_t3 + s2_t4 = 1
s3_class_count: s3_t1 + s3_t2 + s3_t3 + s3_t4 = 1
s4_class_count: s4_t1 + s4_t2 + s4_t3 + s4_t4 = 1
s5_class_count: s5_t2 + s5_t4 = 1
s6_class_count: s6_t1 + s6_t2 + s6_t3 + s6_t4 = 1
s7_class_count: s7_t1 + s7_t2 + s7_t3 + s7_t4 = 1
dev_t1: load_dev_t1 - s1_t1 - s2_t1 - s3_t1 - s4_t1 - s6_t1 - s7_t1 >= -1.75
dev_t2: load_dev_t2 - s1_t2 - s2_t2 - s3_t2 - s4_t2 - s5_t2 - s6_t2 - s7_t2
>= -1.75
dev_t3: load_dev_t3 - s1_t3 - s2_t3 - s3_t3 - s4_t3 - s6_t3 - s7_t3 >= -1.75
dev_t4: load_dev_t4 - s1_t4 - s2_t4 - s3_t4 - s4_t4 - s5_t4 - s6_t4 - s7_t4
>= -1.75
VARIABLES
load_dev_t1 Continuous
load_dev_t2 Continuous
load_dev_t3 Continuous
load_dev_t4 Continuous
0 <= s1_t1 <= 1 Integer
0 <= s1_t2 <= 1 Integer
0 <= s1_t3 <= 1 Integer
0 <= s1_t4 <= 1 Integer
0 <= s2_t1 <= 1 Integer
0 <= s2_t2 <= 1 Integer
0 <= s2_t3 <= 1 Integer
0 <= s2_t4 <= 1 Integer
0 <= s3_t1 <= 1 Integer
0 <= s3_t2 <= 1 Integer
0 <= s3_t3 <= 1 Integer
0 <= s3_t4 <= 1 Integer
0 <= s4_t1 <= 1 Integer
0 <= s4_t2 <= 1 Integer
0 <= s4_t3 <= 1 Integer
0 <= s4_t4 <= 1 Integer
0 <= s5_t2 <= 1 Integer
0 <= s5_t4 <= 1 Integer
0 <= s6_t1 <= 1 Integer
0 <= s6_t2 <= 1 Integer
0 <= s6_t3 <= 1 Integer
0 <= s6_t4 <= 1 Integer
0 <= s7_t1 <= 1 Integer
0 <= s7_t2 <= 1 Integer
0 <= s7_t3 <= 1 Integer
0 <= s7_t4 <= 1 Integer
Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019
At line 2 NAME MODEL
At line 3 ROWS
At line 16 COLUMNS
At line 129 RHS
At line 141 BOUNDS
At line 168 ENDATA
Problem MODEL has 11 rows, 30 columns and 56 elements
...
Result - Optimal solution found
Objective value: 0.75000000
Enumerated nodes: 0
Total iterations: 0
Time (CPU seconds): 0.02
Time (Wallclock seconds): 0.02
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.02 (Wallclock seconds): 0.02
Student-teacher assignments for this week:
teacher t1 t2 t3 t4
student
s1 1 0 0 0
s2 0 1 0 0
s3 1 0 0 0
s4 0 0 1 0
s5 0 1 0 0
s6 0 0 1 0
s7 0 0 0 1
Teacher loads:
load deviation above mean
teacher
t1 2.0 0.25
t2 2.0 0.25
t3 2.0 0.25
t4 1.0 0.00
History penalty instead of history constraint
This does better, and is also simpler:
import pandas as pd
import pulp
teachers = pd.Index(name='teacher', data=("t1", "t2", "t3", "t4"))
students = pd.Index(name='student', data=("s1", "s2", "s3", "s4", "s5", "s6", "s7"))
class_history_data = {
"s1": ["t1", "t4", "t3", "t1"],
"s2": ["t3", "t2", "t4"],
"s3": ["t4"],
"s4": ["t1", "t4", "t3", "t2", "t4"],
"s5": ["t1", "t3"],
"s6": ["t3", "t4", "t1", "t2"],
"s7": ["t2", "t1", "t3"],
}
# 22 long, for every teacher-student history pair
class_history = pd.DataFrame(
[
[student, teacher]
for student, teachers in class_history_data.items()
for teacher in teachers
],
columns=['student', 'teacher'],
)
# 20-long, teacher-student history pair counts, multi-indexed
class_history['n'] = 1
class_history = (
class_history.set_index(['student', 'teacher'])
.groupby(level=['student', 'teacher'])
.sum()
)
# 28-long: for every possible teacher-student pair
combos = pd.DataFrame(
index=pd.MultiIndex.from_product((students, teachers))
)
# 28-long names by student and teacher, e.g. s1_t2
var_names = (
combos.index.get_level_values('student')
+ '_'
+ combos.index.get_level_values('teacher')
).to_series(index=combos.index)
combos['asn'] = var_names.apply(pulp.LpVariable, cat=pulp.LpBinary)
combos['history'] = class_history.n
combos['history'] = combos.history.fillna(value=0).astype(int)
# 4-long: deviation above mean load per teacher (deviation below mean load is ignored).
mean_load = len(students) / len(teachers)
load_dev = (
(
'load_dev_' +
teachers.astype(str)
.to_series(index=teachers)
)
.apply(pulp.LpVariable, cat=pulp.LpContinuous, lowBound=0)
)
# minimize total deviation above mean load, and minimize prior history
load_weight = 0.75
hist_weight = 1.00
prob = pulp.LpProblem(name='class_distribution', sense=pulp.LpMinimize)
prob.objective = load_weight*load_dev.sum() + hist_weight*combos.asn.dot(combos.history)
# Every student must have one class
for student, group in combos.asn.groupby(level='student'):
prob.addConstraint(name=student + '_class_count', constraint=group.sum() == 1)
# The deviation above mean load has a lower bound of sum(load) - mean(load)
for teacher, load in combos.asn.groupby('teacher').sum().items():
prob.addConstraint(
name='dev_' + teacher,
constraint=load_dev.loc[teacher] >= load - mean_load,
)
print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal
print('Student-teacher assignments for this week:')
combos['asn'] = combos.asn.apply(pulp.LpVariable.value).astype(int)
print(combos.asn.unstack(level='teacher', fill_value=0), end='\n\n')
print('Teacher loads:')
loads = pd.DataFrame({
'load': combos.asn.groupby('teacher').sum(),
'deviation above mean': load_dev.apply(pulp.LpVariable.value),
})
print(loads)
Student-teacher assignments for this week:
teacher t1 t2 t3 t4
student
s1 0 1 0 0
s2 1 0 0 0
s3 0 0 1 0
s4 1 0 0 0
s5 0 0 0 1
s6 0 1 0 0
s7 0 0 0 1
Teacher loads:
load deviation above mean
teacher
t1 2 0.25
t2 2 0.25
t3 1 0.00
t4 2 0.25
classes_data
variable stores the history of sessions that each student already had. For instance,"s7": ["t2", "t1", "t3"]
indicates thats7
had a class with teachert2
3 weeks ago, then had a class witht1
2 weeks ago, and witht3
last week. Now I need to attribute what would be his teacher this week. Ideally it should bet4
, but without compromising the overall distribution as stated above. $\endgroup$