Python做Reliability分析，weibull模型为例，计算单侧/双侧可靠性值（含置信度）

主要用到Matthew Reid开发的Reliability包

第一步：导入或创建F和S数据，本案例中假设有4种失效类型，同时有F和S的数据。

import numpy as np
right_censored_2 = np.repeat(5130,3000) #有3000个样本在5130个循环时发生右删失
print(right_censored)

#假设有4种失效类型
failures1=[500,600,xxx] #输出失效样本的失效时间
failures2=[500,600,xxx] #输出失效样本的失效时间
failures3=[500,600,xxx] #输出失效样本的失效时间
failures4=[500,600,xxx] #输出失效样本的失效时间

failures = np.hstack([failures4,failures3,failures2,failures1])
print(failures)

第二步：拟合2p-weibull分布（若样本的失效服从其它分布，则需要选择其它统计分布）并绘图

from reliability.Fitters import Fit_Weibull_2P
import matplotlib.pyplot as plt
from reliability.Other_functions import crosshairs
fit=Fit_Weibull_2P(failures=failures, right_censored=right_censored, show_probability_plot=True,CI=0.90,CI_type="reliability",
                  quantiles=[0.01,0.03,0.05,0.07,0.09,0.1])
plt.xlabel('Cycle')

生成拟合结果效果图

如果需要查看SF曲线

#SF曲线
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator
fit.distribution.SF(label='Fitted Distribution',color='steelblue')
plt.xlabel('Cycle')
plt.legend()
plt.grid()
plt.xlim(0,8000)
plt.ylim(0.96,1.01)
x_major_locator=MultipleLocator(1000)
y_major_locator=MultipleLocator(0.02)
plt.axvline(x=5000,linestyle="--",color="r")
plt.show()

第三步：计算可靠性（区分双侧和单侧）

置信度C50的情况，以下代码计算出在50%置信度下，5000个循环的可靠性为97.76%

dist=fit.distribution
sf_5000 = dist.SF(5000)  #计算置信度50%时，cycle5000的可靠性
print('The value of the SF at cycle 5000 with 50% confidence is', round(sf_5000 * 100, 2), '%') 

双侧和单侧

Theory: R is reliability, T is time/cycle
u estimat =ln(−ln(R))
u estimate =β.(ln(T)−ln(α))

Theory: formula of var(u estimate)

#u_estimate
u_estimate=np.log(-np.log(0.9776)) #填入C50的可靠性值
print('u estimate=',u_estimate)

#var(u_estimate)
one=(np.log(5000)-np.log(fit3.alpha))**2 * fit3.beta_SE ** 2
two=(-fit3.beta/fit3.alpha)**2 * fit3.alpha_SE ** 2
three=2*(np.log(5000)-np.log(fit3.alpha))*(-fit3.beta/fit3.alpha)*fit3.Cov_alpha_beta
var_v=one+two+three
print('var(u estimate)=',var_v)

计算结果：

u estimate= -3.7873883595391376

var(u estimate)= 0.02000064345713326

Theory: the calculation of u

Theory: the calculation of r

双侧C90的计算代码：

#u-two-sided-bound
z=1.645 #置信度90%的双侧
u_lower=u_estimate- z * (var_v**0.5)
u_upper=u_estimate+ z * (var_v**0.5)
print('u lower=',u_lower)
print('u upper=',u_upper)

#reliability-two-sided
r_lower=np.exp(-np.exp(u_lower))
r_upper=np.exp(-np.exp(u_upper))
print('The two-sided reliability upper value with 90% confidence at 5000 cycle is',round(r_lower,4))
print('The two-sided reliability lower value with 90% confidence at 5000 cycle is',round(r_upper,4))

计算结果：

u lower= -4.020030232836033
u upper= -3.5547464862422413
The two-sided reliability upper value with 90% confidence at 5000 cycle is 0.9822
The two-sided reliability lower value with 90% confidence at 5000 cycle is 0.9718

单侧C90的计算代码：

#reliability-one-sided-lower
z=1.28 #C90单侧
u_upper=u_estimate+ z * (var_v**0.5)
r_upper=np.exp(-np.exp(u_upper))
print('u upper=',u_upper)
print('The one-sided reliability lower value with 90% confidence at 5000 cycle is',round(r_upper,4))

计算结果：

u upper= -3.6063661116242276
The one-sided reliability lower value with 90% confidence at 5000 cycle is 0.9732