'''

Load this code into your Python editor, or just open the https://www.onlinegdb.com/ online Python
compiler. Remember to select Language = Python 3

'''
import math
# 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("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): # Objective function
    return c_PM/tau + c_U*(math.gamma(1+1/alpha)/MTTF)**alpha*tau**(alpha-1)

from scipy.optimize import minimize    
MTTF = MTTF/8760 # Numerical routine is more robust for smaler value of the argument
tau0 = MTTF/5 # Initial guess
result = minimize(C, tau0, args = (MTTF,alpha,c_PM,c_U))
print ("tau, scipy (years)=",result.x[0])

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 ()