PFDを求めてみた (シングル構成版)

fs-calculation-1

機能安全規格IEC-61508で登場するPFDは、システムの危険側不動作確率を表します。リスクを軽減するための「万一の為の仕組み」が、「イザ!というときにちゃんと動かないかもしれない確率=PFD」を求め、対策するリスクのレベルに応じて十分な信頼性を持つかを確認することが目的です。ここでは、このPFDの式がどのような背景で求められたのかを見てゆきます。

PFDを求めるながれ

設計したシステムが必要なSILレベルを満たすのかどうか、そのシステムのPFDを算出して証明することが必要です。そのシステムは何らかのリスクを軽減するための装置です。仕組みや保護装置で危険そのものは排除できませんが、そのような危険な事態になった時に、本当にそれが災害につながりにくくすることが装置の目的です。

そんな「万一!」「イザ!」という時にちゃんと働かない装置では、安心してリスクが軽減されたとは言えません。なので、機能安全IEC-61508では必要なSILのレベルに応じて、それに対応するシステムの信頼性をPFDという確率の数字で制限を設けているのです。

当然ながら、SILのレベルの高い=リスクの高い部分に使う装置ほど、求める信頼性=PFDの数値は小さく厳しくなってきます。

 そんなPFDは、システムの部品レベルからの故障率を積み上げることから求めます。内部の部品がどのような壊れ方をした場合に、どのような事態になるのか、それが「イザ!」という時に影響を与えるのか(つまり危険な故障なのか)そうではないのか、また故障が検知できる種類のものか、そうではないのかを全部品にわたって分類分析し、その結果を合計し、次の式に代入することで求めます

ここではシングルのシステムのPFDの式を記載します。

PFD=({ \lambda }_{ DU }+{ \lambda }_{ DD })\times \left[ \frac { { \lambda }_{ DU } }{ { \lambda }_{ D } } \times \left( \frac { T }{ 2 } +MRT \right) +\frac { { \lambda }_{ DD } }{ { \lambda }_{ D } } \times MTTR \right]

 流れとしては下記の図のようになります。算出したPFDの値が見合わなければ、システムの設計を変更するなどして再度トライです。

fs-sil-pfd-flow

 PFDの計算式はどこから来たのか

IEC-61508にはこの式をいろいろな方法で導出できると書いてありますが、イマイチ具体的な計算が書いてあるわけではないです。Markovモデルの記載では、いろいろ書いてありますがMarkovモデルは最終的に連立の一階微分方程式になるということが書いてあって、そらそうだろうと思うだけです。

そこでここからは実際にMarkovモデルを用いて微分方程式を解いてみます。PFDの式はλDUの前半とλDDの後半がありますが、説明を簡単にするために、前半と後半をそれぞれ別に調べてみます。λDU+λDDは結局λDのことですから、式は次のように表現できます。

fs-pfd-formula-exp2

1.前半部分(λDUパート)

時刻 tにおいて、健全な状態の確率をP(t)とします。つぎに、『実は壊れていて適切に動作できないけれど検知できないため修理することが出来ずそのままにされている状態』、これを潜在的な故障状態とし、この状態の確率をF(t)とします。定期検査のおかげで故障が判明し、修理中の状態にある確率をR(t)とします。

fs-pfd-fomura-markov-lambda-du

 

毎回の定期検査周期である時刻 T において、潜在故障状態のものは全て点検され、故障が判明するとすれば、F(T-)は、R(T+)に加えられ、F(T+)は、ゼロリセットします。同様に、R(T-)は、R(T+)においてF(T-)だけ増加します。(T-とは時刻Tのぎりぎり直前。T+とは時刻Tのぎりぎり直後を表すものとします。)このように、F(t)からR(t)への遷移は、図の点線のように間欠的に繰返し発生するデルタ関数のようなδ(t-nT)のように移動するというイメージです。

時刻 t=nT 毎にこのようなFからRへの移動が発生しますが、検査サイクル Tを何度も何度も繰り返した定常状態になると、次のような関係式が成り立つと考えられます。

(1)   \begin{equation*}R\left( \left( n-1 \right) T+ \right) =R\left( nT+ \right) =F(nT-)+R(nT-)\end{equation*}

すなわちインターバル Tの間に、健全な状態から潜在故障状態にゼロから累積していった量は、同じ期間 Tの間に、修理中状態から健全状態に復帰した量に等しくなる、の意味です。この関係は後ほど使用します。

さて、先ほどの図を微分方程式で表現してみます。

(2)   \begin{equation*}\frac { d }{ dt } R(t)=-\mu R(t)\end{equation*}

(3)   \begin{equation*}\frac { d }{ dt } P(t)=-{ \lambda  }_{ DU }P(t)+\mu R(t)\end{equation*}

(2)式を解くと次のようになります。

(4)   \begin{equation*}R(t)=Ro\times { e }^{ -\mu  }\end{equation*}

(4)を(3)に代入し、P(t)を求めます。

(5)   \begin{equation*}\frac { d }{ dt } P(t)=-{ \lambda }_{ DU }P(t)+\mu { R }_{ 0 }\times { e }^{ -\mu }\end{equation*}

(6)   \begin{equation*}\frac { d }{ dt } P(t)+{ \lambda }_{ DU }P(t)=\mu { R }_{ 0 }\times { e }^{ -\mu }\end{equation*}

(7)   \begin{equation*}\frac { d }{ dt } \left( P(t)\ast { e }^{ { \lambda  }_{ DU }t } \right) =\mu { R }_{ 0 }\times { e }^{ \left( { \lambda  }_{ DU }-\mu  \right) t }\end{equation*}

(8)   \begin{equation*}P(t)\times { e }^{ { \lambda  }_{ DU }t }=\frac { \mu { R }_{ 0 } }{ { \lambda  }_{ DU }-\mu  } \times { e }^{ \left( { \lambda  }_{ DU }-\mu  \right) t }+C\end{equation*}

(9)   \begin{equation*}P(t)=\frac { \mu { R }_{ 0 } }{ { \lambda  }_{ DU }-\mu  } \times { e }^{ -\mu  }+C\times { e }^{ -{ \lambda  }_{ DU }t }\end{equation*}

 定期検査をしたぎりぎり直後  t=(nT+)では、潜在故障状態 F(nT+)=0ですから、状態はPとRの2つしかありません。全体の状態の合計は「1」ですから、

(10)   \begin{equation*}P(nT+)+R(nT+)=1\end{equation*}

これは周期Tごとの初期値とも言えますから、

(11)   \begin{equation*}P(0)+R(0)=1\end{equation*}

これを(4)(9)の式で、t=0の時の関係として利用します。

(12)   \begin{eqnarray*} P(0)&=&\frac { \mu { R }_{ 0 } }{ { \lambda  }_{ DU }-\mu  } \times { e }^{ -\mu \times 0 }+C\times { e }^{ -{ \lambda  }_{ DU }\times 0 }\nonumber \\ &=&\frac { \mu { R }_{ 0 } }{ { \lambda  }_{ DU }-\mu  } +C \end{eqnarray*}

R(0)=Ro
P(0)+R(0)=1 なので

(13)   \begin{equation*}C=1-\frac { { \lambda  }_{ DU }{ R }_{ 0 } }{ { \lambda  }_{ DU }-\mu  } \end{equation*}

これでCが求まりましたので、以下のようになります。

(14)   \begin{equation*}P(t)=\frac { \mu { R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \times { e }^{ -\mu t }+\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) { e }^{ -{ \lambda }_{ DU }t }\end{equation*}

(4)と(14)の式は、ともにRoが求められていません。これは後で求めます。 先にF(t)がどのような式になるかを考えます。 最初のMarkovモデルを見ていただくと、健全なP(t)から故障率λDUだけF(t)は変化しますから、次のようにな式になります。

(15)   \begin{eqnarray*} \frac { d }{ dt } F(t)&=&{ \lambda }_{ DU }P(t)\\ F(t)&=&\int { { \lambda }_{ DU } } P(t)dt\\ &=&{ \lambda }_{ DU }\int { \frac { \mu { R }_{ 0 } }{ { \lambda }_{ DU }-\mu } } \times { e }^{ -\mu t }+\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) \times { e }^{ -{ \lambda }_{ DU }t }dt\\ &=&{ \lambda }_{ DU }\left( \left( \frac { \mu { R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) \ast \frac { 1 }{ \mu } \ast \left( 1-{ e }^{ -\mu } \right) +\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) \ast \frac { 1 }{ { \lambda }_{ DU } } \ast \left( 1-{ e }^{ -{ \lambda }_{ DU }t } \right) \right) \\ &=&\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu t} \left( 1-{ e }^{ -\mu } \right) +\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) \left( 1-{ e }^{ -{ \lambda }_{ DU }t } \right) \end{eqnarray*}

これで次のように求まりました。

(16)   \begin{equation*}P(t)&=&\frac { \mu { R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \times { e }^{ -\mu t }+\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) { e }^{ -{ \lambda }_{ DU }t }\end{equation*}

(17)   \begin{equation*}R(t)&=&Ro\times { e }^{ -\mu }\end{equation*}

(18)   \begin{equation*}F(t)&=&\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \left( 1-{ e }^{ -\mu } \right) +\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) \left( 1-{ e }^{ -{ \lambda }_{ DU }t } \right)\end{equation*}

ここで下記の「a」が十分に小さいときの関係を利用します。

(19)   \begin{equation*}{ e }^{ at }=1+a\ast t\end{equation*}

今回でてくる検知不可能な危険側故障率λDU、および修理率μは、十分に小さな値ですので、前述の(17)と(18)式をこの関係を使って簡略化します。

(20)   \begin{equation*}R(t)&=&{ R }_{ o }\ast \left( 1-\mu t \right) \end{equation*}

(21)   \begin{equation*}F(t)&=&\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \ast \mu t +\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) \ast { \lambda }_{ DU } t\end{equation*}

ここで、最初に登場した(1)の関係を使用します。

(22)   \begin{equation*}R\left( \left( n-1 \right) T+ \right)&=&R\left( nT+ \right) =F(nT-)+R(nT-)\end{equation*}

(23)   \begin{equation*}R\left(0\right)&=&F\left(T\right)+R\left(T\right)\end{equation*}

これを(20)、(21)に代入します。

(24)   \begin{equation*}R\left(T \right)&=&{ R }_{ 0 } \left(1-\mu t \right)\end{equation*}

(25)   \begin{equation*}F \left(T \right)&=&\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \ast \mu T +\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) \ast { \lambda }_{ DU } T\end{equation*}

(26)   \begin{equation*}R\left(0\right)&=&F\left(T\right)+R\left(T\right)\end{equation*}

を用います。

(27)   \begin{equation*}R\left(0\right)&=&{ R }_{ 0 } \left(1-\mu T \right)+\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \ast \mu T +\left( 1-\frac { { \lambda }_{ DU }{ R }_{ 0 } }{ { \lambda }_{ DU }-\mu } \right) \ast { \lambda }_{ DU } \end{equation*}

これを整理すると

(28)   \begin{equation*}{ R }_{ 0 }\left( \left( 1-\mu T \right) +\frac { { \lambda  }_{ DU }\mu T }{ { \lambda  }_{ DU }-\mu  } -\frac { { { \lambda  }_{ DU } }^{ 2 }T }{ { \lambda  }_{ DU }-\mu  } -1 \right) =-{ \lambda  }_{ DU }T\end{equation*}

(29)   \begin{equation*}{ R }_{ 0 }\left( \frac { -{ \lambda  }_{ DU }\mu T+{ \mu  }^{ 2 }T }{ { \lambda  }_{ DU }-\mu  } +\frac { { \lambda  }_{ DU }\mu T }{ { \lambda  }_{ DU }-\mu  } -\frac { { { \lambda  }_{ DU } }^{ 2 }T }{ { \lambda  }_{ DU }-\mu  }  \right) =-{ \lambda  }_{ DU }T\end{equation*}

(30)   \begin{equation*}{ R }_{ 0 }\left({ \lambda }_{ DU } + \mu \right)T={ \lambda }_{ DU } T\end{equation*}

(31)   \begin{equation*}{ R }_{ 0 }=\frac {{ \lambda }_{ DU }}{{ \lambda }_{ DU }+ \mu}\end{equation*}

これを(21)に代入してF(t)を求めます。

(32)   \begin{eqnarray*}F(t)&=&\frac { { { \lambda  }_{ DU } }^{ 2 } }{ { { \lambda  }_{ DU } }^{ 2 }-{ \mu  }^{ 2 } } \mu t-\frac { { { \mu  } }^{ 2 } }{ { { \lambda  }_{ DU } }^{ 2 }-{ \mu  }^{ 2 } } { \lambda  }_{ DU }t\\&=&\frac { { { \lambda  }_{ DU } }^{ 2 }\mu -{ \lambda  }_{ DU }{ \mu  }^{ 2 } }{ { { \lambda  }_{ DU } }^{ 2 }-{ \mu  }^{ 2 } } t\\&=&\frac {{ \lambda }_{ DU }\mu}{{ \lambda }_{ DU }+ \mu}t\end{eqnarray*}

PFDというのは危険側不動作確率つまりF(t)の状態とR(t)の状態の平均値の合計です。そこでF(t)とR(t)を検査周期T時間で平均をとって合計します。

(33)   \begin{equation*}{ PFD }_{ F }=\frac { 1 }{ T } \int _{ 0 }^{ T }{ F(t)dt }\\&=&\frac { { \lambda  }_{ DU }\mu  }{ { \lambda  }_{ DU }+\mu  }   \ast \frac { T }{ 2 } \end{equation*}

(34)   \begin{eqnarray*}{ PFD }_{ R }&=&\frac { 1 }{ T } \int _{ 0 }^{ T }{ R(t)dt }\\&=&\frac { 1 }{ T } \int _{ 0 }^{ T }{\frac { { \lambda  }_{ DU }  }{ { \lambda  }_{ DU }+\mu  } \ast { e }^{ -\mu t }dt\\&=&\frac{1}{T} \ast \frac { { \lambda  }_{ DU }  }{ { \lambda  }_{ DU }+\mu  } \ast \frac {1-{e}^{\mu T}}{\mu}\approx \frac { { \lambda  }_{ DU }\mu  }{ { \lambda  }_{ DU }+\mu  } \ast \frac { 1 }{ \mu  } \end{eqnarray*}

まとめると

(35)   \begin{eqnarray*}{ PFD }_{ DU }&=&{ PFD }_{ F }+{ PFD }_{ R }=\frac { { \lambda  }_{ DU }\mu  }{ { \lambda  }_{ DU }+\mu  } \left( \frac { T }{ 2 } +\frac { 1 }{ \mu  }  \right) \\&=&\frac { { \lambda  }_{ DU } }{ \frac { { \lambda  }_{ DU } }{ \mu  } +1 } \left( \frac { T }{ 2 } +\frac { 1 }{ \mu  }  \right) \approx { \lambda  }_{ DU }\left( \frac { T }{ 2 } +\frac { 1 }{ \mu  }  \right) \end{eqnarray*}

 λは一定とするとMTBFの逆数、µも一定とすると修理時間MRT(Mean Repair Time)の逆数になります。MTBFは年単位(10年以上の単位)の一方、修理にかかる時間MRTはインフラ設備では8時間または24時間程度を想定しますから、

(36)   \begin{equation*}{ \lambda  }_{ DU }\ll \mu \end{equation*}

として近似しています。MRTを用いて書き直すと

(37)   \begin{equation*}{ PFD }_{ DU }={\lambda}_{DU}\left(\frac{T}{2}+\frac{1}{\mu}\right)=}={\lambda}_{DU}\left(\frac{T}{2}+MRT \right)\end{equation*}

これで元のPFDの定義式の前半部分が導出できました。

2.PFD定義式でのMRTとMTTRの使い分けのナゼ

IEC-61508:2000年度版では、この部分のMRTがMTTRになっています。規格にある略語の定義ページをみても、MRTとMTTRの違いは書いていません。どちらも修理にかかる時間と書いてあります。おそらくMTTRにはごく短い「異常を検知する時間」が入っていると考えたためでしょう。異常になってから気づいて修理にかかるまでの時間です。MTBF+MTTRは全時間を表しますから、検知する時間はMTTRに含まれると考えます。

一方今回のR(t)は周期的な検査を行った結果判明した故障をまさに修理している正味の状態ですので検知時間を含まないと考えることが適当です。よってMTTRではなく、修理時間そのものとしてMRTを用いることとしたのでしょう。

下記から始まる式の後半部分の議論では、検知可能な故障を検知するわずかな時間も含めて修理時間としています。このため後半の議論ではMTTRを用います。

3.後半部分(λDDパート)

 今度は前半ほど複雑ではありません。時刻tにおいて健全な状態をP(t)、修辞中の時間も含めた故障状態をR(t)と表現します。今回は検知可能な故障ですから、故障が発生したらすぐに修理できると考えて、この2つの状態で表せるとします。よって R(t)=1-P(t)です。

fs-pfd-fomura-markov-lambda-dd

つまり時刻 t から t+dtの間に、P(t)×λDD×dtだけ故障し修理中のR(t)に向かい、その一方、R(t)×μ×dtが修理を終えて健全な状態に戻ってきます。これを式に表します。

(38)   \begin{eqnarray*}\frac{d}{dt}P(t)&=&-{\lambda}_{DD}P(t)+\mu R(t)\\&=&-{\lambda}_{DD}P(t)+\mu \left(1-P(t)\right)\\&=&-\left({\lambda}_{DD}+\mu \right)P(t) +\mu \end{eqnarray*}

(39)   \begin{equation*}\frac{d}{dt}\left(P(t)\ast {e}^{\left({\lambda}_{DD}+\mu}\right)t}\right)=\mu \ast {e}^{\left({\lambda}_{DD}+\mu \right)t}\end{equation*}

(40)   \begin{equation*}R(t)\ast {e}^{\left({\lambda}_{DD}+\mu}\right)t}=\frac{\mu}{{\lambda}_{DD}+\mu}\ast {e}^{\left({\lambda}_{DD}+\mu}\right)t}+C\end{equation*}

(41)   \begin{equation*}R(t)=\frac{\mu}{{\lambda}_{DD}+\mu}+C \ast {e}^{-\left({\lambda}_{DD}+\mu}\right)t}\end{equation*}

ここで最初はすべて健全な状態ですから、P(0)=1です。これを用いてCを求めます。すると以下のようになります。

(42)   \begin{equation*}P(t)=\frac{\mu}{{\lambda}_{DD}+\mu}+\frac{{\lambda}_{DD}}{{\lambda}_{DD}+\mu}\ast {e}^{-\left({\lambda}_{DD}+\mu}\right)t}\end{equation*}

時間tがずっと長く経過した場合、上式からP(t)は

(43)   \begin{equation*}\frac{\mu}{{\lambda}_{DD}+\mu}\end{equation*}

に収斂してゆきます。この値をPとします。

λは一定とするとMTBFの逆数で表されます。また同様にμも一定とするとMTTRの逆数です。これを利用すると次のようになります。

(44)   \begin{equation*}P=\frac{\mu}{{\lambda}_{DD}+\mu}=\frac{\frac{1}{MTTR}}{\frac{1}{MTBF}+\frac{1}{MTTR}}=\frac{MTBF}{MTBF+MTTR}\end{equation*}

PはAvailabilityですね。その時のPFDは、1-Pとあらわされますから、

(45)   \begin{eqnarray*}{PFD}_{DD}&=&1-\frac{MTBF}{MTBF+MTTR}=\frac{MTTR}{MTBF+MTTR}\\&=&\frac{1}{MTBF}\ast \frac {MTTR}{1+\frac {MTTR}{MTBF}}}\end{eqnarray*}

ここでMTBF>>MTTRであることを踏まえると

(46)   \begin{equation*}{PFD}_{DD}=\frac{1}{MTBF}\ast \frac {MTTR}{1+\frac {MTTR}{MTBF}}}\approx \frac{1}{MTBF}\ast MTTR\end{equation*}

MTBFはλの逆数でしたから元に戻すと結局

(47)   \begin{equation*}{PFD}_{DD}={\lambda}_{DD}\ast MTTR\end{equation*}

となり、もとのPFD算出式の後半部分が導出できました。

.前半部分と後半部分を合計し、全体のPFDを求める

 これまでの(37)と(47)から、

故障が検知できない部分と、検知できる部分のPFDが計算できました。これを合計すると、IEC-61508に掲載されている計算式となります。

(48)   \begin{eqnarray*}{PFD}_{G}&=&{PFD}_{DU}+{PFD}_{DD}\\&=&{ \lambda  }_{ DU }\times \left( \frac { T }{ 2 } +MRT \right) +{ \lambda  }_{ DD }\times MTTR\\&=&\left( { \lambda  }_{ DU }+{ \lambda  }_{ DD } \right) \times \left( \frac { { \lambda  }_{ DU } }{ { \lambda  }_{ D } } \times \left( \frac { T }{ 2 } +MRT \right) +\frac { { \lambda  }_{ DD } }{ { \lambda  }_{ D } } \times MTTR \right) \end{eqnarray*}