Pulp staffing optimization problem: shift combination constraint

Question

I'm trying to solve a staffing optimization problem and this is the model I have.

First constraint: a minimum number of employee needs to be scheduled for each shift.

For exemple:

shift	requirement
1	1
2	1
3	2
4	2
5	1

Second constraint: An employee can work 0 shift, 1 shift, or 2 shift and it can't be whatever shift.

Let me explain, if the employee works 2 shifts, they need exactly 1 shift break between each.

This is the list of the possible combination:

shift	shift
1	3
2	4
3	5
1
2
3
4
5

Desired outcome:

Let's assume these are my employees:

1. Robin
2. Mathew
3. George
4. Elisa

This is what a valid solution could look like

shift	employee
1	Robin
2	George
3	Robin, Mathew
4	George, Elisa
5	Mathew

Robin 1, 3

George 2, 4

Mathew 3, 5

Elisa 4

I have no problem with the first constraint. It's more the second constraint I'm struggling to implement in Pulp. I would appreciate it if you could give me pointers.

Thanks!

Edit: This is the code I have right now, trying to go the way @airsquid suggested

# constraint 1: must schedule someone at shift 1, 2, 5
for shift, patterns in vars_by_shift.items():
    if shift == 1 or shift == 2 or shift == 5:
        prob += sum(vars_by_shift[shift]) >= 1
    elif shift == 4 or shift == 5:  # and 2 persons at shift 4 and 5
        prob += sum(vars_by_shift[shift]) >= 2

# constraint 2: agent can at most be scheduled 1 time among all patterns
for agent in agent_names:
    prob += sum(vars_by_agent[agent]) <= 1

Where vars_by_shift looks like this

{
    1: [1_3,1,robin, 1,1,robin, 1_3,1,mathew, 1,1,mathew, 1_3,1,george, 1,1,george, 1_3,1,elisa, 1,1,elisa], 
    2: [2_4,2,robin, 2,2,robin, 2_4,2,mathew, 2,2,mathew, 2_4,2,george, 2,2,george, 2_4,2,elisa, 2,2,elisa], 
    3: [3_5,3,robin, 3,3,robin, 3_5,3,mathew, 3,3,mathew, 3_5,3,george, 3,3,george, 3_5,3,elisa, 3,3,elisa], 
    4: [2_4,4,robin, 4,4,robin, 2_4,4,mathew, 4,4,mathew, 2_4,4,george, 4,4,george, 2_4,4,elisa, 4,4,elisa], 
    5: [3_5,5,robin, 5,5,robin, 3_5,5,mathew, 5,5,mathew, 3_5,5,george, 5,5,george, 3_5,5,elisa, 5,5,elisa]
}

and vars_by_agent like this:

{
    'robin': [1_3,1,robin, 1,1,robin, 2_4,2,robin, 2,2,robin, 3_5,3,robin, 3,3,robin, 2_4,4,robin, 4,4,robin, 3_5,5,robin, 5,5,robin], 
    'mathew': [1_3,1,mathew, 1,1,mathew, 2_4,2,mathew, 2,2,mathew, 3_5,3,mathew, 3,3,mathew, 2_4,4,mathew, 4,4,mathew, 3_5,5,mathew, 5,5,mathew], 
    'george': [1_3,1,george, 1,1,george, 2_4,2,george, 2,2,george, 3_5,3,george, 3,3,george, 2_4,4,george, 4,4,george, 3_5,5,george, 5,5,george], 
    'elisa': [1_3,1,elisa, 1,1,elisa, 2_4,2,elisa, 2,2,elisa, 3_5,3,elisa, 3,3,elisa, 2_4,4,elisa, 4,4,elisa, 3_5,5,elisa, 5,5,elisa]
}

I get an infeasible solution when calling solve() and I'm not sure why.

Final edit:

Thank you so much @AirSquid

Finally got a valid solution.

The key was in making sure I don't duplicate decision variables when storing them in the different dictionaries.

def shift_to_be_sent_to(pattern):
    shifts = []
    for shift in range(1, 6):
        if coverage_map[pattern, shift] == 1:
            shifts.append(shift)
    return shifts

Final code is available here

Answer 1

Discrete-optimization is all about careful formalization of the task which makes it possible to exploit various assumptions formalized. There is always a trade-off involved balancing generalization and specialization.

Your example would allow the following specialization, which might not be possible if we omitted things later introduced. Nonetheless let's work with problem as presented.

As you already modelled the (generalized) assignment-problem parts which implies the cardinalities (0, 1, or 2 shifts) we only need to model the "forbidden assignments" in the case of 2 shifts for some worker.

This is simple: Just use "implications".

FOR EACH WORKER W

    FOR EACH SHIFT S in [1, N_SHIFTS)

        assignment(w, s) -> NOT assignment(w, s+1)

We already assumed a binary/logic-based domain and the following holds:

  the logical implication a -> b

is equivalent to:

  NOT a OR b

This can be linearized as following:

(1-a) + b >= 1

So all we need is:

FOR EACH WORKER W

    FOR EACH SHIFT S in [1, N_SHIFTS)

        (1 - assignment(w, s) ) + (1 - assignment(w, s+1)) >= 1

This might look like a lot of constraints, but every competitive MILP-solver (PULPs solvers too; at least Cbc, less sure about Glpk) has explicit exploitation of these constraints = clique-tables and co. The solver will use powerful reasoning here!

EDIT

Reading the task again, i omitted "they need exactly 1 shift break between each".

Well... we could start from scratch or just use the approach again:

FOR EACH WORKER W

    FOR EACH SHIFT S in [1, N_SHIFTS)

        assignment(w, s) -> assignment(w, s+2) 
        IFF sum(assignment(w, *)) == 2

This works as above (note the non-negative right-hand side), but we somehow must deactivate this implication when there are 0 or 1 assignments!

This we can achieve by a big-M like formulation:

We adapt the right-hand side from:

>= 1

basically meaning: ALWAYS

to:

>= sum() - 1

basically meaning: Only when sum() >= 2

Leading to:

FOR EACH WORKER W

    FOR EACH SHIFT S in [1, N_SHIFTS-1)

        (1 - assignment(w, s) ) + assignment(w, s+2) >= sum(w, *) - 1

which is:

*activated* if sum == 2
*deactivated* if sum < 2
breaks when we got 3 or more assignments -> infeasible

Opposed to above, this one is much more aggressive and might not work as well with future additions! I_ think it's clear, that this exploits a lot!

In more complex pattern-like models i'm a huge fan of automata and would always recommend those! This howewer has some tough initial development-costs as you basically need some graph-library at modelling-time and be able to work with both (the other being your MILP-model) in sync.

Motivation and background (although focused on constraint-programming) is explained here:

"Flexible Optimization: Nurse Scheduling with Constraint Programming and Automata").

Some MILP-focused related work would be available in:

Côté, Marie-Claude, Bernard Gendron, and Louis-Martin Rousseau. "Modeling the regular constraint with integer programming." International Conference on Integration of Artificial Intelligence (AI) and Operations Research (OR) Techniques in Constraint Programming. Springer, Berlin, Heidelberg, 2007.

(Although it's easy to work out the theory when having seen the idea of layered-graphs -> expansion of the automaton over the sequence-length and a network-flow like model / flow-constraints).

I might add, that or-tools cp-sat is an amazing hybrid-solver (CP/SAT/MILP/more) which has support for automata-constraints out of the box.

Answer 2

I really like the other solution presented, which is great for larger generalizations of this type of problem. For the complexity of the problem presented, I think you could get by with a simpler approach that uses enumeration of the possible shifts in a simple way.

In your problem description, each employee may be assigned only 1 of 8 patterns of work (or no assignment at all). That gives rise to a 2-indexed binary variable:

assign[employee, pattern]

The relationship of the "pattern" to the requirement in each shift is known, so you then have a model parameter that pairs these 2 things:

fulfills[pattern, shift] = {('p1', 's1'): 1, ... }

with a little elbow grease, and some summation constraints, you have a model with that variable, parameter, 3 sets (employees, patterns, shifts) and 3 or 4 constraints

This works fine in the case described because it is possible to enumerate the acceptable patterns fairly easily and the relationship of pattern:outcome is a 1-to-1 pairing.

Edit: augmented regarding your posted code...

You are on the right track... kinda... lol. I think your approach can be successful, I had something a little different in mind that is not triple indexed, but yours can work.

2 probs you have: you have typo in your constraint for shifts. You have shift '5' instead of '3' for 2 requirement. Also, if you look at your vars_by_shift you'll note that for shift 3 you are not crediting the 1-3 combo!! You might want to restructure that part a bit. Here is an idea on how you can get the shifts associated with patterns more easily. (Realize there are a bunch of ways to do this):

Example code:

cov_file = 'coverages.csv'

coverages = dict()
with open(cov_file, 'r') as src:
    for line in src:
        data = line.strip().split(',')
        pattern = data[0]
        shifts = [int(t) for t in data[1:]]
        coverages[pattern] = shifts

shifts = [1, 2, 3, 4]
patterns = coverages.keys()

coverage_map = dict()
for (p, s) in [(p, s) for p in patterns for s in shifts]:
    coverage_map[p, s] = 1 if s in coverages[p] else 0

print(coverage_map)

coverages.csv

A,1
B,1,3
C,2,4
D,2,3,4

Yields:

{('A', 1): 1, ('A', 2): 0, ('A', 3): 0, ('A', 4): 0, ('B', 1): 1, ('B', 2): 0, ('B', 3): 1, ('B', 4): 0, ('C', 1): 0, ('C', 2): 1, ('C', 3): 0, ('C', 4): 1, ('D', 1): 0, ('D', 2): 1, ('D', 3): 1, ('D', 4): 1}