'''

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_EP = 0 # Do not consider "punctuality" now, see Problem 11.2 
C_ES = 0 # Do not consider "safety" now, see Problem 11.2
C_F = C_CM+C_EP+C_ES
alpha=3
MTTF=175000

# Calculate the expected cost for various tau values, the "range()" function 
# takes "low", "high" and "step" as arguments
x_i=[]
c_PM_i=[]
c_F_i=[]
c_T_i=[]
best=1e30
for tau in range(80000,110000,1000):
    c=C_PM/tau + lambdaE(tau,MTTF,alpha)*C_F
    if c < best:
        best = c
        best_tau = tau
    print("tau = ",tau,"C(tau) = ",c)
    x_i.append(tau)
    c_PM_i.append(C_PM/tau)
    c_F_i.append(lambdaE(tau,MTTF,alpha)*C_F)
    c_T_i.append(c)
print("\ntau*= ", best_tau)
#
# 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

# The minimize() function takes the following arguments: 
#    <function name>, i.e., C_tau
#    <initial value>, i.e., tau0
#    <arguments>, i.e., arguments/parameters to pass to the function to minimize
#
#    The syntax for the arguments are, args = (arg1, arg2, ...)
res=minimize(C_tau, tau0, args = (MTTF,alpha,C_PM,C_F))    

# The minimize function returns a "result"-variable, comprising several attributes
# The most important ones are:
#
# res.x = the argument that minimizes the objective (cost) function. 
# x could be a vector, i.e., in principle we can minimize wrt several x-values
#
# res.fun = the value of the objective function for the returned x-value

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)

import matplotlib.pyplot as plt


 

# Create a figure with two vertically stacked subplots

plt.plot(x_i, c_PM_i,label= 'PM')
plt.plot(x_i, c_F_i,label= 'CM')
plt.plot(x_i, c_T_i,label= 'Total')
plt.legend()
plt.show()