import math
def lambdaE(tau, MTTF, alpha):
    return   math.gamma(1 / alpha + 1) ** alpha * tau ** (alpha - 1) / \
                    (MTTF ** alpha) *  \
                    (1 - 0.1 * alpha * (tau / MTTF) ** 2 + (0.09 * alpha - 0.2) * (tau / MTTF))

def pEmpty(s , rho):
    from scipy.stats import poisson
    if s == 0: 
        return  1
    else:
        avg = 0
        n = 20
        for i in range(n):
            avg +=  1/n * (1 - poisson.cdf(s - 1,  rho * (i+0.5)/n))
        return  avg
    

# Main program

MTTF = 8760  # MTTF without preventive maintenance
alpha = 3    # Medium ageing
c_PM = 10000 # Cost per preventive maintenance action
c_CM = 20000 # Cost per corrective maintenance action
c_U = 50000  # Unavailability cost per hour downtime
MDT_W = 4    # Mean downtime with a spare part available
MDT_WO = 72  # Mean downtime without a spare part available
rho = 5      # Rate of withdraws from the stock
c_S = 1      # Cost to keep one spare in stock per hour
best = 1e30
for s in range(11):
    for tau in range (1000, 2500, 25):
        p = pEmpty(s, rho) # Average probability of empty stock
        # calculate cost
        cost = c_PM / tau + lambdaE(tau, MTTF, alpha) * (c_CM + c_U * (p * MDT_WO + (1 - p) * MDT_W)) + s * c_S
        if cost < best:  
            best = cost
            bestTau = tau
            bestS = s

print("tau* =",bestTau)
print("s* =", bestS)
print("C(tau*,s*) =",best)