import math # Import standard library for math functions


def Correction(tau, MTTF, alpha):
    return (1 - 0.1 * alpha * (tau / MTTF) ** 2 + \
        (0.09 * alpha - 0.2) * (tau / MTTF))
    

def lambdaE(tau,MTTF,alpha,isSimple = False): # Default is to use correction formula
    l = (math.gamma(1 / alpha + 1)/MTTF) ** alpha * tau** (alpha - 1)
    if not isSimple:
        l *= Correction(tau, MTTF, alpha)
    return l

# Define parameters
MTTF  =	5*8760	
c_PM  =	15000	
c_CM  =	30000	
MDT	  = 12	
kw	  = 6000	
price = 0.5	
alpha = 4	
c_EP  = price*kw*MDT
c_U   = c_CM+c_EP
tau   = (MTTF/math.gamma(1+1/alpha))*((c_PM/((alpha-1)*c_U))**(1/alpha))
print("\nResults using analytical formula:")
print("tau (hours) =",tau)
print("tau (months) =",tau/730)
print("tau (years) =",tau/8760)

# C() is the objective function to minimize by scipy.optimize
def C(tau,MTTF,alpha,c_PM,c_U, isSimple = False): # Objective function
    return c_PM/tau + c_U*lambdaE(tau, MTTF, alpha, isSimple)

from scipy.optimize import minimize    
MTTF = MTTF/8760 # Numerical routine is more robust for smaler value of the argument
tau0 = MTTF/5 # Initial guess

#First, minimize without correction:
result = minimize(C, tau0, args = (MTTF,alpha,c_PM,c_U,True))
print("\nNumerical minimization with scipy, without correction term in lambdaE()")
print ("tau* (years) =",result.x[0], ", C(tau*) = ", result.fun)

# Then with correction:
result = minimize(C, tau0, args = (MTTF,alpha,c_PM,c_U))
print("\nNumerical minimization with scipy, with correction term in lambdaE()")
print ("tau* (years) =",result.x[0], ", C(tau*) = ", result.fun)



from numpy import arange , zeros
import matplotlib.pyplot as plt
# Create empty vectors to hold values to plot
xlist  = zeros([21])
y1list = zeros([21])
y2list = zeros([21])
y3list = zeros([21])
# Calculate values to plot
i=0
for tau in arange (2,4.1,0.1):
    xlist[i]= tau
    y1list[i]= C(tau ,MTTF ,alpha ,c_PM,c_U)
    y2list[i]= c_PM/tau
    y3list[i]= y1list[i]-y2list[i]
    i+=1
plt.plot (xlist , y1list, label= 'Total')
plt.plot (xlist , y2list, label= 'PM')
plt.plot (xlist , y3list, label= 'Failure')
plt.xlabel (r'$\tau$ (years)')
plt.ylabel (r'$C (\tau )$')
plt.title ('Cost as function of maintenance interval ')
plt.legend()
plt.savefig ('output.svg')
plt.show ()