【文献】
大森 敏明,使える!統計検定・機械学習−V,システム/制御/情報,システム制御情報学会,2015年04月15日, 第59巻 第4号,33〜38
(58)【調査した分野】(Int.Cl.,DB名)
前記多変数多項式回帰において、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて前記説明変数の選択を行うことを特徴とする、請求項2に記載の解析パラメータの推定方法。
前記ベイズ推定において、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法を用いることを特徴とする、請求項1から3のいずれか1項に記載の解析パラメータの推定方法。
【背景技術】
【0002】
近年の計算機技術の進展により、物理現象を模擬するための数値解析あるいは数値シミュレーションが、様々な分野で広く用いられるようになっている。工業的には、ある装置や機械が、その目的に適った性能を発揮できるかを事前に評価することを目的として、装置や機械の設計段階において事前に数値シミュレーションを実施し、その結果に基づいて装置や機械の設計を最適化することが広く行われている。
【0003】
こうした工業的な装置や機械における物理現象を対象とした数値シミュレーションを精度よく実行するためには、解析プログラムに与える初期条件や境界条件、ないしは解析プログラムの中で物理現象を記述するために用いられる物理モデルが含有するモデルパラメータなどの解析パラメータを適切に与える必要がある。
【0004】
例えば、ある燃焼装置内の温度分布を評価するために、ナビエ・ストークス方程式や熱伝導方程式、燃焼反応のモデル式などに基づいた燃焼・流体シミュレーションを行うことができるが、その場合、境界条件としての燃焼装置入口での燃料ガスの流速および化学組成の分布や、燃焼反応のモデル式が含有するモデルパラメータ、例えば反応速度定数などの解析パラメータを適切に与えない限り、実用に足る解析精度を得ることはできない。
【0005】
しかしながら、一般に実際の工業的な装置や機械を対象とする場合、これらの解析パラメータを事前に全て正確に知ることは困難であることが多い。その場合、事前に正確に知ることができない解析パラメータについては、何らかの仮定値を与えて数値シミュレーションを実行することになる。そのため、数値シミュレーション自体を実行することは可能であっても、得られた結果が目的に適う十分な精度を有しているとは限らない。
【0006】
すなわち、装置や機器の設計に足る十分な精度において物理現象を模擬するための数値シミュレーションを実行するためには、未知の解析パラメータの値を的確に推定する必要がある。
【0007】
上記の課題に対しては、逆問題(inverse problem)あるいは逆解析(inverse analysis)の有用性が近年注目されている。逆問題あるいは逆解析とは、入力から出力を、または原因から結果を直接的に求める順問題(direct problem)あるいは順解析(direct analysis)とは逆の概念であり、出力から入力を、または結果から原因を求めようとするものであり、その概要については例えば非特許文献1に記述されている。
【0008】
前記の燃焼装置内の温度分布を評価する問題についても、こうした逆解析手法を適用することにより、出力である燃焼装置の出口境界での燃焼ガスの温度および化学組成などの測定結果から、入力である燃焼装置の入口境界での燃料ガスの流速および化学組成の分布や、燃焼反応のモデルパラメータなどを求めることができると期待される。
【0009】
実際、逆解析手法は様々な工業分野への応用が行われている。例えば非特許文献2では、自動車用ディスクブレーキに用いられるブレーキパッド摩擦材に関して、直接的に計測することが困難である摩擦材の弾性定数を、逆解析手法を適用することにより同定する手法について記述されている。具体的には、まず打撃試験によりブレーキパッドの固有振動数とモード形状を測定した上で、弾性定数を適当に仮定した数値シミュレーション(有限
要素解析)を多数実行して、計算結果を入力(推定対象パラメータ)である弾性定数と出力(計測可能な観測量)であるブレーキパッドの固有振動数およびモード形状の関係性として定式化し、最後に最尤推定法を用いて測定データと計算結果の誤差が最小となるよう弾性定数を同定する手法について記述されている。
【0010】
また非特許文献3では、逆解析手法によるパラメータ推定問題における実験方法あるいは計測条件、例えば計測点の位置や計測する物理量などを最適化する手法について記述されている。具体的には、事後推定誤差の共分散行列の最大固有値に着目し、その値を最小化する実験方法あるいは計測条件を採用することで、効果的なパラメータ推定が実現できることが記述されている。
【発明の概要】
【発明が解決しようとする課題】
【0012】
しかしながら、一般に実際の工業的な装置や機械を対象とする場合、これらの解析パラメータを事前に全て正確に知ることは困難であることが多い。例えば前記の例では、装置や機械の形状によって、境界条件としての流路入口の流速および圧力分布を正確に測定することが物理的に困難な場合もあるし、また物性の良くわかっていない新たな流体などを対象とする場合のように、密度や粘性係数などのパラメータを事前に把握することが困難な場合もある。
【0013】
こうした場合、事前に正確に知ることができない解析パラメータについては、何らかの仮定値を与えて数値シミュレーションを実行することになる。そのため、数値シミュレーション自体を実行することは可能であっても、得られた結果が目的に適う十分な精度を有しているとは限らない、という問題があった。すなわち、装置や機器の設計に足る十分な精度において物理現象を模擬するための数値シミュレーションを実行するためには、未知の解析パラメータの値を的確に推定する必要がある、という課題があった。
【0014】
これらの非特許文献に記述されている概念や手法は、様々な工業的な装置や機械を対象とした数値シミュレーションにおいて、目的に適う十分な精度が得られるように、未知の解析パラメータの値を的確に推定するために適用することが可能と考えられるが、実用上はいくつかの課題がある。
【0015】
まず、前記の非特許文献2においては、有限要素解析の計算コスト削減のために、入力(推定対象パラメータ)と出力(計測可能な観測量)の関係性を、予め実行しておいた数値シミュレーションの結果に基づいて、線形あるいは非線形の多項式を用いて応答曲面として定式化しておき、それを実際の観測結果と照らし合わせることで、未知の(推定対象の)パラメータの同定が行われている。
【0016】
ここで、非特許文献2にも記載されているように、応答曲面の精度はパラメータの同定
精度に大きく影響するため、問題の性質に適した定式化の手法を選択することが重要となる。具体的には、問題における入力と出力の関係性を適切に表現できるよう、例えば応答曲面の多項式の次数を適切に選択することが必要となる。
【0017】
しかしながら、推定対象であるパラメータの数(次元)が大きい場合、応答曲面の多項式の次数を増やしていくと、応答曲面を構成する多項式の係数の数が爆発的に増大してしまう。例えば、推定対象であるパラメータの数が20個ある場合、それらのパラメータと出力(観測量)の関係性を多項式で表すとすると、1次多項式の場合は係数の数は21個となるが、2次多項式の場合は231個、3次多項式の場合は1581個となる。
【0018】
応答曲面の定式化にあたって、単純にこれらの多項式の係数を数値シミュレーションの計算結果に基づいて、例えば最小二乗回帰などの手法により決定するためには、係数の数に応じた多数の数値シミュレーションを、入力条件を変化させた上で実行する必要がある。すなわち、1次多項式を採用した場合には最低21回の数値シミュレーションを実行すれば係数を決定することができるが、2次多項式を採用した場合は最低でも231回、3次多項式を採用した場合には最低でも1581回以上の数値シミュレーションを実行しない限り、応答曲面の係数を決定することができない。
【0019】
また、非特許文献3においては、こうした逆解析手法によるパラメータ推定問題において、パラメータの推定精度を高めるために観測量(計測点の位置や計測する物理量など)を最適化する手法について記述されている。しかしながら、「物理現象を模擬するための数値シミュレーションを所望の解析精度において実行する」という本来の目的に対しては、パラメータの推定精度を過度に高めることを志向することは必ずしも現実的ではない。つまり、パラメータの推定精度を高めるためには、より高い精度で前述の応答曲面を作成する必要が生じるため、その作成のために必要な数値シミュレーションの実行回数も増やす必要がある。
【0020】
したがって、現実的に実行可能な数値シミュレーションの実行回数のもとで、できるだけ高次の非線形性を表現できる応答曲面の定式化手法と、得られた応答曲面が、目的とする数値シミュレーションにおいて所望の解析精度を得る上で十分な精度を有しているか否かを適切に評価できる、応答曲面および観測量の最適化手法とが必要となるが、そのような手法はこれまでに実用化されていなかった。
【0021】
そのため、数値シミュレーションの計算コストが大きく、かつ問題を支配する物理現象の非線形性が強い問題、例えば乱流解析や燃焼解析などの問題に対して、このような応答曲面法に基づく逆解析手法を適用して、目的とする数値シミュレーションを所望の解析精度において実行するために最低限必要な精度において、未知の解析パラメータの値を的確に推定することは、これまでは実用的に困難であった。
【0022】
そこで本発明では、物理現象を模擬するための数値シミュレーションを所望の解析精度において実行するために必要な精度において、未知の解析パラメータの値を的確に推定できるとともに、数値シミュレーションの計算コストが大きく、かつ問題を支配する物理現象の非線形性が強い問題に対しても適用可能な、解析パラメータの推定方法を提供することを課題とする。
【課題を解決するための手段】
【0023】
上記課題を解決するため、本発明の第1の態様である解析パラメータの推定方法は、物理現象を模擬するための数値シミュレーションを所望の解析精度において実行するために必要な、1つ以上の未知の解析パラメータX
i(i=1,2,・・・,N)の推定値X
*iを推定する手法であって、対象とする物理現象における、1つ以上の観測点における1
つ以上の観測変数Y
j(j=1,2,・・・,M)について、前記1つ以上の解析パラメータX
iの値を変化させた多数の数値シミュレーションを実行して得られた結果に基づいて、前記X
iと前記Y
jとの関係性を定式化し、さらに前記物理現象について実験的計測を行うことにより、前記Y
jの観測値Y
*j(j=1,2,・・・,M)を得て、前記Y
*jの統計情報と、前記定式化された前記X
iと前記Y
jの関係性に基づいて、ベイズ推定を行うことにより、前記観測値Y
*jが得られたもとでの、前記パラメータX
iの事後確率分布P(X
*i)=P(X
i|Y
*1,Y
*2,・・・
,Y
*M)として、前記パラメータX
iの推定値X
*iの確率分布を得ることを特徴とする。
【0024】
このように構成すると、解析パラメータX
iの値を変化させた多数の数値シミュレーションを実行して得られた結果に基づいて定式化された、解析パラメータX
iと観測変数Y
jの関係性に基づいて、実際に得られた観測値Y
*jの統計情報を適切に説明できる解析パラメータX
iの推定値X
*iの確率分布P(X
*i)を得ることができる。
【0025】
また、本発明の第2の態様である解析パラメータの推定方法では、前記X
iと前記Y
jとの関係性は、Y
j(j=1,2,・・・,M)を目的変数とし、X
i(i=1,2,・・・,N)を説明変数とした多変数多項式回帰により定式化されることを特徴とする。
【0026】
このように構成すると、解析パラメータX
iと観測変数Y
jの関係性を、その関係性に応じて適切な次数および/または項数の多変数多項式を選択することにより、精度よく表現することができる。
【0027】
また、本発明の第3の態様である解析パラメータの推定方法は、前記多変数多項式回帰において、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて変数選択を行うことを特徴とする。
【0028】
このように構成すると、前記多変数多項式を構成する変数および/または項の候補となる数が多い場合においても、解析パラメータX
iと観測変数Y
jの関係性を、変数の数および/または項の数を最小限に抑制しつつ、かつ精度よく表現することができるため、計算負荷を抑制しつつ所望の解析精度を得ることができる。
【0029】
また、本発明の第4の態様である解析パラメータの推定方法は、前記ベイズ推定において、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法を用いることを特徴とする。
【0030】
このように構成すると、前記多変数多項式における変数の数および/または項の数が多い場合や、あるいは解析パラメータX
iと観測変数Y
jの関係性が強い非線形性を有する場合、さらには観測値Y
*jの統計分布が正規分布に従わないような場合など、観測値Y
*jの統計情報から直接的に前記パラメータX
iの推定値X
*iを推定することが困難である場合であっても、観測値Y
*jの統計情報に基づいて、前記パラメータX
iの事後確率分布P(X
*i)=P(X
i|Y
*1,Y
*2,・・・
,Y
*M)を数値的に求めることが可能となる。
【0031】
また、本発明の第5の態様である解析パラメータの推定方法は、前記1つ以上の観測点は、前記パラメータX
iと相関の高い観測変数Y
j1(j
1=1,2,・・・,M
1)が得られるように選ばれた観測点と、前記数値シミュレーションにおいて、計算結果の評価に用いられる評価変数Y
j2(j
2=M
1+1,M
1+2,・・・,M)が得られるように選ばれた観測点とからなることを特徴とする。
【0032】
このように構成すると、前記未知の解析パラメータX
iを精度よく推定するために選ば
れた観測変数Y
j1と、前記数値シミュレーションにおける計算結果の精度を評価するために選ばれた評価変数Y
j2に基づいて、前記解析パラメータX
iの推定値X
*iの確率分布P(X
*i)を効率的に得ることができる。
【0033】
さらに、本発明の第6の態様である解析パラメータの推定方法は、前記事後確率分布P(X
i|Y
*1,Y
*2,・・・
,Y
*M)に基づく、前記パラメータX
iの推定値X
*iの事後確率分布P(X
*i)を、前記定式化された前記X
iと前記Y
jとの関係性に入力することにより得られる、前記評価変数Y
j2の事後確率分布P(Y
**j2)=P(Y
j2|X
*i)=P(Y
j2|Y
*1,Y
*2,・・・
,Y
*M)について、前記事後確率分布P(Y
**j2)の信用区間が、いずれかのj
2について所定の閾値を上回る場合に、前記観測変数の数M
1、M
2および/ないしMを増やすことにより、前記事後確率分布P(Y
**j2)の信用区間が、いずれのj
2についても所定の閾値を下回るように構成したことを特徴とする。
【0034】
このように構成すると、前記評価変数Y
j2の事後確率分布P(Y
**j2)の信用区間が所定の閾値を下回るように、すなわち目的とする数値シミュレーションの実行結果において前記評価変数Y
j2を所望の解析精度において得るために必要な精度において前記の解析パラメータX
iの推定値X
*iの確率分布P(X
*i)を得ることができる。
【発明を実施するための形態】
【0036】
以下、図面を参照して、本発明の好適な実施形態について詳細に説明する。
【0037】
図1は、本発明の一実施形態による未知の解析パラメータX
i(i=1,2,・・・,N)の推定方法のステップを示している。
【0038】
はじめに、ステップ1として、解析対象に観測点および評価点を設定し、観測点における観測変数Y
j1(j
1=1,2,・・・,M
1)と、評価点における評価変数Y
j2(j
2=M
1+1,M
1+2,・・・,M)を定義する。解析対象は、例えば、ごみ焼却プラントやバイオマス燃焼プラントなどを含むが、これらに限られず、装置や機器の設計の過程において数値シミュレーションが用いられるあらゆる装置や機器が、解析対象となりうる。観測点の数と観測変数の数との関係は任意である。例えば、M
1個の観測点に各1つずつ合計M
1個の観測変数を定義してもよいし、1個の観測点にM
1個の観測変数を定義してもよいし、5個の観測点に各M
1/5個ずつ合計M
1個の観測変数を定義してもよい。
【0039】
後述するように、本実施形態における未知の解析パラメータX
i(i=1,2,・・・,N)の推定は、観測点において得られる観測変数Y
j(j=1,2,・・・,M)の統計情報に基づいてなされるため、原理的には観測変数Y
jの数Mが大きいほど、精度よく解析パラメータX
iを推定することができるが、一方で解析パラメータX
iを推定するための計算量が増大するという問題が生じる。
【0040】
そのため本実施形態においては、観測変数Y
jの数Mをできるだけ少なくした上で精度よく解析パラメータX
iを推定するために、以下に記すように、観測変数Y
j(j=1,2,・・・,M)を、観測変数Y
j1(j
1=1,2,・・・,M
1)と、評価変数Y
j2(j
2=M
1+1,M
1+2,・・・,M)の2つに分けて定義する。
【0041】
観測変数Y
j1(j
1=1,2,・・・,M
1)については、対象とする物理現象を考慮して、推定すべき未知の解析パラメータX
i(i=1,2,・・・,N)との相関が高い観測量として定義する。観測変数Y
j1の数M
1を適切に決定する方法については、後述する。
【0042】
評価変数Y
j2(j
2=M
1+1,M
1+2,・・・,M)については、対象とする数値シミュレーションにおいて、計算結果の評価を行う上で着目する変数、すなわち計算結果の精度を評価する対象となる変数として定義する。評価変数Y
j2の数M
2は、対象とする数値シミュレーションにおいて計算結果の評価を行う上で必要となる数とすればよい。
【0043】
例えば前記の燃焼装置内の温度分布を評価する問題において、燃焼装置の入口境界での燃料ガスの流速および化学組成の分布が未知の解析パラメータである場合、観測変数Y
j1としては燃焼装置の入口境界にできるだけ近く、物理的に観測が可能な位置における装置内のガス流速および/または化学組成の計測値を採用することができる。また、評価変数Y
j2としては、装置内の代表的な箇所であり、かつ物理的に観測が可能である箇所における温度の計測値や、装置出口境界でのガス流速および/または化学組成の計測値を採用することができる。
【0044】
次に、ステップ2として、解析パラメータX
i(i=1,2,・・・,N)を変化させた多数の数値シミュレーションを実行する。以下では、ステップ2において、解析パラメータX
i(i=1,2,・・・,N)をK通りに変化させた数値シミュレーションを実行することとし、k番目の計算条件に対応する解析パラメータX
iの値をX
ik(k=1,2,・・・,K)と記述する。また、k番目の計算条件での数値シミュレーション結果における観測変数Y
j1および評価変数Y
j2の値を、それぞれY
j1kおよびY
j2kと記述する。計算条件の数Kを適切に決定する方法については、後述する。
【0045】
次に、ステップ3として、ステップ2で実行した数値シミュレーションの実行結果に基づいて、解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性を定式化する。
【0046】
本実施形態においては、解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性は、解析パラメータX
iについての多変数多項式回帰により定式化される。ここで、多変数多項式回帰式の次数については、模擬対象となる物理現象の複雑さや、解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性の複雑さに応じて決めるのがよい。
【0047】
例えば、多変数多項式回帰式の次数を2とした場合において、解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性は次式のように定式化することができ、数値シミュレーションの計算条件であるX
ikと、その実行結果であるY
j1kおよびY
j2kに基づいて、解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性を最もうまく説明できる次式の係数a、b、cおよびdを、最小二乗回帰や最尤推定法などの手法によって求める回帰問題に帰着する。多変数多項式回帰式の次数を1とした場合や、3以上とした場合についても、同様に定
式化することができる。
【0049】
しかしながら、解析パラメータX
iの数が大きいと、回帰係数a、b、cおよびdの数が増大する。あるいは、多変数多項式の次数が大きくなった場合にも、求めるべき回帰係数の数が増大する。その場合、後述するステップ6において、観測変数Y
j1および評価変数Y
j2の観測値Y
*j1およびY
*j2の統計情報に基づいて、解析パラメータX
iの推定値X
*iの確率分布を得る処理の計算負荷が増大してしまうという問題が生じる。
【0050】
この問題に対処するため、本実施形態においては、ステップ3における解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性の定式化に際して、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて説明変数の選択が行われる。
【0051】
スパース推定手法とは、回帰モデルの損失関数に回帰係数のノルムに基づく正則化項を加えた正則化損失関数を最小化することにより、変数の数が多い場合において、推定の安定化とともに、真に意味のある変数を選択することが可能となる手法であり、リッジ回帰、ラッソ回帰またはエラスティックネット等の手法が広く用いられている。これらのスパース推定手法については公知であるのでその詳細は省略する。
【0052】
次に、ステップ4として、スパース推定手法による多変数多項式回帰により定式化された解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の回帰式の精度を評価する。具体的には、ステップ2において実行した数値シミュレーションの全ての計算条件X
ikについて、ステップ3で得られた回帰式に各X
ikを代入することにより得られる予測値Y
〜j1kおよびY
〜j2kの値と、シミュレーション結果として得られたY
j1kおよびY
j2kの値とを比較する。
【0053】
回帰式の精度の評価指標としては、例えば上記のY
j1kおよびY
j2kの値と、対応するY
〜j1kおよびY
〜j2kの値との二乗平均誤差を用いることができる。二乗平均誤差の値が所定の閾値よりも大きい場合は、シミュレーションの実行ケース数Kを増やすか、または多変量多項式における次数を増やした上で、再度ステップ2からステップ4を実行することにより、回帰式の精度を高めるのがよい。
【0054】
次に、ステップ5として、観測点における観測変数Y
j1の観測値Y
*j1および評価点における評価変数Y
j2の観測値Y
*j2を実験的に得る。
【0055】
次に、ステップ6として、ステップ5で得られた観測値Y
*j1およびY
*j2の統計情報と、ステップ3で定式化された解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性とに基づいて、ベイズ推定を行うことにより、観測値Y
*j1およびY
*j2が得られたもとでの、解析パラメータX
iの事後確率分布P(X
*i)=P(X
i|Y
*1,Y
*2,・・・
,Y
*M)を求める。
【0056】
本ステップ6の目的は、実験的に得られた観測値Y
*j1およびY
*j2の情報から、それらを与える解析パラメータX
*iを逆推定することである。
【0057】
ここで、ステップ5において得られる観測値Y
*j1およびY
*j2は、実際の装置または機器において実験的に計測される量であるため、実験上あるいは計測上の様々な要因を反映した観測誤差を含んでいる。そのため、観測値Y
*j1およびY
*j2は確定値として得られるのではなく、バラつきをもった統計分布として得られることが通常である。
【0058】
観測値Y
*j1およびY
*j2の分布におけるバラつきが小さい場合や、その分布が正規分布に従う場合は、観測値Y
*j1およびY
*j2の各々の分布の平均値ないし中央値をそれらの代表値として確定的に用いることができる。この場合、解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性が線形関係にある場合は、観測値Y
*j1およびY
*j2を与える解析パラメータX
*iを解析的に算出することができる。
【0059】
しかしながら、解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性がより複雑である場合、例えばこれらが2次以上の多項式で表されるような非線形関係にある場合は、前記のように観測値Y
*j1およびY
*j2の各々の分布の平均値ないし中央値をそれらの代表値として確定的に用いることができる場合であっても、それらの観測値Y
*j1およびY
*j2を与える解析パラメータX
*iを解析的に算出することは不可能であることが多い。
【0060】
あるいは、観測値Y
*j1およびY
*j2の分布におけるバラつきが大きい場合や、その分布が正規分布に従わない場合には、そもそも観測値Y
*j1およびY
*j2の代表値を確定的に定義することができないため、たとえ解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性が線形関係にあったとしても、そのような観測値Y
*j1およびY
*j2に対応する解析パラメータX
*iを解析的に算出することは不可能である。
【0061】
そこで本実施形態による解析パラメータの推定方法は、ベイズ推定を用いることにより、観測値Y
*j1およびY
*j2の統計情報に基づいて、それらを与える解析パラメータX
*iの確率分布を、観測値Y
*j1およびY
*j2が得られたもとでの、解析パラメータX
iの事後確率分布P(X
*i)=P(X
i|Y
*1,Y
*2,・・・
,Y
*M)として求めるように構成されている。
【0062】
具体的には、本実施形態においては、ステップ6におけるベイズ推定手法として、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法が用いられている。
【0063】
マルコフ連鎖モンテカルロ法は、ベイズ推定の枠組みにおいて、事後確率分布の計算を数値的に行うことができる手法として広く用いられており、メトロポリス−ヘイスティングアルゴリズム、ギブスサンプラーなどの手法の総称である。またハミルトニアンモンテカルロ法は、マルコフ連鎖モンテカルロ法における目標分布への収束性の改善を目的として近年開発され広く利用されている。これらの手法については公知であるのでその詳細は省略する。
【0064】
これらのベイズ推定手法を用いることにより、前記多変数多項式における変数の数および/または項の数が多い場合や、あるいは解析パラメータX
iと観測変数Y
j1および/または評価変数Y
j2との関係性が強い非線形性を有する場合や、さらには観測変数Y
j1および/または評価変数Y
j2の統計分布が正規分布に従わないような場合など、観測変数Y
j1および評価変数Y
j2の観測値Y
*j1およびY
*j2の統計情報から直接的に前記解析パラメータX
iの推定値X
*iを推定することが困難な場合であっても、観測値Y
*j1およびY
*j2の統計情報に基づいて、解析パラメータX
iの事後確率分布P
(X
*i)=P(X
i|Y
*1,Y
*2,・・・
,Y
*M)を数値的に求めることが可能となる。
【0065】
次に、ステップ7として、ステップ6で得られた解析パラメータX
iの推定値X
*iの事後確率分布P(X
*i)を、ステップ3で定式化された解析パラメータX
iと評価変数Y
j2の関係性に入力することにより、評価変数Y
j2の事後確率分布P(Y
**j2)=P(Y
j2|X
*i)=P(Y
j2|Y
*1,Y
*2,・・・
,Y
*M)を得る。
【0066】
次に、ステップ8として、ステップ7で得られた評価変数Y
j2の事後確率分布P(Y
**j2)の信用区間の評価を行う。
【0067】
前述のように、本実施形態における評価変数Y
j2は、計算結果の評価を行う上で着目する変数、すなわち計算結果の精度を評価する対象となる変数として定義されている。したがって、観測変数Y
j1および評価変数Y
j2の観測値Y
*j1およびY
*j2が得られたもとでの、評価変数Y
j2の事後確率分布P(Y
**j2)の信用区間が、所定の閾値よりも小さくなっていれば、対象とする数値シミュレーションにおいて、所望の解析精度を得るために必要な精度において未知の解析パラメータX
iの推定値X
*iの事後確率分布P(X
*i)が求まったといえることになる。
【0068】
ここで、評価変数Y
j2の事後確率分布P(Y
**j2)の信用区間が、所定の閾値よりも小さくならなかった場合は、対象とする数値シミュレーションにおいて、解析パラメータX
iの推定値X
*iの事後確率分布P(X
*i)を、所望の解析精度を得るために必要な精度で求めるためには、情報が不足していることを意味している。したがって、そのような場合は、観測変数Y
j1の数M
1を増やした上で、再度ステップ1からステップ8までの手順を実行する。これにより、観測変数Y
j1から得られる情報を増やすことができ、その結果として、数値シミュレーションにおいて所望の解析精度を得るために必要な精度で、未知の解析パラメータX
iの推定値X
*iの事後確率分布P(X
*i)を求めることができる。
【実施例】
【0069】
以下、図面を参照して、本発明の実施例について説明する。
【0070】
図2は、ある燃焼装置内の温度分布やガス組成分布を評価するために行う燃焼シミュレーションにおいて、境界条件としての燃焼装置入口境界での燃料ガスの化学組成の分布が未知である場合の一実施例を示している。
【0071】
ここでは、
図2に示されているように4つに分割されている燃焼装置の入口境界の各々の領域F
1、F
2、F
3、F
4における燃料ガスG
1、G
2、G
3、G
4の組成を表現するためのパラメータとして、合計20個の解析パラメータX
i(i=1〜20)が定義されている。すなわち、fをある関数として、G
n=f(X
i)(n=1〜4、i=1〜20)であり、これらX
iの値を、上述した実施形態による解析パラメータの推定方法によって推定するものとする。また、本実施例における数値シミュレーションの要求精度としては、燃焼装置出口のガス温度を±10℃の精度で求めることが要求されているものとする。
【0072】
はじめに、ステップ1として、観測点および観測変数、ならびに評価点および評価変数を設定する。具体的には、燃焼装置の入口境界の各領域F
1、F
2、F
3、F
4の直上であって、かつ燃焼装置内のガス温度およびガス組成の計測が実際に可能な箇所に、4つの観測点R
1、R
2、R
3、R
4を定義し、それらの位置における燃焼装置内のガス温度およびガス組成(未燃ガス濃度、酸素濃度、水分濃度の3つ)を観測変数として定義する。
また、燃焼装置の出口境界に1つの評価点R
5を定義し、その位置におけるガス温度およびガス組成(同じく未燃ガス濃度、酸素濃度、水分濃度の3つ)を評価変数として定義する。
【0073】
つまり、1つの観測点ないし評価点において、4つずつの観測変数ないし評価変数が定義される。以下では本実施例における観測変数Y
j1(j
1=1〜16)および評価変数Y
j2(j
2=17〜20)を、それぞれ次の表1のように定義する。
【0074】
【表1】
【0075】
次に、ステップ2として、解析パラメータX
i(i=1〜20)を変化させた多数の数値シミュレーションを実行する。具体的には、ステップ3で各観測変数Y
j1(j
1=1〜16)および評価変数Y
j2(j
2=17〜20)を各解析パラメータX
iの2次多変数多項式として表現することで応答曲面の定式化を行うことを前提に、2次多変数多項式の係数の数(231個)を上回るように合計250点(k=250)の数値シミュレーションを実行した。
【0076】
次に、ステップ3として、ステップ2で実行した数値シミュレーションの実行結果に基づいて、解析パラメータX
i(i=1〜20)と、観測変数Y
j1(j
1=1〜16)および評価変数Y
j2(j
2=17〜20)との関係性を定式化する。ここでは、2次の多変数多項式を用い、さらにラッソ回帰によるスパース推定手法を適用することで、説明変数の選択を行った。変数選択の結果の詳細は省略するが、もとの多変量多項式の係数の数(231個)に対して、選択された変数の数は各観測変数および評価変数に対して、10個から90個程度となり、効果的に変数選択が行われたことが分かる。
【0077】
次に、ステップ4として、上記により定式化された解析パラメータX
i(i=1〜20)と観測変数Y
j1(j
1=1〜16)および評価変数Y
j2(j
2=17〜20)の回帰式の精度を評価する。具体的には、ステップ2において実行した数値シミュレーションの全ての計算条件X
ik(i=1〜20、k=1〜250)について、ステップ3で得られた回帰式に各X
ikを代入することにより得られる予測値Y
〜j1kおよびY
〜j2kの値と、シミュレーション結果として得られたY
j1kおよびY
j2kの値との二乗平均誤差を比較し、十分な精度が得られていることを確認した。より高い精度が必要である場合には、
図1に示したフローチャートに従って、適宜数値シミュレーションの実行回数な
いし多変量多項式の次数を増やした上で、再度ステップ2からステップ4までの手順を実行すればよい。
【0078】
次に、ステップ5として、観測点における観測変数Y
j1の観測値Y
*j1および評価点における評価変数Y
j2の観測値Y
*j2を実験的に得た。
【0079】
次に、ステップ6として、ステップ5で得られた観測値Y
*j1およびY
*j2の統計情報と、ステップ3で定式化された解析パラメータX
iと観測変数Y
j1、および解析パラメータX
iと評価変数Y
j2の関係性とに基づいて、ベイズ推定を行うことにより、観測値Y
*j1およびY
*j2が得られたもとでの、解析パラメータX
iの事後確率分布P(X
*i)=P(X
i|Y
*j)をベイズ推定手法(ハミルトニアンモンテカルロ法)により求めた。
【0080】
その結果の一例として、解析パラメータX
1からX
10までの事後確率分布P(X
*i)の推定結果を
図3に示す。これらの図に示されるように、実験的に得られる観測値の統計情報に対応する解析パラメータの値の事後確率分布は、一般に幅広い分布を示すことが分かる。
【0081】
次に、ステップ7として、ステップ6で得られた解析パラメータX
iの推定値X
*iの事後確率分布P(X
*i)を、ステップ3で定式化された解析パラメータX
iと評価変数Y
j2の関係性に入力することにより、評価変数Y
j2の事後確率分布P(Y
**j2)=P(Y
j2|X
*i)=P(Y
j2|Y
*j)を得た。その結果の一例として、各観測点および評価点におけるガス温度の事後確率分布P(Y
**1)、P(Y
**5)、P(Y
**9)、P(Y
**13)、P(Y
**17)を
図4に示す。
【0082】
前述のように、本手法により同定された未知の解析パラメータX
i(i=1〜20)の事後確率分布P(X
*i)は、
図3に示すように幅広い分布を示すものの、各観測点および評価点におけるガス温度の事後確率分布は、
図4に示すように観測変数および評価変数の実際の観測値(図中の太線)を含むシャープな分布となっており、それらの実際の観測値Y
*j1およびY
*j2をうまく説明できていることが分かる。
【0083】
次に、ステップ8として、ステップ7で得られた評価変数Y
j2の事後確率分布P(Y
**j2)の信用区間の評価を行う。前述したように本実施例では、評価変数である燃焼装置出口のガス温度を±10℃の精度で求めることが要求されていたが、
図4で示される燃焼装置出口のガス温度Y
17の事後確率分布は、真値(実際の測定値)836℃に対して、約834〜841℃の範囲にシャープに分布しており、十分な精度が得られていることが分かる。より高い精度が必要である場合には、
図1に示したフローチャートに従って、適宜観測変数の数を増やした上で、再度ステップ1からステップ8までの手順を実行すればよい。
【0084】
以上の実施例では、本発明を燃焼装置における未知のガス組成の同定に適用しているが、本発明はこの例に限定されず、その他の様々な工業的装置あるいは機械における未知の解析パラメータの推定、例えば機器や装置における温度分布や濃度分布の同定、境界表面力の同定、電流密度の同定などにも用いることができる。
【0085】
以上、本発明の実施形態および実施例についてその詳細を説明したが、本発明は上記の実施形態および実施例に限定されるものではなく、特許請求の範囲、及び明細書と図面に記載された技術的思想を逸脱しない範囲において種々変更又は修正を行って実施することが可能である。