The Expected Total Number of Tardy JobsΒΆ

We compute the expected total number of tardy jobs by means of the tools in Section 4.3 of [sarin10:_stoch_sched].

import numpy as np
from scipy.stats import norm
from scipy.integrate import quad

class Job(object):
    def __init__(self, mean, variance, due):
        self.mean = mean
        self.variance = variance
        self.due = due
        self.C = 0 # completion time
        self.varC = 0
        self.rv = None

The example of Table 4.1 of Sarin. The schedule below eq. 4.28 is contained in s.

job1 = Job(40., 15., 120)
job2 = Job(20., 15., 80)
job3 = Job(60., 20., 90)
job4 = Job(25., 5., 60)

s = [job4, job2, job1, job3]

Now compute the completion times.

j1 = s[0]
j1.C = j1.mean
j1.varC = j1.variance
for j1, j2 in zip(s[:-1], s[1:]):
    j2.C = j1.C + j2.mean
    j2.varC = j1.varC + j2.variance

# F_{[C[j]]}(d_{[j]})

for j in s:
    j.rv = norm(loc = j.C, scale = np.sqrt(j.varC) )

The number of tardy jobs, see eq 4.24 and 4.26:

E = sum( j.rv.sf(j.due) for j in s)
print E
1.00000000165

The variance

V = sum( j.rv.sf(j.due)*j.rv.cdf(j.due) for j in s )

And the final result.

for i in range(len(s)-1):
    j1 = s[i]
    muX = j1.C
    varX = j1.varC
    sigmaX = np.sqrt(varX)
    for j in range(i+1, len(s)):
        j2 = s[j]
        muY = sum(k.mean for k in s[i+1:j+1] )
        varY = sum(k.variance for k in s[i+1:j+1] )
        sigmaY = np.sqrt(varY)
        res = quad(lambda x: norm.pdf( (x-muX)/sigmaX )*\
                   norm.cdf( (x-j2.due+muY)/sigmaY ), j1.due, np.inf)
        EUU = res[0]/sigmaX
        V += EUU
        V -= j1.rv.sf(j1.due)*j2.rv.sf(j2.due)
print V
1.64856166217e-09