Column generation with employee shift scheduling

Hi, I am trying to apply column generation on an employee shift scheduling problem. Please find below my code. What I see with my code is that solving for the sub-problem always creates the same shift pattern (column), and after a few iterations the total cost remains same.

Can you please let me know what am I doing incorrect with the below program.

import gurobipy as gp
from gurobipy import GRB
import random
import copy

# data
# every employee is available for 24 hours
# e.g. SMITH has to work for a minimum of 6 hours and maximum of 8 hours
# payrate for SMITH is 30 dollars per hour
employees, min_hours, max_hours, cost, availability = gp.multidict(
    {
        "SMITH": [6, 8, 30, range(24)],
        "JOHNSON": [6, 8, 50, range(24)],
        "WILLIAMS": [6, 8, 30, range(24)],
        "JONES": [6, 8, 30, range(24)],
        "BROWN": [6, 8, 40, range(24)],
        "DAVIS": [6, 8, 50, range(24)],
        "MILLER": [6, 8, 45, range(24)],
        "WILSON": [6, 8, 30, range(24)],
        "MOORE": [6, 8, 35, range(24)],
        "TAYLOR": [6, 8, 40, range(24)],
        "ANDERSON": [2, 3, 60, range(24)],
        "THOMAS": [2, 4, 40, range(24)],
        "JACKSON": [2, 4, 60, range(24)],
        "WHITE": [2, 6, 55, range(24)],
        "HARRIS": [2, 6, 45, range(24)],
        "MARTIN": [2, 3, 40, range(24)],
        "THOMPSON": [2, 5, 50, range(24)],
        "GARCIA": [2, 4, 50, range(24)],
        "MARTINEZ": [2, 4, 40, range(24)],
        "ROBINSON": [2, 5, 50, range(24)],
    }
)
employees_copy = copy.deepcopy(employees)

# number of team members required for each hour of the day
requirements = [1, 1, 2, 3, 6, 6, 7, 8, 9, 8, 8, 8, 7, 6, 6, 5, 5, 4, 4, 3, 2, 2, 2, 2]
requirements = {i: requirements[i] for i in range(24)}




# hours in the shift where each team member is working for
# to create an initial feasible solution (dumb solution) we say none of the team member 
# is working
# but we create slack variable to meet the demand. 

# so we have "('SMITH', 0)" zero says, we have only one column for "SMITH"
# but in future if we get more columns for "SMITH" we will add them under
# "('SMITH', 1)"

empl_shifts = gp.tupledict({('SMITH', 0) : [0] * 24,
                ('JOHNSON', 0) : [0] * 24,
                ('WILLIAMS', 0) : [0] * 24,
                ('JONES', 0) : [0] * 24,
                ('BROWN', 0) : [0] * 24,
                ('DAVIS', 0) : [0] * 24,
                ('MILLER', 0) : [0] * 24,
                ('WILSON', 0) : [0] * 24,
                ('MOORE', 0) : [0] * 24,
                ('TAYLOR', 0) : [0] * 24,
                ('ANDERSON', 0) : [0] * 24,
                ('THOMAS', 0) : [0] * 24,
                ('JACKSON', 0) : [0] * 24,
                ('WHITE', 0) : [0] * 24,
                ('HARRIS', 0) : [0] * 24,
                ('MARTIN', 0) : [0] * 24,
                ('THOMPSON', 0) : [0] * 24,
                ('GARCIA', 0) : [0] * 24,
                ('MARTINEZ', 0) : [0] * 24,
                ('ROBINSON', 0) : [0] * 24})

# solve the master problem

model = gp.Model("employee_scheduling")

# decision variable for each of the shift we enumerate for each team member
x = model.addVars(empl_shifts, vtype=GRB.CONTINUOUS, lb=0, ub=1)

# since initially none of the team member is working, to meet the demand
# we introduce slack variable and the cost attached to these slack variables
# cost attached is higher than the team member pay rate.
# so usage of slack variables is discouraged
c_t = {i: 100 for i in range(24)}
s_t = model.addVars(range(24), vtype=GRB.CONTINUOUS)

# in "lnexp" we hold the total cost of putting team members in the
# right shift to meet the demand for each hour
lnexp = gp.LinExpr()

for e in employees:
    num = [k[1] for k in empl_shifts if k[0] == e]
    lnexp.add(
        gp.quicksum(sum(empl_shifts[(e, n)]) * x[(e, n)] * cost[e] for n in num)
    )

# obj is the team member total cost + cost of meeting demand from slack
obj = lnexp + gp.quicksum(c_t[i] * s_t[i] for i in range(24))

# constraint : select at max only one column (shift schedule) for a team member
model.addConstrs(x.sum(e, "*") <= 1 for e in employees)

# constraint : meet the hourly demand either through team members or through slack
constr = model.addConstrs(
    gp.quicksum(empl_shifts[k][i] * x[k] for k in empl_shifts) + s_t[i]
    >= requirements[i]
    for i in range(24)
)

model.setObjective(obj, sense=GRB.MINIMIZE)

model.optimize()


TOL = 1e-6
iter_count = 0
while iter_count <= 10:
    price = [constr[i].pi for i in range(24)]
    print(f"Price = {price}")

    # sub-problem
    sp = gp.Model("subproblem")
    sp.ModelSense = GRB.MAXIMIZE

    new_pattern = sp.addVars(range(24), obj=price, vtype=GRB.BINARY)

    # idea behind "diff" and diff_max" is to ensure consecutiveness of shift
    diff = sp.addVars(range(1, 24), vtype=GRB.BINARY)
    diff_max = sp.addVars(range(1, 24), vtype=GRB.BINARY)
    sp.addConstrs(diff[i] == new_pattern[i] - new_pattern[i - 1] for i in range(1, 24))
    for i in range(1, 23):
        sp.addGenConstrMax(diff_max[i], [0, diff[i]])

    sp.addConstr(gp.quicksum(diff_max) == 1)

    # min number of working hours = 2 and max hours = 8 from the data
    sp.addRange(gp.quicksum(new_pattern), 2, 8)
    sp.optimize()
    # reduced cost
    min_rc = 1 - sp.objVal
    if min_rc < -TOL:
        sol = [int(abs(new_pattern[i].X)) for i in range(24)]

        # now we have got a new shift pattern, or a column,
        # we need to assign this shift pattern to an employee
        # so we shuffle the employeed and pick an employee who would be
        # compatible with this new pattern (or column)
        random.shuffle(employees_copy)
        for e in employees_copy:
            if (min_hours[e] <= sum(sol)) and (max_hours[e] >= sum(sol)):
                ind = max([k[1] for k, v in empl_shifts.items() if k[0] == e]) + 1
                # update empl_shifts with the new pattern
                empl_shifts[(e, ind)] = sol
                break

    # with the new column added solve the master problem again
    # below is the exact same code as above
    model = gp.Model("employee_scheduling")
    
    x = model.addVars(empl_shifts, vtype=GRB.CONTINUOUS, lb=0, ub=1)
    
    c_t = {i: 100 for i in range(24)}
    s_t = model.addVars(range(24), vtype=GRB.CONTINUOUS)
    
    lnexp = gp.LinExpr()
    
    for e in employees:
        num = [k[1] for k in empl_shifts if k[0] == e]
        lnexp.add(
            gp.quicksum(sum(empl_shifts[(e, n)]) * x[(e, n)] * cost[e] for n in num)
        )
    
    obj = lnexp + gp.quicksum(c_t[i] * s_t[i] for i in range(24))
    
    model.addConstrs(x.sum(e, "*") <= 1 for e in employees)
    
    constr = model.addConstrs(
        gp.quicksum(empl_shifts[k][i] * x[k] for k in empl_shifts) + s_t[i]
        >= requirements[i]
        for i in range(24)
    )
    
    model.setObjective(obj, sense=GRB.MINIMIZE)
    
    model.optimize()