(57)【要約】      （修正有）【課題】デジタルコンピュータを用いたモンテカルロシミュレーションに適した手法でシミュレーションを行うシミュレーション装置、シミュレーション方法、及びコンピュータプログラムを提供する。【解決手段】ｎ個のパラメータを表すベクトルｘに基づいて、着目する事象が生じる場合と生じない場合とで、０又は１の値をとる状態指標関数Ｉ（ｘ）、前記ベクトルｘに基づき、着目する事象の生じやすさを示す確率密度関数ｆＸ（ｘ）が関数規定部１４０により規定される。定積分変換部１１０は、状態指標関数、確率密度関数を定積分の形式に変換する。乱数発生部１５０からの乱数を用いてモンテカルロシミュレーションを行い、注目事象の発生確率を予測する。【選択図】図１
2011215686.tif 000002
【特許請求の範囲】
【請求項１】
ｎ個のパラメータを表すベクトルｘに基づいて、着目する事象が生じる場合と生じない場合とで、０又は１の値をとる状態指標関数Ｉ（ｘ）、前記ベクトルｘに基づき、着目する事象の生じやすさを示す確率密度関数ｆＸ（ｘ）に従って、前記事象が生じる確率Ｐｆを求めるシミュレーション装置であって、
確率Ｐｆを、下記（数１）に示される多重定積分形式を用いて算出する
ことを特徴とするシミュレーション装置。
【数１】
2011215686.tif 000009
【請求項２】
前記（数１）を、下記（数２）に変換して、確率Ｐｆの推定値を求める
ことを特徴とする請求項１に記載のシミュレーション装置。
【数２】
2011215686.tif 000010
【請求項３】
積分の上下限値ｂｊ及びａｊを、各基本確率変数Ｘｊの平均値と標準偏差から求めてシミュレーションを行う
ことを特徴とする請求項１又は２に記載のシミュレーション装置。
【請求項４】
前記（数２）に示されるξが一様乱数である
ことを特徴とする請求項２又は３に記載のシミュレーション装置。
【請求項５】
ｎ個のパラメータを表すベクトルｘに基づいて、着目する事象が生じる場合と生じない場合とで、０又は１の値をとる状態指標関数Ｉ（ｘ）、前記ベクトルｘに基づき、着目する事象の生じやすさを示す確率密度関数ｆＸ（ｘ）に従って、前記事象が生じる確率Ｐｆを求めるシミュレーション装置であって、
確率Ｐｆを、下記（数３）に示される多重定積分形式を用いて算出する
ことを特徴とするシミュレーション方法。
【数３】
2011215686.tif 000011
【請求項６】
ｎ個のパラメータを表すベクトルｘに基づいて、着目する事象が生じる場合と生じない場合とで、０又は１の値をとる状態指標関数Ｉ（ｘ）、前記ベクトルｘに基づき、着目する事象の生じやすさを示す確率密度関数ｆＸ（ｘ）に従って、前記事象が生じる確率Ｐｆを求めるコンピュータプログラムであって、
確率Ｐｆを、下記（数４）に示される多重定積分形式を用いて算出する
ことを特徴とするコンピュータプログラム。
【数４】
2011215686.tif 000012
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション装置、シミュレーション方法、及びシミュレーションを行うコンピュータプログラムに関する。
【背景技術】
【０００２】
モンテカルロ法（以下、「ＭＣ法」という。）は、一般的には乱数を用いたシミュレーション法として、極めて広範囲に適用が検討されている。例えば、実際に実験を行うことが容易でないような大規模なものが対象であるような場合（例えば構造物の信頼性予測）の他、目的とするものを見出すために膨大な実験量を必要とすることがある場合（例えば化学反応のコンピュータシミュレーション）、最近では、コンピュータグラフィックなど、その適用可能分野は極めて広い。
【０００３】
本願発明者は、ＭＣ法を、領域積分型ＭＣ法と事象再現型ＭＣ法とに大別した（例えば非特許文献１、非特許文献２）。領域積分型ＭＣ法は効率化が比較的容易であり、その効率化の方法は、Importance Sampling 法、Latin Hypercube Sampling 法等の分散逓減法などが広く研究されている。
【０００４】
しかし、多くのパラメータが複雑に関係するような問題や解析解が得られないような問題では、積分領域の記述に充分な検討が必要な場合がある。
【０００５】
また、構造物システムの信頼性評価を行うような場合には、バイリニア型の応力・歪関係を仮定して領域積分型ＭＣ法を考える場合が多いが、実際には、それ以外の応力・歪関係（あるいは復元力特性）を仮定しなければならない場合（座屈問題、建物の地震応答解析問題、等）があり、領域積分型ＭＣ法の適用は困難である。
【０００６】
このような場合、注目事象を正確に再現し、所要量を注目観測量の平均値として求める方法が便利である。このような方法は応用数学上、Hit-or-miss ＭＣ法と呼ばれることがあるが、Hit-or-miss ＭＣ法は、領域積分型ＭＣ計算にも使われるので、領域積分型ＭＣ計算と区別して、事象再現型ＭＣ法と呼ぶ。事象再現型ＭＣ法は、殆ど全ての場合、安定的に解を得ることができ、広範囲の問題の解決法として有効である。しかし、事象再現型ＭＣ法は効率が悪く、一部の研究を除いて、その効率化法は余り研究されていない。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】S.A. Matsuho, “Efficient hit-or-miss Monte Carlo methods simulating relevantevents for structural reliability evaluation”, Journal of the Society of Materials Science, Japan, Vol.54, No.3, pp.314-319 (2005).
【非特許文献２】S. A. Matsuho, “Basic study on structuralreliability evaluation based on Hit-or-miss Monte-Carlo method”, Journal ofConstructional Steel, Vol.16, No.1, pp.147-152 (2008).
【発明の概要】
【発明が解決しようとする課題】
【０００８】
本発明は、以上のような状況に鑑みてなされたものであって、事象再現型ＭＣ法の効率化を図ることができるシミュレーション装置、シミュレーション方法、及びコンピュータプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００９】
上記の問題点を解決するために、本発明に係るシミュレーション装置は、ｎ個のパラメータを表すベクトルｘに基づいて、着目する事象が生じる場合と生じない場合とで、０又は１の値をとる状態指標関数Ｉ（ｘ）、前記ベクトルｘに基づき、着目する事象の生じやすさを示す確率密度関数ｆＸ（ｘ）に従って、前記事象が生じる確率Ｐｆを求めるシミュレーション装置であって、確率Ｐｆを、下記（数３）に示される多重定積分形式を用いて算出することを特徴としている。
【発明の効果】
【００１０】
本発明に係るシミュレーション装置等によれば、特に事象再現型ＭＣ法において、限界状態関数の定式化が不要となり、デジタルコンピュータによるシミュレーションに好適であるという効果を奏する。
【図面の簡単な説明】
【００１１】
【図１】本実施の形態におけるシミュレーション装置の構成の一例について説明するためのブロック図である。
【図２】数値計算に用いた下路橋タイプのワーレントラス橋について説明するための模式図である。
【図３】図２のトラス橋において、各部材の断面積、及び細長比ｌ／ｒを示す図である。
【図４】事象再現型ＭＣ法に基づいて、破損確率Ｐｆの計算結果を、試行回数とともに示すテーブル２、定積分変換法の式にCrude  ＭＣ法を適用する手法で信頼性評価を行った結果を示すテーブル３を示す図である。
【図５】層別化セルのサンプリングを行わない単純な層別サンプリングの適用を試みた結果（テーブル４）、参考までに、事象再現型ＭＣ法に対する擬似乱数割当（Assigning−Pseudorandom−Numbers）法の適用を試みた結果（テーブル５）を示す図である。
【図６】定積分形式に、準ＭＣ法の一種である準乱数を適用した場合について検討した結果（テーブル６）を示す図である。
【発明を実施するための形態】
【００１２】
以下、本発明の実施の形態について、図面を参照しながら説明する。
図１は、本実施の形態におけるシミュレーション装置の構成の一例について説明するためのブロック図である。
【００１３】
シミュレーション装置１００は、定積分変換部１１０、確率推定部１２０、推定結果出力部１３０を含む。本発明のシミュレーション装置は、デジタルコンピュータに、以下に説明するような処理を行うコンピュータプログラムをインストールして実現することもできる。乱数発生部１５０は、ＭＣ法シミュレーションのための乱数を発生し、確率推定部１２０に出力する。乱数発生部１５０をコンピュータ内に設けても良いことは勿論である。なお、乱数の種類に特に限定はなく、例えばメルセンヌ・ツイスタのような一様乱数の他、必ずしも「乱数」に当たらないような準乱数、低くい違い量列等を用いるような実施形態も考えられる。
【００１４】
関数規定部１４０は、後に詳細に説明する状態指標関数Ｉ（ｘ）、確率密度関数ｆＸ（ｘ）を規定し、定積分変換部１１０に送る。状態指標関数、確率密度関数は、予めコンピュータに入力し、ハードディスク等の記憶手段に格納しておくことができるが、シミュレーションを行うたびに入力するような構成も可能である。また、シミュレーションの結果として得られた状態指標関数Ｉ（ｘ）、確率密度関数ｆＸ（ｘ）の情報を用いる構成も可能である。定積分変換部１１０は、状態指標関数、確率密度関数を定積分の形式に変換する。
【００１５】
事象再現型ＭＣ法では、基本確率変数Ｘ＝（Ｘ１，Ｘ２，・・・，Ｘｎ）Ｔの確率密度関数ｆＸ（ｘ）に従って発生させたｎ個の乱数η（ｉ）＝（η１（ｉ），η２（ｉ），・・・，ηｎ（ｉ））Ｔを用いて、１回の試行を行い、これをＮ回繰り返して、所要量を平均値として推定する。本願発明者の分類する事象再現型ＭＣ法では、ｎ個の乱数η（ｉ）は一様乱数ではない。
【００１６】
後述する構造信頼性評価の場合には、各試行で、対象となるシステムが安全か否かを観測し、安全の場合には０、破損可能性有りの場合には１の値をとる状態指標関数Ｉ（ｘ）を用いて定式化する。即ち、破損確率Ｐｆの推定値Ｐ＾ｆは、下記（数１）で与えられる。
【００１７】
【数１】
2011215686.tif 000003
【００１８】
限界状態関数ｇ（ｘ）は、危険の場合ｇ（ｘ）≦０、安全の場合ｇ（ｘ）＞０となるように定義されるので、上記（数１）の状態指標関数Ｉ（ｘ）の値は、限界状態関数ｇ（ｘ）の値から決めることも可能ではあるが、本実施の形態では、着目事象を計算機内で正確に再現し、各試行で着目事象が生起（破損）するか否かを観測するようにした。
【００１９】
一般に、破損確率Ｐｆの推定には、上記（数１）に従い、事象再現型ＭＣ計算を行うか、あるいは、下記（数２）のように限界状態関数ｇ（ｘ）を用いた形式に対して領域積分型ＭＣ計算を行う。
【００２０】
【数２】
2011215686.tif 000004
【００２１】
しかしながら、（数１）を用いる方法では、稀な事象の生起に膨大な試行が必要となって計算効率が悪く、（数２）を用いる方法では、限界状態関数ｇ（ｘ）の関数形、あるいは、その関数値が与えられる必要があり、その定式化の負担を考えるとデジタルコンピュータでのシミュレーションには適さない。
【００２２】
本実施の形態では、上記（数１）が、状態指標関数Ｉ（ｘ）の標本平均であることに鑑み、いわゆる期待値の定義式より、下記（数３）のような定積分の形式に変換することを試みた。具体的には、定積分変換部１１０は、関数規定部１４０から与えられた状態指標関数、確率密度関数を、下記（数３）のような定積分の形式に変換する。
【００２３】
【数３】
2011215686.tif 000005
【００２４】
前記したように、状態指標関数Ｉ（ｘ）は、各試行で、対象となるシステムが安全の場合には０、破損可能性有りの場合には１の値をとる。このように定積分の形式に変換することにより、限界状態関数ｇ（ｘ）の値の計算、定式化を行う必要がなく、デジタルコンピュータを用いたシミュレーションに適したシミュレーション方法を提供することができる。また、ＭＣ法を用いた積分に関し、これまでに種々検討されてきた有効な手法を適用し、さらなる効率化を図ることが可能になるという利点がある。
【００２５】
例えば、Crude（入門的）ＭＣ法で上記（数３）を評価するとすれば、１回の試行に必要なｎ個の基本確率変数に対して、例えば一様乱数ξ（ｉ）＝（ξ１（ｉ），ξ２（ｉ），・・・，ξｎ（ｉ））Ｔを用いて着目事象を再現し、状態指標関数Ｉ（ｘ）の値を観測し、これをＮ回試行することにより、下記（数４）により、破損確率Ｐ＾ｆを推定する。推定結果出力部１３０は、推定結果を出力する。
【００２６】
【数４】
2011215686.tif 000006
【００２７】
（数４）は、（数３）に示した定積分形式をデジタルコンピュータでの計算に適するように離散化したものであり、{Π（ｂｊ−ａｊ）}は、基本確率変数のｎ次元空間での定義域で定まるHypercube（超立体）の体積である。なお、前記したように、（数３）のように定積分の形式に変換した場合、種々の効率化の手法を適用することが可能となるが、本願発明者は、限界状態関数ｇ（ｘ）の関数形を気にしなくてもよい、という本実施の形態の手法の特徴を活かすためには、被積分関数の形状に余り依存しないLatin Hypercube Sampling 法の適用が好適と考えている。
【００２８】
（数値計算の例）
以下、本実施の形態のシミュレーション装置による実際の数値計算の例について説明する。図２は、数値計算に用いた下路橋タイプのワーレントラス橋について説明するための模式図である。このトラス橋は左右対称であり、図３は、各部材の断面積、及び細長比ｌ／ｒ（ｌは部材長、ｒは断面２次半径）を示す図である。なお、幅員は６．０ｍと仮定する。使用鋼材はＳＭ４００である。
【００２９】
信頼性評価においては、作用荷重として死荷重と活荷重とを仮定する。死荷重は確定量とする。また、本道路橋の最も危険な活荷重載荷状態として、上下線それぞれの中央に大型トレーラ１台ずつ（合計２台）が同時に作用する場合を仮定する。大型トレーラ１台当たりの重量は、過去の交通流実態調査結果を参考にして、平均１５４．３５ｋＮ、分散１０９１２．５ｋＮ２の正規分布に従うものとし、衝撃を考慮して作用応力を算定する。
【００３０】
衝撃は、道路橋示方書で規定される衝撃係数を確定値とした。また、部材強度σＲは、各部材間で完全相関を仮定する。引張強度の平均値はσｔａ＝１４０Ｎ／ｍｍ２（許容引張応力度）、変動係数０．０２を有する正規分布に従うものとする。圧縮部材の場合は、
【００３１】
σｃａ＝１４０    （ｌ／ｒ≦１８の時）
σｃａ＝１４０−０．８２（ｌ／ｒ−１８）    （１８＜ｌ／ｒ≦９２の時）
σｃａ＝１２０００００／（６７００＋（ｌ／ｒ）２）    （ｌ／ｒ≧９２の時）
【００３２】
で与えられる許容圧縮応力度σｃａを平均値として持ち、変動係数０．０２を有する正規分布を仮定する。解析を簡単にするため、１部材の作用応力σｓが、その許容応力度σＲを超過する場合（実際に破損していなくても、少なくとも不利益な状態であるとして）その部材の破損と仮定する。この時、本例題は最弱リンクシステムとして考えることができるので、構成部材の一つでも破損状態となれば構造システム全体が破損状態に至るものとして、構造システムの破損確率Ｐｆを評価する。
【００３３】
以上の条件下で、トラス構造物の信頼性評価を行った。計算環境は、Windows（登録商標）  XP  on  Pentium（登録商標）  ＩＶ  ＣＰＵ  ２．８０ＧＨｚで、富士通Fortran＆Ｃ  Academic Package V4.0 のＣ言語を用い、一様乱数にはメルセンヌ・ツイスタを用いた。
【００３４】
図４は、事象再現型ＭＣ法（数１参照）に基づいて、破損確率Ｐｆの計算結果を、試行回数とともに示すテーブル２、定積分変換法の式（数４参照）にCrude  ＭＣ法を適用する手法で信頼性評価を行った結果を示すテーブル３を示す図である。
【００３５】
各テーブル中、Ｐ^ｆは破損確率Ｐｆの推定値、Time は所要計算時間（秒）を示している。計算処理系の時間計測の精度上、１秒未満は０秒と表記したが、例えばテーブル２において、１０万回試行の０秒は、一千万回試行の所要時間から、約０．５８秒と推定することができる。
【００３６】
テーブル中「―――」と表記されているのは、試行回数が少なく、計算結果が得られていないことを示す。さらに、ρは、推定値Ｐ^ｆの変動係数を示し、下記（数５）及び（数６）で表される。
【００３７】
【数５】
2011215686.tif 000007
【００３８】
【数６】
2011215686.tif 000008
【００３９】
テーブル２から、試行回数の増加とともに、推定値Ｐ^ｆの変動係数ρは小さくなり、推定精度が向上していくのがわかる。
【００４０】
次にテーブル３の結果について考察する。テーブル３の試行において、積分の上下限値ｂｊ、ａｊとしては、ｂｊ＝ｍｊ＋５σｊ、ａｊ＝ｍｊ−５σｊとした。もっとも、積分の上下限値は、これに限定されず、用途、目的等により適宜設定することができる。なお、ここで、ｍｊとσｊとは、各基本確率変数Ｘｊの平均値と標準偏差である。テーブル３より、百万回試行以降は、Ｐ^ｆ＝７．４×１０−４（有効桁数２桁）が得られているのがわかる。即ち、本実施の形態の定積分変換法を用いることにより、より少ない試行回数で解に収束していることが分かった。
【００４１】
ただし、変動係数ρについて明確な傾向は確認できなかった。これは、例えば誤差が大きい推定値であっても、破損確率Ｐｆの推定値が大きければρは小さめの値となる（Ｐｆが小さい場合より、Ｐｆが大きい場合の方が推定しやすい）等の性質を有しているためであると考えられる。従って、ρを参照する場合には、単純にρの大小比較のみで精度の良否を判断することは適切ではないと考えられえる。
【００４２】
本実施の形態の定積分変換法では、基礎式（数３）が定積分で表されているので、種々の有効な積分法を応用することができる。上記には、Latin Hypercube Sampling 法の適用を示唆したが、本計算例は２次元問題であり、次元数が少ないので、層別化セルのサンプリングを行わない単純な層別サンプリングの（数３）への適用を試みた。また、参考までに、（数１）の事象再現型ＭＣ法に対する擬似乱数割当（Assigning−Pseudorandom−Numbers）法の適用を試みた。図５は、それらの結果を示すものであり、それぞれの結果が、テーブル４、テーブル５に示されている。
【００４３】
単純な層別サンプリングについては、例えば、津田孝夫著、「モンテカルロ法とシミュレーション」、改訂版、株式会社培風館、昭和６２年１０月１０日、ｐｐ．８４−１０６に説明がなされている。単純な層別サンプリングの（数３）への適用（テーブル４）において、各基本確率変数の定義域（積分の上下限値）は、簡単のため１０分割に等分した。
【００４４】
テーブル４に示されるように、一千万回以上の試行でＰ^ｆ＝７．４５×１０−４（有効桁数３桁）が得られている。即ち、テーブル３より高精度な計算が実現できているのがわかる。もっとも、テーブル３等と比較して所要計算時間は若干長めとなった。
【００４５】
テーブル５の場合、テーブル２の結果と比較して、同一試行回数での所要計算時間が短く、また、より少ない試行回数で解に収束する傾向が見られた。しかし、本実施の形態の定積分変換法（テーブル３、テーブル４）と比較すると、解への収束状況で見劣りすることが確認できた。
【００４６】
最後に、上記（数３）の定積分形式に、準ＭＣ法の一種である準乱数を適用した場合について検討した。図６は、その計算結果（テーブル６）を示す図である。準乱数は擬似乱数ではないので、公式適用上の打ち切り誤差等が伴う場合があるが、所要量の推定に伴う統計的なバラツキがないので、テーブル６にはρを記載していない。テーブル６の計算結果から、百万回以上の試行でＰ^ｆ＝７．４４９×１０−４（有効桁数４桁）が得られているのがわかる。以上に説明したように、本実施の形態の定積分変換方式は、一様乱数のみならず、準乱数を用いることも可能であり、また、必ずしも「乱数」に限られないことも予測可能である。
【００４７】
本発明に係るシミュレーション手法は、構造物の信頼性評価に限らず、他のシミュレーション、コンピュータグラフィック、さらには金融工学など、適用可能な範囲は極めて多岐にわたり、その実用的効果は極めて大である。
【産業上の利用可能性】
【００４８】
本発明は、例えば、デジタルコンピュータによるシミュレーションに適用することができる。
【符号の説明】
【００４９】
１００                  シミュレーション装置
１１０                  定積分変換部
１２０                  確率推定部
１３０                  推定結果出力部
１４０                  関数規定部
１５０                  乱数発生部
【図１】
2011215686.tif 000013
【図２】
2011215686.tif 000014
【図３】
2011215686.tif 000015
【図４】
2011215686.tif 000016
【図５】
2011215686.tif 000017
【図６】
2011215686.tif 000018