I'm doing a model for scheduling personnel through VRP formulation and I would like to know how can I integrate precedence binary matrix, using the CSR to accelerate the data processing.
My precedence matrix will variate per day, so I don't know how to read it within the model.
Also I would to confirm if i am not forcing the utilization of all staff, because the idea is to use as less as possible.
Here is my code:
# Parameters
N = 5 Personnel
P = 2 Package
J = 21 Tasks
S = 3 Shifts
D = 7 Days
Prec = [[0,1,0,0,0,0,0],
[0,0,0,0,0,1,0],
[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0],
[0,1,0,0,0,0,0],
[0,0,0,0,0,0,0],
[0,1,0,0,0,0,0]]
random.seed(8)
A = {}
for n in range(1, N+1):
A[n] = {}
for j in range(1, J+1):
A[n][j] = random.randint(30,50)
for j in range(J+1, 2*N+J+1): # Dummy tasks
A[n][j] = 0
SL = {1: 480, 2: 960, 3: 1440} # Shift time limits
# Package id per task
ID = {
1: {1: 1, 2: 2, 3: 1, 4: 2, 5: 1, 6: 2, 7: 1},
2: {1: 2, 2: 1, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
3: {1: 1, 2: 1, 3: 1, 4: 2, 5: 2, 6: 1, 7: 2},
4: {1: 2, 2: 2, 3: 1, 4: 1, 5: 2, 6: 1, 7: 2},
5: {1: 1, 2: 2, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
6: {1: 2, 2: 1, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
7: {1: 1, 2: 2, 3: 1, 4: 2, 5: 1, 6: 2, 7: 1},
8: {1: 1, 2: 2, 3: 1, 4: 2, 5: 1, 6: 2, 7: 1},
9: {1: 2, 2: 1, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
10: {1: 1, 2: 1, 3: 1, 4: 2, 5: 2, 6: 1, 7: 2},
11: {1: 2, 2: 2, 3: 1, 4: 1, 5: 2, 6: 1, 7: 2},
12: {1: 1, 2: 2, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
13: {1: 2, 2: 1, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
14: {1: 1, 2: 2, 3: 1, 4: 2, 5: 1, 6: 2, 7: 1},
15: {1: 1, 2: 2, 3: 1, 4: 2, 5: 1, 6: 2, 7: 1},
16: {1: 2, 2: 1, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
17: {1: 1, 2: 1, 3: 1, 4: 2, 5: 2, 6: 1, 7: 2},
18: {1: 2, 2: 2, 3: 1, 4: 1, 5: 2, 6: 1, 7: 2},
19: {1: 1, 2: 2, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
20: {1: 2, 2: 1, 3: 2, 4: 1, 5: 2, 6: 1, 7: 2},
21: {1: 1, 2: 2, 3: 1, 4: 2, 5: 1, 6: 2, 7: 1}
}
IDS = {
1: {1: 1, 2: 2, 3: 1, 4: 3, 5: 1, 6: 2, 7: 3},
2: {1: 2, 2: 2, 3: 1, 4: 3, 5: 3, 6: 1, 7: 2},
3: {1: 3, 2: 3, 3: 3, 4: 2, 5: 1, 6: 3, 7: 2},
4: {1: 3, 2: 1, 3: 2, 4: 3, 5: 2, 6: 1, 7: 3},
5: {1: 2, 2: 1, 3: 3, 4: 1, 5: 2, 6: 3, 7: 1},
6: {1: 1, 2: 3, 3: 1, 4: 2, 5: 3, 6: 2, 7: 1},
7: {1: 1, 2: 2, 3: 2, 4: 3, 5: 1, 6: 3, 7: 2},
8: {1: 1, 2: 2, 3: 1, 4: 3, 5: 1, 6: 2, 7: 3},
9: {1: 2, 2: 2, 3: 1, 4: 3, 5: 3, 6: 1, 7: 2},
10: {1: 3, 2: 3, 3: 3, 4: 2, 5: 1, 6: 3, 7: 2},
11: {1: 3, 2: 1, 3: 2, 4: 3, 5: 2, 6: 1, 7: 3},
12: {1: 2, 2: 1, 3: 3, 4: 1, 5: 2, 6: 3, 7: 1},
13: {1: 1, 2: 3, 3: 1, 4: 2, 5: 3, 6: 2, 7: 1},
14: {1: 1, 2: 2, 3: 2, 4: 3, 5: 1, 6: 3, 7: 2},
15: {1: 1, 2: 2, 3: 1, 4: 3, 5: 1, 6: 2, 7: 3},
16: {1: 2, 2: 2, 3: 1, 4: 3, 5: 3, 6: 1, 7: 2},
17: {1: 3, 2: 3, 3: 3, 4: 2, 5: 1, 6: 3, 7: 2},
18: {1: 3, 2: 1, 3: 2, 4: 3, 5: 2, 6: 1, 7: 3},
19: {1: 2, 2: 1, 3: 3, 4: 1, 5: 2, 6: 3, 7: 1},
20: {1: 1, 2: 3, 3: 1, 4: 2, 5: 3, 6: 2, 7: 1},
21: {1: 1, 2: 2, 3: 2, 4: 3, 5: 1, 6: 3, 7: 2},
}
E = {
1: {1: 0, 2: 481, 3: 0, 4: 961, 5: 0, 6: 481, 7: 961},
2: {1: 481, 2: 481, 3: 0, 4: 961, 5: 961, 6: 0, 7: 481},
3: {1: 961, 2: 961, 3: 961, 4: 481, 5: 0, 6: 961, 7: 481},
4: {1: 961, 2: 0, 3: 481, 4: 961, 5: 481, 6: 0, 7: 961},
5: {1: 481, 2: 0, 3: 961, 4: 0, 5: 481, 6: 961, 7: 0},
6: {1: 0, 2: 961, 3: 0, 4: 481, 5: 961, 6: 481, 7: 0},
7: {1: 0, 2: 481, 3: 481, 4: 961, 5: 0, 6: 961, 7: 481},
8: {1: 0, 2: 481, 3: 0, 4: 961, 5: 0, 6: 481, 7: 961},
9: {1: 481, 2: 481, 3: 0, 4: 961, 5: 961, 6: 0, 7: 481},
10: {1: 961, 2: 961, 3: 961, 4: 481, 5: 0, 6: 961, 7: 481},
11: {1: 961, 2: 0, 3: 481, 4: 961, 5: 481, 6: 0, 7: 961},
12: {1: 481, 2: 0, 3: 961, 4: 0, 5: 481, 6: 961, 7: 0},
13: {1: 0, 2: 961, 3: 0, 4: 481, 5: 961, 6: 481, 7: 0},
14: {1: 0, 2: 481, 3: 481, 4: 961, 5: 0, 6: 961, 7: 481},
15: {1: 0, 2: 481, 3: 0, 4: 961, 5: 0, 6: 481, 7: 961},
16: {1: 481, 2: 481, 3: 0, 4: 961, 5: 961, 6: 0, 7: 481},
17: {1: 961, 2: 961, 3: 961, 4: 481, 5: 0, 6: 961, 7: 481},
18: {1: 961, 2: 0, 3: 481, 4: 961, 5: 481, 6: 0, 7: 961},
19: {1: 481, 2: 0, 3: 961, 4: 0, 5: 481, 6: 961, 7: 0},
20: {1: 0, 2: 961, 3: 0, 4: 481, 5: 961, 6: 481, 7: 0},
21: {1: 0, 2: 481, 3: 481, 4: 961, 5: 0, 6: 961, 7: 481}
}
L = {
1: {1: 480, 2: 960, 3: 480, 4: 1440, 5: 480, 6: 960, 7: 1440},
2: {1: 960, 2: 960, 3: 480, 4: 1440, 5: 1440, 6: 480, 7: 960},
3: {1: 1440, 2: 1440, 3: 1440, 4: 960, 5: 480, 6: 1440, 7: 960},
4: {1: 1440, 2: 480, 3: 960, 4: 1440, 5: 960, 6: 480, 7: 1440},
5: {1: 960, 2: 480, 3: 1440, 4: 480, 5: 960, 6: 1440, 7: 480},
6: {1: 480, 2: 1440, 3: 480, 4: 960, 5: 1440, 6: 960, 7: 480},
7: {1: 480, 2: 960, 3: 960, 4: 1440, 5: 480, 6: 1440, 7: 960},
8: {1: 480, 2: 960, 3: 480, 4: 1440, 5: 480, 6: 960, 7: 1440},
9: {1: 960, 2: 960, 3: 480, 4: 1440, 5: 1440, 6: 480, 7: 960},
10: {1: 1440, 2: 1440, 3: 1440, 4: 960, 5: 480, 6: 1440, 7: 960},
11: {1: 1440, 2: 480, 3: 960, 4: 1440, 5: 960, 6: 480, 7: 1440},
12: {1: 960, 2: 480, 3: 1440, 4: 480, 5: 960, 6: 1440, 7: 480},
13: {1: 480, 2: 1440, 3: 480, 4: 960, 5: 1440, 6: 960, 7: 480},
14: {1: 480, 2: 960, 3: 960, 4: 1440, 5: 480, 6: 1440, 7: 960},
15: {1: 480, 2: 960, 3: 480, 4: 1440, 5: 480, 6: 960, 7: 1440},
16: {1: 960, 2: 960, 3: 480, 4: 1440, 5: 1440, 6: 480, 7: 960},
17: {1: 1440, 2: 1440, 3: 1440, 4: 960, 5: 480, 6: 1440, 7: 960},
18: {1: 1440, 2: 480, 3: 960, 4: 1440, 5: 960, 6: 480, 7: 1440},
19: {1: 960, 2: 480, 3: 1440, 4: 480, 5: 960, 6: 1440, 7: 480},
20: {1: 480, 2: 1440, 3: 480, 4: 960, 5: 1440, 6: 960, 7: 480},
21: {1: 480, 2: 960, 3: 960, 4: 1440, 5: 480, 6: 1440, 7: 960}
}
# Define problem
prob = LpProblem("Nurse_Scheduling", LpMinimize)
# Decision variables
X = LpVariable.dicts("X", [(n, j1, j2, d) for n in range(1, N+1) for j1 in range(1, J+N+1) for j2 in range(1, J+2*N+1) for d in range(1, D+1)], cat=LpBinary) #If nurse i performs j1 before j2
Y = LpVariable.dicts("Y", [(n, s, d) for n in range(1, N+1) for s in range(1, S+1) for d in range(1, D+1)], cat=LpBinary) # if nurse i is allocated in shift k
ST = LpVariable.dicts("ST",[(j, d) for j in range(1, J+2*N+1) for d in range(1, D+1)], lowBound=0, cat=LpContinuous) # Start time for job j
W = LpVariable.dicts("W",[(n, d) for n in range(1, N+1) for d in range(1, D+1)], lowBound=0, upBound=SL[1], cat=LpContinuous) #workload of the nurse i
SP = LpVariable.dicts("SP",[(n, z) for n in range(1, N+1) for z in range(1, Z+1)], lowBound=0, cat=LpBinary) #workload of the nurse i
LateSlack = LpVariable.dicts("LateSlack", range(1, J+1), lowBound=0, cat=LpContinuous)
MaxWorkload = LpVariable.dicts("MaxWorkload", range(1, D+1), lowBound=0, cat=LpContinuous)
MinWorkload = LpVariable.dicts("MinWorkload", range(1, D+1), lowBound=0, cat=LpContinuous)
# Objective function
prob += 10000 * lpSum(Y[n, s, d] for n in range(1, N+1) for s in range(1, S+1) for d in range(1, D+1))
# Constraints
for n in range(1, N+1):
for d in range(1, D+1):
prob += MaxWorkload[d] >= W[n,d]
prob += MinWorkload[d] <= W[n,d]
for n in range(1, N+1):
for d in range(1, D+1):
prob += lpSum(Y[n, s, d] for s in range(1, S+1)) <= 1
for s in range(1, S+1):
for d in range(1, D+1):
prob += lpSum(Y[n, s, d] for n in range(1, N+1)) >= 1
for j in range(1, J+1):
for d in range(1, D+1):
prob += lpSum(X[n, j1, j, d] for n in range(1, N+1) for j1 in list(range(1, J+1)) + [n+J]) == 1 # Each job has just one precedessor
prob += lpSum(X[n, j, j1, d] for n in range(1, N+1) for j1 in list(range(1, J+1)) + [n+J+N]) == 1 # Each job has just one successor
for n in range(1, N+1):
for j in range(1, J+1):
for d in range(1, D+1):
prob += lpSum(X[n, j, j1, d] for j1 in list(range(1, J+1)) + [n+J+N]) == lpSum(X[n, j1, j, d] for j1 in list(range(1, J+1)) + [n+J]) #Route continuoity
for n in range(1, N+1):
for d in range(1, D+1):
prob += lpSum(X[n, n+J, j, d] for j in range(1, J+1)) <= 1 # Each dummy job is attached to one nurse a has just one successor
prob += lpSum(X[n, j, n+J+N, d] for j in range(1, J+1)) <= 1 # Each dummy job is attached to one nurse a has just one precedessor
prob += lpSum(A[n][j] * X[n, j, j1, d] for j in range(1, J+1) for j1 in list(range(1, J+1)) + [n+J+N]) == W[n,d] # Workload equivalence
for n in range(1, N+1):
for s in range(1, S+1):
for d in range(1, D+1):
prob += lpSum(X[n, j, j1, d] for j in range(1, J+1) for j1 in list(range(1, J+1)) + [n+J+N] if IDS[j][d] == s) >= Y[n, s, d]
prob += lpSum(X[n, j, j1, d] for j in range(1, J+1) for j1 in list(range(1, J+1)) + [n+J+N] if IDS[j][d] == s) <= J * Y[n, s, d]
prob += lpSum(X[n, j, j1, d] for j in list(range(1, J+1)) + [n+J] for j1 in list(range(1, J+1)) if IDS[j1][d] == s) <= J * Y[n, s, d]
for j1 in range(1, J+1):
for j2 in range(1, J+1):
for d in range(1, D+1):
if IDS[j1][d] == IDS[j2][d] and ID[j1][d] == ID[j2][d]:
prob += lpSum(X[n, j1, k, d] for n in range(1, N+1) for k in list(range(1, J+1)) + [n+J+N]) == \
lpSum(X[n, j2, k, d] for n in range(1, N+1) for k in list(range(1, J+1)) + [n+J+N])
for n in range(1, N+1):
for j1 in list(range(1, J+1)) + [n+J]:
for j2 in list(range(1, J+1)) + [n+J+N]:
for d in range(1, D+1):
prob += ST[j1, d] + A[n][j1] <= ST[j2, d] + 10000 * (1 - X[n, j1, j2, d]) # Time accumulation
for n in range(1, N+1):
for s in range(1, S+1):
for d in range(1, D+1):
prob += ST[n+J+N, d] <= SL[s] + SL[S] * (1 - Y[n, s, d])
prob += ST[n+J, d] >= (SL[s] - SL[1]) * Y[n, s, d]
# Solve the problem
solver = getSolver('CPLEX_CMD', timeLimit=30)
prob.solve(solver)
print("Status:", LpStatus[prob.status])
print("Objective:", value(prob.objective))