import math # Import standard library for math functions
from scipy.optimize import minimize
from scipy.integrate import quad as numInt
from maintLib import cBRP, tauBRP, lambdaE
from statLib import RWeibul, pdfWeibul
# Define parameters
MTTF  =	5*8760	
c_PM  =	15000	
c_CM  =	30000	
kw	  = 6000	
price = 0.5	
alpha = 4	
p = 0.1
t_W = 14*24
MRT = 4
MLD = 3
MDT = MRT + MLD + p*t_W/2
c_EP  = price*kw*MDT  # Expected production loss
c_F   = c_CM+c_EP
x = 20000
lmbda = math.gamma(1+1/alpha)/MTTF

tau = tauBRP(MTTF,alpha,c_PM,c_F)
C_avg = cBRP(MTTF,alpha,c_PM,c_F)
print("\ntau* = ",tau)
print("C* =",C_avg)

#
# Integrand for \int_0,t_w f(t)*(t_W-t) dt, i.e., expected # of hours down:
#
def integrand(t, alpha,lmbda, t_W, x): # x = age at start of interval
    return pdfWeibul(t+x,alpha,lmbda)*(t_W-t) / RWeibul(x,alpha,lmbda)


def DeltaC(t,tau,MTTF,alpha,c_PM,c_F,C_avg):
    # Expected cost of an "old" component, i.e., the component has been used for t time units
    W = lambdaE(tau,MTTF,alpha)*tau-lambdaE(t,MTTF,alpha)*t # E(# of failures in [t,tau])
    return W*c_F + c_PM - C_avg*(tau-t)


print("\nPM at time t_0:") # First present results if we do maintenance at time t_0
I,err = numInt(integrand,0,t_W,args=(alpha,lmbda,t_W,0))
q_0 = (1- RWeibul(t_W,alpha,lmbda)) # Probability of failure in [t_0,t_1]
C_t0 = c_PM + q_0*(c_CM + price*kw*(MLD+MRT)) + price*kw*I # PM + Failure + Lost production 
print("Expected cost in [t_0,t_1] =", C_t0)
#
# If a failure during the closed maint. window, the component is as good as new afterwords,
# else we have "used" t_w of life-time. Add cost of "used" life-time, "Delta C"
#
C_t0 +=  (1-q_0)* DeltaC(t_W,tau,MTTF,alpha,c_PM,c_F,C_avg)
print("Expected cost in [t_0,t_1] + expected cost of 'used component' =", C_t0)


print("\nNo PM at time t_0:") # Repeat if we do not maintain before the weather window closes
I, err = numInt(integrand,0,t_W,args=(alpha,lmbda,t_W,x))
q_1 = (1- RWeibul(t_W+x,alpha,lmbda)/RWeibul(x,alpha,lmbda))
C_t1 =  q_1*(c_CM + price*kw*(MLD+MRT)) + price*kw*I
print("Expected cost in [t_0,t_1] =", C_t1)

C_t1 +=  (1-q_1)* DeltaC(x+t_W,tau,MTTF,alpha,c_PM,c_F,C_avg)
print("Expected cost in [t_0,t_1] + expected cost of 'used component' =", C_t1)