Seeking help in pulp scheduling problem

Question

I have successfully implemented a program where I allocate N truck drivers to M gathering hubs for each one of the days of the week. The constraints I have implemented are:

1. • A driver cannot work more than 6 days, i.e. 1 day to rest
2. • A driver cannot be allocated in more than 1 hubs for each day
3. • Each hub must satisfy its driver requirements for each day of the week
4. • A driver must work his days at one hub instead of many.

The program runs smoothly, satisfies the overall objective and outputs a schedule for each hub-driver pair.

However, the output schedule allocates 30 drivers who all work in one-hub, but unfortunately their working days are way less than 6 / week. How could that be amended? The optimal solution would be that each driver works 6 days (one day off) and only at one hub, but this, unfortunately is not possible from what I understand. Small tweaks in the constraints or demand are acceptable. Any ideas?

Code below.

import pulp
import pandas as pd
import numpy as np

day_reqs = [[1, 1, 1, 1, 1, 1, 1],
            [3, 4, 3, 4, 5, 3, 3],
            [3, 4, 4, 3, 4, 5, 5],
            [1, 1, 1, 1, 1, 1, 2],
            [2, 2, 2, 2, 2, 2, 2],
            [2, 4, 3, 2, 2, 2, 2],
            [6, 5, 3, 3, 3, 5, 4],
            [2, 3, 2, 2, 2, 1, 2]]

total_day_requirements = ([sum(x) for x in zip(*day_reqs)])

hub_names = {0: 'Hub 1',
             1: 'Hub 2',
             2: 'Hub 3',
             3: 'Hub 4',
             4: 'Hub 5',
             5: 'Hub 6',
             6: 'Hub 7',
             7: 'Hub 8'}

total_drivers = max(total_day_requirements)  # minimum number of drivers
total_hubs = len(day_reqs)  # number of hubs
days = 7


def crosshubbers(dashboard, driver_names):
    test = dashboard.reset_index()
    counter = 0
    for name in driver_names:
        driv = test[test['Driver'] == name]

        temp = list(driv.sum(axis=1).values)
        cnt = 0
        for val in temp:
            if val > 0:
                cnt += 1

        if cnt > 1:
            #     print(f'{cnt} for driver {name}')
            counter += 1
    return counter

def schedule(drivers, hubs, day_requirement):
    driver_names = ['Driver_{}'.format(i) for i in range(drivers)]
    var = pulp.LpVariable.dicts('VAR', (range(hubs), range(drivers), range(days)), 0, 1, 'Binary')
    problem = pulp.LpProblem('shift', pulp.LpMinimize)

    obj = None
    for h in range(hubs):
        for driver in range(drivers):
            for day in range(days):
                obj += var[h][driver][day]
    problem += obj

    #  we add binary variables z indicating if a driver is active on a hub:
    z = pulp.LpVariable.dicts('Z', (range(hubs), range(drivers)), 0, 1, 'Binary')
    # minimize the sum of drivers active on hubs
    for h in range(hubs):
        for driver in range(drivers):
            obj += z[h][driver]
    problem += obj

    # Add constraints to connect z with VAR:
    for driver in range(drivers):
        for h in range(hubs):
            problem += z[h][driver] <= pulp.lpSum(var[h][driver][day] for day in range(days))
            problem += days * z[h][driver] >= pulp.lpSum(var[h][driver][day] for day in range(days))

    for driver in range(drivers):
        problem += pulp.lpSum(z[h][driver] for h in range(hubs)) <= 1

    # if a driver works one day at a hub, he cannot work that day in a different hub obviously
    for driver in range(drivers):
        for day in range(days):
            problem += pulp.lpSum([var[h][driver][day] for h in range(hubs)]) <= 1

    # schedule must satisfy daily requirements of each hub
    for day in range(days):
        for h in range(hubs):
            problem += pulp.lpSum(var[h][driver][day] for driver in range(drivers)) == \
                       day_requirement[h][day]

    # a driver cannot work more than 6 days
    for driver in range(drivers):
        problem += pulp.lpSum([var[h][driver][day] for day in range(days) for h in range(hubs)]) \
                   <= 6

    # Solve problem. We have a very complex solution so we set a timeout at 10secs.
    status = problem.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=30))

    idx = pd.MultiIndex.from_product([hub_names.values(), driver_names], names=['Hub', 'Driver'])

    # col = ['Sunday', 'Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday']
    col = ['Day_{}'.format(i) for i in range(days)]
    dashboard = pd.DataFrame(0, idx, col)

    for h in range(hubs):
        for driver in range(drivers):
            for day in range(days):
                if var[h][driver][day].value() > 0.001:
                    dashboard.loc[hub_names[h], driver_names[driver]][col[day]] = 1

    # print(dashboard)

    driver_table = dashboard.groupby('Driver').sum()
    driver_sums = driver_table.sum(axis=1)
    # print(driver_sums)

    day_sums = driver_table.sum(axis=0)
    # print(day_sums)

    print("Status", pulp.LpStatus[status])
    return driver_sums, dashboard, status


driver_sums, dashboard, status = schedule(total_drivers, total_hubs, day_reqs)
counter = crosshubbers(dashboard, ['Driver_{}'.format(i) for i in range(total_drivers)])
while (driver_sums > 6).any() or status == -1 or counter > 0:
    print('Staus infeasible or cross-hubbers or one or more drivers have been allocated more than 6 '
          'days of: {}->{}'.format(total_drivers, total_drivers + 1))
    print(f'Status: {status}')
    print(f'Cross - Hubbers: {counter}')
    if counter == 0:
        print(f'Driver over limit: {driver_sums[driver_sums > 6].count()}')
        print(driver_sums[driver_sums > 6])
        print('\n')
        print(f'Driver under limit: {driver_sums[driver_sums < 6].count()}')
        print(driver_sums[driver_sums < 6])
    print('\n')
    total_drivers += 1
    driver_sums, dashboard, status = schedule(total_drivers, total_hubs, day_reqs)
    counter = crosshubbers(dashboard, ['Driver_{}'.format(i) for i in range(total_drivers)])
print('########################################################################')
print('Found solution')

```