'''


'''

from maintLib import lambdaE, tauBRP
import numpy as np
import math
from scipy.optimize import minimize
MTTF_A = 8000
q_B = 0.1
MTTF_B = 16000
alpha = 3
c_PM = 4000
c_CM = 9000
c_FT = 1000
S = 750
MDT = 8
c_U = 25000
tau_A_O = 1000
tau_B_O = 14000
tau_B_FT = 7000
c_F = c_CM + q_B * c_U * MDT

print("\na) Since we have the Weibull distribution, replacing / overhauling the pump will reduce")
print("the effective failure rate. In adition, since pump B has a hidden function, it is important")
print("to reveal a hidden failure as soon as possilbe, hence we also do a proof test")
print("\nb) Calculation of q_B with current maint. intervall for pump B:")
print("\nlambda_E,B =",lambdaE(tau_B_O,MTTF_B,alpha))
print("q_B =",lambdaE(tau_B_O,MTTF_B,alpha)*tau_B_FT/2)
print("q_B =",lambdaE(tau_B_O,MTTF_B,alpha,True)*tau_B_FT/2, "(with simple formula for lambda_E")
print("\nThe result is closed to the assumption in the beginning where q_B = 10%")
print("\nThe effective failure rate for pump B is found by the standard formula,")
print("implemented by e.g., lambdaE(tau_B_O,MTTF_B,alpha)")
print("\nq_B is then found by lambdaE(tau_B_O,MTTF_B,alpha)*tau_B_FT/2")

def C_2(tau,MTTF_A,MTTF_B,alpha,c_PM,c_FT,c_CM,c_U,MDT):
    # tau = vector = [tau_A_O,tau_B_O,tau_B_FT]
    lambda_E_B = lambdaE(tau[1],MTTF_B,alpha)
    q = lambda_E_B * tau[2]/2
    c_F = c_CM+q*MDT*c_U
    return c_PM/tau[0] +c_PM/tau[1] + c_FT/tau[2]+ \
        lambdaE(tau[0],MTTF_A,alpha)*c_F + \
            lambda_E_B * c_CM
print("\nc) The effective failure rate is given by:")
print("\nC(tau_A_O, tau_B_O, tau_B_FT) = c_PM/tau_A_O + c_PM/tau_B_O + c_FT/tau_B_FT")
print("    + lambdaE(tau_A_O,)*c_F + lambda_E(tau_B_O)*c_CM")
print("\nwhere c_F = c_CM + q_B*c_U*MDT now uses the calculated value of q_B")


    

print("\nd) Results for scipy minimize:")
tau=[0.5*MTTF_A,0.5*MTTF_B, 0.1*MTTF_B] # Initial guess for tau-vector

res = minimize(C_2, tau, args = (MTTF_A,MTTF_B,alpha,c_PM,c_FT,c_CM,c_U,MDT))
print("tau_A_O =",round(res.x[0]))
print("tau_B_O =",round(res.x[1]))
print("tau_B_FT =",round(res.x[2]))
print("C(tau_A_O, tau_B_O, tau_B_FT)=",round(res.fun,2), "(per hour) =", round(res.fun*8760),"(per year)")    


    
# Calculate new q_B and cost of failure:
lambda_B = lambdaE(tau[1],MTTF_B,alpha)
q = lambda_B * tau[2]/2
c_F = c_CM+q*MDT*c_U    
print("\nNew values: q_B =",round(q,4),",c_F =",round(c_F))
tau= [3674,tau_B_O,tau_B_FT] 
print("Cost given constraints on q_B=10%, =",round(C_2(tau,MTTF_A,MTTF_B,alpha,c_PM,c_FT,c_CM,c_U,MDT),2))
print("\ni.e., the constraints represent a sub-optimal solution")