# how the number of employees affects the married employee scheduling problem

My scheduling problem has a specific pattern of 2 days of work followed by 2 days off. Additionally, there is a week off after every 4 weeks of work. The scheduling matrix (P) has five different schedules, where 0 represents not working, and 1 represents working. The constraints of my problem include pattern assignment constraints and worker efficiency constraints.

P = [[0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0],
[1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1],
[0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1],
[0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0],
[1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]]


My objective function aims to maximize the days spent together by married couples while minimizing the days they both work, to make it easier to take care of the children. I am trying to find a compromise between the needs of the employees and the company by assuring at least a certain percentage of the workforce (e.g., 40%) is present every day throughout the year while providing generous rest days for the employees. The goal of my project is to increase the workforce proportion to its highest while satisfying the objective functions for 100 employees. However, my current algorithm works well only for up to 30 people. I have noticed that for some reason, certain values of N (such as 14 and 19) result in the highest availability (K can increase up to 4), but I am not sure why. Currently I believe it might be related to how I distribute the worker’s efficiency.

NOTATION:

*𝑐_𝑚𝑛: 𝑀𝑎𝑡𝑟𝑖𝑥 𝑤ℎ𝑒𝑟𝑒 𝑝𝑎𝑡𝑡𝑒𝑟𝑛 𝑚 𝑎𝑛𝑑 𝑛 𝑏𝑜𝑡ℎ 𝑎𝑟𝑒 𝑜𝑓𝑓 𝑤𝑜𝑟𝑘

𝑤_𝑚𝑛: 𝑀𝑎𝑡𝑟𝑖𝑥 𝑤ℎ𝑒𝑟𝑒 𝑝𝑎𝑡𝑡𝑒𝑟𝑛 𝑚 𝑎𝑛𝑑 𝑛 𝑏𝑜𝑡ℎ 𝑎𝑟𝑒 𝑎𝑡 𝑤𝑜𝑟𝑘

𝑝𝑎_𝑖𝑗: 𝐵𝑖𝑛𝑎𝑟𝑦 𝑝𝑎𝑟𝑎𝑚𝑒𝑡𝑒𝑟 𝑤ℎ𝑖𝑐ℎ 𝑠ℎ𝑜𝑤𝑠 𝑤ℎ𝑒𝑡ℎ𝑒𝑟 𝑝𝑒𝑟𝑠𝑜𝑛 𝑖 𝑎𝑛𝑑 𝑗 𝑎𝑟𝑒 𝑐𝑜𝑢𝑝𝑙𝑒𝑠 (Assumed half of the workforce was married)

𝑎_𝑖: 𝐸𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑐𝑦 𝑜𝑓 𝑝𝑒𝑟𝑠𝑜𝑛 𝑖(𝐶𝑢𝑟𝑟𝑒𝑛𝑡𝑙𝑦 𝑚𝑜𝑑𝑒𝑙𝑙𝑒𝑑 𝑎𝑠 𝑛𝑜𝑟𝑚𝑎𝑙 𝑑𝑖𝑠𝑡𝑟𝑖𝑏𝑢𝑡𝑖𝑜𝑛)

𝑥_𝑖𝑚: 𝐵𝑖𝑛𝑎𝑟𝑦 𝐷𝑒𝑐𝑖𝑠𝑖𝑜𝑛 𝑉𝑎𝑟𝑖𝑎𝑏𝑙𝑒 𝑤ℎ𝑒𝑡ℎ𝑒𝑟 𝑝𝑒𝑟𝑠𝑜𝑛 𝑖 𝑔𝑒𝑡𝑠 𝑎𝑠𝑠𝑖𝑔𝑛𝑒𝑑 𝑝𝑎𝑡𝑡𝑒𝑟𝑛 𝑚

𝑦_𝑖𝑑: 𝑊ℎ𝑒𝑡ℎ𝑒𝑟 𝑝𝑒𝑟𝑠𝑜𝑛 𝑖 𝑤𝑜𝑟𝑘𝑠 𝑜𝑛 𝑑𝑎𝑦 𝑑

𝑧_𝑚𝑑: Binary parameter whether pattern m works on day d

N: number of employees *

MODEL

Max Σ_𝑖 Σ_𝑗 Σ_𝑚 Σ_𝑛 (𝛼(𝑝𝑎_𝑖𝑗 𝑐_𝑚𝑛 𝑥_𝑖𝑚 𝑥_𝑗𝑛)-(1−𝛼)(𝑝𝑎_𝑖𝑗 𝑤_𝑚𝑛 𝑥_𝑖𝑚 𝑥_𝑗𝑛 ))

s.t.

1. sum(𝑥_𝑖𝑚)==1≫Person can only have one schedule

2. 𝑦_𝑖𝑑 == 𝑠𝑢𝑚(𝑧_𝑚𝑑 𝑥_𝑖𝑚 )≫ 𝐹𝑜𝑟 𝑝𝑒𝑟𝑠𝑜𝑛 𝑖 𝑡𝑜 𝑤𝑜𝑟𝑘 𝑜𝑛 𝑑𝑎𝑦 𝑑

3. sum(𝑦_𝑖𝑑 𝑎_𝑖) ≥𝑘×𝑠𝑢𝑚𝑠𝑢𝑚 𝑊ℎ𝑒𝑟𝑒 𝑘 𝑖𝑠 𝑡ℎ𝑒 𝑝𝑟𝑜𝑝𝑜𝑟𝑡𝑖𝑜𝑛 𝑎𝑛𝑑 𝑠𝑢𝑚𝑠𝑢𝑚 𝑖𝑠 𝑡ℎ𝑒 𝑠𝑢𝑚 𝑜𝑓 𝑒𝑓𝑓𝑖𝑐𝑖𝑛𝑐𝑖𝑒𝑠 𝑜𝑓 𝑒𝑚𝑝𝑙𝑜𝑦𝑒𝑒𝑠 k initially set as 0.3 and improved

4. 𝑢_𝑖𝑗𝑚𝑛≤𝑥_𝑖𝑚 Linearise 𝑢_𝑖𝑗𝑚𝑛≤𝑥_𝑗𝑛 𝑢_𝑖𝑗𝑚𝑛≥𝑥_𝑖𝑚+ 𝑥_𝑗𝑛−1

5. 𝛼 is currently set as 0.7

My questions are:

1. Is there a better way to model the distribution of worker efficiency? Currently, I am using a normal distribution

a =  * N
smart = round(N * 0.021)
ok = round(smart + N * 0.136)
notbad = round(ok + N * 0.682)
dumb = round(bad + N * 0.021)
for j in range(N):
if j <= smart:
a[j] = 5
elif smart < j <= ok:
a[j] = 4
elif ok < j <= notbad:
a[j] = 3
a[j] = 2
elif bad < j <= dumb:
a[j] = 1


2. Why do certain values of N, such as 14 and 19, result in better outcomes in terms of availability (K)? >> the maximum for 14 and 19 was 0.4 while it was only 0.34 for 11 employees

Here is a full version of the code.

from ortools.linear_solver import pywraplp
import numpy as np
import time
import random
import signal
import multiprocessing as mp

class Me:

def __init__(self, kkk, N):

self.hihi(kkk, N)

def hihi(self, kkk, N):

startTime = time.time()

print("When N is ", N)
# D = 35 # five weeks

# while()
D = 35

I = [i for i in range(N)]  # set of workers : 0~39, 40명
# k = 0.3
k = kkk
print("k is ", k)
# Patterns
P = []

#5 DIFFERENT PATTERNS EMPLOYEES CAN BE ASSIGNED TO
P = [[0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0],
[1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1],
[0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1],
[0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0],
[1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

z = P.copy()

# pattern-pattern (c, w)
c = [ * len(P)] * len(P)  # When both patterns are off work
w = [ * len(P)] * len(P)  # When both patterns are on work

for m in range(len(P)):  # mth pattern
c[m] = [sum([(1 - a) * (1 - b) for a, b in zip(P[m], P[n])]) for n in range(len(P))]
w[m] = [sum([a * b for a, b in zip(P[m], P[n])]) for n in range(len(P))]

#OMIT C MATRIX LESS THAN 9 AND W MATRIX BIGGER THAN 5;
#AS EVEN WITH A CONVENTIONAL 5 DAY WEEK, THERE WOULD BE 10 DAYS WHEN THEY ARE BOTH OFF WORK
for m in range(len(P)):
for n in range(len(c[m])):
if c[m][n] < 9 or w[m][n] > 5:

c[m][n] = 0
w[m][n] = 0

print("C matrix", c)
print("W matrix", w)

#EFFICIENCY OF WORKERS OF NORMAL DIST
a =  * N
smart = round(N * 0.021)
ok = round(smart + N * 0.136)
notbad = round(ok + N * 0.682)
dumb = round(bad + N * 0.021)
for j in range(N):
if j <= smart:
a[j] = 5
elif smart < j <= ok:
a[j] = 4
elif ok < j <= notbad:
a[j] = 3
a[j] = 2
elif bad < j <= dumb:
a[j] = 1
print("original a", a)
for kk in range(N):
b = random.randint(0, N - 1)
tmp = a[kk]
a[kk] = a[b]
a[b] = tmp
print("a", a)
sumsum = 0
for aa in a:
sumsum = sumsum + aa
print(sumsum)

# ASSUMED HALF OF THE WORKFORCE IS MARRIED
pa = [ * len(I)] * len(I)
for i in [2 * j for j in range(int(len(I) / 4))]:
pa[i] =  * len(I)
pa[i + 1] =  * len(I)
pa[i][i + 1] = 1
pa[i + 1][i] = 1
# p[2*i+1][2*i] = 1
print("pa", pa)

# Solver
# Create the mip solver with the CBC backend.
solver = pywraplp.Solver.CreateSolver('CBC')

x = []
y = []
u = []

for i in I:
t1 = []
t2 = []
for m in range(len(P)):
# for m in range(20):
t1.append(solver.IntVar(0, 1, ''))
for d in range(len(P)):
# for d in range(35):
t2.append(solver.IntVar(0, 1, ''))
x.append(t1)
y.append(t2)

for i in I:
t1 = []
for j in I:
t2 = []
for m in range(len(P)):
# for m in range(20):
t3 = []
for n in range(len(P)):
# for n in range(20):
t3.append(solver.IntVar(0, 1, ''))
t2.append(t3)
t1.append(t2)
u.append(t1)

#print('Number of variables =', model.NumVari)

# CONSTRAINTS
# pattern assignment constraint
for i in I:

# worker efficiency
for d in range(len(P)):

solver.Add(sum(y[i][d] * a[i] for i in I) >= k * sumsum)

# defining y
for i in I:
for d in range(len(P)):

solver.Add(y[i][d] == sum(z[m][d] * x[i][m] for m in range(len(P))))
#THIS SHOULD BE A SOFT CONSTRAINT INSTEAD OF A HARD CONSTRAINT BUT NEED MULTIPLE FEASIBLE SOLUTIONS

# defining u : linearize!
for i in I:
for j in I:
for m in range(len(P)):
# for m in range(20):
for n in range(len(P)):
# for n in range(20):
solver.Add(u[i][j][m][n] >= x[i][m] + x[j][n] - 1)

# OBJECTIVE FUNCTION
objective_terms = []
alpha = 0.7
for i in I:
for j in I:
for m in range(len(P)):
# for m in range(20):
for n in range(len(P)):
# for n in range(20):
# objective_terms.append(0)
objective_terms.append(alpha * pa[i][j] * c[m][n] * u[i][j][m][n])
objective_terms.append(- (1 - alpha) * pa[i][j] * w[m][n] * u[i][j][m][n])

solver.Maximize(solver.Sum(objective_terms))  # MAXIMISING THE OBJECTIVE FUNCTION

# Solve
solver.SetTimeLimit(100000)#SETTING A TIME LIMIT OF 100 SECONDS
status = solver.Solve()

print("hae")

# Print solution.
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:

print('Optimal objective function = ', solver.Objective().Value(), '\n')
print('Optimal assignment\n')
for i in I:
print('Schedule to worker ', i, ': ', int(sum(m * x[i][m].solution_value() for m in range(len(P)))))
for m in range(len(P)):
if x[i][m].solution_value() == 1:
print(P[m])
runTime = time.time() - startTime
print("--- %s seconds ---" % (runTime))
print("finding better solution")
# if N<11 or k<0.33:
# if N!=11 and N!=12:
if N != 0:
kiki = k + 0.01
kfork = round(kiki, 3)
self.hihi(kfork, N)

else:
print("Sorry we couldn't find any feasible or optimal solution.")

if __name__ == '__main__':

for i in range(10):
for N in range(5, 45):
a = Me(0.3, N)

• It might help if you explained what $K$ and $N$ are and add tags for the language and solver you are using.
– prubin
Jul 1 at 22:00
• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Jul 1 at 22:59
• @prubin Thanks for the suggestion! I have changed it according to your suggestion Jul 3 at 4:41

Just confirming if you are using a binary variable $$u_{ij}^d$$ to count married couple assignment with constraints
$$p_dx_{i,p}+p'_{d}x_{j,p'} \le u_{ij}^d + 1 \ \ \forall p,p' \in P$$
$$u_{ij}^d \le p_dx_{i,p}$$
$$u_{ij}^d \le p_dx_{j,p} \ \ \forall p \ \ \forall d \ \ \forall (i,j) \in$$ Couples
And then minimizing $$\sum_d \sum_{(ij)}u_{ij}^d$$