################################################################### # # # This program plots I_B^{(i)}(t), i = 1, ..., 5 as a function # # of t for a threshold system using a ROBDD representation of # # the system. # # # # Contrary to script_B_3_7 which uses Monte Carlo simulation, # # this script calculates importance analytically. Thus, the # # results are exact I_B^{(i)}(t)-curves. # # # ################################################################### from math import exp import matplotlib.pyplot as plt import numpy as np from c_robdd import ROBDDSystem from c_threshold import ROBDDThreshold save_plots = True # Time parameters max_t = 120 steps = 120 # Weibull-parameters for the components alpha = [2.0, 2.5, 3.0, 1.0, 1.5] beta = [50.0, 60.0, 70.0, 40.0, 50.0] # Number of components in the system num_comps = len(alpha) # Threshold system weights and threshold a = [2.0, 4.0, 2.3, 3.5, 1.9] b = 5.0 # Calculate the weight sum asum = 0 for i in range(num_comps): asum += a[i] # The ROBDD representation of the threshold system sys = ROBDDSystem(ROBDDThreshold(a, asum, b)) # The reliability function of the system def reliability(pp): result = sys.calculateReliability(pp) return result[1] # The survival function of a Weibull-distribution def weibull(aa, bb, t): if t >= 0: return exp(- (t / bb) ** aa) else: return 1.0 # Component reliability vector p = np.zeros(num_comps) # Time values and Birnbaum measures time = np.linspace(0.0, max_t, steps) IB = np.zeros((num_comps, steps)) # Calculate the Birnbaum measures for step in range(steps): for i in range(num_comps): p[i] = weibull(alpha[i], beta[i], time[step]) for i in range(num_comps): temp = p[i] p[i] = 1 h1 = reliability(p) p[i] = 0 h0 = reliability(p) p[i] = temp IB[i][step] = h1 - h0 # Plot the results fig = plt.figure(figsize = (7, 4)) for i in range(num_comps): plt.plot(time, IB[i], label='IB(' + str(i+1) + ", t)") plt.xlabel('Time') plt.ylabel('Importance') plt.title("Importance as a function of time") plt.legend() if save_plots: plt.savefig("birnbaum_03b/birnbaum_03b.pdf") plt.show()