ベイズ統計を用いた故障分布関数の更新
公開日:
カテゴリ: 第2回
1.緒言
プラントの保全合理化の1つのテーマとして、機器 /部品類の取替周期の最適化がある[1][2][3]。こうした 最適化を図るためには、故障率や故障分布関数のデー タが不可欠である。一方、日本のプラントの場合、ト ラブルを恐れるあまり、寿命に達するかなり前に交換 されてしまうため、安全に関わる機器はともかく、一 般機器のデータまで十分整備されている訳ではない。 特に、故障分布関数の整備は遅れている。こうした状 況では、また、故障分布関数を推定するための一貫し たデータセットを収集する事は難しいと考えられ、社 内試験や異なるプラントあるいは更に異なる種類のプ ラント等におけるデータをベースに、少ない実プラン トでのデータと合わせて推定しなければならない。そ のための手法の1つとしてベイズ推定がある[4]。しかし、こうした推定には誤差が付きのもで、リス クを考慮するためには、この誤差も考慮することが必 要である。ベイズ推定では、事後分布によりこういっ た誤差を考慮することができる[4]。実プラントでのデ ータそのもので無くても、参考データがあり、ある程 度誤差が推測されているならば、それを事前分布とし て、実プラントでのデータから事後分布を求めること により、誤差をリスク評価に利用することができる。 ここでは、従来の最尤推定法[5]とベイズ推定法によ
L-11-4““ newp-clr““, expre-ry] のEL LTTり故障分布関数のパラメータを推定し比較検討した。 2. 故障分布関数のパラメータ決定 2.1 最尤推定法 * n個の機器あるいは部品を一定時間(以降、試験期間) 使用し、r個の機器/部品が時刻 1,,... , に故障した とする。一定期間(打切時間)を T とすると、例えばワ イブル分布の場合、尤度関数 Lは次式で表される。_L - 117-45““ - exper-day) | capt-ry] ()i=I LIIワイブルパラメータ m、T は、(1)式で与えられる尤度 関数 L が最大となるように決められる。(1)式の代わり に(1)式の対数 log L を考え、これを最大にする m、tは 次式で与えられる。““ + m_logt, r)-““ log(0, 1z) + (n - r)log(T / r)} = 0__ +mZlogc, r)rum)-2-{““ log(t, It) + (n - r)T““ log(T / r)} = 0[立x““ +(n-1). T)T＝=くこれらの(2)式および(3)式を連立して解くことにより、 ワイブルパラメータ m、t が求められる。399故障分布関数を用いて機器/部品等の取替周期の最 適化を図る場合、リスクを的確に評価するためには、 分布関数のパラメータの点推定値のみでなく区間推定 値も必要である。ワイブル分布の場合、2つのパラメ ータ m、tの両方とも未知の場合、両者の区間推定値 を求めるのは難しい。ここでは、m の点推定値が分か っているとして、T の区間推定を求めた場合の結果を 照会する[5]。 (4)式の統計量を考えると、2S, Ir““ が自 由度 2F の x 分布に従う。即ち、信頼係数1-a の両側 信頼区間は(5)式で表される。* s = +(n-17““( 12S28 (x (2ryu/2)““ x (2r;1- al2))片側区間も同様に(6)式で表される。-20* L1 ,, tr;m, )p(m, T) '(m,t) ==SSL(t1,.., tr; m, 7)(m, t)dmdr 2.2 ベイズ推定法 ・ ワイブル分布のパラメータm、Tは統計的に推定さ れており、その事前分布はゆ(m,)と書けるとする。 - 2.1 で述べた故障時刻データより、ワイブル分布の パラメータ m、t の事後分布の(m, r)はベイズ推定法 を用いると次のようになる。L(tys.,ltr; m, t)(m,t) 小(m,t) =““So So Lify, tr; n, 7)(m, )dmdr * ここで、L(1,, ,; m, )は(1)で与えられる尤度関数 である。ワイブル分布のパラメータ m、t が独立な場合、 小(m,n)は次のように表される。 ゆ(m, t) = pn(m). p. ()(8) この場合、(m) あるいは中.(r)は、ゆ(m, )より次 のように求めることが出来る。-9(m) = [ $(m,)dr spy(r) = [ p(m,t)dm但し、(m) および中.(r)はパラメータ m およびT に関する事前/事後分布である。てが独立な場合、-8不信賴度は、ゆ(m,)より次、0_50100-タ m およびベイズ推定の場合、こうして事後分布を求めること により、点推定値および区間推定値が、難しい数学的 取扱い無しに求める事が出来る。3. 故障分布関数のパラメータ決定試計算3.1 試計算に用いた故障分布関数ここでは、故障分布関数がワイブル分布の場合につ いて、最尤推定法とベイズ推定法による結果を比較す る。手順は、先ずパラメータ m、t を仮定し、その故 障分布関数(以降、基準関数)から、乱数により故障時刻 を求め、そのデータを基に最尤推定法とベイズ推定法 により m と を求め、それらの解の妥当性をみること とした。 今回の試計算の条件を Table 1 にまとめる。Table 1 試計算条件 Im (1,000h)|T(1,000h) |機器数n|故障数r | 2.5 150| 400 1000 -m = 2.5、T = 150,000hのワイブル分布の場合、機器 数 100台のうち、試験期間 T = 40,000h 内に故障する期 待数は約 3.6台となるため、ここでは、故障機器数を4 台とし、4台の機器の故障時刻を乱数により決定した た。 Fig.1 に m = 2.5、T = 150,000 h のワイブル分布に 4台の機器の故障データをプロットした図を示す。但し、機器の故障時刻から故障分布関数の値 F を算 出する際には、下式で示すメジアンランクを用いた。i-0.3 11n +0.4ここで、i は当該時刻までに故障した機器数を表し、 nは総機器数である。F%3D-10不信賴度10_50250300100 150 200 作動時間(1,000h)Fig.1 基準関数、 本図より、=150,000 hに対し、試験期間 T を 40,000 h とした場合、故障分布関数の極初期の段階に400おける故障しか見ていない事が分かる。こうした設定 は、機器の寿命の遙か以前に機器を交換してしまう日 本の原子力プラントの実情を考慮した結果である。 より分かりやすく、作動時間の範囲 T = 40,000 h に狭め て表示した図を Fig. 2 に示す。0.050.04。不信賴度。。。10_5_10 15 20 25 30 3540作動時間(1,000h) - Fig.2 故障データに基づく不信頼度(case 1) [ 本図に示した故障時刻データは、基準関数に良く合 致しているケースである。3.2 最尤推定法とベイズ推定法の比較 - 3.1に示したデータに対して、最尤推定法にて求め たm、T は、m = 4.0、T =90,000 h となった。m、t が この値の時の故障分布関数を Fig.3 に示す。0.050.040.03世VK* 0_5_101900/01/142025301900/02/0340IL EN AAT, anni!5_10152025303540には日(10001,不信賴度10_5_1030354015 20 25 作動時間(1,000 h)Fig.3 最尤推定法により得た故障分布関数(case 1) 作動時間が時間スケールに比べて極限られた範囲で あるため、m、t の推定値は、故障を生み出した本来 の故障分布関数とは異なったものとなっているが、デ ータの範囲では妥当と判断される結果を与えている。これに対し、同じ故障のデータについてベイズ推定 を行った結果を Fig. 4 に示す。なお、m、T の事前分布■30-35 ロ 25-30 ■20-25 015-20 ロ 10-15 ■5-10 00-5mは一様分布を仮定している。■30-35 025-30 ■20-25 015-20 010-15 ■ 5-10 00-595125)1551953215円mFig. 4(a) ベイズ推定による m、t 事後分布(case 1)11111111111HAH24311TTTTTTTTTB9876543210 5555555555HTTFig.4(6) ベイズ推定による m、事後分布(case 1) 事後分布の最大値は最尤推定法で得られた結果とほ ぼ同じ m = 4.0、T =90,000h で与えられている。しかし 事後分布を見ると、特定の方向ではあるが、かなり広 い分布を示しており、m=2.5、t=150,000h である可 能性も依然高い事が分かる。また、パラメータ m とて の間には強い相関がある事も理解される。こうした強 い相関の情報は、最尤推定法からは得る事が出来ない。 m = 2.5 として、(5)式より区間推定を行うと、95%信頼 区間として、105(x1,000h) < < 242(x1,000h) を得る。 また、m = 4.0 とすると、73(x1,000h)pp33-39 [4] 繁桝算男著、”ベイズ統計入門”、東京大学出版会 [5] 三根 久,河合 一著,信頼性・保全性の基礎数理,日科技連, (1984) pp.151-159.403“ “ベイズ統計を用いた故障分布関数の更新“ “笠井 雅夫,Masao KASAI
プラントの保全合理化の1つのテーマとして、機器 /部品類の取替周期の最適化がある[1][2][3]。こうした 最適化を図るためには、故障率や故障分布関数のデー タが不可欠である。一方、日本のプラントの場合、ト ラブルを恐れるあまり、寿命に達するかなり前に交換 されてしまうため、安全に関わる機器はともかく、一 般機器のデータまで十分整備されている訳ではない。 特に、故障分布関数の整備は遅れている。こうした状 況では、また、故障分布関数を推定するための一貫し たデータセットを収集する事は難しいと考えられ、社 内試験や異なるプラントあるいは更に異なる種類のプ ラント等におけるデータをベースに、少ない実プラン トでのデータと合わせて推定しなければならない。そ のための手法の1つとしてベイズ推定がある[4]。しかし、こうした推定には誤差が付きのもで、リス クを考慮するためには、この誤差も考慮することが必 要である。ベイズ推定では、事後分布によりこういっ た誤差を考慮することができる[4]。実プラントでのデ ータそのもので無くても、参考データがあり、ある程 度誤差が推測されているならば、それを事前分布とし て、実プラントでのデータから事後分布を求めること により、誤差をリスク評価に利用することができる。 ここでは、従来の最尤推定法[5]とベイズ推定法によ
L-11-4““ newp-clr““, expre-ry] のEL LTTり故障分布関数のパラメータを推定し比較検討した。 2. 故障分布関数のパラメータ決定 2.1 最尤推定法 * n個の機器あるいは部品を一定時間(以降、試験期間) 使用し、r個の機器/部品が時刻 1,,... , に故障した とする。一定期間(打切時間)を T とすると、例えばワ イブル分布の場合、尤度関数 Lは次式で表される。_L - 117-45““ - exper-day) | capt-ry] ()i=I LIIワイブルパラメータ m、T は、(1)式で与えられる尤度 関数 L が最大となるように決められる。(1)式の代わり に(1)式の対数 log L を考え、これを最大にする m、tは 次式で与えられる。““ + m_logt, r)-““ log(0, 1z) + (n - r)log(T / r)} = 0__ +mZlogc, r)rum)-2-{““ log(t, It) + (n - r)T““ log(T / r)} = 0[立x““ +(n-1). T)T＝=くこれらの(2)式および(3)式を連立して解くことにより、 ワイブルパラメータ m、t が求められる。399故障分布関数を用いて機器/部品等の取替周期の最 適化を図る場合、リスクを的確に評価するためには、 分布関数のパラメータの点推定値のみでなく区間推定 値も必要である。ワイブル分布の場合、2つのパラメ ータ m、tの両方とも未知の場合、両者の区間推定値 を求めるのは難しい。ここでは、m の点推定値が分か っているとして、T の区間推定を求めた場合の結果を 照会する[5]。 (4)式の統計量を考えると、2S, Ir““ が自 由度 2F の x 分布に従う。即ち、信頼係数1-a の両側 信頼区間は(5)式で表される。* s = +(n-17““( 12S28 (x (2ryu/2)““ x (2r;1- al2))片側区間も同様に(6)式で表される。-20* L1 ,, tr;m, )p(m, T) '(m,t) ==SSL(t1,.., tr; m, 7)(m, t)dmdr 2.2 ベイズ推定法 ・ ワイブル分布のパラメータm、Tは統計的に推定さ れており、その事前分布はゆ(m,)と書けるとする。 - 2.1 で述べた故障時刻データより、ワイブル分布の パラメータ m、t の事後分布の(m, r)はベイズ推定法 を用いると次のようになる。L(tys.,ltr; m, t)(m,t) 小(m,t) =““So So Lify, tr; n, 7)(m, )dmdr * ここで、L(1,, ,; m, )は(1)で与えられる尤度関数 である。ワイブル分布のパラメータ m、t が独立な場合、 小(m,n)は次のように表される。 ゆ(m, t) = pn(m). p. ()(8) この場合、(m) あるいは中.(r)は、ゆ(m, )より次 のように求めることが出来る。-9(m) = [ $(m,)dr spy(r) = [ p(m,t)dm但し、(m) および中.(r)はパラメータ m およびT に関する事前/事後分布である。てが独立な場合、-8不信賴度は、ゆ(m,)より次、0_50100-タ m およびベイズ推定の場合、こうして事後分布を求めること により、点推定値および区間推定値が、難しい数学的 取扱い無しに求める事が出来る。3. 故障分布関数のパラメータ決定試計算3.1 試計算に用いた故障分布関数ここでは、故障分布関数がワイブル分布の場合につ いて、最尤推定法とベイズ推定法による結果を比較す る。手順は、先ずパラメータ m、t を仮定し、その故 障分布関数(以降、基準関数)から、乱数により故障時刻 を求め、そのデータを基に最尤推定法とベイズ推定法 により m と を求め、それらの解の妥当性をみること とした。 今回の試計算の条件を Table 1 にまとめる。Table 1 試計算条件 Im (1,000h)|T(1,000h) |機器数n|故障数r | 2.5 150| 400 1000 -m = 2.5、T = 150,000hのワイブル分布の場合、機器 数 100台のうち、試験期間 T = 40,000h 内に故障する期 待数は約 3.6台となるため、ここでは、故障機器数を4 台とし、4台の機器の故障時刻を乱数により決定した た。 Fig.1 に m = 2.5、T = 150,000 h のワイブル分布に 4台の機器の故障データをプロットした図を示す。但し、機器の故障時刻から故障分布関数の値 F を算 出する際には、下式で示すメジアンランクを用いた。i-0.3 11n +0.4ここで、i は当該時刻までに故障した機器数を表し、 nは総機器数である。F%3D-10不信賴度10_50250300100 150 200 作動時間(1,000h)Fig.1 基準関数、 本図より、=150,000 hに対し、試験期間 T を 40,000 h とした場合、故障分布関数の極初期の段階に400おける故障しか見ていない事が分かる。こうした設定 は、機器の寿命の遙か以前に機器を交換してしまう日 本の原子力プラントの実情を考慮した結果である。 より分かりやすく、作動時間の範囲 T = 40,000 h に狭め て表示した図を Fig. 2 に示す。0.050.04。不信賴度。。。10_5_10 15 20 25 30 3540作動時間(1,000h) - Fig.2 故障データに基づく不信頼度(case 1) [ 本図に示した故障時刻データは、基準関数に良く合 致しているケースである。3.2 最尤推定法とベイズ推定法の比較 - 3.1に示したデータに対して、最尤推定法にて求め たm、T は、m = 4.0、T =90,000 h となった。m、t が この値の時の故障分布関数を Fig.3 に示す。0.050.040.03世VK* 0_5_101900/01/142025301900/02/0340IL EN AAT, anni!5_10152025303540には日(10001,不信賴度10_5_1030354015 20 25 作動時間(1,000 h)Fig.3 最尤推定法により得た故障分布関数(case 1) 作動時間が時間スケールに比べて極限られた範囲で あるため、m、t の推定値は、故障を生み出した本来 の故障分布関数とは異なったものとなっているが、デ ータの範囲では妥当と判断される結果を与えている。これに対し、同じ故障のデータについてベイズ推定 を行った結果を Fig. 4 に示す。なお、m、T の事前分布■30-35 ロ 25-30 ■20-25 015-20 ロ 10-15 ■5-10 00-5mは一様分布を仮定している。■30-35 025-30 ■20-25 015-20 010-15 ■ 5-10 00-595125)1551953215円mFig. 4(a) ベイズ推定による m、t 事後分布(case 1)11111111111HAH24311TTTTTTTTTB9876543210 5555555555HTTFig.4(6) ベイズ推定による m、事後分布(case 1) 事後分布の最大値は最尤推定法で得られた結果とほ ぼ同じ m = 4.0、T =90,000h で与えられている。しかし 事後分布を見ると、特定の方向ではあるが、かなり広 い分布を示しており、m=2.5、t=150,000h である可 能性も依然高い事が分かる。また、パラメータ m とて の間には強い相関がある事も理解される。こうした強 い相関の情報は、最尤推定法からは得る事が出来ない。 m = 2.5 として、(5)式より区間推定を行うと、95%信頼 区間として、105(x1,000h) < < 242(x1,000h) を得る。 また、m = 4.0 とすると、73(x1,000h)