First, we load the required libraries to help us.
library(lpSolve)
library(lpSolveAPI)Next, we specify the objective function direction as well as a vector for the costs. Remember, if we have 10 bucks per hour every day for 40 hours each week, that is 400 bucks. This works for the Monday shift only, though, because that shift doesn’t work the weekends. Tuesday and Sunday shifts have only one weekend day (32 hours x 10 + 8 hours x 15 ) for 440 bucks. All other shift starts cost 480 bucks (24 hours x 10 + 16 hours x 15).
dir="min" #This will be used to tell the solver to minimize rather than maximize.
z=c(400, 440, 480, 480, 480, 480, 440) #Obj Fn: The cost per FTE starting on M, T, W...SIn optimization, the solution is always found in the constraint matrix. The objective function just tells us how much the solution is worth or costs. Remember, we have shifts. So the constraint matrix is just an indicator of whether a shift starting on Day X is working. For example, the Monday shift works Monday through Friday with Saturday and Sunday off. The constraint here is m=c(1,1,1,1,1,0,0), which just tells us who is working when.
m=c(1,1,1,1,1,0,0) #working days for Monday
t=c(0,1,1,1,1,1,0) #working days for Tuesday
w=c(0,0,1,1,1,1,1) #working days for Wednesday
r=c(1,0,0,1,1,1,1) #...
f=c(1,1,0,0,1,1,1) #...
s=c(1,1,1,0,0,1,1) #...
y=c(1,1,1,1,0,0,1) #working days for Sunday
#We can then bind these together into a matrix. I could have done this in one step, but I think it is less confusing to do so step-by-step.
A=matrix(rbind(m,t,w,r,f,s,y), nrow=7)
rownames(A)=c("M","T","W","R","F","S","Sy")
colnames(A)=c("M","T","W","R","F","S","Sy")
A## M T W R F S Sy
## M 1 1 1 1 1 0 0
## T 0 1 1 1 1 1 0
## W 0 0 1 1 1 1 1
## R 1 0 0 1 1 1 1
## F 1 1 0 0 1 1 1
## S 1 1 1 0 0 1 1
## Sy 1 1 1 1 0 0 1To build the right-hand side of the constraints, I simply provide the direction of the constraint and the associated value.
const.dir=c(rep('>=',7)) #Constraint Direction
b=c(17,13,15,19,14,16,11) #The minimum number of workers by dayThe last thing to do is put the pieces together and solve.
#the only tricky thing here is that lpSolveAPI wants the constraints in column order.
x=lp(dir,z,A,const.dir,b,transpose.constraints="F",compute.sen=1)
#To run this as an integer program, do the following.
#y=lp(dir,z,A,const.dir,b,transpose.constraints="F",int.vec=1:7,compute.sen=1) x #Objective Function## Success: the objective function is 10013.33x$solution #Number of Workers by Shift## [1] 6.0000000 5.3333333 0.0000000 7.3333333 0.0000000 3.3333333 0.3333333Shadow prices tell us what a one unit increase in the left-hand side of the constraints would do to the objective function. There is a range for where this analysis is valid. The example below shows that if you increased required FTEs on Monday by 1, that you would have an associated increase of 133 dollars in the objective function. It also shows that the unit increase / decrease is valid for between 11 and 17.5 FTEs.
#
options(scipen=1)
duals=round(x$duals[1:7],2)
duals[1]## [1] 133.33lower=round(x$duals.from[1:7],2)
upper=round(x$duals.to[1:7],2)
lower[1]## [1] 11upper[1]## [1] 17.5Reduced costs tell us how much we could change the cost of a decision variable before it changes the solution set. The example below shows that the solution set does not change on Monday if the cost is between 400 and 440.
round(x$sens.coef.from,2)## [1] 400.00 40.00 480.00 80.00 346.67 80.00 400.00round(x$sens.coef.to,2)## [1] 4.4e+02 4.4e+02 1.0e+30 6.8e+02 1.0e+30 6.8e+02 4.4e+02