Column generation for workforce scheduling

A client of mine showed some interest in developing some MIP based workforce planning models. Column generation has been used successfully in such applications. To show how such an approach could look like, I have added a small example here: http://www.amsterdamoptimization.com/pdf/colgen.pdf.

The master problem is simply:

I.e. make sure you have enough shifts to handle all forecasted calls. Inside the algorithm we use a relaxed version of this.

The subproblems are a little bit more complicated: they form allowed shifts, and deal with minimum and maximum work time requirements and the introduction of appropriate breaks in a shift. A complete description of the conditions is:

• A day shift cannot exceed 9 hours, but is not allowed to be shorter than 4 hours.
• If a shift exceeds 8 hours, 0.75 break time is required of which 0.5 hour contiguous. If a shift exceeds 5.5 hours, a single break of 0.5 hours or two breaks of 0.25 hours is required. For shorter shifts a single break of 0.25 hours applies.
• Maximum contiguous time an operator can work without a break is 3.25 hours.

However as the subproblem only needs to form a single candidate shift, the subproblem is a very small MIP and solves fast. As a result the above schedule only takes 18 seconds to generate.

The complete model can look like:

$ontext    Column generation for work force scheduling    Erwin Kalvelagen    erwin@amsterdamoptimization.com    References:        Annemieke van Dongen, "Personeelsplanning en kolomgeneratie",           Vrije Universiteit Amsterdam, 2005 $offtext set t 'periods' /t1*t48/; parameter req(t) 'required number of employees' /   t1  2,  t2  5,  t3  6,  t4  5,  t5  7,  t6  7,  t7  4,  t8  6,  t9  7,  t10 7   t11 7,  t12 7,  t13 6,  t14 5,  t15 6,  t16 6,  t17 5,  t18 4,  t19 7,  t20 6   t21 6,  t22 4,  t23 5,  t24 5,  t25 6,  t26 6,  t27 5,  t28 7,  t29 7,  t30 6   t31 9,  t32 6,  t33 7,  t34 6,  t35 6,  t36 7,  t37 7,  t38 4,  t39 5,  t40 5   t41 4,  t42 5,  t43 3,  t44 4,  t45 5,  t46 3,  t47 5,  t48 4 /; set   slen 'shift lengths' /short, medium, long/ ; table shiftlen(slen,*)  'lengths of shifts in 0.25 hours'            minwork maxwork  sumbreak  maxnumbreak   short     16      22        1          1   medium    23      32        2          2   long      33      36        3          2 ; scalar maxwlen 'max periods of contiguous work' /13/; *----------------------------------------------------- * Gilmore-Gomory column generation algorithm *----------------------------------------------------- sets    s    'possible shifts' /shift1*shift1000/    iter 'maximum iterations' /iter1*iter1000/ ; *----------------------------------------------------- * Master model *----------------------------------------------------- parameters    a(t,s)   'matrix growing in dimension s'    c(s)     'cost vector growing in dimension s' ; variables    x(s)    'shifts used'    z       'objective variable' ; integer variable x; x.up(s) = sum(t, req(t)); set ds(s) 'dynamic subset'; equations   masterobj  'cost: total number of periods worked'   edemand(t) 'meet required manpower' ; masterobj..  z =e= sum(ds,c(ds)*x(ds)); edemand(t).. sum(ds, a(t,ds)*x(ds)) =g= req(t); model master /masterobj,edemand/; * reduce amount of information written to the listing file master.solprint = 2; master.limrow = 0; master.limcol = 0; master.solvelink = 2; *----------------------------------------------------- * Sub model *----------------------------------------------------- variables   sigma       'reduced cost'   work(t)     'new shift: quarters worked'   shift(t)    'operator is present'   break(t)    'operator is on break'   totwork     'number of periods worked'   totbreak    'break time'   startshift(t) 'shift starts'   startbreak(t) 'break starts'   tshift(slen) 'type of shift: short, medium, long' ; binary variables work, present, break, startshift, startbreak, tshift; positive variables totwork, totbreak, totpresent; * relax: these variables are automatically integer break.prior(t) = INF; startbreak.prior(t) = INF; startshift.prior(t) = INF; equations   subobj        'objective function'   eshift(t)     'shift = work or break'   etotwork      'calculate total quarter hours worked'   etotbreak     'calculate total quarter hours used for breaks'   estshift(t)   'start of shift'   estbreak(t)   'start of break'   estartsh1     'start shift once'   etshift1      'limit work for short and long shifts'   etshift2      'limit work for short and long shifts'   etshift3      'limit breaks for short and long shifts'   ebreak1       'required break time'   ebreak2       'limits on number of breaks (forces longer breaks)'   emaxcont(t)   'max contiguous working time' ; subobj..       sigma =e= sum(t, (1-edemand.m(t))*work(t)); eshift(t)..    shift(t) =e= work(t) + break(t); etotwork..     totwork =e= sum(t, work(t)); etotbreak..    totbreak =e= sum(t, break(t)); estshift(t)..  startshift(t) =g= shift(t) - shift(t-1); estbreak(t)..  startbreak(t) =g= break(t) - break(t-1); estartsh1..    sum(t, startshift(t)) =e= 1; etshift1..     totwork =g= sum(slen, shiftlen(slen,'minwork')*tshift(slen)); etshift2..     totwork =l= sum(slen, shiftlen(slen,'maxwork')*tshift(slen)); etshift3..     sum(slen, tshift(slen)) =e= 1; ebreak1..      totbreak =e= sum(slen, shiftlen(slen,'sumbreak')*tshift(slen)); ebreak2..      sum(t, startbreak(t)) =l= sum(slen, shiftlen(slen,'maxnumbreak')*tshift(slen)); alias (t,tt); set tmaxw(t,tt); tmaxw(t,tt)$(ord(tt)>=ord(t) and ord(tt)<=ord(t)+maxwlen and ord(t)<=card(t)-maxwlen) = yes; display tmaxw; emaxcont(t)$sum(tmaxw(t,tt),1).. sum(tmaxw(t,tt), work(tt)) =l= maxwlen; model sub /   subobj, eshift, etotwork, etotbreak, estshift, estbreak,   estartsh1, etshift1, etshift2, etshift3, ebreak1, ebreak2,   emaxcont /; equation edummy; variable dummy; edummy.. dummy =e= 0; model subcheck /   edummy, eshift, etotwork, etotbreak, estshift, estbreak,   estartsh1, etshift1, etshift2, etshift3, ebreak1, ebreak2,   emaxcont /; * reduce amount of information written to the listing file sub.solprint = 2; sub.limrow = 0; sub.limcol = 0; sub.solvelink = 2; sub.optcr=0; subcheck.solprint = 2; subcheck.limrow = 0; subcheck.limcol = 0; subcheck.solvelink = 2; *----------------------------------------------------- * initial feasible solution * invent some candidate shifts *----------------------------------------------------- table init(t,s) '1=work,2=break'         shift1   t1      1   t2      1   t3      1   t4      1   t5      1   t6      1   t7      1   t8      1   t9      1   t10     1   t11     1   t12     2   t13     1   t14     1   t15     1   t16     1   t17     1   t18     1   t19     1   t20     1   t21     1   t22     1   t23     1 ; set s2(s) /shift2*shift26/; loop(s2(s),    init(t,s) = init(t-1,s-1); ); display init; * * check initial solution. * loop(s$sum(t$init(t,s),1),    shift.fx(t)$init(t,s) = 1;    work.fx(t)$(init(t,s)=1) = 1;    break.fx(t)$(init(t,s)=2) = 1;    solve subcheck minimizing dummy using mip;    abort$(subcheck.modelstat<>1 and subcheck.modelstat<>8) "Initial solution not correct";    shift.lo(t) = 0;   shift.up(t) = 1;    work.lo(t) = 0;    work.up(t) = 1;    break.lo(t) = 0;   break.up(t) = 1; ); ds(s)$sum(t$init(t,s),1) = 1; a(t,ds)$(init(t,ds)=1) = 1; c(ds) = sum(t$(init(t,ds)=1),1); set nexts(s); nexts(s)$(ord(s)=card(ds)+1) = yes; option a:0; option c:0; display   "------------------------------------------------"   "initial solution",   "------------------------------------------------",   a,c,ds,nexts; *----------------------------------------------------- * column generation algorithm starts here *----------------------------------------------------- scalar done /0/; scalar iteration; loop(iter$(not done), * * solve master problem *    solve master using rmip minimizing z; * * solve knapsack problem *    solve sub using mip minimizing sigma; * * new column found? *    if(sigma.l < -0.001,       a(t,nexts) = work.l(t);       c(nexts) = sum(t, work.l(t));       ds(nexts) = yes;       iteration = ord(iter);       display "--------------------------------------------------------------",               iteration,               "--------------------------------------------------------------",               ds,a;       nexts(s) = nexts(s-1);    else       done = 1;    ); ); abort$(not done) "Too many iterations."; * * solve final mip * display "--------------------------------------------------------------",         "Final MIP",         "--------------------------------------------------------------"; master.optcr=0; solve master using mip minimizing z; display z.l,x.l; *-------------------------------------------------------------- * reporting *-------------------------------------------------------------- scalars    totreq 'total required periods'    totusd 'total allocated periods' ; totreq = sum(t, req(t)); totusd = z.l; display totreq, totusd; alias(s,s3); set curs(s3) /shift1/; scalar k; parameter schedule(s,t); loop(s$(x.l(s)>0.5),    for(k=1 to round(x.l(s)),        schedule(curs,t) = a(t,s);        curs(s3) = curs(s3-1);    ); ); option schedule:0; display schedule;

If speed is everything, often a more brute force can be employed for these type of models: generate all possible shifts and solve just that single set covering problem. Even if done relatively unintelligently with pieces of code like:

for k1  := 1 to 13 do      for k2 := 1 to 13 do         for k3 := 1 to 13 do begin            k := k1+k2+k3;            if (k>=33) and (k<=36) then               writeshift;         end;

the generation part goes very fast. Especially after stepping through the code with a debugger, it is always surprising how fast such code executes. For this problem we only needed to generate 20579 shifts (using my interpretation of the conditions – slightly more strict than the referenced paper). Also the resulting master problem solves very fast. An additional advantage is that this approach provides proven global optimal solutions. In this case we saw that the column generation heuristic gave the same objective as the complete enumeration model.