JP6945492B2 - 解析パラメータの推定方法 - Google Patents

Description

本発明は、数値シミュレーションにおける解析パラメータの推定方法に関する。

近年の計算機技術の進展により、物理現象を模擬するための数値解析あるいは数値シミュレーションが、様々な分野で広く用いられるようになっている。工業的には、ある装置や機械が、その目的に適った性能を発揮できるかを事前に評価することを目的として、装置や機械の設計段階において事前に数値シミュレーションを実施し、その結果に基づいて装置や機械の設計を最適化することが広く行われている。

こうした工業的な装置や機械における物理現象を対象とした数値シミュレーションを精度よく実行するためには、解析プログラムに与える初期条件や境界条件、ないしは解析プログラムの中で物理現象を記述するために用いられる物理モデルが含有するモデルパラメータなどの解析パラメータを適切に与える必要がある。

例えば、ある燃焼装置内の温度分布を評価するために、ナビエ・ストークス方程式や熱伝導方程式、燃焼反応のモデル式などに基づいた燃焼・流体シミュレーションを行うことができるが、その場合、境界条件としての燃焼装置入口での燃料ガスの流速および化学組成の分布や、燃焼反応のモデル式が含有するモデルパラメータ、例えば反応速度定数などの解析パラメータを適切に与えない限り、実用に足る解析精度を得ることはできない。

しかしながら、一般に実際の工業的な装置や機械を対象とする場合、これらの解析パラメータを事前に全て正確に知ることは困難であることが多い。その場合、事前に正確に知ることができない解析パラメータについては、何らかの仮定値を与えて数値シミュレーションを実行することになる。そのため、数値シミュレーション自体を実行することは可能であっても、得られた結果が目的に適う十分な精度を有しているとは限らない。

すなわち、装置や機器の設計に足る十分な精度において物理現象を模擬するための数値シミュレーションを実行するためには、未知の解析パラメータの値を的確に推定する必要がある。

上記の課題に対しては、逆問題（inverse problem）あるいは逆解析（inverse analysis）の有用性が近年注目されている。逆問題あるいは逆解析とは、入力から出力を、または原因から結果を直接的に求める順問題（direct problem）あるいは順解析（direct analysis）とは逆の概念であり、出力から入力を、または結果から原因を求めようとするものであり、その概要については例えば非特許文献１に記述されている。

前記の燃焼装置内の温度分布を評価する問題についても、こうした逆解析手法を適用することにより、出力である燃焼装置の出口境界での燃焼ガスの温度および化学組成などの測定結果から、入力である燃焼装置の入口境界での燃料ガスの流速および化学組成の分布や、燃焼反応のモデルパラメータなどを求めることができると期待される。

実際、逆解析手法は様々な工業分野への応用が行われている。例えば非特許文献２では、自動車用ディスクブレーキに用いられるブレーキパッド摩擦材に関して、直接的に計測することが困難である摩擦材の弾性定数を、逆解析手法を適用することにより同定する手法について記述されている。具体的には、まず打撃試験によりブレーキパッドの固有振動数とモード形状を測定した上で、弾性定数を適当に仮定した数値シミュレーション（有限
要素解析）を多数実行して、計算結果を入力（推定対象パラメータ）である弾性定数と出力（計測可能な観測量）であるブレーキパッドの固有振動数およびモード形状の関係性として定式化し、最後に最尤推定法を用いて測定データと計算結果の誤差が最小となるよう弾性定数を同定する手法について記述されている。

また非特許文献３では、逆解析手法によるパラメータ推定問題における実験方法あるいは計測条件、例えば計測点の位置や計測する物理量などを最適化する手法について記述されている。具体的には、事後推定誤差の共分散行列の最大固有値に着目し、その値を最小化する実験方法あるいは計測条件を採用することで、効果的なパラメータ推定が実現できることが記述されている。

久保司郎、「逆問題の解析手法」、材料、日本材料学会、１９９２年、第４１巻、第４７０号、ｐ．１５９５−１６０４ 井上正則ほか、「逆解析を用いたディスクブレーキパッド摩擦材の異方性弾性定数の非侵襲同定」、日本機械学会論文集、２０１６年、第８２巻、第８４０号、ｐ．１−１５ Ameyaほか、「Optimization of Measurements For Inverse Problem」、3rd International Conference on Inverse Problems in Engineering、The American Society of Mechanical Engineers、１９９９年

しかしながら、一般に実際の工業的な装置や機械を対象とする場合、これらの解析パラメータを事前に全て正確に知ることは困難であることが多い。例えば前記の例では、装置や機械の形状によって、境界条件としての流路入口の流速および圧力分布を正確に測定することが物理的に困難な場合もあるし、また物性の良くわかっていない新たな流体などを対象とする場合のように、密度や粘性係数などのパラメータを事前に把握することが困難な場合もある。

こうした場合、事前に正確に知ることができない解析パラメータについては、何らかの仮定値を与えて数値シミュレーションを実行することになる。そのため、数値シミュレーション自体を実行することは可能であっても、得られた結果が目的に適う十分な精度を有しているとは限らない、という問題があった。すなわち、装置や機器の設計に足る十分な精度において物理現象を模擬するための数値シミュレーションを実行するためには、未知の解析パラメータの値を的確に推定する必要がある、という課題があった。

これらの非特許文献に記述されている概念や手法は、様々な工業的な装置や機械を対象とした数値シミュレーションにおいて、目的に適う十分な精度が得られるように、未知の解析パラメータの値を的確に推定するために適用することが可能と考えられるが、実用上はいくつかの課題がある。

まず、前記の非特許文献２においては、有限要素解析の計算コスト削減のために、入力（推定対象パラメータ）と出力（計測可能な観測量）の関係性を、予め実行しておいた数値シミュレーションの結果に基づいて、線形あるいは非線形の多項式を用いて応答曲面として定式化しておき、それを実際の観測結果と照らし合わせることで、未知の（推定対象の）パラメータの同定が行われている。

ここで、非特許文献２にも記載されているように、応答曲面の精度はパラメータの同定
精度に大きく影響するため、問題の性質に適した定式化の手法を選択することが重要となる。具体的には、問題における入力と出力の関係性を適切に表現できるよう、例えば応答曲面の多項式の次数を適切に選択することが必要となる。

しかしながら、推定対象であるパラメータの数（次元）が大きい場合、応答曲面の多項式の次数を増やしていくと、応答曲面を構成する多項式の係数の数が爆発的に増大してしまう。例えば、推定対象であるパラメータの数が２０個ある場合、それらのパラメータと出力（観測量）の関係性を多項式で表すとすると、１次多項式の場合は係数の数は２１個となるが、２次多項式の場合は２３１個、３次多項式の場合は１５８１個となる。

応答曲面の定式化にあたって、単純にこれらの多項式の係数を数値シミュレーションの計算結果に基づいて、例えば最小二乗回帰などの手法により決定するためには、係数の数に応じた多数の数値シミュレーションを、入力条件を変化させた上で実行する必要がある。すなわち、１次多項式を採用した場合には最低２１回の数値シミュレーションを実行すれば係数を決定することができるが、２次多項式を採用した場合は最低でも２３１回、３次多項式を採用した場合には最低でも１５８１回以上の数値シミュレーションを実行しない限り、応答曲面の係数を決定することができない。

また、非特許文献３においては、こうした逆解析手法によるパラメータ推定問題において、パラメータの推定精度を高めるために観測量（計測点の位置や計測する物理量など）を最適化する手法について記述されている。しかしながら、「物理現象を模擬するための数値シミュレーションを所望の解析精度において実行する」という本来の目的に対しては、パラメータの推定精度を過度に高めることを志向することは必ずしも現実的ではない。つまり、パラメータの推定精度を高めるためには、より高い精度で前述の応答曲面を作成する必要が生じるため、その作成のために必要な数値シミュレーションの実行回数も増やす必要がある。

したがって、現実的に実行可能な数値シミュレーションの実行回数のもとで、できるだけ高次の非線形性を表現できる応答曲面の定式化手法と、得られた応答曲面が、目的とする数値シミュレーションにおいて所望の解析精度を得る上で十分な精度を有しているか否かを適切に評価できる、応答曲面および観測量の最適化手法とが必要となるが、そのような手法はこれまでに実用化されていなかった。

そのため、数値シミュレーションの計算コストが大きく、かつ問題を支配する物理現象の非線形性が強い問題、例えば乱流解析や燃焼解析などの問題に対して、このような応答曲面法に基づく逆解析手法を適用して、目的とする数値シミュレーションを所望の解析精度において実行するために最低限必要な精度において、未知の解析パラメータの値を的確に推定することは、これまでは実用的に困難であった。

そこで本発明では、物理現象を模擬するための数値シミュレーションを所望の解析精度において実行するために必要な精度において、未知の解析パラメータの値を的確に推定できるとともに、数値シミュレーションの計算コストが大きく、かつ問題を支配する物理現象の非線形性が強い問題に対しても適用可能な、解析パラメータの推定方法を提供することを課題とする。

上記課題を解決するため、本発明の第１の態様である解析パラメータの推定方法は、物理現象を模擬するための数値シミュレーションを所望の解析精度において実行するために必要な、１つ以上の未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）の推定値Ｘ^＊ _ｉを推定する手法であって、対象とする物理現象における、１つ以上の観測点における１
つ以上の観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）について、前記１つ以上の解析パラメータＸ_ｉの値を変化させた多数の数値シミュレーションを実行して得られた結果に基づいて、前記Ｘ_ｉと前記Ｙ_ｊとの関係性を定式化し、さらに前記物理現象について実験的計測を行うことにより、前記Ｙ_ｊの観測値Ｙ^＊ _ｊ（ｊ＝１，２，・・・，Ｍ）を得て、前記Ｙ^＊ _ｊの統計情報と、前記定式化された前記Ｘ_ｉと前記Ｙ_ｊの関係性に基づいて、ベイズ推定を行うことにより、前記観測値Ｙ^＊ _ｊが得られたもとでの、前記パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）として、前記パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布を得ることを特徴とする。

このように構成すると、解析パラメータＸ_ｉの値を変化させた多数の数値シミュレーションを実行して得られた結果に基づいて定式化された、解析パラメータＸ_ｉと観測変数Ｙ_ｊの関係性に基づいて、実際に得られた観測値Ｙ^＊ _ｊの統計情報を適切に説明できる解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布Ｐ（Ｘ^＊ _ｉ）を得ることができる。

また、本発明の第２の態様である解析パラメータの推定方法では、前記Ｘ_ｉと前記Ｙ_ｊとの関係性は、Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）を目的変数とし、Ｘ_ｉ（ｉ＝１，２，・・・，Ｎ）を説明変数とした多変数多項式回帰により定式化されることを特徴とする。

このように構成すると、解析パラメータＸ_ｉと観測変数Ｙ_ｊの関係性を、その関係性に応じて適切な次数および／または項数の多変数多項式を選択することにより、精度よく表現することができる。

また、本発明の第３の態様である解析パラメータの推定方法は、前記多変数多項式回帰において、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて変数選択を行うことを特徴とする。

このように構成すると、前記多変数多項式を構成する変数および／または項の候補となる数が多い場合においても、解析パラメータＸ_ｉと観測変数Ｙ_ｊの関係性を、変数の数および／または項の数を最小限に抑制しつつ、かつ精度よく表現することができるため、計算負荷を抑制しつつ所望の解析精度を得ることができる。

また、本発明の第４の態様である解析パラメータの推定方法は、前記ベイズ推定において、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法を用いることを特徴とする。

このように構成すると、前記多変数多項式における変数の数および／または項の数が多い場合や、あるいは解析パラメータＸ_ｉと観測変数Ｙ_ｊの関係性が強い非線形性を有する場合、さらには観測値Ｙ^＊ _ｊの統計分布が正規分布に従わないような場合など、観測値Ｙ^＊ _ｊの統計情報から直接的に前記パラメータＸ_ｉの推定値Ｘ^＊ _ｉを推定することが困難である場合であっても、観測値Ｙ^＊ _ｊの統計情報に基づいて、前記パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を数値的に求めることが可能となる。

また、本発明の第５の態様である解析パラメータの推定方法は、前記１つ以上の観測点は、前記パラメータＸ_ｉと相関の高い観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）が得られるように選ばれた観測点と、前記数値シミュレーションにおいて、計算結果の評価に用いられる評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）が得られるように選ばれた観測点とからなることを特徴とする。

このように構成すると、前記未知の解析パラメータＸ_ｉを精度よく推定するために選ば
れた観測変数Ｙ_ｊ１と、前記数値シミュレーションにおける計算結果の精度を評価するために選ばれた評価変数Ｙ_ｊ２に基づいて、前記解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布Ｐ（Ｘ^＊ _ｉ）を効率的に得ることができる。

さらに、本発明の第６の態様である解析パラメータの推定方法は、前記事後確率分布Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）に基づく、前記パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、前記定式化された前記Ｘ_ｉと前記Ｙ_ｊとの関係性に入力することにより得られる、前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）＝Ｐ（Ｙ_ｊ２｜Ｘ^＊ _ｉ）＝Ｐ（Ｙ_ｊ２｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）について、前記事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、いずれかのｊ_２について所定の閾値を上回る場合に、前記観測変数の数Ｍ_１、Ｍ_２および／ないしＭを増やすことにより、前記事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、いずれのｊ_２についても所定の閾値を下回るように構成したことを特徴とする。

このように構成すると、前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が所定の閾値を下回るように、すなわち目的とする数値シミュレーションの実行結果において前記評価変数Ｙ_ｊ２を所望の解析精度において得るために必要な精度において前記の解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布Ｐ（Ｘ^＊ _ｉ）を得ることができる。

一実施形態による未知の解析パラメータの推定方法のステップを示す図である。 境界条件としての燃焼装置入口境界での燃料ガスの化学組成の分布が未知である場合における燃焼シミュレーションの一実施例を示す。 解析パラメータＸ_１からＸ_１０までの事後確率分布Ｐ（Ｘ^＊ _ｉ）の推定結果の一例を示す。 各観測点および評価点におけるガス温度の事後確率分布Ｐ（Ｙ^＊＊ _１）、Ｐ（Ｙ^＊＊ _５）、Ｐ（Ｙ^＊＊ _９）、Ｐ（Ｙ^＊＊ _１３）、Ｐ（Ｙ^＊＊ _１７）の一例を示す。

以下、図面を参照して、本発明の好適な実施形態について詳細に説明する。

図１は、本発明の一実施形態による未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）の推定方法のステップを示している。

はじめに、ステップ１として、解析対象に観測点および評価点を設定し、観測点における観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）と、評価点における評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）を定義する。解析対象は、例えば、ごみ焼却プラントやバイオマス燃焼プラントなどを含むが、これらに限られず、装置や機器の設計の過程において数値シミュレーションが用いられるあらゆる装置や機器が、解析対象となりうる。観測点の数と観測変数の数との関係は任意である。例えば、Ｍ_１個の観測点に各１つずつ合計Ｍ_１個の観測変数を定義してもよいし、１個の観測点にＭ_１個の観測変数を定義してもよいし、５個の観測点に各Ｍ_１／５個ずつ合計Ｍ_１個の観測変数を定義してもよい。

後述するように、本実施形態における未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）の推定は、観測点において得られる観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）の統計情報に基づいてなされるため、原理的には観測変数Ｙ_ｊの数Ｍが大きいほど、精度よく解析パラメータＸ_ｉを推定することができるが、一方で解析パラメータＸ_ｉを推定するための計算量が増大するという問題が生じる。

そのため本実施形態においては、観測変数Ｙ_ｊの数Ｍをできるだけ少なくした上で精度よく解析パラメータＸ_ｉを推定するために、以下に記すように、観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）を、観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）と、評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）の２つに分けて定義する。

観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）については、対象とする物理現象を考慮して、推定すべき未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）との相関が高い観測量として定義する。観測変数Ｙ_ｊ１の数Ｍ_１を適切に決定する方法については、後述する。

評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）については、対象とする数値シミュレーションにおいて、計算結果の評価を行う上で着目する変数、すなわち計算結果の精度を評価する対象となる変数として定義する。評価変数Ｙ_ｊ２の数Ｍ_２は、対象とする数値シミュレーションにおいて計算結果の評価を行う上で必要となる数とすればよい。

例えば前記の燃焼装置内の温度分布を評価する問題において、燃焼装置の入口境界での燃料ガスの流速および化学組成の分布が未知の解析パラメータである場合、観測変数Ｙ_ｊ１としては燃焼装置の入口境界にできるだけ近く、物理的に観測が可能な位置における装置内のガス流速および／または化学組成の計測値を採用することができる。また、評価変数Ｙ_ｊ２としては、装置内の代表的な箇所であり、かつ物理的に観測が可能である箇所における温度の計測値や、装置出口境界でのガス流速および／または化学組成の計測値を採用することができる。

次に、ステップ２として、解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）を変化させた多数の数値シミュレーションを実行する。以下では、ステップ２において、解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）をＫ通りに変化させた数値シミュレーションを実行することとし、ｋ番目の計算条件に対応する解析パラメータＸ_ｉの値をＸ_ｉ ^ｋ（ｋ＝１，２，・・・，Ｋ）と記述する。また、ｋ番目の計算条件での数値シミュレーション結果における観測変数Ｙ_ｊ１および評価変数Ｙ_ｊ２の値を、それぞれＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋと記述する。計算条件の数Ｋを適切に決定する方法については、後述する。

次に、ステップ３として、ステップ２で実行した数値シミュレーションの実行結果に基づいて、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性を定式化する。

本実施形態においては、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性は、解析パラメータＸ_ｉについての多変数多項式回帰により定式化される。ここで、多変数多項式回帰式の次数については、模擬対象となる物理現象の複雑さや、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性の複雑さに応じて決めるのがよい。

例えば、多変数多項式回帰式の次数を２とした場合において、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性は次式のように定式化することができ、数値シミュレーションの計算条件であるＸ_ｉ ^ｋと、その実行結果であるＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋに基づいて、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性を最もうまく説明できる次式の係数ａ、ｂ、ｃおよびｄを、最小二乗回帰や最尤推定法などの手法によって求める回帰問題に帰着する。多変数多項式回帰式の次数を１とした場合や、３以上とした場合についても、同様に定
式化することができる。

しかしながら、解析パラメータＸ_ｉの数が大きいと、回帰係数ａ、ｂ、ｃおよびｄの数が増大する。あるいは、多変数多項式の次数が大きくなった場合にも、求めるべき回帰係数の数が増大する。その場合、後述するステップ６において、観測変数Ｙ_ｊ１および評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報に基づいて、解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布を得る処理の計算負荷が増大してしまうという問題が生じる。

この問題に対処するため、本実施形態においては、ステップ３における解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性の定式化に際して、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて説明変数の選択が行われる。

スパース推定手法とは、回帰モデルの損失関数に回帰係数のノルムに基づく正則化項を加えた正則化損失関数を最小化することにより、変数の数が多い場合において、推定の安定化とともに、真に意味のある変数を選択することが可能となる手法であり、リッジ回帰、ラッソ回帰またはエラスティックネット等の手法が広く用いられている。これらのスパース推定手法については公知であるのでその詳細は省略する。

次に、ステップ４として、スパース推定手法による多変数多項式回帰により定式化された解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の回帰式の精度を評価する。具体的には、ステップ２において実行した数値シミュレーションの全ての計算条件Ｘ_ｉ ^ｋについて、ステップ３で得られた回帰式に各Ｘ_ｉ ^ｋを代入することにより得られる予測値Ｙ^〜 _ｊ１ ^ｋおよびＹ^〜 _ｊ２ ^ｋの値と、シミュレーション結果として得られたＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋの値とを比較する。

回帰式の精度の評価指標としては、例えば上記のＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋの値と、対応するＹ^〜 _ｊ１ ^ｋおよびＹ^〜 _ｊ２ ^ｋの値との二乗平均誤差を用いることができる。二乗平均誤差の値が所定の閾値よりも大きい場合は、シミュレーションの実行ケース数Ｋを増やすか、または多変量多項式における次数を増やした上で、再度ステップ２からステップ４を実行することにより、回帰式の精度を高めるのがよい。

次に、ステップ５として、観測点における観測変数Ｙ_ｊ１の観測値Ｙ^＊ _ｊ１および評価点における評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ２を実験的に得る。

次に、ステップ６として、ステップ５で得られた観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報と、ステップ３で定式化された解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性とに基づいて、ベイズ推定を行うことにより、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２が得られたもとでの、解析パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を求める。

本ステップ６の目的は、実験的に得られた観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の情報から、それらを与える解析パラメータＸ^＊ _ｉを逆推定することである。

ここで、ステップ５において得られる観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２は、実際の装置または機器において実験的に計測される量であるため、実験上あるいは計測上の様々な要因を反映した観測誤差を含んでいる。そのため、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２は確定値として得られるのではなく、バラつきをもった統計分布として得られることが通常である。

観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の分布におけるバラつきが小さい場合や、その分布が正規分布に従う場合は、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の各々の分布の平均値ないし中央値をそれらの代表値として確定的に用いることができる。この場合、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性が線形関係にある場合は、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２を与える解析パラメータＸ^＊ _ｉを解析的に算出することができる。

しかしながら、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性がより複雑である場合、例えばこれらが２次以上の多項式で表されるような非線形関係にある場合は、前記のように観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の各々の分布の平均値ないし中央値をそれらの代表値として確定的に用いることができる場合であっても、それらの観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２を与える解析パラメータＸ^＊ _ｉを解析的に算出することは不可能であることが多い。

あるいは、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の分布におけるバラつきが大きい場合や、その分布が正規分布に従わない場合には、そもそも観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の代表値を確定的に定義することができないため、たとえ解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性が線形関係にあったとしても、そのような観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２に対応する解析パラメータＸ^＊ _ｉを解析的に算出することは不可能である。

そこで本実施形態による解析パラメータの推定方法は、ベイズ推定を用いることにより、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報に基づいて、それらを与える解析パラメータＸ^＊ _ｉの確率分布を、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２が得られたもとでの、解析パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）として求めるように構成されている。

具体的には、本実施形態においては、ステップ６におけるベイズ推定手法として、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法が用いられている。

マルコフ連鎖モンテカルロ法は、ベイズ推定の枠組みにおいて、事後確率分布の計算を数値的に行うことができる手法として広く用いられており、メトロポリス−ヘイスティングアルゴリズム、ギブスサンプラーなどの手法の総称である。またハミルトニアンモンテカルロ法は、マルコフ連鎖モンテカルロ法における目標分布への収束性の改善を目的として近年開発され広く利用されている。これらの手法については公知であるのでその詳細は省略する。

これらのベイズ推定手法を用いることにより、前記多変数多項式における変数の数および／または項の数が多い場合や、あるいは解析パラメータＸ_ｉと観測変数Ｙ_ｊ１および／または評価変数Ｙ_ｊ２との関係性が強い非線形性を有する場合や、さらには観測変数Ｙ_ｊ１および／または評価変数Ｙ_ｊ２の統計分布が正規分布に従わないような場合など、観測変数Ｙ_ｊ１および評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報から直接的に前記解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉを推定することが困難な場合であっても、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報に基づいて、解析パラメータＸ_ｉの事後確率分布Ｐ
（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を数値的に求めることが可能となる。

次に、ステップ７として、ステップ６で得られた解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、ステップ３で定式化された解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性に入力することにより、評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）＝Ｐ（Ｙ_ｊ２｜Ｘ^＊ _ｉ）＝Ｐ（Ｙ_ｊ２｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を得る。

次に、ステップ８として、ステップ７で得られた評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間の評価を行う。

前述のように、本実施形態における評価変数Ｙ_ｊ２は、計算結果の評価を行う上で着目する変数、すなわち計算結果の精度を評価する対象となる変数として定義されている。したがって、観測変数Ｙ_ｊ１および評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２が得られたもとでの、評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、所定の閾値よりも小さくなっていれば、対象とする数値シミュレーションにおいて、所望の解析精度を得るために必要な精度において未知の解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）が求まったといえることになる。

ここで、評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、所定の閾値よりも小さくならなかった場合は、対象とする数値シミュレーションにおいて、解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、所望の解析精度を得るために必要な精度で求めるためには、情報が不足していることを意味している。したがって、そのような場合は、観測変数Ｙ_ｊ１の数Ｍ_１を増やした上で、再度ステップ１からステップ８までの手順を実行する。これにより、観測変数Ｙ_ｊ１から得られる情報を増やすことができ、その結果として、数値シミュレーションにおいて所望の解析精度を得るために必要な精度で、未知の解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を求めることができる。

以下、図面を参照して、本発明の実施例について説明する。

図２は、ある燃焼装置内の温度分布やガス組成分布を評価するために行う燃焼シミュレーションにおいて、境界条件としての燃焼装置入口境界での燃料ガスの化学組成の分布が未知である場合の一実施例を示している。

ここでは、図２に示されているように４つに分割されている燃焼装置の入口境界の各々の領域Ｆ_１、Ｆ_２、Ｆ_３、Ｆ_４における燃料ガスＧ_１、Ｇ_２、Ｇ_３、Ｇ_４の組成を表現するためのパラメータとして、合計２０個の解析パラメータＸ_ｉ（ｉ＝１〜２０）が定義されている。すなわち、ｆをある関数として、Ｇ_ｎ＝ｆ（Ｘ_ｉ）（ｎ＝１〜４、ｉ＝１〜２０）であり、これらＸ_ｉの値を、上述した実施形態による解析パラメータの推定方法によって推定するものとする。また、本実施例における数値シミュレーションの要求精度としては、燃焼装置出口のガス温度を±１０℃の精度で求めることが要求されているものとする。

はじめに、ステップ１として、観測点および観測変数、ならびに評価点および評価変数を設定する。具体的には、燃焼装置の入口境界の各領域Ｆ_１、Ｆ_２、Ｆ_３、Ｆ_４の直上であって、かつ燃焼装置内のガス温度およびガス組成の計測が実際に可能な箇所に、４つの観測点Ｒ_１、Ｒ_２、Ｒ_３、Ｒ_４を定義し、それらの位置における燃焼装置内のガス温度およびガス組成（未燃ガス濃度、酸素濃度、水分濃度の３つ）を観測変数として定義する。
また、燃焼装置の出口境界に１つの評価点Ｒ_５を定義し、その位置におけるガス温度およびガス組成（同じく未燃ガス濃度、酸素濃度、水分濃度の３つ）を評価変数として定義する。

つまり、１つの観測点ないし評価点において、４つずつの観測変数ないし評価変数が定義される。以下では本実施例における観測変数Ｙ_ｊ１（ｊ_１＝１〜１６）および評価変数Ｙ_ｊ２（ｊ_２＝１７〜２０）を、それぞれ次の表１のように定義する。

次に、ステップ２として、解析パラメータＸ_ｉ（ｉ＝１〜２０）を変化させた多数の数値シミュレーションを実行する。具体的には、ステップ３で各観測変数Ｙ_ｊ１（ｊ_１＝１〜１６）および評価変数Ｙ_ｊ２（ｊ_２＝１７〜２０）を各解析パラメータＸ_ｉの２次多変数多項式として表現することで応答曲面の定式化を行うことを前提に、２次多変数多項式の係数の数（２３１個）を上回るように合計２５０点（ｋ＝２５０）の数値シミュレーションを実行した。

次に、ステップ３として、ステップ２で実行した数値シミュレーションの実行結果に基づいて、解析パラメータＸ_ｉ（ｉ＝１〜２０）と、観測変数Ｙ_ｊ１（ｊ_１＝１〜１６）および評価変数Ｙ_ｊ２（ｊ_２＝１７〜２０）との関係性を定式化する。ここでは、２次の多変数多項式を用い、さらにラッソ回帰によるスパース推定手法を適用することで、説明変数の選択を行った。変数選択の結果の詳細は省略するが、もとの多変量多項式の係数の数（２３１個）に対して、選択された変数の数は各観測変数および評価変数に対して、１０個から９０個程度となり、効果的に変数選択が行われたことが分かる。

次に、ステップ４として、上記により定式化された解析パラメータＸ_ｉ（ｉ＝１〜２０）と観測変数Ｙ_ｊ１（ｊ_１＝１〜１６）および評価変数Ｙ_ｊ２（ｊ_２＝１７〜２０）の回帰式の精度を評価する。具体的には、ステップ２において実行した数値シミュレーションの全ての計算条件Ｘ_ｉ ^ｋ（ｉ＝１〜２０、ｋ＝１〜２５０）について、ステップ３で得られた回帰式に各Ｘ_ｉ ^ｋを代入することにより得られる予測値Ｙ^〜 _ｊ１ ^ｋおよびＹ^〜 _ｊ２ ^ｋの値と、シミュレーション結果として得られたＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋの値との二乗平均誤差を比較し、十分な精度が得られていることを確認した。より高い精度が必要である場合には、図１に示したフローチャートに従って、適宜数値シミュレーションの実行回数な
いし多変量多項式の次数を増やした上で、再度ステップ２からステップ４までの手順を実行すればよい。

次に、ステップ５として、観測点における観測変数Ｙ_ｊ１の観測値Ｙ^＊ _ｊ１および評価点における評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ２を実験的に得た。

次に、ステップ６として、ステップ５で得られた観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報と、ステップ３で定式化された解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性とに基づいて、ベイズ推定を行うことにより、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２が得られたもとでの、解析パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _ｊ）をベイズ推定手法（ハミルトニアンモンテカルロ法）により求めた。

その結果の一例として、解析パラメータＸ_１からＸ_１０までの事後確率分布Ｐ（Ｘ^＊ _ｉ）の推定結果を図３に示す。これらの図に示されるように、実験的に得られる観測値の統計情報に対応する解析パラメータの値の事後確率分布は、一般に幅広い分布を示すことが分かる。

次に、ステップ７として、ステップ６で得られた解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、ステップ３で定式化された解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性に入力することにより、評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）＝Ｐ（Ｙ_ｊ２｜Ｘ^＊ _ｉ）＝Ｐ（Ｙ_ｊ２｜Ｙ^＊ _ｊ）を得た。その結果の一例として、各観測点および評価点におけるガス温度の事後確率分布Ｐ（Ｙ^＊＊ _１）、Ｐ（Ｙ^＊＊ _５）、Ｐ（Ｙ^＊＊ _９）、Ｐ（Ｙ^＊＊ _１３）、Ｐ（Ｙ^＊＊ _１７）を図４に示す。

前述のように、本手法により同定された未知の解析パラメータＸ_ｉ（ｉ＝１〜２０）の事後確率分布Ｐ（Ｘ^＊ _ｉ）は、図３に示すように幅広い分布を示すものの、各観測点および評価点におけるガス温度の事後確率分布は、図４に示すように観測変数および評価変数の実際の観測値（図中の太線）を含むシャープな分布となっており、それらの実際の観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２をうまく説明できていることが分かる。

次に、ステップ８として、ステップ７で得られた評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間の評価を行う。前述したように本実施例では、評価変数である燃焼装置出口のガス温度を±１０℃の精度で求めることが要求されていたが、図４で示される燃焼装置出口のガス温度Ｙ_１７の事後確率分布は、真値（実際の測定値）８３６℃に対して、約８３４〜８４１℃の範囲にシャープに分布しており、十分な精度が得られていることが分かる。より高い精度が必要である場合には、図１に示したフローチャートに従って、適宜観測変数の数を増やした上で、再度ステップ１からステップ８までの手順を実行すればよい。

以上の実施例では、本発明を燃焼装置における未知のガス組成の同定に適用しているが、本発明はこの例に限定されず、その他の様々な工業的装置あるいは機械における未知の解析パラメータの推定、例えば機器や装置における温度分布や濃度分布の同定、境界表面力の同定、電流密度の同定などにも用いることができる。

以上、本発明の実施形態および実施例についてその詳細を説明したが、本発明は上記の実施形態および実施例に限定されるものではなく、特許請求の範囲、及び明細書と図面に記載された技術的思想を逸脱しない範囲において種々変更又は修正を行って実施することが可能である。

Claims (6)

1. 物理現象を模擬するための数値シミュレーションを所望の解析精度において実行するために必要な、１つ以上の未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）の推定値Ｘ^＊ _ｉを推定する方法であって、
   対象とする物理現象における、１つ以上の観測点における１つ以上の観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）について、前記１つ以上の解析パラメータＸ_ｉの値を変化させた複数の数値シミュレーションを実行するステップと、
   前記複数の数値シミュレーションを実行して得られた結果に基づいて、前記解析パラメータＸ_ｉと前記観測変数Ｙ_ｊとの関係性を定式化するステップと、
   前記物理現象について実験的計測を行うことにより、前記観測変数Ｙ_ｊの観測値Ｙ^＊ _ｊ（ｊ＝１，２，・・・，Ｍ）を得るステップと、
   前記観測変数Ｙ_ｊの観測値Ｙ^＊ _ｊの統計情報と、前記定式化された前記解析パラメータＸ_ｉと前記観測変数Ｙ_ｊの関係性とに基づいて、ベイズ推定を行うことにより、前記観測値Ｙ^＊ _ｊが得られたもとでの、前記解析パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）として、前記解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布を得るステップと、
   を含むことを特徴とする、解析パラメータの推定方法。
2. 前記解析パラメータＸ_ｉと前記観測変数Ｙ_ｊとの関係性は、前記観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）を目的変数とし、前記解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）を説明変数とした多変数多項式回帰により定式化されることを特徴とする、請求項１に記載の解析パラメータの推定方法。
3. 前記多変数多項式回帰において、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて前記説明変数の選択を行うことを特徴とする、請求項２に記載の解析パラメータの推定方法。
4. 前記ベイズ推定において、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法を用いることを特徴とする、請求項１から３のいずれか１項に記載の解析パラメータの推定方法。
5. 前記１つ以上の観測点は、前記解析パラメータＸ_ｉとの相関性の高い観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）が得られるように選ばれた観測点と、前記数値シミュレーションにおいて計算結果の評価に用いられる評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）が得られるように選ばれた観測点とからなることを特徴とする、請求項１から４のいずれか１項に記載の解析パラメータの推定方法。
6. 前記解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、前記定式化された前記解析パラメータＸ_ｉと前記観測変数Ｙ_ｊとの関係性に入力することにより、前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）＝Ｐ（Ｙ_ｊ２｜Ｘ^＊ _ｉ）＝Ｐ（Ｙ_ｊ２｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を得るステップと、
   前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、いずれかのｊ_２について所定の閾値を上回る場合に、前記観測変数Ｙ_ｊ１の数Ｍ_１を増やすことにより、前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、いずれのｊ_２についても前記所定の閾値を下回るようにするステップと、
   をさらに含むことを特徴とする、請求項５に記載の解析パラメータの推定方法。
   

Images