「信頼性工学のはなし―信頼度99.9999…%をめざして」1～2章
大村平「信頼性工学のはなし―信頼度99.9999…%をめざして」日科技連 の読後行間補充メモ
同書籍は、印刷数表による信頼性工学の啓蒙書。
本稿では、同書籍の設例について、Python ［Google Colaboratory］による演算処理例を示す。
１章 信頼性を支える３本の柱
p5 双発と４発、どちらを信頼するか。
efp = float(input('エンジンの故障確率',))
print(f'双発機の墜落確率：{efp**2}')
print(f'４発機の墜落確率：{efp**4}')
print(f'双発機が遅れる[=エンジンのいずれか又は全部が故障する]確率：{1-(1-efp)**2}')
print(f'４発機が遅れる[=エンジンのいずれか又は全部が故障する]確率：{1-(1-efp)**4}')
エンジンの故障確率0.001 双発機の墜落確率：1e-06 ４発機の墜落確率：1.0000000000000002e-12 双発機が遅れる[=エンジンのいずれか又は全部が故障する]確率：0.001998999999999973 ４発機が遅れる[=エンジンのいずれか又は全部が故障する]確率：0.003994003998999962
- \(1e-06\) とあるのは、\(1\times 10^{-6}\) すなわち、0.000001 を指す。
- \(1.00.......002e-12\) とあるのは、\(1.00.......002\times 10^{-12}\) すなわち、0.00000000000100.......002 を指す。
p9 ２進法で表す 表1.1
import pandas as pd decimal = [] binary = [] for i in range(10): decimal.append(i) # 10進数 bin = format(i,'04b') #４桁の２進数 binary.append(bin) del i # 変数を廃棄 DF = pd.DataFrame([decimal,binary]).T DF.columns = ['10進法','２進法'] DF # 表1.1
10進法 ２進法 0 0 0000 1 1 0001 2 2 0010 3 3 0011 4 4 0100 5 5 0101 6 6 0110 7 7 0111 8 8 1000 9 9 1001
p10 パリティ・ビットを付ける 表1.2
（処理工程）
- ２進法の数値における１の個数をカウントし、奇数であれば０を末尾に加える
- ２進法の数値における１の個数をカウントし、偶数であれば１を末尾に加える
- 上記処理結果を、PB付き列として生成して、一覧表化する
# コードは上記の続き
import re
def parity_odd(string):
new_str = ""
string = re.sub('[^0-1]', '', string)
if string.count("1") % 2 == 1: #１の個数を２で除した余りが１[奇数]か否か[偶数]
new_str = string +"0" # 奇数の場合、０を末尾に
else:
new_str = string +"1" # 偶数の場合、１を末尾に
return new_str
# 新列（PB付）を加える
PB_col = []
for n in range(DF.shape[0]): # 行数分、繰り返す
PB_col.append(parity_odd(DF['２進法'][n]))
del n # 変数を廃棄
DF['PB付き'] = PB_col
DF
10進法 ２進法 PB付き 0 0 0000 00001 1 1 0001 00010 2 2 0010 00100 3 3 0011 00111 4 4 0100 01000 5 5 0101 01011 6 6 0110 01101 7 7 0111 01110 8 8 1000 10000 9 9 1001 10011
p11 ５桁のうちのどれか１つを間違える確率
import math
# 関数 ｎ個のうちｒ個のみで間違いが発生する確率
def error_p(n,r,p):
return math.comb(n,r)*(p**r)*(1-p)**(n-r)
n = int(input('全体n 個：',))
r = int(input('間違い r個：',))
p = float(input('0と1を間違える確率：',))
print(f'{n}個のうち{r}個だけ間違える確率：{round(error_p(n,r,p),7)}')
全体n 個：5 間違い r個：1 0と1を間違える確率：0.001 5個のうち1個だけ間違える確率：0.00498
p11 ５桁のうちのどれか２つを間違える確率
コードは上記と同じ。
全体n 個：5 間違い r個：2 0と1を間違える確率：0.001 5個のうち2個だけ間違える確率：1e-05
２章 バスタブが語る故障のパターン
p26 信頼性と故障の続柄
- 信頼度\(R\)
- 故障確率\(F\)
\(R=1-F\)
p27 人間の死亡率をグラフに描く 表2.1
（処理工程）
- 生存者データを設定
- 生存者データの差分 →死亡者数
- 死亡者数を加算 →死亡者累計数
- 死亡者数／生存者 →死亡率
- 一覧表化
# 表２．１ import pandas as pd import numpy as np # 年齢の列 age_col = [] ages = np.array([0,10]) for i in range(11): age_col.append(ages+10*i) # 生存者データ live_col = [100,95,94,93,91,85,69,42,17,4,0] DFal = pd.DataFrame([age_col,live_col]).T DFal.columns = ['Ages', 'Living'] # 死亡者数 Dead = [] for i in range(10): Dead.append(DFal.loc[:,'Living'][i] -DFal.loc[:,'Living'][i+1]) Dead_column = pd.Series(Dead) DFald = pd.concat([DFal, Dead_column],axis=1) DFald.columns = ['Ages', 'Living','Dead'] # 死亡者累計数 Dead_Accum = [] for j in range(10): Dead_Accum.append(DFald.loc[:,'Living'][0] -DFald.loc[:,'Living'][j+1]) DA_column = pd.Series(Dead_Accum) DFalda = pd.concat([DFald, DA_column],axis=1) DFalda.columns = ['Ages', 'Living','Dead','Dead_Accum'] # 死亡率％ Dead_percent = [] for k in range(10): Dead_percent.append(DFalda.loc[:,'Dead'][k] /DFalda.loc[:,'Living'][k]) DP_column = pd.Series(Dead_percent) DF = pd.concat([DFalda, round(DP_column*100,1)],axis=1) DF.columns = ['Ages', 'Living','Dead','Dead_Accum', 'Dead_Percent'] DF
Ages Living Dead Dead_Accum Dead_Percent 0 [0, 10] 100 5.0 5.0 5.0 1 [10, 20] 95 1.0 6.0 1.1 2 [20, 30] 94 1.0 7.0 1.1 3 [30, 40] 93 2.0 9.0 2.2 4 [40, 50] 91 6.0 15.0 6.6 5 [50, 60] 85 16.0 31.0 18.8 6 [60, 70] 69 27.0 58.0 39.1 7 [70, 80] 42 25.0 83.0 59.5 8 [80, 90] 17 13.0 96.0 76.5 9 [90, 100] 4 4.0 100.0 100.0 10 [100, 110] 0 NaN NaN NaN
p30 図21 死亡者数のグラフ
# コードは上記の続き
import matplotlib.pyplot as plt
# 年齢をリスト化
Ages_list = []
for m in DF['Ages']:
Ages_list.append(m[0])
# グラフ描画
plt.bar(Ages_list, DF['Dead'],width=9,align='edge') #alignをedgeとするとx座標は棒の左端に変更
plt.title('Dead')
plt.xlabel('Ages')
p30 図21 死亡者累計のグラフ
# コードは上記の続き
plt.bar(Ages_list, DF['Dead_Accum'],width=9,align='edge', color='orange')
plt.title('Dead_Accumulated')
plt.xlabel('Ages')
p30 図2.1 死亡率のグラフ
# コードは上記の続き
plt.bar(Ages_list, DF['Dead_Percent'],width=9,align='edge', color='lightgreen')
plt.title('Dead_Percent')
plt.grid()
plt.xlabel('Ages')
p31 図2.2 人間の死亡率
- 厚生労働省「第23回生命表（完全生命表）の概況」 …グラフあり
- 死亡率： 幼年時は高い。５～１０歳で極小化。その後、漸増する。
p32～37 図2.3、図2.4 故障率\(k\) のバスタブ曲線
ワイブル故障分布（p151）
\begin{eqnarray*} f(t)&=&k(t)\times R(t)\\&=&\cfrac{m}{\eta}\left(\cfrac{t-\gamma}{\eta}\right)^{m-1} \times e^{-\left( \cfrac{t-\gamma}{\eta}\right)^m}\end{eqnarray*}
を構成する故障率関数
\(k(t)=\cfrac{m}{\eta} \left(\cfrac{t-\gamma}{\eta}\right)^{m-1}\)
を使用する。
式中の各パラメータの趣旨は、以下のとおり。
- \(m\) …形状パラメータ
- \(\eta\) …尺度パラメータ
- \(\gamma\) …位置パラメータ
（処理工程）
- 形状パラメータを期間別に３つ用意する（初期、偶発、摩耗期間）
- 尺度パラメータ、位置パラメータを設定
- パラメータ値を故障率関数に入力して、各故障率を演算する
- ３つの演算結果を合計 →観測故障率
- グラフ化
# 故障確率k
import numpy as np
import math
from matplotlib import pyplot as plt
import japanize_matplotlib
# 設定値
m_list = [0.5, 1, 10] # 形状パラメータ（初期、偶発、摩耗）
eta_value = 1 # 尺度パラメータ(半減期、時定数)
gamma_value = 0 # 位置パラメータ
# 各時期（初期、偶発、摩耗）のパラメータ設定
g_wide =1 # グラフ横幅
step = 0.0001 # グラフ描画刻み
p1 = [0, g_wide, m_list[0], eta_value, gamma_value] # 始点、終点、形状、尺度、位置
p2 = [0, g_wide, m_list[1], eta_value, gamma_value]
p3 = [0, g_wide, m_list[2], eta_value, gamma_value]
graph_list = [p1,p2,p3] # リスト化
# 関数 故障確率 k (ワイブル)
def k_weibull(m_value,eta_value,gamma_value,time): # 形状、尺度、位置、時間
khead = m_value/eta_value
ktail = (time - gamma_value)/eta_value
kpower = m_value -1
y_value = khead * ktail ** kpower
return y_value
# 関数 グラフ用データ作成
def xylist(start,end,m_value,eta_value,gamma_value): # 始点、終点、形状、尺度、位置
x_list = []
y_list = []
for time in np.arange(start,end,step):
x_list.append(time) # 時間リスト
y_list.append(k_weibull(m_value,eta_value,gamma_value,time)) # 演算結果リスト
return x_list, y_list
# 各グラフのパラメータを関数（グラフ用データ）に渡して、描画
ky_lists = []
for graph in graph_list:
x_data = xylist(graph[0],graph[1],graph[2],graph[3],graph[4])[0]
y_data = xylist(graph[0],graph[1],graph[2],graph[3],graph[4])[1]
plt.plot(x_data, y_data, linestyle="--")
ky_lists.append(y_data)
# バスタブ曲線
ky_lists[2] = np.nan_to_num(ky_lists[2], nan=0) # nan を0に置換
ky_sums = [x+y+z for (x,y,z) in zip(ky_lists[0], ky_lists[1],ky_lists[2])] # ３データを合計
plt.plot(x_data, ky_sums, label="観測故障率", color='red', linewidth=3)
# グラフ凡例等
plt.xlim(0,g_wide)
plt.ylim(0,15)
plt.yticks([])
plt.xlabel('time')
plt.ylabel('k [/t]')
plt.legend()
plt.title('バスタブ曲線')
plt.show()
