Rostering problem with column generation
Answered
Hello, I would like to solve the following scheduling model with column generation.
The indices are \(i\in I: Nurse\), \(s \in S: Shifts\) and \(t \in T: Days\). The decision variables are:
- \(x_{its}\)= 1, if nurse \(i\) works shift \(d\) on day \(t\)
- \(slack_{ts}\): slack for shift \(s\) on day \(t\)
- \(motivation_{its}\): for nurse \(i\) in shift \(s\) on day \(t\)
- \(mood_{it}\): for nurse \(i\) on day \(t\)
The parameters / inputs are:
- \(demand_{ts}: \text{ demand per shift $s$ on day $t$}\)
- \(\alpha\): random number \(\in \left[ 0,1 \right]\)
- \(M\): Big number
- \(Max_W\): Number of maximum allowed consecutive working days
Coming to the column generation part:
- \(r\): Index for schedules/roster (a feasible assignment of shifts to a nurse for the planning horizon)
- \(R(i)\): Set of alternative schedules/roster for nurse \(i\)
- \(motivation^r_{its}\): motivation of nurse \(i\) is assigned to shift \(s\) in schedule \(r\) on day \(d\).
This leads to the following master problem:
\[\begin{align}\mathcal{MP}&\quad \min \sum_{t\in T}^{}\sum_{s\in S}^{}slack_{ts}&\\&s.t.\\&\sum_{i\in I}^{}\sum_{r\in R(i)}^{}\lambda _{ir}motivation^r_{its}+slack_{ts} = demand_{ts}\quad &\forall s\in S, t\in T\\&\sum_{r\in R(i)}^{}\lambda _{ir}=1\quad &\forall i\in I\\&slack_{ts}\ge 0\quad &\forall t\in T, s\in S\\&\lambda_{ir}\in \mathbb{Z}^{+}_0\quad &\forall i\in I, r\in R(i)\end{align}\]
and the individual sub problems for every nurse \(i\) go as follows. Note that within the code, the index i is dropped, apart from the objective function:
\[\begin{align}\mathcal{SP(i)}\quad&\min -\sum_{t,s} \pi_{ts} \text{motivation}_{its} - \mu_i \\
&s.t.\\&-M(1-x_{its}) \le \text{motivation}_{its} - \text{mood}_{it} \le M(1-x_{its}) &\forall i\in I,t\in T,s\in S \\&\text{motivation}_{its} \le x_{its} &\forall i\in I,t\in T,s\in S \\&\alpha \sum_s x_{its} + \text{mood}_{it} = 1 &\forall i\in I,t\in T\\&\sum_{l=t}^{t+Max_W}\sum_{s}^{}x_{ils}\le Max_W&\quad\forall i\in I,~t\in \{1,\ldots,T - Max_{W}\}\\&mood_{i1}=1&\forall i\in I\\&\text{motivation}_{its} \in[0,1] &\forall i\in I,t\in T,s\in S \\&\text{mood}_{it} \in[0,1] &\forall i\in I,t\in T \\&x_{its}\in \{0,1\} &\forall t\in T,s\in S\end{align}\]
This is my implementation so far:
from gurobipy import *
import gurobipy as gu
import pandas as pd
import itertools
import time
import matplotlib.pyplot as plt
# Create DF out of Sets
I_list = [1,2,3]
T_list = [1,2,3,4,5,6,7]
K_list = [1,2,3]
I_list1 = pd.DataFrame(I_list, columns=['I'])
T_list1 = pd.DataFrame(T_list, columns=['T'])
K_list1 = pd.DataFrame(K_list, columns=['K'])
DataDF = pd.concat([I_list1, T_list1, K_list1], axis=1)
Demand_Dict = {(1, 1): 2, (1, 2): 1, (1, 3): 0, (2, 1): 1, (2, 2): 2, (2, 3): 0, (3, 1): 1, (3, 2): 1, (3, 3): 1,
(4, 1): 1, (4, 2): 2, (4, 3): 0, (5, 1): 2, (5, 2): 0, (5, 3): 1, (6, 1): 1, (6, 2): 1, (6, 3): 1,
(7, 1): 0, (7, 2): 3, (7, 3): 0}
class MasterProblem:
def __init__(self, dfData, DemandDF, iteration):
self.iteration = iteration
self.physicians = dfData['I'].dropna().astype(int).unique().tolist()
self.days = dfData['T'].dropna().astype(int).unique().tolist()
self.shifts = dfData['K'].dropna().astype(int).unique().tolist()
self.roster = list(range(1,self.iteration+2))
self.demand = DemandDF
self.model = gu.Model("MasterProblem")
self.cons_demand = {}
self.cons_lmbda = {}
def buildModel(self):
self.generateVariables()
self.generateConstraints()
self.generateObjective()
self.setStartSolution()
self.modelFlags()
self.model.update()
def generateVariables(self):
self.slack = self.model.addVars(self.days, self.shifts, vtype=gu.GRB.CONTINUOUS, lb=0, name='slack')
self.motivation_i = self.model.addVars(self.physicians, self.days, self.shifts, self.roster, vtype=gu.GRB.CONTINUOUS, lb= 0, ub= 1, name = 'motivation_i')
self.lmbda = self.model.addVars(self.physicians, self.roster, vtype=gu.GRB.INTEGER, lb = 0, name = 'lmbda')
def generateConstraints(self):
for i in self.physicians:
self.cons_lmbda[i] = self.model.addLConstr(gu.quicksum(self.lmbda[i, r] for r in self.roster) == 1)
for t in self.days:
for s in self.shifts:
self.cons_demand[t,s] = self.model.addConstr(gu.quicksum(self.motivation_i[i, t, s, r] for i in self.physicians for r in self.roster) + self.slack[t, s] >= self.demand[t, s])
return self.cons_lmbda, self.cons_demand
def generateObjective(self):
self.model.setObjective(gu.quicksum(self.slack[t, s] for t in self.days for s in self.shifts), sense = gu.GRB.MINIMIZE)
def solveRelaxModel(self):
for v in self.model.getVars():
v.setAttr('vtype', 'C')
self.model.optimize()
def getDuals_i(self):
Pi_cons_lmbda = self.model.getAttr("Pi", self.cons_lmbda)
return Pi_cons_lmbda
def getDuals_ts(self):
Pi_cons_demand = self.model.getAttr("Pi", self.cons_demand)
return Pi_cons_demand
def getObjValues(self):
obj = self.model.objVal
return obj
def updateModel(self):
self.model.update()
def addColumn(self, newSchedule, iter, i):
ctName = f"ScheduleUsed[{i},{iter}]"
newColumn = gu.Column(newSchedule, self.model.getConstrs())
self.model.addVar(vtype=gu.GBR.CONTINOUS, lb=0, obj=1.0, column=newColumn, name=ctName)
self.model.update()
def setStartSolution(self):
startValues = {}
for i, t, s, r in itertools.product(self.physicians, self.days, self.shifts, self.roster):
startValues[(i, t, s, r)] = 0
for i, t, s, r in startValues:
self.motivation_i[i, t, s, r].Start = startValues[i, t, s, r]
def modelFlags(self):
self.model.Params.OutputFlag = 0
def solveModel(self, timeLimit, EPS):
self.model.setParam('TimeLimit', timeLimit)
self.model.setParam('MIPGap', EPS)
self.model.optimize()
class Subproblem:
def __init__(self, duals_i, duals_ts, dfData, i):
self.days = dfData['T'].dropna().astype(int).unique().tolist()
self.shifts = dfData['K'].dropna().astype(int).unique().tolist()
self.duals_i = duals_i
self.duals_ts = duals_ts
self.Max = 5
self.M = 100
self.alpha = 0.5
self.model = gu.Model("Subproblem")
self.i = i
def buildModel(self):
self.generateVariables()
self.generateConstraints()
self.generateObjective()
self.modelFlags()
self.model.update()
def generateVariables(self):
self.x = self.model.addVars(self.days, self.shifts, vtype=GRB.BINARY, name='x')
self.mood = self.model.addVars(self.days, vtype=GRB.CONTINUOUS, lb = 0, ub = 1, name='mood')
self.motivation = self.model.addVars(self.days, self.shifts, vtype=GRB.CONTINUOUS, lb = 0, ub = 1, name='motivation')
def generateConstraints(self):
for t in self.days:
self.model.addConstr(self.mood[t] == 1 - self.alpha * quicksum(self.x[t, s] for s in self.shifts))
for t in self.days:
self.model.addConstr(gu.quicksum(self.x[t, s] for s in self.shifts) <= 1)
for t in range(1, len(self.days) - self.Max + 1):
self.model.addConstr(gu.quicksum(self.x[u, s] for s in self.shifts for u in range(t, t + 1 + self.Max)) <= self.Max)
for t in self.days:
for s in self.shifts:
self.model.addConstr(self.mood[t] + self.M*(1-self.x[t, s]) >= self.motivation[t, s])
self.model.addConstr(self.motivation[t, s] >= self.mood[t] - self.M * (1 - self.x[t, s]))
self.model.addConstr(self.motivation[t, s] <= self.x[t, s])
def generateObjective(self):
self.model.setObjective(0-gu.quicksum(self.motivation[t,s]*self.duals_ts[t,s] for t in self.days for s in self.shifts)-self.duals_i[self.i], sense = gu.GRB.MINIMIZE)
def getNewSchedule(self):
return self.model.getAttr("X", self.model.getVars())
def getObjVal(self):
obj = self.model.getObjective()
value = obj.getValue()
return value
def getStatus(self):
return self.model.status
def solveModel(self, timeLimit, EPS):
self.model.setParam('TimeLimit', timeLimit)
self.model.setParam('MIPGap', EPS)
self.model.optimize()
def modelFlags(self):
self.model.Params.OutputFlag = 0
#### Column Generation
# Prerequisites
modelImprovable = True
objValHist = []
t0 = time.time()
max_itr = 3000
itr = 0
## Start
print(' *****Column Generation Iteration***** \n')
while (modelImprovable) and itr < max_itr:
# Build MP
master = MasterProblem(DataDF, Demand_Dict, itr)
master.buildModel()
# Start
itr += 1
print('current iteration time: ', itr)
# Solve RMP
master.updateModel()
master.solveRelaxModel()
objValHist.append(master.getObjValues)
print('Current rmp objval: ', objValHist)
# Get Duals
duals_i = master.getDuals_i()
duals_ts = master.getDuals_ts()
# Solve SPs
for i in I_list:
subproblem = Subproblem(duals_i, duals_ts, DataDF, i)
subproblem.buildModel()
subproblem.solveModel(3600, 1e-6)
if subproblem.getStatus != GRB.status.OPTIMAL:
raise Exception("Pricing-Problem can not reach optimal!")
reducedCost = subproblem.getObjVal()
print('reduced cost', reducedCost)
if reducedCost < -1e-6:
ScheduleCuts = subproblem.getNewSchedule()
master.addColumn(ScheduleCuts, itr)
else:
modelImprovable = False
# Solve MP
master.solveModel(3600, 0.01)
# Results
print('*** Results ***')
print('Total iteration: ', itr)
t1 = time.time()
print('Total elapsed time: ', t1 - t0)
print('Exact solution cost:', master.getAttr("ObjVal"))
# Plot
plt.scatter(list(range(len(rmp_objvals))), rmp_objvals, c='r')
plt.xlabel('history')
plt.ylabel('objective function value')
title = 'solution: ' + str(rmp_objvals[-1])
plt.title(title)
plt.show()
Obviously, the code does not produce a solution, leading to the following problems / questions:
- How do i fix the duals error?
- How do i model the dynamic index r?
- Are the columns, i.e. the optimal values of motivation of the subproblem, properly added as columns to the MP?
- How do i assign an physician i in the MP to the newly added column (roster)?
- How to retrieve to optimal values for x from the final SP iterations?
Thank you so much for your help!
-
Hi Carl,
My suggestion would be to try and implement this step-by-step. Looking at your code, I get the impression some parts are not yet valid Python (some syntax errors, and variables referenced without creating them) while there's already a lot of complexity... start small, test, then iterate. I would start with the RMP (only!) and a trivial solution as initial column(s). First make sure that part is solved correctly before working on the subproblems.
Kind regards,
Ronald0 -
Ronald van der Velden Sorry, i accidentally added some code with errors. I fixed some of them; some of them still remain, however. What exactly do you mean by solving the RMP only? Just initially solving the objective with the demand constraint? The problem is, i dont really know how to create a trivial (initial) solution. Excuse me, but I am a real beginner, but I have to solve this problem somehow. I am therefore grateful for any help.
0 -
Exactly, without "real" rosters produced by a subproblem. I believe since you have slack variables, there is always a solution to the master problem right? So you could generate the obvious empty roster (no shifts) for each nurse and include that in your initial master problem. Or any heuristic approach where you fill a roster left-to-right using a few shifts.
0 -
Ronald van der Velden I see. How would the code look to implement this?
0 -
Hi Carl,
Please give it a try yourself and feel free to ask a more specific question once you need help with a particular step.
Kind regards,
Ronald0 -
Ronald van der Velden I have just updated my code and adapted the structure. Unfortunately, I still can't manage to link the dynamic index r in the master problem with the variable \(motivation_{its}\) and to add the columns from the individual subproblems. Are the individual subproblems called correctly?
I would be infinitely grateful if you could help me solve these problems.0 -
Just to be sure - I think the current issue is that you solve the relaxed model, but then request "Pi" for the original model (which has not been solved, hence Pi is not available there). Are you observing the same on your side?
0 -
Ronald van der Velden Yes, i updated the code, the error is now gone, but now the Constr cant be found.
0 -
Please have a look at (the end of) this question about the same problem.
0 -
Ronald van der Velden Thanks, that did help. I also fixed a few other errors, but now i run into this error
self.model.setObjective(-gu.quicksum(self.motivation[t,s]*self.duals_ts[t,s] for t in self.days for s in self.shifts)-self.duals_i[i], sense = gu.GRB.MINIMIZE)
~~~~~~~~~~~~~^^^^^
KeyError: (1, 1)I guess the problem is, that the duals are not available yet? However, at this point, i dont really know where i should look at into the code, to fix this problem. I also updated my current problems and where i still need help!
0 -
Ronald van der Velden And i also guess that the column added is the wrong on, as i am "giving" over the whole solution matrix, but obvioulsy i only want to hand over the \(motivation_{ts}\) values for all nurses \(i\) to the suproblem and assign the index \(r\). Maybe replace
def getNewSchedule(self): return self.model.getAttr("X", self.model.getVars())
with
def getNewSchedule(self):
sol = self.motivation_i.values()
sol_dict = dict((var.VarName, var.X) for var in sol)
return sol_dict0 -
Hi Carl,
Unfortunately I can't help you debug this step-by-step. I don't believe the error you see is optimization-specific: you should be able to add breakpoints to your code and inspect variable values, see where you set and retrieve duals_ts etcetera and this will help you solve the problem.
Kind regards,
Ronald0 -
Ronald van der Velden Of course it just was a stupid double return in the variable function. Now the optimization runs through but obviously there is now optimal solution. Unfortunately, I'm at my wit's end now, as I'm really not good at programming (see this post) and don't know how the individual functions in Gruobi work. Should I create a new post for the sake of clarity?
0 -
don't know how the individual functions in Gruobi work
Could this help?
0
Please sign in to leave a comment.
Comments
14 comments