Schedule Optimization Algorithm Walkthrough

Schedule Optimization Algorithm Walkthrough provides a step-by-step guide to building a contact center schedule optimizer from scratch. This article walks through the set-covering formulation, constraint encoding, solver execution, result interpretation, and post-optimization processing — with working Python code using PuLP.

Problem Formulation

Schedule optimization LP: shift-interval coverage matrix

The core scheduling problem: given demand across time intervals, select how many of each shift type to assign so that demand is covered at minimum cost.

This is a set-covering integer program:

$\min \sum_{j \in J} c_{j} x_{j}$

subject to:

$\sum_{j \in J} a_{i j} x_{j} \geq d_{i} \forall i \in I$

$x_{j} \in ℤ_{\geq 0} \forall j \in J$

where:

$J$ = set of shift types (the shift catalog)
$I$ = set of intervals in the planning day
$c_{j}$ = cost of one agent working shift $j$
$a_{i j}$ = 1 if shift $j$ covers interval $i$ (agent is productive), 0 otherwise
$d_{i}$ = required staffing in interval $i$
$x_{j}$ = decision variable: how many agents assigned to shift $j$

Step 1: Define the Shift Catalog

The shift catalog lists every possible shift configuration. Each shift has a start time, end time, break placement, and cost.

Example shift catalog for a center operating 7:00-23:00:

Shift ID	Start	End	Duration	Break	Cost/Shift
S1	07:00	15:00	8h	11:00-11:30	$160
S2	07:00	15:00	8h	12:00-12:30	$160
S3	08:00	16:00	8h	12:00-12:30	$160
S4	08:00	16:00	8h	13:00-13:30	$160
S5	09:00	17:00	8h	12:30-13:00	$160
S6	09:00	17:00	8h	13:30-14:00	$160
S7	10:00	18:00	8h	13:00-13:30	$160
S8	11:00	19:00	8h	14:00-14:30	$160
S9	12:00	20:00	8h	15:00-15:30	$160
S10	14:00	22:00	8h	17:00-17:30	$168

Design considerations:

More shift types = better coverage flexibility but larger problem size
Differentiate break times within same start/end — this is how break optimization integrates with shift optimization
Cost differentials: evening/weekend premiums, overtime rates
Part-time shifts (4h, 6h) should be included if workforce mix allows

Step 2: Build the Coverage Matrix

The coverage matrix $A$ is binary: $a_{i j} = 1$ if shift $j$ has an agent productive (not on break) during interval $i$ .

import numpy as np

def build_coverage_matrix(shifts, intervals):
    """
    Build binary coverage matrix.

    shifts: list of dicts with 'start', 'end', 'break_start', 'break_end'
            (all as interval indices, e.g., 0=07:00, 1=07:30, ...)
    intervals: list of interval indices
    """
    n_intervals = len(intervals)
    n_shifts = len(shifts)
    A = np.zeros((n_intervals, n_shifts), dtype=int)

    for j, shift in enumerate(shifts):
        for i in range(n_intervals):
            # Agent is productive if interval is within shift AND not on break
            in_shift = shift['start'] <= i < shift['end']
            on_break = shift['break_start'] <= i < shift['break_end']
            A[i, j] = 1 if (in_shift and not on_break) else 0

    return A

# Example: 32 half-hour intervals (07:00-23:00)
intervals = list(range(32))

shifts = [
    {'id': 'S1', 'start': 0,  'end': 16, 'break_start': 8,  'break_end': 9,  'cost': 160},
    {'id': 'S2', 'start': 0,  'end': 16, 'break_start': 10, 'break_end': 11, 'cost': 160},
    {'id': 'S3', 'start': 2,  'end': 18, 'break_start': 10, 'break_end': 11, 'cost': 160},
    {'id': 'S4', 'start': 2,  'end': 18, 'break_start': 12, 'break_end': 13, 'cost': 160},
    {'id': 'S5', 'start': 4,  'end': 20, 'break_start': 11, 'break_end': 12, 'cost': 160},
    {'id': 'S6', 'start': 4,  'end': 20, 'break_start': 13, 'break_end': 14, 'cost': 160},
    {'id': 'S7', 'start': 6,  'end': 22, 'break_start': 12, 'break_end': 13, 'cost': 160},
    {'id': 'S8', 'start': 8,  'end': 24, 'break_start': 14, 'break_end': 15, 'cost': 160},
    {'id': 'S9', 'start': 10, 'end': 26, 'break_start': 16, 'break_end': 17, 'cost': 160},
    {'id': 'S10','start': 14, 'end': 30, 'break_start': 20, 'break_end': 21, 'cost': 168},
]

A = build_coverage_matrix(shifts, intervals)

Step 3: Define Demand

Demand comes from the forecasting pipeline: Erlang C output converted to required staff per interval.

# Required staff per half-hour interval (07:00-23:00)
demand = [
    12, 15, 22, 28, 35, 38, 40, 42,  # 07:00-11:00 (morning ramp)
    40, 38, 36, 35, 33, 30, 28, 25,  # 11:00-15:00 (midday)
    22, 20, 18, 18, 20, 22, 20, 18,  # 15:00-19:00 (afternoon)
    15, 12, 10, 8,  6,  5,  4,  3,   # 19:00-23:00 (evening wind-down)
]

Step 4: Build and Solve the Integer Program

from pulp import *

# Create problem
prob = LpProblem("Shift_Scheduling", LpMinimize)

# Decision variables: how many agents per shift type
x = {}
for j, shift in enumerate(shifts):
    x[j] = LpVariable(f"x_{shift['id']}", lowBound=0, cat='Integer')

# Objective: minimize total cost
prob += lpSum(shifts[j]['cost'] * x[j] for j in range(len(shifts)))

# Constraints: coverage >= demand for each interval
for i in range(len(intervals)):
    prob += (
        lpSum(A[i, j] * x[j] for j in range(len(shifts))) >= demand[i],
        f"Coverage_interval_{i}"
    )

# Optional: maximum agents per shift type (operational limit)
for j in range(len(shifts)):
    prob += x[j] <= 20, f"Max_shift_{shifts[j]['id']}"

# Solve
solver = PULP_CBC_CMD(msg=1, timeLimit=60)  # CBC is open-source, bundled with PuLP
status = prob.solve(solver)

print(f"Status: {LpStatus[prob.status]}")
print(f"Total Cost: ${value(prob.objective):,.0f}")
print()
for j in range(len(shifts)):
    if value(x[j]) > 0:
        print(f"  {shifts[j]['id']}: {int(value(x[j]))} agents  (${shifts[j]['cost'] * int(value(x[j])):,.0f})")

Expected output:

Status: Optimal
Total Cost: $7,048

  S1: 5 agents  ($800)
  S3: 8 agents  ($1,280)
  S5: 12 agents ($1,920)
  S7: 8 agents  ($1,280)
  S8: 5 agents  ($800)
  S9: 3 agents  ($480)
  S10: 3 agents ($504)

Step 5: Interpret Results

Reading the Solution

The solver output tells you the optimal number of agents per shift type. This is not yet a roster — you don't know which agents work which shifts. That comes in Step 7.

Understanding the Gap

# Check coverage vs demand
total_staff = np.zeros(len(intervals))
for j in range(len(shifts)):
    total_staff += A[:, j] * int(value(x[j]))

print("\nInterval | Required | Scheduled | Delta")
print("-" * 45)
for i in range(len(intervals)):
    time = f"{7 + i//2:02d}:{(i%2)*30:02d}"
    delta = int(total_staff[i] - demand[i])
    flag = " ***" if delta < 0 else ""
    print(f"  {time}   |    {demand[i]:3d}   |    {int(total_staff[i]):3d}    |  {delta:+d}{flag}")

# Overstaffing cost
overstaffing_intervals = sum(max(0, total_staff[i] - demand[i]) for i in range(len(intervals)))
print(f"\nTotal overstaffing: {overstaffing_intervals:.0f} agent-intervals")
print(f"Overstaffing %: {overstaffing_intervals / sum(demand) * 100:.1f}%")

MIP Gap

The MIP gap measures how close the solution is to the theoretical optimum:

$Gap = \frac{Best Integer - Best Bound}{Best Integer} \times 100 %$

Gap = 0%: Provably optimal solution found
Gap < 1%: Effectively optimal for WFM purposes (forecast error is much larger)
Gap > 5%: Solver needs more time, or problem formulation should be reviewed

Step 6: Add Real-World Constraints

The basic formulation above ignores many operational realities. Here is how to encode common constraints:

Maximum Headcount

# Total agents cannot exceed available headcount
prob += lpSum(x[j] for j in range(len(shifts))) <= 55, "Max_headcount"

Minimum Staffing Floor

# Never go below 3 agents in any interval (even if demand says 1)
for i in range(len(intervals)):
    prob += lpSum(A[i,j] * x[j] for j in range(len(shifts))) >= max(demand[i], 3)

Shift Mix Constraints

# At least 30% of shifts must be full-time (8h)
full_time_shifts = [j for j, s in enumerate(shifts) if s['end'] - s['start'] >= 16]
prob += lpSum(x[j] for j in full_time_shifts) >= 0.30 * lpSum(x[j] for j in range(len(shifts)))

Overtime Limit

# Overtime shifts (if defined) cannot exceed 10% of total
overtime_shifts = [j for j, s in enumerate(shifts) if s.get('is_overtime', False)]
prob += lpSum(x[j] for j in overtime_shifts) <= 0.10 * lpSum(x[j] for j in range(len(shifts)))

Skill Coverage

For multi-skill environments, add per-skill coverage constraints:

# For each skill group, demand must be met by agents with that skill
for skill in skill_groups:
    skill_shifts = [j for j, s in enumerate(shifts) if skill in s['skills']]
    for i in range(len(intervals)):
        prob += (
            lpSum(A[i,j] * x[j] for j in skill_shifts) >= demand_by_skill[skill][i],
            f"Skill_{skill}_interval_{i}"
        )

Step 7: Post-Optimization Processing

Rostering (Agent-to-Shift Assignment)

The optimization output says "assign 12 agents to shift S5." Rostering decides which 12 agents. This is a separate optimization problem:

Inputs:

Shift counts from optimizer (Step 4)
Agent attributes: skills, preferences, constraints (max hours/week, day-off patterns)
Fairness rules: equitable weekend distribution, rotation equity

Approach: Formulate as a binary assignment problem:

# y[a,j] = 1 if agent a is assigned to shift j
agents = list(range(55))

y = {}
for a in agents:
    for j in range(len(shifts)):
        y[a, j] = LpVariable(f"assign_{a}_{j}", cat='Binary')

# Each agent gets exactly one shift (or zero if not working today)
for a in agents:
    roster_prob += lpSum(y[a,j] for j in range(len(shifts))) <= 1

# Shift counts must match optimizer output
for j in range(len(shifts)):
    roster_prob += lpSum(y[a,j] for a in agents) == int(value(x[j]))

# Agent-specific constraints
for a in agents:
    # Agent can only work shifts they're skilled for
    for j in range(len(shifts)):
        if not agent_can_work_shift(a, j):
            roster_prob += y[a,j] == 0

Break Insertion

If the shift catalog does not pre-define break times (many production systems don't), break placement is a separate optimization:

Objective: Place breaks to minimize the maximum coverage gap across all intervals
Constraints: Break must fall within legal window (e.g., between hour 3 and hour 6 of shift), each agent gets exactly one break, minimum/maximum break duration

Infeasibility Handling

When the solver returns "Infeasible," the constraints cannot all be satisfied simultaneously. Common causes and fixes:

Cause	Diagnosis	Fix
Insufficient headcount	Max_headcount < peak demand	Relax headcount or accept SL degradation
Conflicting constraints	Shift mix + coverage + headcount incompatible	Rank constraints by priority, relax lowest priority
Too-narrow shift catalog	No shift type covers a critical interval	Add shift types to catalog
Skill mismatch	Not enough skilled agents for a queue	Cross-train or relax skill requirement

Constraint relaxation approach: Add slack variables to soft constraints with penalty costs:

# Soft coverage constraints with penalty
slack = {}
for i in range(len(intervals)):
    slack[i] = LpVariable(f"slack_{i}", lowBound=0)
    prob += lpSum(A[i,j] * x[j] for j in range(len(shifts))) + slack[i] >= demand[i]

# Add penalty to objective
penalty_cost = 500  # $/agent-interval of understaffing
prob += lpSum(shifts[j]['cost'] * x[j] for j in range(len(shifts))) + \
       penalty_cost * lpSum(slack[i] for i in range(len(intervals)))

Solver Selection

Solver	License	Speed	Best For
CBC (COIN-OR)	Open-source, bundled with PuLP	Adequate for < 500 agents	Small-to-medium problems, prototyping
HiGHS	Open-source	Fast, competitive with commercial	Medium problems, budget-constrained teams
GLPK	Open-source	Slower than CBC for MIP	Educational use
Gurobi	Commercial ($$$, free academic)	Fastest for large MIP	Large centers (1000+ agents), tight SLA on solve time
CPLEX (IBM)	Commercial ($$$)	Comparable to Gurobi	Enterprise environments already using IBM
OR-Tools (Google)	Open-source	Good for CP-SAT formulations	Constraint programming approach, rostering

Recommendation: Start with CBC (free, bundled). If solve times exceed 5 minutes for your problem size, switch to HiGHS (free) or Gurobi (commercial, free for academic).

Full Working Example

Complete end-to-end script combining all steps:

"""
Contact center shift optimizer — complete working example.
Requires: pip install pulp numpy
"""
from pulp import *
import numpy as np

# --- Configuration ---
INTERVALS_PER_HOUR = 2  # 30-minute intervals
OPERATING_START = 7     # 07:00
OPERATING_END = 23      # 23:00
N_INTERVALS = (OPERATING_END - OPERATING_START) * INTERVALS_PER_HOUR

def interval_to_time(i):
    hour = OPERATING_START + i // INTERVALS_PER_HOUR
    minute = (i % INTERVALS_PER_HOUR) * (60 // INTERVALS_PER_HOUR)
    return f"{hour:02d}:{minute:02d}"

# --- Demand (from Erlang C output) ---
demand = [
    12, 15, 22, 28, 35, 38, 40, 42,
    40, 38, 36, 35, 33, 30, 28, 25,
    22, 20, 18, 18, 20, 22, 20, 18,
    15, 12, 10, 8,  6,  5,  4,  3,
]

# --- Shift Catalog ---
shifts = [
    {'id': 'Early-8h-A',  'start': 0,  'end': 16, 'break_start': 8,  'break_end': 9,  'cost': 160},
    {'id': 'Early-8h-B',  'start': 0,  'end': 16, 'break_start': 10, 'break_end': 11, 'cost': 160},
    {'id': 'Mid-8h-A',    'start': 2,  'end': 18, 'break_start': 10, 'break_end': 11, 'cost': 160},
    {'id': 'Mid-8h-B',    'start': 2,  'end': 18, 'break_start': 12, 'break_end': 13, 'cost': 160},
    {'id': 'Core-8h-A',   'start': 4,  'end': 20, 'break_start': 11, 'break_end': 12, 'cost': 160},
    {'id': 'Core-8h-B',   'start': 4,  'end': 20, 'break_start': 13, 'break_end': 14, 'cost': 160},
    {'id': 'Late-8h-A',   'start': 6,  'end': 22, 'break_start': 12, 'break_end': 13, 'cost': 160},
    {'id': 'Swing-8h',    'start': 8,  'end': 24, 'break_start': 14, 'break_end': 15, 'cost': 160},
    {'id': 'Close-8h',    'start': 10, 'end': 26, 'break_start': 16, 'break_end': 17, 'cost': 160},
    {'id': 'Evening-8h',  'start': 14, 'end': 30, 'break_start': 20, 'break_end': 21, 'cost': 168},
    {'id': 'AM-Part-4h',  'start': 0,  'end': 8,  'break_start': -1, 'break_end': -1, 'cost': 80},
    {'id': 'PM-Part-4h',  'start': 16, 'end': 24, 'break_start': -1, 'break_end': -1, 'cost': 84},
]

# --- Build Coverage Matrix ---
A = np.zeros((N_INTERVALS, len(shifts)), dtype=int)
for j, shift in enumerate(shifts):
    for i in range(N_INTERVALS):
        in_shift = shift['start'] <= i < min(shift['end'], N_INTERVALS)
        on_break = shift['break_start'] <= i < shift['break_end']
        A[i, j] = 1 if (in_shift and not on_break) else 0

# --- Formulate and Solve ---
prob = LpProblem("WFM_Shift_Optimizer", LpMinimize)
x = {j: LpVariable(f"x_{shifts[j]['id']}", lowBound=0, cat='Integer') for j in range(len(shifts))}

# Objective: minimize cost
prob += lpSum(shifts[j]['cost'] * x[j] for j in range(len(shifts)))

# Coverage constraints
for i in range(N_INTERVALS):
    prob += lpSum(A[i,j] * x[j] for j in range(len(shifts))) >= demand[i], f"Cov_{i}"

# Part-time limit: <= 20% of total
part_time = [j for j, s in enumerate(shifts) if 'Part' in s['id']]
full_time = [j for j, s in enumerate(shifts) if 'Part' not in s['id']]
prob += lpSum(x[j] for j in part_time) <= 0.20 * lpSum(x[j] for j in range(len(shifts)))

# Solve
status = prob.solve(PULP_CBC_CMD(msg=0, timeLimit=120))

# --- Output ---
print(f"Status: {LpStatus[status]}")
print(f"Optimal Cost: ${value(prob.objective):,.0f}")
print(f"\nShift Assignments:")
total_agents = 0
for j in range(len(shifts)):
    count = int(value(x[j]))
    if count > 0:
        print(f"  {shifts[j]['id']:15s}: {count:3d} agents  ${shifts[j]['cost'] * count:>6,}")
        total_agents += count
print(f"  {'TOTAL':15s}: {total_agents:3d} agents")

# Coverage analysis
staff = np.zeros(N_INTERVALS)
for j in range(len(shifts)):
    staff += A[:, j] * int(value(x[j]))

print(f"\nCoverage Analysis:")
print(f"{'Time':>7s} | {'Req':>4s} | {'Sched':>5s} | {'Delta':>5s}")
print("-" * 30)
for i in range(N_INTERVALS):
    delta = int(staff[i] - demand[i])
    print(f"{interval_to_time(i):>7s} | {demand[i]:4d} | {int(staff[i]):5d} | {delta:+5d}")

total_req = sum(demand)
total_sched = int(sum(staff))
overstaff = sum(max(0, staff[i] - demand[i]) for i in range(N_INTERVALS))
print(f"\nTotal required agent-intervals:  {total_req}")
print(f"Total scheduled agent-intervals: {total_sched}")
print(f"Overstaffing: {overstaff:.0f} agent-intervals ({overstaff/total_req*100:.1f}%)")

Extending the Basic Model

Multi-Day Scheduling

The example above optimizes a single day. Production scheduling optimizes an entire week (or pay period) simultaneously, adding constraints that span days:

Weekly hour limits: Each agent works 36-40 hours per week
Consecutive day constraints: Maximum 5 consecutive working days
Day-off patterns: Two consecutive days off per week preferred
Rotation fairness: Weekend shifts distributed equitably across agents

Multi-day formulation multiplies problem size by the number of days. A 7-day problem with 15 shift types and 48 intervals has 7 × 15 = 105 shift-day combinations and 7 × 48 = 336 coverage constraints — still tractable for modern solvers.

Cost Function Refinement

The basic model uses constant cost per shift. Realistic cost functions include:

Cost Component	How to Encode
Base wage	Shift hours × hourly rate
Overtime premium	Hours above 8/day or 40/week at 1.5x rate
Shift differential	Evening/night/weekend premiums added to base
Overstaffing penalty	Cost per excess agent-interval (soft constraint in objective)
Understaffing penalty	Cost per missing agent-interval (typically 3-5x overstaffing cost)
Preference violation	Small penalty for not granting agent shift preferences

By adjusting the relative size of understaffing vs overstaffing penalties, you control the optimizer's tendency to overstaff (conservative) or understaff (aggressive).

Iterative Refinement Workflow

In practice, schedule optimization is not a single solver run. It is an iterative process:

Run initial optimization with base constraints
Review coverage chart — identify intervals with excessive overstaffing or understaffing
Adjust shift catalog — if a consistent gap exists at 7 PM, add a shift type that covers 7 PM
Tighten/relax constraints — if the optimizer cannot find a feasible solution, identify which constraint is binding and discuss with operations whether it can be relaxed
Re-run and compare — measure total cost, overstaffing %, understaffing %, and schedule efficiency against the previous iteration
Present to operations — show the coverage chart, total cost, and tradeoffs. Get approval.
Proceed to rostering — assign specific agents to the approved shift plan

Typical iteration count: 3-5 iterations before a schedule is ready for publication. The first run is rarely the final answer.

When to Build vs When to Use Vendor Optimizer

Factor	Build Custom (PuLP/OR-Tools)	Use Vendor WFM Optimizer
Constraint flexibility	Unlimited — encode any constraint you can formulate	Limited to vendor's supported constraint types
Transparency	Full visibility into formulation and solver behavior	Black box — you see inputs and outputs, not the model
Integration	Requires custom integration with WFM tool	Native integration with vendor ecosystem
Maintenance	You own it — updates, bugs, and infrastructure	Vendor maintains and updates
Skill requirement	Needs OR/optimization expertise on team	WFM analyst can operate with training
Speed to production	Weeks to months for custom build	Days to configure

Recommendation: Use vendor optimizer for daily operations. Build custom optimizer for (a) research and experimentation, (b) validating vendor results, (c) constraints the vendor cannot handle, or (d) environments where no WFM tool is deployed.

Anonymous

Search