JP2021043687A - 故障確率評価装置、及び、故障確率評価方法 - Google Patents

Abstract

【課題】 応力や強度、もしくは、例えば荷重のように、応力や強度と関連のある物理量の発生頻度分布の中央部と裾部が、同一の確率分布に従わない場合であっても、裾部の確率密度関数を精度良く推定することにより故障確率の評価精度を向上させる技術及び高精度の故障確率評価装置を提供する【解決手段】 機械システムの故障確率を演算する故障モデル、及び、故障モデルに用いる確率変数の発生頻度分布を記憶する記憶部と、発生頻度分布の裾部の確率密度関数を極値統計モデルに基づいて推定する裾部推定部と、発生頻度分布の裾部以外の確率密度関数を推定する中央推定部と、裾部の確率密度関数と裾部以外の確率密度関数を用いて、発生頻度分布の全体の確率密度関数を推定する接続部と、を有する確率密度関数推定部と、全体の確率密度関数と故障モデルに基づき、機械システムの故障確率を演算する故障確率演算部と、を有することを特徴とする故障確率評価装置。【選択図】 図２

Description

本発明は、計測データから推定した確率密度関数を用いて、機械システムの故障確率を評価する、故障確率評価装置、及び、故障確率評価方法に関する。

例えば、風力発電機や建設機械などの機械システムの大型化と複雑化に伴い、機械システムが故障した場合の影響が大きくなってきている。機械システムの故障が周囲に及ぼすリスク（故障確率と影響度の積）を許容範囲内に抑えるためには、機械システムの設計段階と運用保守段階の双方において、機械システムの将来の故障確率を正確に評価し、適切な時期に保守作業や部品交換などの適切な対策を取ることが肝要である。故障確率の正確な評価のためには、計測データを活用して、故障確率を計算する際に用いる確率密度関数を、高精度に推定することが有効である。

なお、計測データとは、例えば、機械システムが風力発電機であるときには風速など、機械システムが稼動する環境に関する実計測データが該当する。また、機械システムの稼動状態に関する歪等の実計測データや、例えば降伏応力など、材料試験や抜き取り検査により得られる強度データのような直接測定可能な物理量に加えて、例えば応力など、機械システムが稼動する環境に関するデータと機械システムの稼動状態に関するデータを基に、コンピュータシミュレーションにより得られるシミュレーションデータのように、直接測定可能な物理量から間接的に推定可能な物理量も計測データに含むものとする。

このような計測データを活用して、機械システムの故障確率を評価する技術として、例えば特許文献１に記載の技術がある。この特許文献１の請求項１には、「経年劣化した構造物の実態に関する調査データ並びに専門家へのアンケート調査データを基に特性変数毎の確率密度関数を推定すると共に前記特性変数毎の確率密度から前記特性変数の値の組み合わせ毎の相対的な起こり易さを算出し、前記組み合わせ毎に行った有限要素解析で得られた発生応力の頻度分布を前記組み合わせ毎の相対的な起こり易さを用いて変換した前記発生応力の頻度分布から確率密度関数を推定し、該発生応力の確率密度関数と既存データから設定した強度の確率密度関数とに基づいて破壊確率と安全性指標を算出する経年劣化した構造物の信頼性評価方法」が開示されている。

特開２００７−１１３９６７号公報

機械システムの故障は、機械システムに生じる応力が機械システムの強度を超過する際に発生する。数学的には、応力Ｓの確率密度関数ｆ_Ｓと強度Ｒの確率密度関数ｆ_Ｒを用いて多重積分を行い、応力Ｓが強度Ｒを超過する確率を計算することにより、故障確率を求める。ところで、高信頼で故障確率の低い機械システムでは、特許文献１の図６に示されるように、強度Ｒに対して応力Ｓが十分に小さくなるように設計されている。よって、前述の多重積分において、応力Ｓが強度Ｒを超過する確率に最も寄与するのは、応力Ｓの確率密度関数ｆ_Ｓの右裾部と強度Ｒの確率密度関数ｆ_Ｒの左裾部が交差する領域の近傍となる。そのため、応力と強度の確率密度関数の裾部、もしくは、例えば荷重のように、応力Ｓもしくは強度Ｒと関連のある物理量の確率密度関数の裾部を精度良く推定することが、故障確率の評価精度向上において重要となる。

実際の機械システムには、例えば風力発電機おける風荷重のように、複雑でかつ瞬間的な高負荷が加わることがある。また、負荷が増した場合は、機械システム自体が複雑な非線形挙動を示すこともある。そのため、最も頻繁に発生するオーダーの応力Ｓ（発生応力の頻度分布の中央部）と、稀に生じる高い応力（発生応力の頻度分布の右裾部）が、同一の確率分布に従う保証はない。しかしながら、特許文献１では、発生応力の頻度分布全体を１つの確率密度関数により推定している。そのため、発生応力の頻度分布の中央部と裾部が同一の確率分布に従わない場合には、発生頻度の高い中央部に適合する確率密度関数が選択され、裾部の推定精度が低下するという課題があった。

本発明はこのような事情を鑑みてなされたもので、応力や強度、もしくは、例えば荷重のように、応力や強度と関連のある物理量の発生頻度分布の中央部と裾部が、同一の確率分布に従わない場合であっても、裾部の確率密度関数を精度良く推定することにより故障確率の評価精度を向上させる故障確率評価方法および故障確率評価装置を提供することを目的とする。

本発明に係る故障確率評価装置は、機械システムの故障確率を評価するものであって、前記機械システムの故障確率を演算する故障モデル、及び、前記故障モデルに用いる確率変数の発生頻度分布を記憶する記憶部と、前記発生頻度分布の裾部の確率密度関数を極値統計モデルに基づいて推定する裾部推定部と、前記発生頻度分布の裾部以外の確率密度関数を推定する中央推定部と、前記裾部の確率密度関数と前記裾部以外の確率密度関数を用いて、前記発生頻度分布の全体の確率密度関数を推定する接続部と、を有する確率密度関数推定部と、前記全体の確率密度関数と前記故障モデルに基づき、前記機械システムの故障確率を演算する故障確率演算部と、を有するものとした。

また、本発明に係る故障確率評価方法は、機械システムの故障確率を演算する故障モデルに用いられる確率変数の発生頻度分布の裾部の確率密度関数を極値統計モデルに基づいて推定するステップと、前記発生頻度分布の裾部以外の確率密度関数を推定するステップと、前記裾部の確率密度関数と前記裾部以外の確率密度関数を用いて、前記発生頻度分布の全体の確率密度関数を推定するステップと、前記全体の確率密度関数と前記故障モデルに基づき、前記機械システムの故障確率を演算するステップと、を有するものとした。

本発明によれば、応力や強度の確率密度関数の裾部、もしくは、例えば荷重のように、応力や強度と関連のある物理量の確率密度関数の裾部を精度良く推定することが可能となり、高精度の故障確率評価装置を提供することが可能となる。上記した以外の課題や構成及び効果は、以下の実施形態の説明により明らかにされる。

一実施例の故障確率評価装置の全体概略構成図。 図１の故障確率評価装置の機能ブロック図。 図１の故障確率評価装置により確率密度関数を推定するフローチャート。 取得データの発生頻度分布と推定された確率密度関数の関係を示す模式図。 図１の故障確率評価装置により疲労損傷の故障確率を評価するフローチャート。 累積損傷度と疲労強度の確率密度関数の関係を示す模式図。 表示装置の表示画面例を示す図。 数値実験で発生させた応力の発生頻度分布図。 図６の応力の頻度分布右裾部に対応する、本発明で推定した確率密度関数の形状を示す図。

以下、図１〜７を用いて、本発明の一実施例に係る故障確率評価装置について説明する。なお、本実施例では、機械システムの具体例として風力発電プラントを挙げるが、本発明の適用される機械システムは風力発電プラントに限定されるものではない。

図１は、本実施例の故障確率評価装置１の全体概略構成図である。ここに示すように、故障確率評価装置１は、機械システム２、設計・材料データベース３、計算サーバ４と、信号線、インターネット、イントラネットなどのネットワークを介して接続されており、相互にデータの送受信を行うことができる。

機械システム２には、使用環境や運転状態を計測するためのセンサが多数実装されている。それらのセンサで計測された環境データ、運転データなどの計測データは、故障確率評価装置１に送られ、故障確率の評価に利用される。ここで、環境データとは、機械システム２が曝される環境に関するデータを含むデータであり、風力発電プラントの場合、例えば風車の風速、風向などの風況データは環境データに含まれる。洋上に設置された風力発電プラントの場合、風況データに加えて、波長や波高などの海況データも環境データの範疇である。また、運転データとは、速度、加速度、回転速度、回転角など、機械システム２の稼働状態に関連するデータである。風力発電プラントの場合、風車の発電量、発電機の回転速度、アジマス角、ナセル角などは運転データの範疇である。プラントでは、環境データ、運転データは制御計測（SCADA、SCADA,Supervisory Control And Data Aquisition）として計測されることが多い。ただし、状態監視（CMS、Condition Monitoring System）や構造物監視（SHM、Structural Health Monitoring）を含むプラントに関しては、例えば歪など、これらで計測されるデータでも機械システム２の使用環境や稼働状態に関連するデータであれば、本実施例の環境データや運転データに含まれる。

設計・材料データベース３のデータも故障確率評価装置１に送られ、評価に利用される。風力発電プラントの場合、設計データは、例えば風力発電機の図面など機械システム２の形状に関するデータ、または風況調査や気象庁データベースなどから得られる風速・風向データなどを含んでいる。また、材料データは、風力発電機を構成する材料の特性や、ボルト締結や溶接継手などの構造の特性を含んでいる。

計算サーバ４は、機械システム２から故障確率評価装置１に送られた環境データや運転データ、または、設計・材料データベースから故障確率評価装置１に送られた各種データの一部を、故障確率評価装置１から受け取り、それを基に故障確率評価装置１が要求する物理量を計算した後に、故障確率評価装置１へデータを送り返す。風力発電プラントの場合、例えば、風速や歪などの環境データや運転データ及び空力モデルなどの設計・材料データから、風力発電プラントの各構成部材に発生する応力やモーメントを計算する。なお、計算サーバ４は必ずしも故障確率評価装置１から独立している必要は無く、両者が一体となった構成としても良い。

故障確率評価装置１は、設計・材料データベース３や計算サーバ４と連携して機械システム２の故障確率を評価する装置であり、故障確率評価部１０、入力部１１、表示部１２、記憶部１３を備える。故障確率評価部１０は、確率密度関数推定部１０ａ、故障確率演算部１０ｂを備える。確率密度関数推定部１０ａは、裾部推定部１０ａ１、中央推定部１０ａ２、接続部１０ａ３を備える。入力部１１は、例えば、キーボード又はマウスなどの入力装置であり、利用者が故障確率評価装置１に関して何らかの入力をする際に用いられる。表示部１２は、例えば、液晶ディスプレイ（ＬＣＤ）又は有機ＥＬディスプレイなどの表示装置、及びプリンタなどの印字出力装置であり、故障確率評価部１０による処理の過程や結果、或いは故障確率評価装置１の利用者のための対話的な処理のための画面を表示する。記憶部１３は、故障モデル１３ａ、設計・材料データ１３ｂ、発生頻度分布１３ｃの各データを記憶する。

図２は、図１に示す故障確率評価装置１の機能ブロック図である。ここに示すように、故障確率評価装置１は、図１で示した構成に加え、データ取得部１４、通信Ｉ／Ｆ１５、入力Ｉ／Ｆ１６、表示Ｉ／Ｆ１７を備えており、これらは相互に内部バス１８にてアクセス可能に接続されている。通信Ｉ／Ｆ１５は、機械システム２、設計・材料データベース３、計算サーバ４から送信される各種データを受信する。また、入力Ｉ／Ｆ１６は入力部１１に接続されており、表示Ｉ／Ｆ１７は表示部１２に接続されている。なお、以下では、データの送受信に、内部バス１８、通信Ｉ／Ｆ１５、入力Ｉ／Ｆ１６、表示Ｉ／Ｆ１７等の介在が自明な場合はその介在を適宜省略する。

故障確率評価装置１は実際にはパソコン等の計算機で構成されている。そして、図２に示す故障確率評価部１０、記憶部１３、データ取得部１４は、例えば、図示しないＣＰＵ（Central Processing Unit）などのプロセッサ、各種プログラムを格納するＲＯＭ、演算過程のデータを一時的に格納するＲＡＭ、外部記憶装置などの記憶装置にて実現されると共に、ＣＰＵなどのプロセッサがＲＯＭに格納された各種プログラムを読み出し実行し、実行結果である演算結果をＲＡＭ又は外部記憶装置に格納する。

記憶部１３は、入力部１１と入力Ｉ／Ｆ１６を介して、もしくは、設計・材料データベース３と通信Ｉ／Ｆ１５を介して、利用者により予め入力された、機械システム２の状態が健全か故障かを判断するための故障モデル１３ａ、及び、設計・材料データ１３ｂを、所定の記憶領域に格納している。

データ取得部１４は、機械システム２からの各種計測データを取得し、例えば、Ａ／Ｄ変換処理、平滑化処理（ノイズ除去）、或いは正規化処理などを実行する。また、データ取得部１４は、故障モデル１３ａで用いられる計測データを取得するため、必要に応じて、機械システム２の計測データに各種の処理を施した後に、計算サーバ４にデータを送信し、計算結果を受信する。

故障モデル１３ａは、機械システム２の使用環境や稼動状態に関わる様々な不確定性を考慮するための確率変数を含む。例えば、風力発電プラントの場合、確率変数は風速や風力発電プラントの各構成部材の強度などである。また、確率変数は、歪や風速などの直接測定可能な物理量に限らず、直接測定可能な物理量から間接的に推定される応力などの物理量でも構わない。もちろん、故障モデル１３ａには確率変数以外にも確定的な値をとる変数も含まれていて構わなく、例えば、構造物の板厚や質量など、設計・材料データ１３ｂにおいて不確定性の小さな変数を確定値として組み込むことが可能である。

故障モデル１３ａは、例えば、限界状態式Ｇ（Ｘ_１、Ｘ_２、・・・、Ｘ_ｎ）である。ここで、Ｘは確率変数であり、Ｇ＞０ならば機械システム２は健全、Ｇ＜０ならば故障と判断する。例えば、最も単純な限界状態式Ｇは、Ｇ（Ｓ、Ｒ）＝Ｒ−Ｓとなる（Ｓは応力、Ｒは強度）。なお、限界状態式Ｇには確率変数以外にも、設計・材料データ１３ｂの変数を確定値として組み込んでも構わない。もちろん、限界状態式Ｇは、機械システムの構成部品毎に定義可能であり、その場合は記憶部１３に複数の故障モデル１３ａが記憶される。

故障モデル１３ａにおける確率変数のうち、例えば風力発電プラントの各構成部材の強度など、事前に材料試験などによって確率密度関数が求められている確率変数については、記憶部１３は、設計・材料データ１３ｂとして、所定の記憶領域に格納している。一方で、データ取得部１４が、機械システム２や計算サーバ４から取得した各種データに関しては、機械システム２が稼動する以前には確率密度関数ｆが求められていない。そのため、記憶部１３は、確率密度関数ｆが求められていない計測データをデータ取得部１４から取得し、各データの発生頻度分布１３ｃを所定の記憶領域に格納する。なお、設計・材料データ１３ｂとして確率密度関数ｆが記憶されている確率変数に関しても、確率密度関数ｆの推定に用いられた計測データが入手できるならば、発生頻度分布１３ｃに記憶しても構わない。

確率密度関数推定部１０ａを構成する、裾部推定部１０ａ１、中央推定部１０ａ２、接続部１０ａ３は、記憶部１３から取得した各種計測データの発生頻度分布１３ｃを参照し、確率密度関数ｆを推定する。まず、裾部推定部１０ａ１が、発生頻度分布１３ｃの右及び／又は左の裾部端部を調査し、裾部の確率密度関数ｆ_Ｂが適合する範囲を確定する。裾部の確率密度関数ｆ_Ｂが適合する範囲は、右裾の場合は、裾部推定部１０ａ１が確定した閾値ｔｈ_Ｒ以上の範囲であり、左裾の場合は、裾部推定部１０ａ１が確定した閾値ｔｈ_Ｌ以下の範囲となる。裾部推定部１０ａ１は、中央推定部１０ａ２には確定した閾値の値を、接続部１０ａ３には確定した閾値の値に加えて、推定した裾部に適合する確率密度関数ｆ_Ｂの情報を送る。

中央推定部１０ａ２は、裾部推定部１０ａ１から得た閾値を基に、裾部以外の領域において、発生頻度分布１３ｃに適合する確率密度関数ｆ_ｃを推定する。中央推定部１０ａ２の評価対象となる発生頻度分布１３ｃの領域は、裾部推定部１０ａ１から右裾の閾値ｔｈ_Ｒのみが得られた場合は閾値ｔｈ_Ｒ以下の範囲、裾部推定部１０ａ１から左裾の閾値ｔｈ_Ｌのみが得られた場合は閾値ｔｈ_Ｌ以上の範囲、また、裾部推定部１０ａ１から右裾と左裾の閾値ｔｈ_Ｒ、ｔｈ_Ｌが両方得られた場合は、２つの閾値を上限と下限とする範囲となる。中央推定部１０ａ２は、接続部１０ａ３に、推定した中央部に適合する確率密度関数ｆ_Ｃの情報を送る。

接続部１０ａ３は、裾部推定部１０ａ１から閾値と裾部に適合した確率密度関数ｆ_Ｂを得る。また、中央推定部１０ａ２から中央部に適合した確率密度関数ｆ_Ｃを得る。そして、接続部１０ａ３は、裾部の確率密度関数ｆ_Ｂと中央部の確率密度関数ｆ_Ｃが閾値近傍で滑らかに接続するように、閾値近傍での確率密度関数ｆ_Ｊ形状を微修正することで、発生頻度分布１３ｃの全領域に適合する確率密度関数ｆを推定する。

なお、必ずしも、発生頻度分布１３ｃに記憶された全てのデータに対して、裾部推定部１０ａ１、中央推定部１０ａ２、接続部１０ａ３を用いて確率密度関数を推定する必要はない。すなわち、計測データが、物理的な妥当性を持って１つの確率密度関数で表現できる場合は、確率密度関数推定部１０ａは、従来技術と同様に、計測データの全域に適合する確率密度関数を推定する。例えば、脆性材料の引張強度は、微視的に最も弱い部分が破壊の起点となるため（最弱リンク説）、計測データの大小に関わらず３パラメータのワイブル分布が適合する。また、金属材料の引張強度や降伏応力は、材料を構成する微細な要素の特性の平均値であると考えられるので、中心極限定理より、これらの物理量は正規分布に従う。強度データは、実験コストや時間的な制約上、データ数が少なくなることが一般的であるので、確率密度関数の推定精度を上げるためには、計測データを分割して複数の確率密度関数で推定するよりも、データ全体を１つのパラメトリックな確率密度関数（正規分布、対数正規分布、ワイブル分布など）により推定することが望ましい。一方、荷重や応力と言った負荷に関係する物理量では、最も頻繁に発生するオーダーの負荷（発生頻度分布の中央部）が、稀に生じる高い負荷（発生頻度分布の右裾部）と、同一の確率分布に従う保証はなく、また、強度データと比較して計測可能なデータ数も多いため、裾部と中央部を異なる確率密度関数により推定することが望ましい。

また、故障確率演算部１０ｂは、故障モデル１３ａに用いられる確定値の変数を設計・材料データ１３ｂから取得し、故障モデル１３ａに用いられる確率変数の確率密度関数ｆを、確率密度関数推定部１０ａまたは設計・材料データ１３ｂに記憶された確率密度関数から取得する。故障モデル１３ａには、上述の変数と確率変数から故障確率を演算する方法が定義されており、故障確率演算部１０ｂは、その演算方法に従い、故障確率を計算する。

例えば故障モデル１３ａが限界状態式Ｇ（Ｘ_１、Ｘ_２、・・・、Ｘ_ｎ）で与えられている場合、故障確率Ｐ_ｆは、Ｘ_１、Ｘ_２、・・・Ｘ_ｎの同時確率密度関数ｆ（Ｘ_１，Ｘ_２，・・・、Ｘ_ｎ）を用いて、（式１）で与えられる。

（式１）はモンテカルロシミュレーションなどにより評価できる。また、（式１）の多重積分にかかる計算コストを削減するため、（式１）の近似解を与える手法として知られている一次信頼性理論や二次信頼性理論により、故障確率Ｐ_ｆを求めても良い。例えば、（式１）でＧ（Ｓ、Ｒ）＝Ｒ−Ｓとすれば（Ｓは応力、Ｒは強度）、静的な故障、すなわち過大な応力が１回加わることにより、機械システム２が故障する場合の故障確率が求められる。

表示部１２は、故障確率評価部１０が求めた故障確率Ｐ_ｆを取得し、表示部１２を構成する表示装置の表示画面上に表示する。これにより、故障確率評価装置１の利用者は、容易に機械システム２のリスクを視認することが可能となる。なお、表示部１２の表示画面の詳細は後述する。

＜静的な故障確率の推定方法＞
図３Ａは、故障確率評価装置１の確率密度関数推定部１０ａにおけるデータ処理の詳細を示すフローチャート図である。ここでは、裾部推定部１０ａ１にて、確率密度関数ｆの右裾部の確率密度関数ｆ_Ｂを推定するフローを例にとる。

まず、ステップＳ３１にて、裾部推定部１０ａ１は、記憶部１３から計測データの発生頻度分布１３ｃを取得する。図３Ｂに示すヒストグラムＨが、発生頻度分布１３ｃの具体的な内容の一例である。

次に、ステップＳ３２にて、裾部推定部１０ａ１は、極値統計モデルに基づき、閾値ｕ以上の領域において、計測データのヒストグラムＨの右裾部に適合する確率密度関数ｆ_Ｂを推定する。ここでは、極値統計モデルとして、Peak Over Threshold Method（ＰＯＴ法）を採用する。閾値ｕを十分大きく取れば、閾値ｕ以上の領域では、計測データが従う確率分布の関数形に依らず、右裾部の確率密度関数ｆ_Ｂは、以下の（式２）に示す、一般パレート分布を変形することにより得られる確率密度関数で表される。

ここで、Ｆは計測データの累積分布関数、σは尺度パラメータ、ξは形状パラメータ（極値指数）である。（式２）の、１−Ｆ（ｕ）は、ヒストグラムＨにおいて、閾値ｕを超過する計測データの数と全データの数の比により推定できる。また、最尤推定により、尺度パラメータσおよび形状パラメータξが推定できる。ＰＯＴ法で適切な閾値ｕを決定するための手法は、周知であるので詳細は割愛する。ＰＯＴ法では、計測データの右裾部が従う確率分布が不明な場合でも、推定に使用するのは（式２）のみで良いという利点がある。また閾値ｕを求める手段が提供されているため、ヒストグラムＨにおいて、裾部の領域を定義しやすいという利点がある。

次に、ステップＳ３３にて、中央推定部１０ａ２は、閾値ｕ以下の領域の確率密度関数ｆ_Ｃの推定を行う。確率密度関数ｆ_Ｃを、例えば、正規分布、対数正規分布、ワイブル分布などのパラメトリックな確率密度関数により推定する場合は、適合度を判定し、例えばＡＩＣ（赤池情報量基準）が最小となる確率分布を選択する。または、中央部はデータ数も多くなるため、カーネル法により確率密度関数を推定しても良い。この場合、パラメトリックな確率密度関数の推定と異なり、計測データのヒストグラムＨが従う確率密度関数ｆの関数形が不明でも推定が可能となる。

最後に、ステップＳ３４にて、接続部１０ａ３は、確率密度関数ｆ_Ｂと確率密度関数ｆ_Ｃを滑らかに接続することで、ヒストグラムＨの全体に適合する確率密度関数ｆを推定する。より具体的には、まず、閾値ｕの左側のｖ≦ｘ＜ｕの範囲に３次関数の確率密度関数ｆ_Ｊ：ａｘ^３＋ｂｘ^２＋ｃｘ＋ｄを設定した後、ｘ＜ｖに確率密度関数ｆ_Ｃを、ｖ≦ｘ＜ｕに確率密度関数ｆ_Ｊを、ｕ≦ｘに確率密度関数ｆ_Ｂを割り当てることで、ヒストグラムＨの全体に適合する確率密度関数ｆを推定する。

このステップＳ３４では、（ａ、ｂ、ｃ、ｄ、ｖ）の５つの未知数を定める必要があるが、これらは次の指針により求めることができる。すなわち、まず、ｘ＝ｖにおいては、確率密度関数ｆ_Ｃと確率密度関数ｆ_Ｊが同値であること、及び、両関数の導関数（微分値）が同値であること、の２条件を満たす必要がある。また、ｘ＝ｕにおいては、確率密度関数ｆ_Ｂと確率密度関数ｆ_Ｊが同値であること、及び、両関数の導関数（微分値）が同値であること、の２条件を満たす必要がある。さらに、確率密度関数ｆの全体での積分が１になるという１条件を満たす必要がある。これらの５条件に基づいて、上記した５つの未知数（ａ、ｂ、ｃ、ｄ、ｖ）を決定することができる。

もちろん、確率密度関数ｆ_Ｂと確率密度関数ｆ_Ｃの接続方法は、上記の方法に限定されないことは言うまでもない。なお、確率密度関数推定部１０ａで推定するのが応力Ｓの確率密度関数ｆ_Ｓの場合は、強度の確率密度関数ｆ_Ｒと交差する領域を予め把握しておき、交差領域近傍での適合度が高くなる接続方法を採用することが望ましい。

＜疲労損傷の故障確率の推定方法＞
図４Ａは、故障確率評価装置１の確率密度関数推定部１０ａと故障確率演算部１０ｂで、静的な故障ではなく、疲労損傷の故障確率を評価する際のフローチャート図である。なお、このフローチャートにおいては、応力振幅δが第一の確率変数、累積損傷度Ｄが第二の確率変数である。

まず、ステップＳ４１にて、記憶部１３は、計算サーバ４が計算した応力Ｓの時系列データを取得し、それに対してレインフロー法を用いて、応力振幅δの発生頻度分布１３ｃを作成し、記憶する。また、故障確率演算部１０ｂは、設計・材料データ１３ｂに記憶されている材料の寿命曲線及び、応力振幅δの発生頻度分布１３ｃに対してマイナー則を用い、時刻ｔ_０における累積損傷度Ｄ_０を評価する。

次に、ステップＳ４２にて、確率密度関数推定部１０ａは、応力振幅δの発生頻度分布１３ｃから、応力振幅δの確率密度関数ｆ_δを推定する。材料の寿命曲線は、応力振幅δの冪乗に従うので、高い応力振幅δでは、１回あたりに発生する損傷度が非線形に増大する。よって、応力振幅δの確率密度関数ｆ_δの右裾部を精度良く推定することが、疲労損傷確率の評価において重要であるため、確率密度関数推定部１０ａにより応力振幅δの確率密度関数ｆ_δを推定する。ここでの確率密度関数の推定法は、図３と同様なので、重複する説明は割愛する。

なお、上記のステップＳ４２の説明では、弾性変形範囲内で使用され、高サイクル疲労が生じ得る部品の故障確率評価を念頭において、応力振幅δの発生頻度分布１３ｃを求めたが、発生頻度分布１３ｃは応力振幅δに限定されない。例えば、電子部品の接合部に用いられるはんだは、電子部品の使用中に弾性変形範囲を超して塑性変形を起こすため、材料の寿命曲線は応力振幅δでなく歪振幅により整理される。よって、塑性変形により低サイクル疲労が生じ得る部材の故障確率評価の場合は、歪振幅の発生頻度分布を求めることが望ましい。歪振幅の発生頻度分布に対しても、マイナー則を用いて時刻ｔ_０における累積損傷度Ｄ_０が評価できることは言うまでもない。

次にステップＳ４３にて、故障確率演算部１０ｂは、将来の任意の時刻である時刻ｔ_０＋Δｔにおける累積損傷度Ｄの確率密度関数ｆ_Ｄを推定する。時刻ｔ_０までにＮ_０回の応力振幅δが発生していた場合、機械システム２の利用状況が同等であれば、時刻ｔ_０から時刻ｔ_０＋Δｔまでの間に、さらにＮ_０×Δｔ／ｔ_０回の応力振幅δが発生すると考えられる。よって、故障確率演算部１０ｂにおいて、ステップＳ４２で求めた応力振幅δの確率密度関数ｆ_δに従うランダム変数をＮ_０×Δｔ／ｔ_０回発生させ、それにマイナー則を適用し、時刻ｔ_０における累積損傷度Ｄ_０を加えれば、時刻ｔ_０＋Δｔにおける累積損傷度Ｄ_１が得られる。この試行を繰り返すと、時刻ｔ_０＋Δｔにおける累積損傷度Ｄ_１は、ばらついた値となり、試行回数を十分大きく取ることで、時刻ｔ_０＋Δｔにおける累積損傷度の確率密度関数ｆ_Ｄを得ることができる。

次にステップＳ４４にて、故障確率演算部１０ｂは、設計・材料データ１３ｂに記憶されている材料の疲労強度の確率密度関数ｆ_Ｆを取得し、ステップＳ４３で取得した時刻ｔ_０＋Δｔにおける累積損傷度の確率密度関数ｆ_Ｄと比較する。マイナー則により得られる累積損傷度Ｄ_０と比較するため、材料の疲労強度の確率密度関数ｆ_Ｆは、平均を１とする確率密度関数で与えられる。疲労強度データの取得には多大な時間を要し、データ数が少なくなるため、疲労強度の分布は１つのパラメトリック関数で推定することが望ましい。疲労強度の確率密度関数ｆ_Ｆは、対数正規分布や３パラメータのワイブル分布で良く近似できることが実験的に確認されている。それぞれを材料の疲労強度Ｆ_Ｒ、時刻ｔ_０＋Δｔにおける累積損傷度Ｆ_Ｄとすると、限界状態式はＧ＝Ｆ_Ｒ−Ｆ_Ｄとなる。確率密度関数ｆ_Ｄと確率密度関数ｆ_Ｆを用いて、Ｇ＜０となる確率を、例えばモンテカルロ法などにより評価する。

＜表示部の表示画面＞
図５は、図１に示す表示部１２を構成する表示装置の表示画面１２ａを示す図である。故障確率評価装置１の利用者が、変数選択領域１２ｂのプルダウンから、目的の変数「応力」を選択すると、描画領域１２ｄには、選択された変数「応力」の発生頻度分布（ヒストグラムＨ）と、確率密度関数推定部１０ａにより推定された応力Ｓの確率密度関数ｆ_Ｓが重ねて描画される。このとき、裾部の確率密度関数ｆ_ＳＢと中央部の確率密度関数ｆ_ＳＣは異なる線で表示される。また、裾部と中央部の確率密度関数のパラメータの詳細は、表示領域１２ｅと１２ｆに表示される。そして、評価期間選択領域１２ｃのプルダウンから、今後の評価期間（例えば５０年）を選択すると、評価結果表示領域１２ｇには、静的な負荷による故障確率、今後の評価期間（５０年間）に疲労損傷が発生する確率や、今後の評価期間内（５０年）に平均１回発生するレベルの値（再現レベル）等が表示される。これにより、故障確率評価装置１の利用者は、選択した変数の確率密度関数の適合性と機械システム２のリスクを容易に視認することが可能となる。

＜本実施例の有効性の確認＞
図６および図７は、本実施例の効果を示す数値実験の結果を示す図である。数値実験で発生させた応力Ｓの真の確率密度関数ｆ_ｓは、２７０ＭＰａ未満では正規分布（平均２００ＭＰａ、標準偏差３０ＭＰａ）に、２７０ＭＰａ以上ではワイブル分布（尺度パラメータ１７０ＭＰａ、形状パラメータ２．８）に従うというものである。これにより、機械システム２の非線形性により、応力Ｓの確率密度関数ｆ_ｓの右裾部は中央部と異なることを模擬している。また、ここでは、数値実験に基づく強度の確率密度関数ｆ_Ｒ０は正規分布（平均３８０ＭＰａ、標準偏差２７．５ＭＰａ）とした。

図６では、信頼性の高い機械システム２を模擬するため、応力Ｓの確率密度関数ｆ_ｓの右裾部と強度の確率密度関数ｆ_Ｒ０の左裾の重複領域を小さくしており、この機械システム２が故障しにくいものであることを表現している。確率密度関数ｆ_ｓと確率密度関数ｆ_Ｒ０を用いて１０^７回のモンテカルロシミュレーションを行い、応力Ｓが強度よりも大きくなった回数から故障確率を求めたところ、真の確率密度関数ｆ_ｓを前提とした場合の故障確率は２．０５×１０^−４となった。

確率密度関数ｆ_Ｓに従う乱数を３６５０回発生せた場合の応力Ｓの発生頻度分布が図６に示すヒストグラムＨ_０である。このヒストグラムＨ_０の全域を、従来技術により、１つの確率密度関数ｆ_Ｓ０で推定したところ、最も適合度の高い関数は、発生頻度分布の中央部が合致する確率分布である正規分布（平均２０３ＭＰａ、標準偏差３２ＭＰａ）であった。そのため、ヒストグラムＨ_０の右裾では、ワイブル分布に従うヒストグラムＨ_０と、正規分布に従う確率密度関数ｆ_Ｓ０の差が大きくなっている。

図７は、図６のヒストグラムＨ_０の右裾の拡大図である。また、極値統計モデルの一種であるＰＯＴ法により推定した確率密度関数ｆ_ＳＢを重ねて示す。ＰＯＴ法による確率密度関数ｆ_ＳＢは、データが極めて少ない３４０〜３５０ＭＰａの領域まで含めて、ヒストグラムＨ_０と良く一致している。

従来技術による確率密度関数ｆ_Ｓ０と、本発明の技術による確率密度関数ｆ_ＳＢと、強度の確率密度関数ｆ_Ｒ０の三者を用いて、１０^７回のモンテカルロシミュレーションを行ったところ、従来技術で推定した故障確率は１．４７×１０^−５となったのに対し、本発明で推定した故障確率は２．４１×１０^−４であった。すなわち、真の確率密度関数ｆ_ｓを用いた故障確率（２．０５×１０^−４）と、本実施例で推定した故障確率（２．４１×１０^−４）が良く一致するのに対し、従来技術で推定した故障確率（１．４７×１０^−５）は真値より１桁小さく精度が悪いことから、従来技術に対する本実施例の有効性が確認できた。

以上で説明したように、本実施例によれば、応力や強度、もしくは、例えば荷重のように、応力や強度と関連のある物理量の発生頻度分布の中央部と裾部が、同一の確率分布に従わない場合であっても、裾部の確率密度関数を精度良く推定することにより故障確率の評価精度を向上させることができる。

なお、上述の実施例においては、１基の機械システム２の故障確率評価を行う場合を一例として示したが、複数の機械システム２が稼動する場合にも同様に適用できることは言うまでもない。

また、故障確率評価対象の機械システム２が新たに開発された機械システムの場合に、故障評価精度を向上させるため、同型の機械システム、類似の機械システムで推定された確率密度関数を、新たに開発された機械システムの設計・材料データとして活用することができるのは言うまでもない。

また、本発明は上記した実施例に限定されるものではなく、様々な変形例が含まれる。
例えば、上記した実施例は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。

１ 故障確率評価装置
１０ 故障確率評価部
１０ａ 確率密度関数推定部
１０ａ１ 裾部推定部
１０ａ２ 中央推定部
１０ａ３ 接続部
１０ｂ 故障確率演算部
１１ 入力部
１２ 表示部
１２ａ 表示画面
１２ｂ 変数選択領域
１２ｃ 評価期間選択領域
１２ｄ 選択された変数の確率密度関数表示領域
１２ｅ 裾部の確率密度関数情報の表示領域
１２ｆ 中央部の確率密度関数情報の表示領域
１２ｇ 評価結果表示領域
１３ 記憶部
１３ａ 故障モデル
１３ｂ 設計・材料データ
１３ｃ 発生頻度分布
１４ データ取得部
１５ 通信Ｉ／Ｆ
１６ 入力Ｉ／Ｆ
１７ 表示Ｉ／Ｆ
１８ 内部バス
２ 機械システム
３ 設計・材料データベース
４ 計算サーバ
Ｈ 取得データの発生頻度分布
Ｈ_０ 数値実験で発生させた発生頻度分布
ｔｈ_Ｒ、ｔｈ_Ｌ、ｕ 閾値
ｆ、ｆ_Ｓ 確率密度関数
ｆ_Ｂ、ｆ_ＳＢ 裾部の確率密度関数
ｆ_Ｃ、ｆ_ＳＣ 中央部の確率密度関数
ｆ_Ｊ、ｆ_ＳＪ 接続点近傍の確率密度関数
ｆ_Ｓ０ 従来技術により推定した応力の確率密度関数
Ｄ_０ 時刻ｔ_０における累積損傷度
ｆ_Ｄ 時刻ｔ_０＋Δｔにおける累積損傷度の確率密度関数
ｆ_Ｒ 疲労強度の確率密度関数
ｆ_Ｒ０ 数値実験に基づく強度の確率密度関数

Claims (11)

1. 機械システムの故障確率を評価する故障確率評価装置であって、
   前記機械システムの故障確率を演算する故障モデル、及び、前記故障モデルに用いる確率変数の発生頻度分布を記憶する記憶部と、
   前記発生頻度分布の裾部の確率密度関数を極値統計モデルに基づいて推定する裾部推定部と、
   前記発生頻度分布の裾部以外の確率密度関数を推定する中央推定部と、
   前記裾部の確率密度関数と前記裾部以外の確率密度関数を用いて、前記発生頻度分布の全体の確率密度関数を推定する接続部と、
   を有する確率密度関数推定部と、
   前記全体の確率密度関数と前記故障モデルに基づき、前記機械システムの故障確率を演算する故障確率演算部と、
   を有することを特徴とする故障確率評価装置。
2. 請求項１に記載の故障確率評価装置において、
   前記裾部推定部は、一般パレート分布により確率密度関数を推定することを特徴とする故障確率評価装置。
3. 請求項２に記載の故障確率評価装置において、
   前記中央推定部は、カーネル法により確率密度関数を推定することを特徴とする故障確率評価装置。
4. 請求項１〜３のいずれかに記載の故障確率評価装置において、
   前記故障モデルに用いる確率変数のうち、前記接続部で全体の確率密度関数を推定した確率変数以外の少なくとも１つの確率変数は、１つのパラメトリックな確率分布で与えられることを特徴とする故障確率評価装置。
5. 請求項４に記載の故障確率評価装置において、
   前記接続部で全体の確率密度関数を推定した確率変数は、前記機械システムの応力であり、前記パラメトリックな確率分布で与えられる確率変数は、前記機械システムの強度であることを特徴とする故障確率評価装置。
6. 請求項１〜３のいずれかに記載の故障確率評価装置において、
   前記記憶部には、第一の確率変数を用いて前記故障モデルに用いる第二の確率変数を、演算する方法を記憶しており、
   前記故障確率演算部は、前記第一の確率変数の確率密度関数を用いたモンテカルロシミュレーションと、前記第一の確率変数から前記第二の確率変数を演算する方法を組み合わせることにより、演算時点から所定の時刻が経過した時点での前記第二の確率変数の確率密度関数を求め、演算時点から所定の時刻が経過した時点での故障確率を演算することを特徴とする故障確率評価装置。
7. 請求項６に記載の故障確率評価装置において、
   前記第一の確率変数は応力振幅または歪振幅であり、
   前記第二の確率変数は応力振幅または歪振幅と前記記憶部に記憶された疲労寿命曲線から演算される累積損傷度であることを特徴とする故障確率評価装置。
8. 請求項１〜３に記載のいずれかの故障確率評価装置において、
   前記故障確率演算部は、一次信頼性理論もしくは二次信頼性理論に基づき故障確率を演算することを特徴とする故障確率評価装置。
9. 請求項１〜３に記載のいずれかの故障確率評価装置において、
   さらに、前記接続部により推定された確率変数の確率密度関数の形状及び故障確率を表示する表示部を有することを特徴とする故障確率評価装置。
10. 請求項９に記載の故障確率評価装置において、
    前記表示部には、前記接続部により推定された確率密度関数の極値モデルのパラメータと再現レベルを表示することを特徴とする故障確率評価装置。
11. 機械システムの故障確率を演算する故障モデルに用いられる確率変数の発生頻度分布の裾部の確率密度関数を極値統計モデルに基づいて推定するステップと、
    前記発生頻度分布の裾部以外の確率密度関数を推定するステップと、
    前記裾部の確率密度関数と前記裾部以外の確率密度関数を用いて、前記発生頻度分布の全体の確率密度関数を推定するステップと、
    前記全体の確率密度関数と前記故障モデルに基づき、前記機械システムの故障確率を演算するステップと、
    を有することを特徴とする故障確率評価方法。

Images