'''

Standard model for optimizing preventive maintenance. Weibull model

Timing belt example (Registerreim)
Note: # is used for comments

'''
import math # Import standard library for math functions

# To create functions in Python, we use "def", and NOT "function" as in many languages
# Arguments are specified in the "header", and the def statement is terminated with colon (:)
# Indentation is used to specify the code to execute, and there is no "end function"
# Use the "return" keyword to pass back the function value
# To calculate powers, i.e, x to the power of n use x**n 
def lambdaE(tau,MTTF,alpha):
    return (math.gamma(1 / alpha + 1)/MTTF) ** alpha * tau** (alpha - 1)
    
# Define parameters for the optimization
C_PM=7000
C_CM=35000
C_EP=0.1*5000
C_ES=0.005*0.2*25000000
C_F = C_CM+C_EP+C_ES
alpha=3
MTTF=175000


#
# To run the minimization routine, we need to specify the cost (objective) function:
# The first argument (here: tau) is always the one we would like to optimize wrt
# Then we can add as many arguments as we want
#
def C_tau(tau,MTTF,alpha,C_PM,C_F):
    return C_PM/tau + lambdaE(tau,MTTF,alpha)*C_F

# The minimize function is found in the scipy.optimize library, need to import:
from scipy.optimize import minimize

MTTF = 175 # Optimizing with very large values is not stable, use 1000 km as basis;
tau0 = 80 # Initial guesss
res=minimize(C_tau, tau0, args = (MTTF,alpha,C_PM,C_F))    


tau = res.x[0] # The resutlt is obtained from the first element, i.e., element 0 of x
print("\n******************\nscipy minimize:\n")
print("tau (1000 km) = ",tau,"C(tau) = ",res.fun)
print("tau (km)      = ",tau*1000,"C(tau) = ",res.fun/1000)