# My Professor couldn't complete the model for this optimization problem. how do i model this problem?

(Edit: there was a slight translation error, to be clear, we tried this since the start with Binary variables (IP), even then we couldn't crack it)

I'm an undergraduate in Industrial Engineering and taking an OR class. I've taken an interest in Linear Programming and brought a real problem to my professor to see if it could be solved (This wasn't an assignment)

My Professor is considered a massive reference in OR, with many publications, so i thought he'd be able to help me out with the mathematical model easily, but after two hours in his office we couldn't really finish the model. (Perhaps that is normal) At one point he indicated that this model would be too complex/large/too many variables with our approach, so i wanted to give it another go with you guys and see if it really can't be done.

Context: The problem comes from a student organization that wishes to optimize their member's schedules considering:

1 - (Soft restriction, penalize the max function) Each member's personal preference for working in a given time-slot (and also the time-slots that they CANNOT work in because they have a class)

2 - (Hard restriction, can't overcome this one) The time-slots that they CANNOT work in because they have a class

3 - The degree to which members are working alongside - or, at the same time as - other members belonging to the same department/project team (This would reward the max function)

Example of a member's schedule:

The organization allocates what are called “credits” in specific time slots. A time-slot and a credit each lasts 50 minutes, from monday to friday. We work with a 24-hour clock. The first time slot is at 07:30 and the following ones are at 8:20, 9:10, 10:10, 11:00, 13:30, 14:20. 15:10, 16:20 and 17:10.

These are the constraints:

• Credits can only be allocated in pairs, triplets and quadruplets (that means that for any particular day of the week, if a credit is to be allocated, there MUST be a credit before AND/OR after that one)

• There cannot be 5 credits allocated in a row

• When there are allocated credits in 11:00 AND 13:30, they are NOT considered to be in a row (for example, a worker can have credits allocated in 9:10, 10:10, 11:00, 13:30 AND 14:20 in the same day)

• No credit may be allocated in a timeslot that a worker has classified with the letter F in the preference matrix.

• A worker must have exactly 12 credits allocated in the week.

• The schedule for a workweek cannot be changed after it is defined (The schedule for week 1 is the same as week [x])

The inputs that we are going to submit to the program are one preference matrix for each worker, which has 10 lines (each time slot) and 5 columns (each business day of the week). The worker will specify their preference using the letters A, B, C or F, where A means “It is ideal to allocate a credit in this time slot”, B means “It is not preferable to allocate a credit in this time slot”, C means “it is very much not preferable to allocate a credit in this time slot” and F means “it is impossible to allocate a credit in this time slot”.

Another input will be a list of groups and their respective members. There are two categories of groups: Each worker is part of one "project group" AND one "operational group."

That way, one of the objectives of the program is to build the schedule in a way that workers are the most satisfied according to the preference matrix they submitted and the schedule built by the program. It is ideal that workers who share a group have credits allocated in the same time slots as their fellow group members, so the main objective is to maximize the occurrence of such event.

How would you construct a model (MIP, LP...) for this?

(Edit: there was a slight translation error, to be clear, we tried this since the start with Binary variables (IP), even then we couldn't crack it)

There are at least two approaches that can be used to model this. The more obvious one uses binary variables $$x_{wsd}$$ equal to 1 if and only if worker $$w$$ is scheduled in slot $$s$$ on day $$d.$$ Some of the constraints (such as 12 credits per week) will be straightforward, and some (such as the limits on consecutive credits) will require auxiliary binary variables. If you can express the constraints as logical expressions in conjunction normal form, you can turn them into inequalities. Rewarding team members working side by side would involve finding $$s, d$$ pairs where $$x_{wsd} + x_{w'sd}=2$$ for workers $$w, w'$$ on the same team.

The less obvious approach is to first construct the set of feasible schedules for an individual on a given day. Feasible here would mean satisfying the pairs, triplets, quadruplets, not five in a row stuff. With 10 slots in a day, there would be ~1,000 schedules, but after restricting to feasible ones there would be somewhat fewer. The initial set of binary variables would now be $$x_{wsd}=1$$ if worker $$w$$ is assigned to schedule $$s$$ on day $$d.$$ Workers could only be assigned one schedule in a day, and only a schedule feasible for that worker (i.e., $$x_{wsd}=0$$ if worker $$w$$ has a class that would conflict with schedule $$s$$). You could precompute the overlap between two schedules (number of times both schedules would have someone working) and use that to compute the team overlap reward (which would involve the product of two binary variables, which is straightforward to linearize).

I'm pretty sure the second approach would involve more variables, although it's a bit hard to say without writing out the full models (since both involve a gaggle of auxiliary variables). The ability to solve either one would depend on how many workers are involved and what quality solver you are using.

Without the group score it's a very easy problem. With the group score it becomes (for an open-source solver) a very difficult problem, even though the number of auxiliary variables isn't very high.

This approach makes a compromise where the optimiser attempts to maximise the number of co-working slots with at least two members present.

With

• $$i$$ the worker index;
• $$1 \le j \lt 6$$ the day index;
• $$0 \le k \lt 10$$ the slot index;
• $$l$$ the project index;
• $$0 \le p_{ijk} < 3$$ the slot cost, known;
• $$0 < \mu$$ the group score weight, known;
• $$a_{ijk} \in \lbrace 0, 1 \rbrace$$ the slot assignment, an auxiliary variable;
• $$h_{zijk} \in \lbrace {0, 1 \rbrace}$$ the pattern assignment decision variables;
• $$q_l$$ the project group score, an auxiliary variable;
• $$r_l$$ the operational group score, an auxiliary variable;
• $$s_{jklq}$$ and $$s_{jklr} \in I$$ the floor-sum group slot population auxiliary variable;
• $$t_{zk} \in \lbrace {0, 1 \rbrace}$$ the pattern vectors, known;

minimise

$$\min \sum_i \sum_j \sum_k p_{ijk} a_{ijk} - \mu \sum_l \left( q_l + r_l \right)$$

s.t. pattern constraints

$$a_{ijk} = \sum_z t_{z,k} h_{zijk} \quad\forall i,j; 0 \le k \lt 5$$ $$a_{ijk} = \sum_z t_{z,k-5} h_{zijk} \quad\forall i,j; 5 \le k \lt 10$$ $$1 = \sum_z h_{zijk} \quad\forall i,j,k$$

credit constraints

$$a_{ijk} = 0 \quad\forall \text{ 'F' values of p}$$ $$12 = \sum_j \sum_k a_{ijk} \quad\forall i$$

group constraints

$$q_l \le \frac 1 2 \left( 0.5+\sum_i a_{ijkl} \right) \quad\forall j,k,l$$ $$r_l \le \frac 1 2 \left( 0.5+\sum_i a_{ijkl} \right) \quad\forall j,k,l$$

The implicit floor function on $$q,r$$ has the effect of not contributing to the score until the sum of credits in that slot is at least 2.

import numpy as np
import pandas as pd
import pulp

credits_per_worker = 12

# Problem dimensions by index
project_groups = pd.RangeIndex(3, name='project_group')
operational_groups = pd.RangeIndex(3, name='operational_group')
weekdays = pd.RangeIndex(start=1, stop=6, name='weekday')
slots = pd.TimedeltaIndex(
name='slot', dtype='timedelta64[s]',
data=(
'07:30:00', '08:20:00', '09:10:00', '10:10:00', '11:00:00',
'13:30:00', '14:20:00', '15:10:00', '16:20:00', '17:10:00',
),
)

# Synthetic group memberships per worker
rand = np.random.default_rng(seed=0)
workers = pd.DataFrame(
index=pd.RangeIndex(10, name='worker'),
data={
'project_group': rand.choice(a=project_groups.size, size=10),
'operational_group': rand.choice(a=operational_groups.size, size=10),
},
)

# Cartesian product of worker-slot data, starting with the index
worker_slot_index = pd.MultiIndex.from_product((
workers.index, weekdays, slots,
))
n = len(worker_slot_index)

worker_slots = pd.DataFrame(
index=worker_slot_index,
data={'preference': 'B'},
)

# synthetic preference densities
prefs = worker_slots['preference']
prefs[:n*10//100] = 'F'
prefs[n*10//100: n*30//100] = 'C'
prefs[n*30//100: n*50//100] = 'A'
y = np.arange(n)
rand.shuffle(y)
worker_slots['preference'] = prefs.values[y]

# costs based on worker prefs
worker_slots['pref_cost'] = 6  # placeholder for infeasible
worker_slots.loc[worker_slots['preference'] == 'A', 'pref_cost'] = 0
worker_slots.loc[worker_slots['preference'] == 'B', 'pref_cost'] = 1
worker_slots.loc[worker_slots['preference'] == 'C', 'pref_cost'] = 2

# individual slot assignment variables
worker_slots['assign'] = object()  # to be replaced with LpAffineExpression

prob = pulp.LpProblem(name='schedule', sense=pulp.LpMinimize)

# Within each 5-hour morning or afternoon block, there must be exactly one of the following patterns:
patterns = (
(0,0,0,0,0),
(1,1,0,0,0),
(0,1,1,0,0),
(0,0,1,1,0),
(0,0,0,1,1),
(1,1,0,1,1),
(1,1,1,0,0),
(0,1,1,1,0),
(0,0,1,1,1),
(1,1,1,1,0),
(0,1,1,1,1),
)
for (worker, weekday), _ in worker_slots['assign'].groupby(['worker', 'weekday']):
for offset in (0, 5):
name = f'pat_w{worker}_d{weekday}_o{offset}'
pattern_options = pulp.LpVariable.matrix(
name=name + '_%02d', cat=pulp.LpBinary,
indices=range(len(patterns)),
)
name=name, constraint=pulp.lpSum(pattern_options) == 1,
)
assigns = [
pulp.lpDot(
pattern_options,
[val[islot] for val in patterns],
)
for islot in range(5)
]
worker_slots.loc[
(worker, weekday, slots[offset: offset + 5]),
'assign',
] = assigns

# Exclude 'F' slots
worker_slots['preference'] == 'F',
'assign',
].items():
name=f'unavailable_w{worker}_d{weekday}_s{slot}',
constraint=0 == assign,
)

# Aggregate some summaries to each worker
worker_slots['real_pref_cost'] = worker_slots['assign'] * worker_slots['pref_cost']
workers[['credits', 'pref_cost']] = worker_slots[[
'assign', 'real_pref_cost',
]].groupby('worker').agg(pulp.lpSum)

# Credits per worker
for worker, credits in workers['credits'].items():
name=f'credits_w{worker}', constraint=credits_per_worker == credits,
)

# broadcast group memberships into cartesian-product data
worker_slots[[
'project_group', 'operational_group',
]] = workers.align(worker_slots)[0][['project_group', 'operational_group']]

'''
For each group-slot pair, there is a sum m, the number of assigned workers.
Maximize a new variable c, the minimum of assigned workers for a group over all slots.
We use the minimum and not the maximum because the maximum would optimise one high-population
slot while ignoring all others; the minimum will try hard to pair up every worker with at least one
other worker, and ideally converge on a few large groups.
'''

minima_list = [[] for _ in range(project_groups.size)]
for group_idx in (project_groups, operational_groups):
for group_val, worker_assigns in worker_slots.groupby(group_idx.name)['assign']:

# for each weekday and slot, the total number of credit assignments
credits_series = worker_assigns.groupby(['weekday', 'slot']).agg(pulp.lpSum)
suffix = f'_{group_idx.name}{group_val}'
minimum = 0  # pulp.LpVariable(name='min' + suffix, cat=pulp.LpContinuous)

for (weekday, slot), credits in credits_series.items():
slot_name = f'{suffix}_d{weekday}_s{slots.searchsorted(slot)}'
is_pair = pulp.LpVariable(
name='pair' + slot_name,
cat=pulp.LpInteger,
)
name='pair' + slot_name,
constraint=is_pair <= (credits + 0.5)/2,
)
minimum = is_pair + minimum

# Very difficult problem
# has_any = pulp.LpVariable(name='hasany' + slot_name, cat=pulp.LpBinary)
# M = 2*len(workers)
#     name='hasany' + slot_name,
#     constraint=has_any >= (credits - 0.5)/M,
# )
#
# # Either the minimum is less than credits, or has_any is 0
#     name='min' + slot_name,
#     constraint=minimum <= credits + M*(1 - has_any),
# )

minima_list[group_val].append(minimum)

group_minima = pd.DataFrame(
index=pd.RangeIndex(project_groups.size, name='group_value'),
columns=pd.Index(name='group_kind', data=('project_group', 'operational_group')),
data=minima_list,
).fillna(0)

group_weight = 0.5
prob.setObjective(
pulp.lpDot(
worker_slots['assign'], worker_slots['pref_cost'],
)
- group_weight*group_minima.values.sum()
)

prob.solve()
assert prob.status == pulp.LpStatusOptimal

worker_slots['assign'] = worker_slots['assign'].apply(pulp.value)
workers[['credits', 'pref_cost']] = workers[['credits', 'pref_cost']].map(pulp.value)
schedule = worker_slots['assign'].astype(int).unstack('slot')
group_minima = group_minima.map(pulp.value)

print(workers)
print(group_minima)

for group_idx in (project_groups, operational_groups):
assigns = worker_slots.groupby([
group_idx.name, 'weekday', 'slot',
])['assign'].sum().unstack('weekday')
print(assigns)


With this formulation, for ten students, three project groups and two operational groups, it runs like:

At line 2 NAME          MODEL
At line 3 ROWS
At line 415 COLUMNS
At line 12786 RHS
At line 13197 BOUNDS
At line 14548 ENDATA
Problem MODEL has 410 rows, 1350 columns and 8434 elements
...

Result - Optimal solution found

Objective value:                8.00000000
Enumerated nodes:               1943
Total iterations:               160190
Time (CPU seconds):             29.22
Time (Wallclock seconds):       29.22


and does produce (seemingly) reasonable results:

Workers, their group assignment, total credits assigned, and final assignment cost:

        project_group  operational_group  credits  pref_cost
worker
0                   2                  1     12.0        6.0
1                   1                  2     12.0        5.0
2                   1                  1     12.0        5.0
3                   0                  1     12.0        5.0
4                   0                  2     12.0       10.0
5                   0                  2     12.0        6.0
6                   0                  1     12.0        4.0
7                   0                  1     12.0        4.0
8                   0                  1     12.0        5.0
9                   2                  2     12.0        7.0


Integer scores assigned to each project:

group_kind   project_group  operational_group
group_value
0                     31.0                0.0
1                      5.0               33.0
2                      8.0               21.0


Credits on the project schedule:

weekday                          1    2    3    4    5
project_group slot
0             0 days 07:30:00  0.0  0.0  2.0  2.0  2.0
0 days 08:20:00  1.0  0.0  4.0  2.0  2.0
0 days 09:10:00  2.0  2.0  2.0  0.0  1.0
0 days 10:10:00  1.0  2.0  0.0  0.0  2.0
0 days 11:00:00  0.0  0.0  0.0  0.0  1.0
0 days 13:30:00  2.0  1.0  0.0  4.0  0.0
0 days 14:20:00  2.0  2.0  2.0  4.0  0.0
0 days 15:10:00  2.0  2.0  3.0  4.0  2.0
0 days 16:20:00  3.0  2.0  1.0  3.0  2.0
0 days 17:10:00  2.0  0.0  0.0  1.0  0.0
1             0 days 07:30:00  1.0  0.0  0.0  2.0  0.0
0 days 08:20:00  1.0  0.0  0.0  2.0  0.0
0 days 09:10:00  2.0  0.0  0.0  1.0  1.0
0 days 10:10:00  1.0  0.0  0.0  0.0  1.0
0 days 11:00:00  0.0  0.0  0.0  0.0  0.0
0 days 13:30:00  0.0  0.0  0.0  0.0  0.0
0 days 14:20:00  0.0  0.0  0.0  1.0  0.0
0 days 15:10:00  0.0  0.0  1.0  1.0  2.0
0 days 16:20:00  0.0  1.0  1.0  1.0  2.0
0 days 17:10:00  0.0  1.0  1.0  0.0  0.0
2             0 days 07:30:00  0.0  0.0  0.0  2.0  0.0
0 days 08:20:00  0.0  0.0  0.0  2.0  0.0
0 days 09:10:00  0.0  0.0  0.0  1.0  0.0
0 days 10:10:00  0.0  0.0  0.0  0.0  1.0
0 days 11:00:00  0.0  0.0  0.0  0.0  1.0
0 days 13:30:00  1.0  0.0  1.0  0.0  0.0
0 days 14:20:00  2.0  0.0  1.0  0.0  0.0
0 days 15:10:00  2.0  1.0  0.0  0.0  2.0
0 days 16:20:00  2.0  2.0  0.0  0.0  2.0
0 days 17:10:00  0.0  1.0  0.0  0.0  0.0

weekday                              1    2    3    4    5
operational_group slot
1                 0 days 07:30:00  1.0  0.0  0.0  4.0  0.0
0 days 08:20:00  2.0  0.0  2.0  4.0  0.0
0 days 09:10:00  2.0  2.0  2.0  0.0  2.0
0 days 10:10:00  0.0  2.0  0.0  0.0  4.0
0 days 11:00:00  0.0  0.0  0.0  0.0  2.0
0 days 13:30:00  2.0  1.0  0.0  2.0  0.0
0 days 14:20:00  2.0  2.0  2.0  3.0  0.0
0 days 15:10:00  2.0  3.0  2.0  4.0  4.0
0 days 16:20:00  2.0  3.0  0.0  4.0  4.0
0 days 17:10:00  0.0  0.0  0.0  1.0  0.0
2                 0 days 07:30:00  0.0  0.0  2.0  2.0  2.0
0 days 08:20:00  0.0  0.0  2.0  2.0  2.0
0 days 09:10:00  2.0  0.0  0.0  2.0  0.0
0 days 10:10:00  2.0  0.0  0.0  0.0  0.0
0 days 11:00:00  0.0  0.0  0.0  0.0  0.0
0 days 13:30:00  1.0  0.0  1.0  2.0  0.0
0 days 14:20:00  2.0  0.0  1.0  2.0  0.0
0 days 15:10:00  2.0  0.0  2.0  1.0  2.0
0 days 16:20:00  3.0  2.0  2.0  0.0  2.0
0 days 17:10:00  2.0  2.0  1.0  0.0  0.0