「信頼性工学のはなし―信頼度99.9999…%をめざして」3章
大村平「信頼性工学のはなし―信頼度99.9999…%をめざして」日科技連 の読後行間補充メモ
同書籍は、印刷数表による信頼性工学の啓蒙書。
本稿では、同書籍の設例について、Python ［Google Colaboratory］による演算処理例を示す。
３章 偶発故障と信頼性
p48 表3.1 １時間に10％ずつ減ると
（処理工程）
- 経過時間の列を設定
- 各経過時間に対して残量を計算
- 一覧表化
import pandas as pd import numpy as np # 経過時間 pass_time = np.array(range(0,20,1)) # 残グラム grams = [] for time in pass_time: grams.append(100 * 0.9 ** time) # DF化 time_s = pd.Series(pass_time) grams_s = round(pd.Series(grams),2) DF = pd.DataFrame([time_s,grams_s]).T DF.columns = ['Passed_Time','Residual_Grams'] DF
Passed_Time Residual_Grams 0 0.0 100.00 1 1.0 90.00 2 2.0 81.00 3 3.0 72.90 4 4.0 65.61 5 5.0 59.05 6 6.0 53.14 7 7.0 47.83 8 8.0 43.05 9 9.0 38.74 10 10.0 34.87 11 11.0 31.38 12 12.0 28.24 13 13.0 25.42 14 14.0 22.88 15 15.0 20.59 16 16.0 18.53 17 17.0 16.68 18 18.0 15.01 19 19.0 13.51
p49 図3.1 １時間に10％ずつ減ると （グラフ）
# コードは上記の続き
import matplotlib.pyplot as plt
plt.plot(DF['Passed_Time'],DF['Residual_Grams'], label= r'$y=100 e^{-t/10}$')
# 目盛値
plt.xticks(np.arange(0, 19, step=1))
plt.yticks(np.arange(0, 101, step=10))
# 軸名
plt.xlabel(DF.columns[0])
plt.ylabel(DF.columns[1])
# 表示設定
plt.xlim(0,) # 表示領域
plt.legend() # 凡例
plt.grid() # 目盛線
plt.show()
p50 図3.2 故障率 k 毎のグラフ
（処理工程）
- 複数の故障率\(k\) を設定
- 経過時間\(t\) の列（グラフ横軸情報）を設定
- 故障率\(k\) と経過時間\(t\) から残存率\(R\) を求める関数を設定
- 関数で演算した結果を、グラフ化
import math
import numpy as np
import matplotlib.pyplot as plt
# 関数（時間と残存率）
def resudial_func(k_value):
time_data = [] # 時間_x軸データ用
residual_data = [] # 残存率_y軸データ用
for time in np.arange(0,19, 0.1): # 時間データ
time_data.append(time)
residual = math.e**(-k_value * time)
residual_data.append(residual)
return time_data,residual_data
# 複数グラフ
for k_value in [0.05, 0.1, 0.2, 0.3, 0.5]: # 故障率k のリスト
time_data = resudial_func(k_value)[0] # 関数より時間列を取得
residual_data = resudial_func(k_value)[1] # 関数より残存率列を取得
plt.plot(time_data, residual_data, label= f'k ={k_value}')
# 軸名
plt.xlabel('Time Passed [t]')
plt.ylabel('Residual Rate [R]')
# 目盛値
plt.xticks(np.arange(0, 19, step=1))
plt.yticks(np.arange(0, 1.1, step=0.2))
plt.legend()
plt.xlim(0,)
plt.ylim(0,)
plt.title(r'$R=e^{-kt}$')
plt.plot()
p52 図3.3 減衰曲線
複数の故障率\(k\) を一本化したグラフ
- 横軸（故障率と経過時間の積） \(kt\)
- 縦軸（残存量） \(y\)
import math
import numpy as np
import matplotlib.pyplot as plt
# 設定値
A_value = 5 # 係数A(描画データ作成用)
# 関数（時間と残存数）
def resudial_func():
kt_data = [] # kt(x軸データ用)
residual_data = [] # 残存数(y軸データ用)
for kt in np.arange(0,6, 0.1): # ktデータ
kt_data.append(kt)
residual = A_value * math.e**(-kt)
residual_data.append(residual)
return kt_data,residual_data
# グラフ用データ作成
kt_data = resudial_func()[0] # 関数より kt列を取得
residual_data = resudial_func()[1] # 関数より y列を取得
plt.plot(kt_data, residual_data, label= r'$ y=Ae^{-kt}$')
plt.xlabel(f'failure rate × passed time 'r'$ [kt] $')
plt.yticks([0,1,2,3,4,5], [0,'0.2A','0.4A','0.6A','0.8A','A']) # 目盛値の設定
plt.ylabel('Residual Amount [y]')
plt.legend()
plt.xlim(0,5)
plt.ylim(0,)
plt.title(f'Residual Curve')
plt.plot()
p53 極限値の計算
\(\lim\limits_{n\to\infty}\left(1-\cfrac{0.1}{n}\right)^n\)
import sympy
infn = float('inf') # 無限大
sympy.var('x') # 文字数ｘ
fx = (1-0.1/x)**x # 関数
sympy.limit(fx, x, infn) # 極限値
p54 表3.2 e の -x 乗 の値
import math import pandas as pd x_list = [0, 0.01, 0.02, 0.03, 0.05, 0.07, 0.1, 0.2, 0.3, 0.5, 0.7, 1, 1.5, 2,3,4,5,7,10] y_list_b = [math.e ** -x_value for x_value in x_list] # exp(-x)の計算 y_list = [round(y_value,5) for y_value in y_list_b] # 小数点以下５桁に変換 DF = pd.DataFrame([x_list,y_list]).T DF.columns = ['x','exp(-x)'] DF
x exp(-x)
0 0.00 1.00000
1 0.01 0.99005
2 0.02 0.98020
3 0.03 0.97045
4 0.05 0.95123
5 0.07 0.93239
6 0.10 0.90484
7 0.20 0.81873
8 0.30 0.74082
9 0.50 0.60653
10 0.70 0.49659
11 1.00 0.36788
12 1.50 0.22313
13 2.00 0.13534
14 3.00 0.04979
15 4.00 0.01832
16 5.00 0.00674
17 7.00 0.00091
18 10.00 0.00005
p54 その１ 水槽の中の小魚
経過時間\(t\) と犠牲率\(k\) → 残存率\(R\)
import math
k_value = 0.01 # １分あたりの犠牲率 [ /分]
time = 10 # 経過時間 [分]
residual = math.e ** -(k_value * time)
print(f'残存率：{round(residual*100,1)}％')
残存率：90.5％
p55 その２ 水槽の中の小魚
残存率\(R\) と犠牲率\(k\) →経過時間\(t\)
import math
k_value = 0.01 # １分あたりの犠牲率 [ /分]
e_value = 0.5 # 残存率
time_passed = -math.log(e_value)/k_value
print(f'経過時間：{round(time_passed,1)}分')
経過時間：69.3分
p56 その３前半 デバッキングが終わった製品（偶発故障期間にある製品）
総数、経過時間\(t\)と故障数 →故障率\(k\)
import math
# 設定情報
item_num = 1350 # 製品数
f_num = 40 # 故障数
time_passed = 20 # 経過日数
# 各種数値の演算
residual = item_num - f_num # 残存数
e_value = residual/item_num # 残存率
k_value = -math.log(e_value)/time_passed # 故障率 [/day]
print(f'故障率k：{round(k_value,4)}[/day]')
故障率k：0.0015[/day]
p57 その３後半 偶発故障期間にある製品の信頼度
予定経過時間\(t\)と故障率\(k\) →信頼度（残存率）\(R\)
# コードは上記の続き
time_passed_2 = 200 # 使用予定日数
e_value_2 = math.e ** (-k_value * time_passed_2) # 残存率
print(f'信頼度（残存率）R：{round(e_value_2*100,2)}[％]')
信頼度（残存率）R：74.02[％]
p60 平均寿命（MTBF）
平均寿命の求め方２種
- 減衰曲線\(Ae^{-kt}\) 下の面積（寿命の総和）を総個数で除する
- 故障率\(k\) の逆数
故障率\(k\) が同じ場合、総個数によらず、平均寿命（MTBF）は、\(\cfrac{1}{k}\) で決まる。
import sympy
import math
All = int(input('総個数：',))
k_value = 0.0015 # 故障率 [/day]
time = sympy.Symbol('time')
inf = float('inf')
formula = All * math.e ** (-k_value * time) # 減衰曲線の関数
sum_days = sympy.integrate(formula, (time,0,inf)) # 定積分（寿命の総和）
average = sum_days/All # 平均寿命
print(f'寿命の総和：{int(sum_days)}[day]')
print(f'製品の平均寿命 by 寿命の総和／総個数：{round(average,1)}[day]')
print(f'製品の平均寿命 by 1／故障率：{round(1/k_value,1)}[day]')
総個数：1350 寿命の総和：900000[day] 製品の平均寿命(by 寿命の総和／総個数)： 666.7[day] 製品の平均寿命(by 1／故障率)： 666.7[day]
総個数：555555 寿命の総和：370370000[day] 製品の平均寿命(by 寿命の総和／総個数)： 666.7[day] 製品の平均寿命(by 1／故障率)： 666.7[day]
p65 図3.5
# 図3.5
import math
import numpy as np
from matplotlib import pyplot as plt
# 設定値
A_value = 5 # 係数A(描画データ作成用)
# 関数（時間と残存数）
def resudial_func():
kt_data = [] # kt(x軸データ用)
residual_data = [] # 残存数(y軸データ用)
for kt in np.arange(0,6, 0.1): # ktデータ
kt_data.append(kt)
residual = A_value * math.e**(-kt)
residual_data.append(residual)
return kt_data,residual_data
# グラフ用データ作成
kt_data = resudial_func()[0] # 関数より kt列を取得
residual_data = resudial_func()[1] # 関数より y列を取得
# MTBFとAとの領域の塗りつぶし
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(111)
ax.set_xlim(0, 4)
ax.set_ylim(0, 6)
ax.axvspan(0, 1, 0, 5/6, color="lightgrey")
# 描画
plt.plot(kt_data, residual_data, label= r'$ y=Ae^{-kt}$')
plt.xlabel(f'failure rate × passed time 'r'$ [kt] $')
# 目盛値の設定
plt.xticks([0,1,2,3,4,5], [0,'k\n×\nMTBF',2,3,4,5])
plt.yticks([0,1,2,3,4,5], [0,'0.2A','0.4A','0.6A','0.8A','A'])
# 表示設定
plt.ylabel('Residual Amount [y]')
plt.legend()
plt.xlim(0,5)
plt.ylim(0,)
plt.title(f'Residual Curve')
plt.plot()
p64 半減期
import math
e_value = 0.5 # 残存率
half_life = -math.log(e_value)
print(f'半減期ｔ= {round(half_life,2)}/k')
半減期ｔ= 0.69/k
p64 図3.6＝p68 図3.8 半減期とMTBF
# 図3.6
import math
import numpy as np
from matplotlib import pyplot as plt
# 設定値
A_value = 5 # 係数A(描画データ作成用)
# 関数（時間と残存数）
def resudial_func():
kt_data = [] # kt(x軸データ用)
residual_data = [] # 残存数(y軸データ用)
for kt in np.arange(0,6, 0.1): # ktデータ
kt_data.append(kt)
residual = A_value * math.e**(-kt)
residual_data.append(residual)
return kt_data,residual_data
# 曲線グラフ用データ作成
kt_data = resudial_func()[0] # 関数より kt列を取得
residual_data = resudial_func()[1] # 関数より y列を取得
# 曲線描画
plt.plot(kt_data, residual_data, label= r'$ y=Ae^{-kt}$')
plt.xlabel(f'failure rate × passed time 'r'$ [kt] $')
# 補助線描画
x_HL = half_life
y_MTBF = A_value * math.e**(-1)
plt.hlines(A_value*0.5,0,x_HL, color = 'orange',linestyle="--") # 半減期横
plt.vlines(x_HL,0,A_value*0.5, color = 'orange',linestyle="--") # 半減期縦
plt.hlines(y_MTBF,0,1, color = 'grey',linestyle="--") # MTBF横
plt.vlines(1,0,y_MTBF, color = 'grey',linestyle="--") # MTBF縦
# 目盛値の設定
plt.xticks([0,x_HL,1,2,3,4,5], [0,'k\n×\nHL','k\n×\nMTBF',2,3,4,5])
plt.yticks([0, y_MTBF, A_value*0.5,A_value], [0,'0.37A','0.5A','A'])
# 表示設定
plt.ylabel('Residual Amount [y]')
plt.legend()
plt.xlim(0,5)
plt.ylim(0,)
plt.title(f'Residual Curve')
plt.plot()
p65 図3.7 時定数
\begin{eqnarray*}Average \, Life &=& MTBF\\ &=& Time\, Constant \end{eqnarray*}
import math
import numpy as np
from matplotlib import pyplot as plt
# 設定値
A_value = 5 # 係数A(描画データ作成用)
# 関数（時間と残存数）
def resudial_func():
kt_data = [] # kt(x軸データ用)
residual_data = [] # 残存数(y軸データ用)
for kt in np.arange(0,6, 0.1): # ktデータ
kt_data.append(kt)
residual = A_value * math.e**(-kt)
residual_data.append(residual)
return kt_data,residual_data
# 曲線グラフ用データ作成
kt_data = resudial_func()[0] # 関数より kt列を取得
residual_data = resudial_func()[1] # 関数より y列を取得
# 曲線描画
plt.plot(kt_data, residual_data, label= r'$ y=Ae^{-kt}$')
plt.xlabel(f'failure rate × passed time 'r'$ [kt] $')
# MTBF線描画
y_MTBF = A_value * math.e**(-1)
plt.hlines(y_MTBF,0,1, color = 'grey',linestyle="--") # MTBF横
plt.vlines(1,0,y_MTBF, color = 'grey',linestyle="--") # MTBF縦
# 時定数線 = 0,A接線 描画
plt.plot([0, 1], [A_value, 0], color = 'green', linestyle="--")
# 目盛値の設定
plt.xticks([0,1,2,3,4,5], [0,'k\n×\nMTBF\n(Time Const.)',2,3,4,5])
plt.yticks([0, y_MTBF, A_value], [0,'0.37A','A'])
# 表示設定
plt.ylabel('Residual Amount [y]')
plt.legend()
plt.xlim(0,5)
plt.ylim(0,)
plt.title(f'Residual Curve')
plt.plot()
p71 図3.10 故障率 k が一定の場合（偶発故障期間）における諸グラフ
import numpy as np
import math
from matplotlib import pyplot as plt
t_line = np.arange(0,10,0.1) # 時間
A_value = 5 # 係数
k_value = 0.6 # 故障率［ ／時間］
y1_line = [] # グラフ① ｙ軸データ用
y2_line = [] # グラフ② ｙ軸データ用
y3_line = [] # グラフ③ ｙ軸データ用
y4_line = [] # グラフ④ ｙ軸データ用
y5_line = [] # グラフ⑤ ｙ軸データ用
y6_line = [] # グラフ⑥ ｙ軸データ用
for t in t_line:
y1_value = A_value * math.e**(-k_value*t) # ①の関数
y1_line.append(y1_value)
y2_value = k_value * y1_value # ②の関数
y2_line.append(y2_value)
y3_value = y2_value/A_value # ③の関数
y3_line.append(y3_value)
y5_value = A_value - y1_value # ⑤の関数
y5_line.append(y5_value)
y4_value = y5_value/A_value # ④の関数
y4_line.append(y4_value)
y6_value = y1_value/A_value # ⑥の関数
y6_line.append(y6_value)
# グラフの描画
plt.plot(t_line, y1_line,label=r'$①\quad y=A e^{-kt} $')
plt.plot(t_line, y2_line,label=r'$②\quad y=kA e^{-kt} $')
plt.plot(t_line, y3_line,label=r'$③\quad f=k e^{-kt} $')
plt.plot(t_line, y4_line,label=r'$④\quad F=1 -e^{-kt} $',linestyle="--")
plt.plot(t_line, y5_line,label=r'$⑤\quad y=A -Ae^{-kt} $',linestyle="--")
plt.plot(t_line, y6_line,label=r'$⑥\quad R= e^{-kt} $')
# タイトル
td1 = '① Remaining amount'
td2 = '② Number of failures'
td3 = '③ Distribution of failures'
td4 = '④ Probability of failure'
td5 = '⑤ Cumulative amount of failures'
td6 = '⑥ Reliability'
plt.title(f'{td1}←→{td5}\n{td6}←→{td4}\n{td2}\n{td3}\n A={A_value}, k={k_value}/time', fontsize = '11')
# 表示設定
plt.xlim(0,)
plt.ylim(0,)
plt.xlabel(r'$time$')
plt.ylabel(r'$y,f,F,R$')
plt.legend(loc ='right', fontsize = '12')
plt.grid()
plt.show()
p71 グラフ ① ～ ⑥ の相互関係
残存系①⑥ vs 故障系②③④⑤
- 総数 \(Ａ\) で除する： ①②⑤平面［個数］ \(\Longrightarrow\) ⑥③④平面［割合］
- 故障率 \(k\) を乗ずる： ①⑥軸 \(\Longrightarrow\) ②③軸［瞬間故障］
- 補数をとる： ①⑥軸 \(\Longrightarrow\) ⑤④軸［故障累計］
- 微分する： ⑤④軸 \(\Longrightarrow\) ②③軸
import sympy
import math
A = sympy.Symbol('A')
math.e = sympy.Symbol('e')
k = sympy.Symbol('k')
t = sympy.Symbol('t')
f1 = A*e**(-k*t)
f2 = f1*k
f3 = f2/A
f5 = A-f1
f4 = f5/A
f6 = f1/A
equation_list = [f1,f2,f3,f4,f5,f6]
number = 1
for equation in equation_list:
print(f'\n式{number}')
display(equation)
number += 1
