How can I minimise a rota or roster from a set of days I know somebody can work?

Question

I have a number of workers, in this case 5, who I want to spread over as many days as possible based on the days they are able to work. So, for example, I have this set up:

workers <- data.frame(worker = c("A", "B", "C", "D", "E"),
                      monday = c(1, 1, 1, 0, 0),
                      tuesday = c(0, 0, 1, 0, 1),
                      wednesday = c(0, 0, 1, 1, 0),
                      thursday = c(1, 0, 1, 0, 0),
                      friday = c(0, 1, 0, 1, 0))

The output of this optimisation I'm looking for is something like the following data frame:

workers <- data.frame(worker = c("A", "B", "C", "D", "E"),
                      monday = c(0, 0, 1, 0, 0),
                      tuesday = c(0, 0, 0, 0, 1),
                      wednesday = c(0, 0, 0, 1, 0),
                      thursday = c(1, 0, 0, 0, 0),
                      friday = c(0, 1, 0, 0, 0))

I should mention that there will sometimes be occasions where it's not possible to do put someone on all 5 days and that doesn't matter, I just like the best possible fit. If that means 2 people working on Monday, 1 person on Tuesday, Wednesday and Friday, but nobody on Thursday then that's fine.

I've taken a look at lpSolve, but I couldn't get it to define the answer for me. With a little manipulation of workers above, it will only shows an output of:

> workers_edit
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    1    0
[2,]    1    0    0    0    1
[3,]    0    1    0    0    0
[4,]    1    1    1    1    0
[5,]    0    0    1    0    1
> lp.assign(workers_edit)
Success: the objective function is 0 
> lp.assign(workers_edit)$solution
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    0    0    0
[2,]    0    0    0    1    0
[3,]    0    0    1    0    0
[4,]    0    0    0    0    1
[5,]    1    0    0    0    0

Answer 1

Not completely clear what the model is, so first fill-in the gaps.

1. Each worker can only work exactly 1 day
2. "best fit" is here defined as: the number of workers working at a given day is most evenly distributed.

Now let's make this more precise with developing a mathematical model:

We define a binary variable:

work(w,t) = 1 if worker w works during day t
            0 otherwise

Also, introduce a variable

 wnum(t) = number of workers active at day t

The most evenly distributed schedule is designed by minimizing the difference between the maximum number of workers wnum(t) and the minimum. So, we end up with:

This is a simple Mixed Integer Programming model that can be solved with any MIP solver. I implemented it here using CVXR:

> library(CVXR)
> 
> workers <- data.frame(worker = c("A", "B", "C", "D", "E"),
+                       monday = c(1, 1, 1, 0, 0),
+                       tuesday = c(0, 0, 1, 0, 1),
+                       wednesday = c(0, 0, 1, 1, 0),
+                       thursday = c(1, 0, 1, 0, 0),
+                       friday = c(0, 1, 0, 1, 0))
> workers
  worker monday tuesday wednesday thursday friday
1      A      1       0         0        1      0
2      B      1       0         0        0      1
3      C      1       1         1        1      0
4      D      0       0         1        0      1
5      E      0       1         0        0      0
> 
> # convert to matrix
> available <- as.matrix(workers[,-1])
> available
     monday tuesday wednesday thursday friday
[1,]      1       0         0        1      0
[2,]      1       0         0        0      1
[3,]      1       1         1        1      0
[4,]      0       0         1        0      1
[5,]      0       1         0        0      0
> 
> # sizes
> m <- nrow(available)
> n <- ncol(available)
> 
> # decision variables
> work <- Variable(m,n,boolean=T)
> wnum <- Variable(n)
> wmax <- Variable(1)
> wmin <- Variable(1)
> 
> # optimization model
> problem <- Problem(Minimize(wmax-wmin),
+                    list(sum_entries(work,axis=1)==1,
+                         wnum==sum_entries(work,axis=2),
+                         work <= available,
+                         wmax >= wnum,
+                         wmin <= wnum))
> 
> result <- solve(problem,verbose=T)
GLPK Simplex Optimizer, v4.47
45 rows, 32 columns, 100 non-zeros
      0: obj =  0.000000000e+000  infeas = 5.000e+000 (10)
*    19: obj =  2.000000000e+000  infeas = 0.000e+000 (2)
*    33: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
OPTIMAL SOLUTION FOUND
GLPK Integer Optimizer, v4.47
45 rows, 32 columns, 100 non-zeros
25 integer variables, all of which are binary
Integer optimization begins...
+    33: mip =     not found yet >=              -inf        (1; 0)
+    34: >>>>>  0.000000000e+000 >=  0.000000000e+000   0.0% (1; 0)
+    34: mip =  0.000000000e+000 >=     tree is empty   0.0% (0; 1)
INTEGER OPTIMAL SOLUTION FOUND
> cat("status:",result$status)
status: optimal
> cat("objective:",result$value)
objective: 0
> print(result$getValue(work))
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    1    0
[2,]    1    0    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    0    1
[5,]    0    1    0    0    0
>