Thursday, March 10, 2011

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.

schedule

The master problem is simply:

image

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.