##################################################################### # # # This program calculates the reliability of a threshold system # # by finding the distribution of: # # # # S = (a_{1} X_{1} + ... + a_{n} X_{n}). # # # # # ##################################################################### import numpy as np # Component reliabilities (p[0] is not in use) p = [0.0, 0.75, 0.80, 0.82, 0.65, 0.88, 0.91, 0.92, 0.86] # Component weights (a[0] is not in use) a = [0, 26, 39, 65, 78, 91, 104, 104, 143] # Threshold b = 410 # Number of components n = len(p) - 1 # Cumulative weights (d[0] is not in use) d = np.zeros(n+1, dtype=int) for j in range(1, n + 1): d[j] = a[j] + d[j-1] # After each iteration ps[s] = a_{1} X_{1} + ... + a_{n} X_{n} = s) ps = np.zeros(d[n] + 1) # Temporary storage of ps[s] ph = np.zeros(d[n] + 1) # Calculate ps[0] = P(X_{1} = 0) and ps[1] = P(X_{1} = 1): ps[0] = 1 - p[1] ps[a[1]] = p[1] # Save these values for the next iteration ph[0] = ps[0] ph[a[1]] = ps[a[1]] for j in range(1, n): # Calculate P(a_{1} X_{1} + ... + a_{j} X_{j+1} = s), s = 0, 1, ... , (a[j+1]-1) for s in range(a[j+1]): ps[s] = ph[s] * (1 - p[j+1]) # Calculate P(a_{1} X_{1} + ... + a_{j} X_{j+1} = s), s = a[j+1], ..., d[j]: for s in range(a[j+1], d[j] + 1): ps[s] = ph[s] * (1 - p[j+1]) + ph[s-a[j+1]] * p[j+1] # Calculate P(a_{1} X_{1} + ... + a_{j} X_{j+1} = s), s = d[j]+1 , ... , d[j+1]: for s in range(d[j] + 1, d[j+1] + 1): ps[s] = ph[s-a[j+1]] * p[j+1] # Save these values for the next iteration for s in range(d[j+1] + 1): ph[s] = ps[s] # At this stage ps[s] = P(a_{1} X_{1} + ... + a_{j} X_{j+1} = s), s = 0, 1, ... , d[n] for s in range(d[n] + 1): print("P(S = ", s, ") = ", ps[s]) # Calculate the system reliability h = P(S = b) + ... + P(S = d[n]) h = 0 # for s in range(b, d[n] + 1): # h += ps[s] for s in range(d[n] + 1): if s >= b: h += ps[s] print("h = ", h)